key: cord-0625553-q3l0iik3
authors: Li, Shuolin; Henriquez, Craig; Katul, Gabriel
title: The Role of Geographic Spreaders in Infectious Pattern Formation and Front Propagation Speeds
date: 2021-08-23
journal: nan
DOI: nan
sha: f5c716f415fed9b34c7bb86655e04bd822e96d12
doc_id: 625553
cord_uid: q3l0iik3

The pattern formation and spatial spread of infectious populations are investigated using a kernel-based Susceptible-Infectious-Recovered (SIR) model applicable across a wide range of basic reproduction numbers $R_o$. The focus is on the role of geographic spreaders defined here as a portion of the infected population ($phi$) experiencing high mobility between identical communities. The spatial organization of the infected population and invasive front speeds ($c_{max}$) are determined when the infections are randomly initiated in space within multiple communities. For small but finite $phi$, scaling analysis in 1-dimension and simulation results in 2-dimensions suggest that $c_{max}sim (1-phi) gamma (R_o-1) sigma$, where $gamma$ is the inverse of the infectious duration, and $sigma^2$ is the variance of the spatial kernel describing mobility of long-distance spreaders across communities. Hence, $c_{max}$ is not significantly affected by the small $phi$ though reductions in $phi$ act as retardation factors to the attainment of $c_{max}$. The $sigma$ determines the spatial organization of infections across communities. When $sigma>5dr$ (long-distance mobility, where $dr$ is the minimum spatial extent defining adjacent communities), the infectious population will experience a transient but spatially coherent pattern with a wavelength that can be derived from the spreading kernel properties.

Spatial models of infectious disease transmission remain the optimal viable 'in-silico' tactic where the movement of infected individuals can be combined with a quantitative description of the infection process and disease properties to explore appropriate intervention responses [Riley, 2007 , Riley et al., 2015 , Grassly and Fraser, 2008 , Sanz et al., 2014 , Valdano et al., 2015 , Keeling and Rohani, 2011 . The mobility of infectious individuals in a susceptible population may lead to a recognizable spatial pattern, which when understood, can serve as a basis for concentrating testing capacity or planning control strategies. For this reason, the emergence of coherent spatial patterns in epidemiology, even in the most idealized settings where the pattern is endogenous to the model system describing the spread, continue to draw attention [Murray, 2001 , Eisinger and Thulke, 2008 , Boerlijst and Van Ballegooijen, 2010 , Sun et al., 2016 , Soriano-Paños et al., 2018 .

Pattern formation and deterministic models (e.g. activation-inhibition, reaction-diffusion, networking) in infectious disease have a long tradition [Belik et al., 2011 , Soriano-Paños et al., 2018 , Koher et al., 2019 , and reviewing all their aspects is well beyond the scope of a single study. However, these approaches delineate two broad classes of spatial patterns [Murray, 2001 , Fuentes et al., 2003 , Sun et al., 2016 . The first is stationary patterns that remain unchanged or stable over an extended period of time when referenced to the duration of the outbreak. These patterns often exhibit intermittent areas of a high density of infectious individuals. The locally concentrated areas may favor persistence of a disease and are of logical concern to epidemiology.

The other, often labeled as spatio-temporal patterns, change over time either as (i) quasi-periodic or oscillatory waves [Sherratt and Smith, 2008] or as (ii) spatio-temporal (intermittent) chaos [Fuentes et al., 2003 , Sun et al., 2016 . The work here shows that spatial patterns of infections with dynamical features that share some resemblance to these two categories can emerge. These features depend on the nature of the local competitive interaction (e.g. mass action within a community) as well as on the spatial extent of the competitive range.

Compartment models with mass action allow infectious individuals to access all the susceptible population within an isolated community [Girvan and Newman, 2002] . For clarity, we label the space occupied by such a community as box so as to distinguish it from the term 'compartment' commonly used when classifying population into susceptible, infected, or recovered. The spatial domain is assumed to be much larger than the size of the box so that the box center may be identified by coordinates representing the box location. This assumption allows space to be treated as a continuous variable. Long-distance mobility of infectious individuals from one box to another, however, leads to interaction with another susceptible population that is non-local in space, which we naively label as geographical spreaders. The transport time of infectious individuals is assumed to be much faster than recovery or removal time scale so that time may also be treated as a continuous variable. The geographical spreaders share some similarity with the commonly used term 'superspreader' [Small et al., 2006] . A precise definition of superspreader remains elusive and often defining syndromes are used. In its broadest epidemiological sense, a superspreader has the propensity to infect a larger than an average number of susceptible individuals (where mobility is only one of several factors resulting in elevated infectious rates). Pragmatically, a superspreader is an individual who can infect more individuals than say a number constrained by the basic reproduction number R o [Galvani and May, 2005] . From this pragmatic perspective, geographical spreaders and superspreaders share common attributes but geographical spreaders do not encompasses biological, behavioral, and environmental variables relevant to disease transmission. For this reason, we only focus on infectious individuals with high mobility to distinguish them from the widely used superspeader term. This focus on spatial spread of individuals also ignores numerous other factors such as those associated with particular events where a large congregation of individuals in a small area (colloquially labeled a superspreader event) facilitate infectious spread [Kuebart and Stabler, 2020] . The choice to associate high mobility only with geographical spreaders in compartment models simplifies the mathematical treatment for continuous space-time models. It leads to well-recognized integro-differential equations (IDE) [Lefever and Lejeune, 1997 , Fuentes et al., 2003 , 2004 , Guo et al., 2020 , Keeling and Rohani, 2011 , Murray, 2001 ] that may be viewed as a hybrid between local diffusion-based schemes [Lang et al., 2018] and mobility on networks [Pastor-Satorras et al., 2015] , where long-distance movement is ubiquitous.

In the IDE approach, mobility of geographical spreaders is accommodated using a spatial kernel that describes the probability p(x − x , y − y ) of infectious individuals to travel from point defined by Cartesian coordinates (x, y) to (x , y ). As such, detailed knowledge about long-range connections can, in principle, be accommodated using non-homogeneous spreading kernel functions. By "non-homogeneous", we mean that the spreading kernel is no longer dependent on relative distance between points (x,y) and (x',y') but also on the position or origin (x,y). This situation arises when the statistics of the spreading kernel (e.g. kernel variance) vary with position instead of relative distance r = (x − x ) 2 + ( y − y ) 2 ). Likewise, spreading kernels can be non-isotropic with certain preferential pathways breaking the symmetry of the movement. The

to reflect different mobility patterns during the course of a day (daytime versus night-time), week (weekday versus weekends), or season (summer versus winter) relative to a reference time t o . The IDE approach proved effective in modeling fast propagation fronts of biological systems such as vegetation migration e.g. invasion by extremes, vegetation patterning in dryland ecosystems, the 'chic and hen' pattern in seedlings around adult trees in tropical ecosystems, among others [Lefever and Lejeune, 1997 , Thompson and Katul, 2008 , Thompson et al., 2008 , Meron, 2018 . The IDE approach proposed here is shown to exhibit qualitatively similar invasion thresholds reported in reaction-diffusion processes operating in meta-population models with heterogeneous connectivity patterns [Colizza and Vespignani, 2007] . Moreover, the approach can accommodate (indirectly) bidirectional movements between 'base' and 'destination' locations, which have been shown to generate some saturation attack velocities in mobility networks [Belik et al., 2011 , Soriano-Paños et al., 2018 . More recently, it has been demonstrated that the IDE approach is able to reproduce measured spectral (and multi-fractal) space-time (and intermittent) properties of COVID-19 infections data for a full year across a wide range of scales (from county to continent) [Geng et al., 2021] . However, the mathematics that govern the spread and recession speeds and the concomitant pattern formation during the outbreak have not been investigated and motivate the work here.

The main objective then is to use the IDE approach to establish links between the properties describing the spatial kernel of geographical spreaders and (i) transient spatial pattern formation of infectious individuals during the outbreak, and (ii) invasive and recession speeds for the disease.

These links are tested across a wide range of basic reproduction numbers. The specific questions to be addressed are: (i) What spatial patterns arise when multiple communities are randomly infected resulting in interfering wavefronts and (ii) what is the role of the spreading kernel properties (its variance and spatial extent) and the fraction of geographical spreaders in a community on spatial patterning of infectious individuals. The overall impact of geographical spreaders on compartment models employing only mass action when describing encounters between infectious and susceptible populations is also analyzed and discussed.

A brief review of standard compartment models and their extension to include spatial mobility by diffusion are covered to introduce nomenclature. The use of IDE to expand diffusion schemes along with key mathematical properties for spatial pattern formation and front propagation are then developed for simplified kernel shapes.

A generic mathematical model of disease outbreak, known as the Susceptible-Infectious-Recovered (SIR) model, subdivides the total population N o residing in a closed box into 3 distinct compartments: susceptible (S), infectious (I), and recovered or removed (R). Infectious individuals spread the disease to a susceptible individual and remain in the infectious compartment for a given period before moving into the recovered (or removed) compartment. Individuals in the recovered or removed compartment are assumed to be immune (or removed) from the population.

The dynamical system describing the SIR equations is given as [Kermack and McKendrick, 1927 , Hethcote, 2000 , Murray, 2007 , Keeling and Rohani, 2011 , can be evaluated from the initial condition (t = 0).

The SIR model makes several assumptions that include no large imbalances in natural births or natural deaths during the short-lived outbreak. Direct transmission occurs through individualto-individual contact represented by a coefficient of β. The infection is assumed to have a zero latent period so that an individual becomes infectious as soon as they become infected. However, the most restrictive assumption in SIR dynamics is the mass action of individuals. Mass action law assumes that the rate of an encounter between I and S is proportional to their product. This requires that infected and susceptible individuals be uniformly distributed within the space of the box. It may be argued that for high I and low S densities, the mass action law must be revised to accommodate saturation effects due to competition for space (as infected individuals become reasonably dense, most infected individuals are surrounded by infected individuals thus having lower possibility of interacting with the susceptible individuals). This saturation or second-order effect is ignored here. The γ and β encode the main properties of the epidemics and the population response to it and will be assumed constant in time and space.

For illustration purposes only, we use values for γ based on those used in models of the dynamics of the recent COVID-19 pandemic i.e. the time from the onset of symptoms to natural recovery is, on average, set to γ −1 = 14 d [Katul et al., 2020] . To be clear, the goal is not to explore any particular spreading event of COVID-19 for a problem-specific setting. The choice of a γ −1 being of order 10 days is selected because geographical mobility usually occurs on time scales much faster than γ −1 . The remaining model parameter β must be determined empirically or from separate studies, which can pose difficulties [Keeling et al., 2009 , Turkyilmazoglu, 2021 . It is evident 

In a closed system which we label as box, the SIR model prohibits mobility of individuals across the box boundaries. Mathematical approaches extending the SIR model from a single box (that may be experiencing mass action mixing) into a Cartesian (x, y) plane represent the spatial domain by a square-lattice with characteristic size L and area A d = L 2 . This domain is then filled with individual boxes of dimensions d x by d y with each box representing a community characterized by an initial population N o . The center of an arbitrary box is defined by its coordinates x, y.

The basic SIR dynamics (i.e. conversion of S to I to R) along with mass action only applies within the box positioned at (x, y). The size of each box is assumed to be much smaller than the domain size to allow continuous treatment of space. Each box is subject to the initial condition

The concern here is the spatial mobility of I across boxes (geographical spreaders) and its subsequent spatial patterns when few infectious individuals randomly occur in few boxes within the domain. Because infectious individuals can be mobile, conventional spatial models revise the original SIR to allow for their mobility using a diffusion term given as

where D o is a diffusion coefficient that must be a priori known. This representation is continuous in space and time, meaning d x and d y are infinitely small -at least when compared to the domain size. Such systems are known to admit travelling fronts whose stationary speed c max is, to a leading order, given by the early times dynamics (when S ≈ N o ) [Fisher, 1937 , Volpert and Petrovskii, 2009 , Murray, 2001 , Belik et al., 2011 

when R o > 1. This relation between c ma x and D o γ(R o − 1) can also be derived from dimensional considerations alone. It is to be noted that equation 3 requires (R o − 1) > 0, which agrees with the standard SIR model without a spatial component where R o > 1 was necessary for an epidemic to form. To allow for non-local mobility of geographical spreaders, we formulate the I budget as

with the convolution function C u (x, y) defined as

where φ is the fraction of infectious individuals (geographical spreaders) that are mobile across boxes within a given time step d t γ −1 and p(x , y ) is a spatial spread kernel satisfying the normalizing condition

The conceptual difference between the two models in equations 2 and 4 is explained in Figure 1 .

When φ = 0, the IDE approach reduces to spatially independent or autonomous SIR models operating in closed communities with no connectivity or spatial interaction between boxes (i.e.

immobile model). When φ = 1, all infectious individuals experience spatial mobility per unit time. The fact that the spread kernel acts directly on I is in keeping with recent findings derived from network theory [Brockmann and Helbing, 2013 ] that suggested the number of infectious individuals is more significant than effective distances between connected areas.

Several objections can be raised when setting φ to a finite constant. The first is that N o is no longer a constant within a given box. The second is that spatial mobility is exclusive to I not S and R, which is not realistic. However, a recent analysis on a similar IDE model has demonstrated that the spatial mobility of I is far more critical to disease spread than the spatial mobility of S and R [Geng et al., 2021] . To arrive at key analytical results, the analysis is restricted to small φ( 1) values to illustrate how few geographical spreaders impact pattern formation and infectious fronts.

There are a number of choices that can be made about p(x , y ) [Jo et al., 2014 , Gleeson et al., 2016 , Keeling and Rohani, 2011 . For simplicity and to recover c max predicted by diffusion, we selected

where

representing the box centers, σ is a measure of the spread of the spatial kernel (to be related to D o ), and α N is a normalizing constant. Since the interest here is in spatial spread kernels with Fuentes et al., 2003 , 2004 , Da Cunha et al., 2009 

This condition yields the normalizing constant

To maintain a near-Gaussian shape for the spreading kernel, we selected R a /(2σ) > 1. This IDE approach along with a Gaussian spreading kernel was recently employed to explain the multifractal properties of COVID-19 infections from scales spanning 10 km to 2800 km across the U.S. [Geng et al., 2021] and from 1 week to a 1 year period when using γ −1 = 14 d. Other spatial kernels can also be specified and subjected to the same normalizing conditions thereby making the IDE approach flexible in terms of choices about spatial spread of infectious individuals across communities.

The IDE system exhibits two properties relevant to the objectives here and are derived in 1-D chosen along x for illustration. The first is the existence of travelling waves and the determination of their speed. The second considers necessary (but not sufficient) conditions for pattern formation using stability analysis around the only equilibrium point (I = 0). Both properties are analyzed in a 1-D case to show-case expected connections between these properties, φ, the spatial kernel shape, and epidemiological constants R o and γ. The derivation is also linearized around the early times (i.e. infections are increasing exponentially) when S(x, t) ≈ S(x, 0) so that the IDE can be Page 8 of 28

The convolution between I and p in Equation 10 is now simplified as

For establishing links between the diffusion representation leading to equation 3 and the 1-D IDE system, the travelling wave speed is now considered. A variable transformation z = x − c t along with d x = dz is used to convert the 1-D IDE in equation 10 to the new coordinate system, given

Assuming the solution is of the form I(z) = Aexp(λz) and using a zero-mean Gaussian kernel p(z), two Fourier functions can be used to aid in the evaluation of C u (z) and are defined aŝ

where A is a constant, k is the wavenumber, i 2 = −1, δ(.) is the Dirac-delta function, andŝ denotes the corresponding Fourier transform of function s(z). To evaluate C u (z), the product ofÎ(k) andp(k) is determined in the Fourier domain followed by an inverse Fourier transform

defines the moment-generating function of a zero-mean Gaussian kernel with variance σ 2 . Substituting this estimate of C u (z) into equation 12 and upon further simplification, a characteristic equation for λ is derived as

The interest here is in monotonic (rather than oscillatory) wave fronts and real (rather than complex) roots. Real roots emerge as multiple roots at the nth-order contact [Kot et al., 1996 , Thompson and Katul, 2008 , Volpert and Petrovskii, 2009 , Tian and Yuan, 2017 given by differentiating equation 15 with respect to λ. Noting that d M g /dλ = λσ 2 M g and differentiating equation

15 with respect to λ yields

or as

When R o > 1 (outbreak), equation 16 requires that c/λ < 0 (i.e. c and λ must have opposite signs reflecting expansion and contraction phases of the disease spread). Combining equations 15 and 16 yields a c given by

whose solution is

where c k = λσ is a control factor that varies with the initial conditions. This solution admits two traveling wave velocities (i.e. two λ values) depending on the initial and boundary conditions imposed on the IDE: a forward velocity (when λ < 0) and a backward velocity (when λ > 0).

Conventional approaches to by-pass c k and determine a single c max assume that the discriminant of equation 18 with respect to λ is zero to ensure double roots [Kot et al., 1996] . In this application, the discriminant ∆ is given by

and demonstrates that no double root exists (i.e. disease expansion and contraction fronts must both occur).

When φ 1, as is the case here, c ma x appears to be insensitive to φ though 1 − φ does act as a retardation factor to the attainment of c max (discussed later). Unlike the diffusion model, the kernel-based c ma x varies linearly with γ(R o − 1) instead of γ(R o − 1). The c max also scales linearly with the kernel spread measure σ. Equating the two speeds from equations 3 and 19

yields an effective D o ,

Thus, diffusion models can be made analogous to the IDE approach only for the purposes of estimating stationary front speeds. The variables in equation 19 are derived here mainly to guide the numerical model analysis for invasion front speeds (in 2D) and for differing initial conditions for disease attack. The analytical solution in 1-D here must only be viewed as informing a scaling relation between travelling wave speeds, epidemiological parameters, and kernel spread properties for more complex initial conditions and in a higher dimension.

As common to pattern formation studies, linear stability analysis is conducted by perturbing the disease free equilibrium state I e (x, t). Near this equilibrium, a small perturbation is introduced

where ε is the amplitude of the small perturbation Page 10 of 28 manuscript submitted to Physica D [Fuentes et al., 2004 , Kefi et al., 2008 , Murray, 2001 . Inserting I(x, t) into equation 10 and noting that at equilibrium ∂ I e (x, t)/∂ t = 0 yields

where h(x) = g p(x) is the convolution of g = cos(ωx) and p(x). To evaluate h(x), recall that the Fourier transforms of g (in k or wave space) iŝ

Upon computing the product ofĝ(k) andp(k) in Fourier space and inverse Fourier transforming

Inserting this outcome from equation 24 into equation 22 yields

A necessary condition for pattern formation is that ϕ > 0 [Murray, 2001 , Fuentes et al., 2004 , Segal et al., 2013 . Since R o > 1 (outbreak), ϕ is always positive. It may be stated that the IDE is de-stabilized by small-amplitude perturbations at any wavenumber (in 1-D), which is a necessary (but not sufficient) condition for pattern formation.

The simulation results are presented to address the two main questions: geographical spreader fraction φ on spatial patterns or any switching between them in time. Other disease-related parameters are also varied within certain ranges including R o =2.2-8.0 reported for COVID-19 [Katul et al., 2020 , Kucharski et al., 2020 , and the fractional mobility of the infected population (φ=0.01 and 0.1). The latter is needed to assess whether coherent pattern shapes in infections are driven by finite bounds of the spatial domain. Time evolution of spatiallyaveraged quantities are also required to highlight a number of features about the domain-averaged SIR system. The domain-averaged quantity for an arbitrary variable Ψ(x, y, t) is indicated by

For notational simplicity, the time dependence is dropped and 〈Ψ〉 is assumed to represent the time-dependent quantity 〈Ψ〉(t).

Before presenting the pattern formation results for multiple attack points as initial conditions, the c ma x are first discussed for different kernel properties and R o with idealized initial conditions. gradients in I to ascertain that the results derived from a spatially integrated measure such as A s is comparable to standard gradient based on estimates of fronts.

In Figure 2 , the expansion area calculated with A s (t) captures the front movement indicated by the comparison of white and red circles lending confidence to the use of A s (t) for c max detection.

In Figure 3 , the numerical model results show that c max is indeed proportional to (1 − φ)(R o − 1) when setting σ constant and to σ(1 − φ) when setting R o to a constant, which agrees with the 1-D analysis. More significant, the c ma x values do not vary much even when φ is increased from 0.01 to 0.10, implying that the fraction of geographical spreaders in the infectious population has a minor impact on the magnitude of the maximal spread velocity attained. However, Figure   2 also demonstrates that φ dictates the time at which c max is reached. Thus, φ is pertinent for identifying 'early-warning' signals about impending outbreaks. action assumption representing interactions between I, S, and R within the box. Therefore, a constraint on the kernel standard deviation σ/d r > 2 is imposed in all later calculations. This constraint is above and beyond the prior constraint R a /σ > 2 to ensure a near-Gaussian shape.

When 3 < σ/d r < 30, a linear relation between c max and σ is confirmed by the simulations as conjectured from equation 19. Hence, we imposed R o > 2, σ/d r > 2, and R a /σ > 2 in later numerical simulations where the domain is initially infected at multiple locations and pattern formation is tracked.

In all the following simulations, 5% of the communities delineated by boxes in the domain are randomly selected and simultaneously infected with 10 individuals as initial conditions. Based 

where m(= N /2) is the number of Fourier modes used in the computation of the Fourier energy spectrum, E * (k) is the squared amplitude at the scalar wavenumber k, and the normalization is the Shannon entropy for a white-noise spectrum determined at m wavenumbers. When En s = 1, the spatial energy spectrum of I(x, y, t) is white and no patterning is expected (at t=0) whereas En s = 0 implies that the entire spatial variability is explained by a single Fourier mode (periodic function). Thus, En s (t) is used here as a scalar measure whose time evolution determines how energy becomes concentrated among Fourier modes starting from an initial white-noise spectrum.

It is to be noted that the 1-D linearized analysis predicts that the IDE is unstable at all modes.

However, such instability does not guarantee pattern formations. The calculations of En s (t) are complemented by the constraint that energetic modes bounded in the interval [5d r, L/2] may be used in identifying coherent spatial patterns if this pattern arises (i.e. low En s ). This constraint is to ensure that any deterministic spatial patterns do not have wavelengths too close to the domain size or grid resolution. The temporal variation of the Shannon entropy is also featured in Figure   4 for illustration. Figure 4 shows that En s (t) has distinct features at four times: (i) initial state (t o ), (ii) initial exponential rise in infections (t i ), (iii) peak infections (t p ) and (iv) the final quasi-steady state stage where infections are slowly approaching their zero-equilibrium value (t s ). When the spatial pattern of I(x, y, 0) is random in space, En s (t) is close to unity as expected at t o . The structure at t i is transient due to the persistent impact from the random initial condition, while the emerging organized structure at t p can show different spatial patterns characterized by different En s (t)

values. These differences frame much of the discussion about pattern formation.

Before delving into those differences, several points can be made about the dynamics of 〈I〉 first.

First, it is clear that increases in σ/d r has minor (i.e. graphically not discern-able) impact on the temporal patterns of 〈I〉 (Figure 4) . However, exploration of the phase space of d 〈I〉 /d t versus 〈I〉 reveals a more nuanced perspective in Figure 5 .

When σ/d r increases, there is a clear shift in the phase-space until saturation at a large σ/d r is attained. This finding hints that there may be a critical σ for which disease spread becomes quasi-independent of σ when measured by a spatially aggregated quantity such as d 〈I〉 /d t. For early-times SIR, all three σ/d r collapse on an approximate line with a (R o − 1) slope as discussed elsewhere [Katul et al., 2020] . However, as t p is gradually approached, d 〈I〉 /d t begins to decline from its maximum value, and the effects of σ/d r become more evident. During the expansion phase, higher σ leads to higher d 〈I〉 /d t at a given 〈I〉. Moreover, at t p where d 〈I〉 /d t = 0, the maximum 〈I〉 that can be supported mildly increases with increasing σ but saturates to a (x, y, t) at key time instances t o , t i , t p , and t s with different kernel parameters, which correspond to the patterns shown in Figure  4 . The second row shows the corresponding spectral density E * (k), while the third row shows the energy content R a kE * (k). The wavenumbers are normalized by R a , the kernel bound. As evidenced by the third row, the so-called pre-multiplied spectral representation kE * (k) better delineates key scale transitions associated with changes in the spectrum.

〈I〉 /N o = 0.45. It is to be noted that increasing σ/d r by a factor of 8 only leads to an increase in the maximum 〈I〉 /N o at t p by some 12%.

Moving beyond the dynamics of 〈I〉, and depending on the kernel properties (Ra and σ), the spatial patterns of I(x, y, t p ) vary from disordered to organized (defined below). This spatial pattern shift at t p does not have an appreciable impact on 〈I〉 as earlier noted. The analysis is now focused on the pattern formation of I(x, y, t p ), namely, the variation of pattern topology exemplified in panels (h), (g), and (i). Corresponding to Figure 4 , the following descriptors are used on spatial patterns for clarity. These descriptors are aimed at qualitatively identifying key spatial patterning features in physical and spectral spaces ( Figure 6 ). In all calculations, the dimensionless energy E * (k) at dimensionless time γt here is determined using a 2-D FFT of I(x, y, t), computing the squared amplitudes (in k x and k y directions), normalizing by the overall variance, and setting the energy content as representative of wavenumber k = k 2 x + k 2 y , k = 1/r o is the wavelength and d r < r o < L is the spatial scale. Figure 6 shows the initial condition (at t o ) that then set a white noise spectral pattern identical for all simulations. Likewise, the steady-state patterns (at t s ) are similar for all simulations. At t i , when En s begins its fastest decline in time (i.e. energy is being concentrated across scales), two possible spatial patterns emerge: diffusive (sensitive to domain-size) and disordered. In both cases, energy content begins to concentrate at low wavenumbers. At t p , the pattern forming structures are either disordered or coherent (both not sensitive to domain size). For both pattern types, disordered and coherent patterns have dominant energetic amplitudes at an intermediate scale that is smaller than the domain size but much larger than the grid resolution. The kernel properties that decide on which pattern persists at t p are to be investigated and discussed.

Deterministic spatial patterns occur when energetic modes in the Fourier spectra are restricted to a narrow band of wavenumbers. They are stationary when the aforementioned narrow range of energetic modes do not vary in time. Such energetic modes are determined here from the maxima of the squared Fourier amplitudes during the course of disease spread in dimensionless time (γt). For an arbitrary kernel cut off of the 1-D system, as shown elsewhere [Fuentes et al., 2003] , pattern separation occurs at the scale σ/L = 2/9 for a 1-D system. A different pattern separation result is also observed in Figure 7 for the 2-D SIR system with R o > 2 at t p .

As shown in Figure 7 , kernel size can influence the type of spatial pattern setting at t p . When R a /d r < 5, disordered spatial patterning always forms initially (Phase I) for a wide range of R a /σ.

For R a /σ > 3 values, the separation between the disordered (Phase I) and coherent (Phase II) patterns follows an approximate linear relation with σ/d r ≈ 5. Hence, the transition between disordered and coherent pattern at t p depends on the shape (i.e. σ/d r, normalized standard deviation) of the kernel for such kernel-based IDE [Fuentes et al., 2003 ].

In Figure 2 , all infected individuals are placed at the center of the domain thus allowing a single wavefront to develop, expand, then contract at long times. Here, random initial conditions in space are used to explore whether the 1-D scaling analysis holds even when the expanding waves interfere. The expansion speeds are examined in Figure 8 . Figure 8 shows that when random initial conditions are imposed spatially, the magnitudes of the expansion speeds increase with increasing R o and σ, as expected from the 1-D analysis. The linear relation still holds, which indicates that even when several traveling waves form and later interfere with each other, the overall expansion speed remains linear in R o . The c max also occurs earlier when compared to a single infection initial condition at the center of the domain (as expected). Moreover, φ has a small impact on the magnitude of the maximum velocity but acts as a retardation factor to the attainment of c ma x as before. Unsurprisingly, the time to peak infection also decreases accordingly with increasing c ma x and increasing φ. The coefficient c k is higher than those featured in Figure   2 , indicating dependency on initial conditions. Figure 9 shows that the transition between the disordered pattern and the coherent structure remains persistent at σ/d r ≈ 5, regardless of the domain size enlargement from N=512

to 1024 and 2048.

A kernel-based model is used to investigate the impact of geographical spreaders (i.e., highly mobile infected individuals) across identical communities -each governed by identical SIR dynamics.

Hence, mass action between susceptible and infectious individuals is still assumed to occur locally within communities whereas mobility across communities is restricted to a small fraction of infectious individuals per unity time (φ). This representation leads to a kernel-based SIR, where the infectious population is governed by an integro-differential equation (IDE). The IDE allows for the description of disease spread in continuous instead of discrete time and in space when the size of the communities is much smaller than the domain size. Moreover, the proposed IDE accommodates non-local spatial processes when compared to diffusive transfer. The resulting IDE is shown to be unstable to small perturbations at any scale around the equilibrium state for early times, which is a necessary (but not sufficient) condition for pattern formation in space.

Scaling analysis (in 1-D) within an unbounded domain and numerical simulations (in 2-D) for a bounded domain with periodic boundary conditions indicate that the maximum spreading speed of the disease is determined by conventional epidemiological parameters γ and R o , and on the spatial mobility kernel properties of geographical spreaders, including its standard deviation σ and finite-size R a . The small φ does not alter the maximal spreading speed, but can retard the time needed to attain this maximum. Numerical simulations confirm that when the standard deviation of the spreading kernel is larger than a critical value, spatial patterns can change from disordered to coherent structures at the time of peak infections in the domain. This change resembles phase-transitions observed in other IDE systems. Future work will consider (i) the role of the kernel shape (exponential, Laplace, and Wald) itself, (ii) the initial spatial distribution of the susceptible population within the lattice. The latter is known to be reasonably represented by a multi-fractal process.

The results reported here may be heuristic in guiding the early prevention of the epidemics that can be described by the SIR model. Reducing the kernel size, equivalent to constraining the movement of the geographical spreaders, is beneficial for reducing the maximal epidemic spreading also changes, confirming that the pseudo-diffusive pattern is sensitive to the domain size. This is consistent with the spectral features observed in Figure A1 . However, the more pertinent finding is that for both N = 512 and N = 1024, the pattern morphology (i.e. pseudo-diffusive classification) still holds despite domain size doubling.

Natural human mobility patterns and spatial spread of infectious diseases

Spatial pattern switching enables cyclic evolution in spatial epidemics

The hidden geometry of complex, network-driven contagion phenomena

Size and timescale of epidemics in the SIR framework

Invasion threshold in heterogeneous metapopulation networks

Self-organization analysis for a nonlocal convective Fisher equation

Daniel Bernoulli's epidemiological model revisited

Spatial pattern formation facilitates eradication of infectious diseases

The wave of advance of advantageous genes

Nonlocal interaction effects on pattern formation in population dynamics

Analytical considerations in the study of spatial patterns arising from nonlocal interaction effects

Dimensions of superspreading

A kernel-modulated sir model for COVID-19 contagious spread from county to continent

Community structure in social and biological networks

Effects of network structure, competition and memory time on social spreading phenomena

Mathematical models of infectious disease transmission

Spatial dynamics of an epidemic model with nonlocal infection

The mathematics of infectious diseases

Cyclic epidemics and extreme outbreaks induced by hydro-climatic variability and memory

Analytically solvable model of spreading dynamics with non-Poissonian processes

Global convergence of COVID-19 basic reproduction number and estimation from early-time sir dynamics

Modeling Infectious Diseases in Humans and Animals

Mathematical modelling of infectious diseases

Seasonally forced disease dynamics explored as switching between attractors

Vegetation pattern shift as a result of rising atmospheric CO2 in arid ecosystems

Proceedings of the Royal Society of London. Series A-Containing papers of a mathematical and physical character

Contact-based model for epidemic spreading on temporal networks

Dispersal data and the spread of invading organisms

Early dynamics of transmission and control of COVID-19: a mathematical modelling study. The Lancet Infectious Diseases

Infectious diseases as socio-spatial processes: The COVID-19 outbreak in Germany

Analytic models for SIR disease spread on random spatial networks

On the origin of Tiger Bush

From patterns to function in living systems: Dryland ecosystems as a case study

Mathematical Biology II: Spatial Models and Biomedical Applications

Mathematical Biology: I. An Introduction

Epidemic processes in complex networks

Large-scale spatial-transmission models of infectious disease

Five challenges for spatial epidemic models

Dynamics of interacting diseases

Pattern formation in a model of competing populations with nonlocal interactions

Periodic travelling waves in cyclic populations: field studies and reaction-diffusion models

Super-spreaders and the rate of transmission of the SARS virus

Spreading processes in multiplex metapopulations containing different mobility networks

Pattern transitions in spatial epidemics: Mechanisms and emergent properties

Plant propagation fronts and wind dispersal: an analytical model to upscale from seconds to decades using superstatistics

Secondary seed dispersal and its role in landscape organization

Role of biomass spread in vegetation pattern formation within arid ecosystems

Spatial organization of vegetation arising from non-local excitation with local inhibition in tropical rainforests

Traveling waves for a diffusive seir epidemic model with non-local reaction

Explicit formulae for the peak time of an epidemic from the sir model

Analytical computation of the epidemic threshold on temporal networks

Reaction-diffusion waves in biology

As shown in Figure 4 , the pseudo-diffusive pattern appears at t i when d〈I〉/d t is maximum. This pattern is (i) transient since the initial condition still dominates the entropy plot and (ii) sensitive to the domain size. The appendix here shows that the transition from disordered to pseudodiffusive is also governed by a phase-transition depending on the kernel properties. The spectral features of this pattern are shown in Figure A1 . The pattern partition map between disordered and pseudo-diffusive is now shown in Figure A2 .The spectral structure of the pseudo-diffusive pattern shown in Figure A1 is similar to the coherent structure in Figure 4 . However, the energetic mode of the pseudo-diffusive pattern of R a kE * (k) emerges in the first few wavenumbers. This suggests that the pattern is dictated by the largest length scale in the domain. To test whether its morphology (i.e. labyrinth-like shape) is influenced by the domain size, further simulations with N = 1024 are conducted and shown in Figure A3 .It is shown in Figure A3 that when the domain size changes, the maximal energetic modes at t i