# The Role of Geographic Spreaders in Infectious Pattern Formation and Front
Propagation Speeds

Shuolin Li  Department of Civil and Environmental Engineering, Duke
University, Durham, NC, USA  Craig Henriquez  Department of Biomedical
Engineering, Duke University, Durham, NC, USA  Gabriel Katul  Nicholas School
of the Environment and Earth Science, Duke University, Durham, NC, USA
Department of Civil and Environmental Engineering, Duke University, Durham,
NC, USA

###### Abstract

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 () 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 , scaling
analysis in 1-dimension and simulation results in 2-dimensions suggest that
$c_{max}\sim(1-\phi)\gamma(R_{o}-1)\sigma$, where is the inverse of the
infectious duration, and 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 though reductions in act as
retardation factors to the attainment of $c_{max}$. The 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.

## Keywords

Geographic Spreaders, Integro-Differential Equations, Spatial Pattern,
Susceptible-Infectious-Recovered model

##  Introduction

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^{\prime},y-y^{\prime})$
of infectious individuals to travel from point defined by Cartesian
coordinates $(x,y)$ to $(x^{\prime},y^{\prime})$. 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=\sqrt{(x-x^{\prime})^{2}+(y-y^{\prime})^{2}}$). Likewise, spreading kernels
can be non-isotropic with certain preferential pathways breaking the symmetry
of the movement. The $p(x-x^{\prime},y-y^{\prime})$ can be made time-dependent
i.e. $p(x-x^{\prime},y-y^{\prime},t-t_{o})$ 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,
Thompson and Katul, 2009, Thompson et al., 2009, 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.

##  Theory

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.

###  Definitions and Nomenclature

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],

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$
$\displaystyle-\left(\beta\frac{S}{N_{o}}\right)I,$ (1)
$\displaystyle\frac{dI}{dt}$ $\displaystyle=$
$\displaystyle\left(\beta\frac{S}{N_{o}}\right)I-\gamma I,$
$\displaystyle\frac{dR}{dt}$ $\displaystyle=$ $\displaystyle\gamma I,$

where $t$ is time, is an average number of contacts per individual per unit
time, and -1 is the average time period required for an infected individual to
get recovered. The $\beta S/N_{o}$ describes the force of infection [Dietz and
Heesterbeek, 2002, Cadoni and Gaeta, 2020] and coefficients and are externally
supplied constants that depend on the pathogen, environmental conditions, and
the behaviour of the infected population (discussed later). This system
requires as initial conditions $I(0)$ to be finite, $S(0)=N_{o}-I(0)$ and
$R(0)=0$ at the initial phase of the outbreak. By summing the three equations
describing the SIR system yields $d(S+I+R)/dt=0$, and then upon integration,
result in $S(t)+I(t)+R(t)=N_{o}$ being constant at all times. The constant
$N_{o}$ 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 individual-to-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 that
$dI/dt$ will be positive (outbreak) or negative (epidemic contained) depending
on the sign of $(\beta S/N_{o}-\gamma)$. Defining the basic reproduction
number $R_{o}=\beta/\gamma$ and assuming $S\approx N_{o}$ at early times (i.e.
$I(0)/N_{o}\ll 1$) results in the necessary condition for an outbreak:
$R_{o}>1$. We track the sensitivity of front speeds and spatial pattern
formation of infections on $R_{o}$ while setting -1 constant to 14 days.
Typical values of $R_{o}$ have been reported for a number of diseases and
range between 2 to 8 (including COVID-19 [Katul et al., 2020]).

Spatially, and $R_{o}$ are assumed constant for each box during an outbreak
though effectiveness of certain treatment can reduce the average time between
onset of symptoms and recovery. Likewise, a constant $R_{o}$ and ignores other
factors such as stay-at home orders, masking, social distancing, hand washing,
etc….In general, these two parameters may also vary with other environmental
or hydro-climatic or seasonal conditions [Keeling et al., 2001, Hooshyar et
al., 2020] not covered and are outside the scope of the work here.

###  A Spatial Model with Non-local Mobility of Infectious Individuals

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 $dx$ by $dy$ 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 $S(0)+I(0)+R(0)=N_{o}$.

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

$\frac{\partial I}{\partial
t}=\left(\beta\frac{S}{N_{o}}-\gamma\right)I-D_{o}\left(\frac{\partial^{2}I}{\partial
x^{2}}+\frac{\partial^{2}I}{\partial y^{2}}\right),$ (2)

where $D_{o}$ is a diffusion coefficient that must be a priori known. This
representation is continuous in space and time, meaning $dx$ and $dy$ 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\approx N_{o}$)
[Fisher, 1937, Volpert and Petrovskii, 2009, Murray, 2001, Belik et al., 2011]

$c_{max}=2\sqrt{D_{o}\gamma(R_{o}-1)},$ (3)

when $R_{o}>1$. This relation between $c_{max}$ and
$\sqrt{D_{o}\gamma(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

$\frac{\partial I}{\partial
t}=\left(\beta\frac{S}{N_{o}}-\gamma\right)\left[(1-\phi)I+\phi
C_{u}(x,y)\right],$ (4)

with the convolution function $C_{u}(x,y)$ defined as

$C_{u}(x,y)=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}I(x^{\prime},y^{\prime})p(x-x^{\prime},y-y^{\prime})dx^{\prime}dy^{\prime}$
(5)

where is the fraction of infectious individuals (geographical spreaders) that
are mobile across boxes within a given time step $dt\ll{}^{-1}$ and
$p(x^{\prime},y^{\prime})$ is a spatial spread kernel satisfying the
normalizing condition

$\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}p(x^{\prime},y^{\prime})dx^{\prime}dy^{\prime}=1.$
(6)

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

Figure 1: A conceptual diagram of the SIR model and its extension in space. In
each box (dashed square), there are three compartments ($S$, $I$, and $R$)
describing the composition of its population. Each box can be understood as an
a community abiding by SIR dynamics when the infectious individual(s) arrive
in this box. To describe the interaction between any two communities, two
different transport models are presented. The left shows spatial spread across
communities (or boxes) represented by a diffusion scheme with a diffusion
coefficient $D_{o}$, while the right plot shows the spread between communities
using the kernel-based model accounting for the non-local effects (i.e.
infectious individuals can access all points in the spatial domain). In the
diffusion case, only spatial gradients in $I$ drive spatial mobility whereas
in the kernel-based case, few infectious individuals can travel far from their
parent community and access susceptibles in other communities or boxes. These
individuals are labelled here as geographical spreaders and constitute a small
fraction () of the infected population.

When $\phi=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 $\phi=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 $\phi(\ll 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^{\prime},y^{\prime})$ [Jo et al., 2014, Gleeson et al., 2016, Keeling and
Rohani, 2011]. For simplicity and to recover $c_{max}$ predicted by diffusion,
we selected

$p(r)=\frac{{}_{N}}{\sqrt{2\pi}\sigma}\exp\left[-\frac{1}{2}\left(\frac{r}{\sigma}\right)^{2}\right],$
(7)

where $r=\sqrt{(x-x_{o})^{2}+(y-y_{o})^{2}}$ is the distance, $(x,y)$ and
$(x_{o},y_{o})$ are two points in the domain 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 finite support $R_{a}$ set by $dr<R_{a}<L$ where
$dr^{2}=(dx)^{2}+(dy)^{2}$, N is determined so that [Fuentes et al., 2003,
2004, Da Cunha et al., 2009]

$\frac{{}_{N}}{\sqrt{2\pi}\sigma}\int_{-R_{a}/\sqrt{2}}^{R_{a}/\sqrt{2}}\exp\left[-\frac{1}{2}\left(\frac{r^{\prime}}{\sigma}\right)^{2}\right]dr^{\prime}=1.$
(8)

This condition yields the normalizing constant

${}_{N}(R_{a},\sigma)=\left[\rm{erf}\left(\frac{R_{a}}{2\sigma}\right)\right]^{-1}.$
(9)

To maintain a near-Gaussian shape for the spreading kernel, we selected
$R_{a}/(2\sigma)>1$. This IDE approach along with a Gaussian spreading kernel
was recently employed to explain the multi-fractal 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.

###  Basic Properties of the IDE Model

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)\approx S(x,0)$ so that
the IDE can be reduced to

$\frac{\partial I}{\partial
t}\approx\frac{\gamma}{\left(R_{o}-1\right)}\left[(1-\phi)I+\phi
C_{u}(x)\right].$ (10)

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

$C_{u}(x)=\int_{-\infty}^{+\infty}I(x^{\prime})p(x-x^{\prime})dx^{\prime}=I(x^{\prime})\star
p(x^{\prime}).$ (11)

####  Travelling waves and front speed

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-ct$ along with
$dx^{\prime}=dz^{\prime}$ is used to convert the 1-D IDE in equation 10 to the
new coordinate system, given as

$-c\frac{\partial I}{\partial
z}\approx\gamma\left(R_{o}-1\right)\left[(1-\phi)I+\phi C_{u}(z)\right].$ (12)

Assuming the solution is of the form $I(z)=A\exp(\lambda 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 as

$\hat{I}(k)=2A\pi\delta(\lambda i+k),\leavevmode\nobreak\
\hat{p}(k)=\exp\left(-\frac{1}{2}{}^{2}k^{2}\right),$ (13)

where $A$ is a constant, $k$ is the wavenumber, $i^{2}=-1$, $\delta(.)$ is the
Dirac-delta function, and $\hat{s}$ denotes the corresponding Fourier
transform of function $s(z)$. To evaluate $C_{u}(z)$, the product of
$\hat{I}(k)$ and $\hat{p}(k)$ is determined in the Fourier domain followed by
an inverse Fourier transform yielding $C_{u}(z)=I(z)M_{g}(\lambda)$ where

$\displaystyle M_{g}(\lambda)=\exp\left(\frac{1}{2}{}^{2}{}^{2}\right)$ (14)

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

$-c\lambda=\gamma\left(R_{o}-1\right)\left[(1-\phi)+\phi
M_{g}(\lambda)\right].$ (15)

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 $dM_{g}/d\lambda=\lambda{}^{2}M_{g}$ and
differentiating equation 15 with respect to yields

$\displaystyle-c=\gamma\left(R_{o}-1\right)\lambda\left[\phi{}^{2}M_{g}(\lambda)\right],$
(16)

or as

$\displaystyle\phi
M_{g}(\lambda)=-\frac{c}{\lambda}\left[\frac{1}{\gamma\left(R_{o}-1\right){}^{2}}\right].$
(17)

When $R_{o}>1$ (outbreak), equation 16 requires that $c/\lambda<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

$-c{}^{2}-\gamma(R_{o}-1)(1-\phi)\lambda+\frac{c}{{}^{2}}=0,$ (18)

whose solution is

$c_{max}=\frac{c_{k}}{1-c_{k}^{2}}\gamma\sigma(1-\phi)(R_{o}-1),$ (19)

where $c_{k}=\lambda\sigma$ 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 $\lambda<0$) and a backward velocity (when
$\lambda>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

$\Delta=\left[\gamma(1-\phi)(R_{o}-1)\right]^{2}+\left[\frac{2c_{max}}{\sigma}\right]^{2}>0,$
(20)

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

When $\phi\ll 1$, as is the case here, $c_{max}$ appears to be insensitive to
though $1-\phi$ does act as a retardation factor to the attainment of
$c_{max}$ (discussed later). Unlike the diffusion model, the kernel-based
$c_{max}$ varies linearly with $\gamma(R_{o}-1)$ instead of
$\sqrt{\gamma(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}$,

$\frac{D_{o}}{\gamma{}^{2}}=\frac{c_{k}^{2}}{4(1-c_{k}^{2})^{2}}{(1-\phi)^{2}(R_{o}-1)}.$
(21)

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.

####  Necessary Conditions for Pattern Formation

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 so that
$I(x,t)=I_{e}(x,t)+\epsilon\cos(\omega x)\exp(\varphi t)$, where is the
amplitude of the small perturbation [Fuentes et al., 2004, Kefi et al., 2008,
Murray, 2001]. Inserting $I(x,t)$ into equation 10 and noting that at
equilibrium $\partial I_{e}(x,t)/\partial t=0$ yields

$\varphi=\left(R_{o}-1\right)\left[(1-\phi)+\phi\frac{h(x)}{\cos(\omega
x)}\right],$ (22)

where $h(x)=g\star p(x)$ is the convolution of $g=\cos(\omega x)$ and $p(x)$.
To evaluate $h(x)$, recall that the Fourier transforms of $g$ (in $k$ or wave
space) is

$\displaystyle\hat{g}(k)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\pi}{2}}\left[\delta(\omega-k)+\delta(\omega+k)\right]$
(23)

Upon computing the product of $\hat{g}(k)$ and $\hat{p}(k)$ in Fourier space
and inverse Fourier transforming yields

$h(x)=\cos(\omega x)\exp\left(-\frac{1}{2}{}^{2}{}^{2}\right).$ (24)

Inserting this outcome from equation 24 into equation 22 yields

$\varphi=\left(R_{o}-1\right)\left[(1-\phi)+\phi\exp\left(-\frac{1}{2}{}^{2}{}^{2}\right)\right].$
(25)

A necessary condition for pattern formation is that $\varphi>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.

##  Simulation Results

The simulation results are presented to address the two main questions: (i) to
what degree $c_{max}$ derived from a linearized analysis in 1-D captures the
key state variables describing propagation fronts in 2-D for an isolated
outbreak in the domain center (set at $x_{o},y_{o}$). (ii) What spatial
patterns are expected when the domain becomes finite and is attacked at
multiple communities leading to interfering wavefronts? and (iii) what is the
role of the kernel properties ($R_{a}$, ) and geographical spreader fraction
on spatial patterns or any switching between them in time.

###  Numerical Simulations and their Presentation

All calculations are conducted on a square domain with size $L\times L$
comprising of $N\times N$ square boxes (or lattices) with $dx=dy=1$ of
arbitrary spatial units. For $N$, three different values ($N$=512, 1024 and
2048) are used and compared to ensure that the patterns are not sensitive to
the domain size. The boundary conditions imposed on equation 10 are periodic.
In such calculations, time $t$ is made dimensionless using . Time-integration
employs forward differences with dimensionless time increments $\gamma
dt=0.1$. Spatial integration is conducted using two-dimensional Fast Fourier
Transform (FFT) convolutions between $I(x^{\prime},y^{\prime},t)$ and
$p(x-x^{\prime},y-y^{\prime})$. Variations in $R_{a}/L$ and $\sigma/R_{a}$ are
explored with certain constraints set by (i) the smallest size
$dr=\sqrt{dx^{2}+dy^{2}}=\sqrt{2}$ (arbitrary spatial units) and (ii) the
largest size determined by the number of grid nodes
$N^{2}=(L/dx)\times(L/dy)=2(L/dr)^{2}$.

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 spatially-averaged quantities are also required to highlight a number of
features about the domain-averaged SIR system. The domain-averaged quantity
for an arbitrary variable $\Psi(x,y,t)$ is indicated by

$\langle\Psi\rangle(t)=\frac{1}{L^{2}}\int_{-L/2}^{L/2}\int_{-L/2}^{L/2}\Psi(x,y,t)dxdy.$
(26)

For notational simplicity, the time dependence is dropped and
$\langle\Psi\rangle$ is assumed to represent the time- dependent quantity
$\langle\Psi\rangle(t)$.

###  Maximal Expansion and Contraction Speeds

Before presenting the pattern formation results for multiple attack points as
initial conditions, the $c_{max}$ are first discussed for different kernel
properties and $R_{o}$ with idealized initial conditions. These conditions
feature a concentrated initial infectious event of 50 individuals positioned
at the center of the domain in one box. The spreading area $A_{s}(t)$ is used
throughout to delineate overall spreading within the entire domain. It is
computed from the fractional area where $I(x,y,t)/N_{o}>0.1\%$. While this
threshold appears arbitrary, it is selected as a compromise between the need
for front detection and a sufficient number of infections to initiate early
dynamics in SIR locally within a box. This allows estimates of a spatially-
integrated front speed $c(t)=d\sqrt{A_{s}(t)}/dt$ (expansion when
$dA_{s}/dt>0$ and contraction when $dA_{s}/dt<0$) when multiple infectious
points produce interfering front waves in later analyses. A simulation result
featuring the idealized initial condition and the computed area expansion and
contraction as well as front speed are shown in Figures 2 and 3. This front
detection is compared with another approach based on maximal spatial 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.

Figure 2: Representative contours of the spreading area $A_{s}(t)$ and the
maximal expansion and contraction speeds for an initial infection at the
center of the domain ($I(256,256,0)=50$). The blue and red circles outline the
epidemiological expansion area directly captured by the maximal spatial
gradient and by $A_{s}(t)$, respectively. In all the plots $\phi=0.01$,
$R_{o}=4.5$, $R_{a}/dr=64$, $\sigma/dr=16$, and $\gamma t$=2, 3, 7 and 18.6
respectively.

Figure 3: Variation of the fractional spreading area $A_{s}(t)$ versus
dimensionless time and the maximal expansion speeds $c_{max}$ for an initial
infection at the center of the domain. The blue and red dashed lines represent
the analyzed velocity using Equation 19. In all plots, $R_{a}/dr=64$,
$\sigma/dr=16$ in the top panel and $R_{o}$=4.5 in the bottom panel.

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-\phi)(R_{o}-1)$ when setting constant and to $\sigma(1-\phi)$ when setting
$R_{o}$ to a constant, which agrees with the 1-D analysis. More significant,
the $c_{max}$ 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.

Figure 3 suggests that for the values of and explored here, a front cannot
form and persist when $R_{o}<3$. This implies that $R_{o}-1>1$ in all model
calculations to be conducted later. Interestingly, the 1-D model reasonably
captures the 2D simulation results. However, a saturation in front speeds
begins to emerge at very large $R_{o}>7$ irrespective of , which cannot be
predicted from the 1-D linearized analysis. This saturation in $c_{max}$ with
increasing $R_{o}$ arises because of a rapid areal expansion resulting in a
jump between no infections to an all-domain infection in a short time interval
difficult to fully resolve with finite $dA_{s}/dt$. When the initial condition
is fixed, the fitted constant $c_{k}$ to modeled $c_{max}$ does not change
suggesting it is not sensitive to the kernel properties or . When
$\sigma/dr<3$, kernel shapes are too narrow where much of the $I$ spreading
occurs within an individual box instead of across boxes. This mode of
spreading violates the mass action assumption representing interactions
between $I$, $S$, and $R$ within the box. Therefore, a constraint on the
kernel standard deviation $\sigma/dr>2$ is imposed in all later calculations.
This constraint is above and beyond the prior constraint $R_{a}/\sigma>2$ to
ensure a near-Gaussian shape. When $3<\sigma/dr<30$, a linear relation between
$c_{max}$ and is confirmed by the simulations as conjectured from equation 19.
Hence, we imposed $R_{o}>2$, $\sigma/dr>2$, and $R_{a}/\sigma>2$ in later
numerical simulations where the domain is initially infected at multiple
locations and pattern formation is tracked.

###  Formation of Infectious Spatial Patterns

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 on the IDE solution for $\langle
I(x,y,t)\rangle$ versus $\gamma t$, the SIR model exhibits transient spatial
patterns depending on choices made on $R_{a}$ and . Typical transient patterns
of three structures are featured and analyzed in Figures 4 and 6. To
interrogate a broad range of parameter space reflecting disease transmission
($R_{o}$) and kernel shape ( and $R_{a}$) on pattern formation, a
dimensionless spectral Shannon entropy is proposed, given by

$En_{s}(t)=\frac{\sum_{i=1}^{i=m}E^{*}(k)\log[E^{*}(k)]}{\log(m)},$ (27)

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 $[5dr,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: The variations of $\langle S\rangle$, $\langle I\rangle$, and
$\langle R\rangle$ normalized by the total initial population ($N_{o}$ as a
function of dimensionless time $\gamma t$). The solution starts at $t_{0}=0$
with spatially random infectious locations. The $t_{i}$ is the time when
$d\langle I\rangle/dt$ reach its maximum, $t_{p}$ is the point when $\langle
I\rangle$ becomes its maximum and $t_{s}$ is the quasi-steady state where
infections are slowly approaching their zero equilibrium. The representative
instances under consideration are labeled with red arrows along with the
typical patterns expected. Plots (d) to (f) show the variation of entropy
$En_{s}$ and spreading area $A_{s}$ over time, where the shapes of the kernels
with different are also sketched to illustrate the kernel spread within the
domain. Plots (g) to (l) show the variation of the spatial structure of
$I(x,y)$ depending on different kernel information ($\sigma/dr$) at the four
characteristic times. In all the plots $R_{a}/dr=128$, $R_{o}=4.5$.

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 $\left\langle I\right\rangle$ first. First, it is clear that
increases in $\sigma/dr$ has minor (i.e. graphically not discern-able) impact
on the temporal patterns of $\left\langle I\right\rangle$ (Figure 4). However,
exploration of the phase space of $d\left\langle I\right\rangle/dt$ versus
$\left\langle I\right\rangle$ reveals a more nuanced perspective in Figure 5.

Figure 5: Phase space of the $N_{o}^{-1}d\left\langle I\right\rangle d(\gamma
t)$ and $\left\langle I\right\rangle N_{o}^{-1}$. The inset is a zoom-in over
the initial period where the black solid line is the expected ($R_{o}-1$)
slope.

When $\sigma/dr$ increases, there is a clear shift in the phase-space until
saturation at a large $\sigma/dr$ 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\left\langle
I\right\rangle/dt$. For early-times SIR, all three $\sigma/dr$ 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\left\langle
I\right\rangle/dt$ begins to decline from its maximum value, and the effects
of $\sigma/dr$ become more evident. During the expansion phase, higher leads
to higher $d\left\langle I\right\rangle/dt$ at a given $\left\langle
I\right\rangle$. Moreover, at $t_{p}$ where $d\left\langle
I\right\rangle/dt=0$, the maximum $\left\langle I\right\rangle$ that can be
supported mildly increases with increasing but saturates to a $\left\langle
I\right\rangle/N_{o}=0.45$. It is to be noted that increasing $\sigma/dr$ by a
factor of $8$ only leads to an increase in the maximum $\left\langle
I\right\rangle/N_{o}$ at $t_{p}$ by some 12$\%$.

Moving beyond the dynamics of $\left\langle I\right\rangle$, 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 $\left\langle I\right\rangle$
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).

Figure 6: The first row shows spatial variations of $I(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) Disordered Structure: In physical space, these patterns are characterized
by sharp interfaces that delineate infected from non-infected regions. In
spectral space, this pattern can be discerned when an energetic mode with
maximal spectral energy content (instead of spectral density) arises. This
mode is detected at an intermediate wavenumber ($k_{e}$) where the spectral
energy content $\sigma kE^{*}(k)$ reaches its maximum and obeys the constraint
$k_{o}\ll k_{e}\ll k_{r}$ where $k_{o}\sim 1/L$ and $k_{r}\sim 1/dr$ are the
limiting wavenumbers in Fourier space dictated by domain size $L$ and spatial
resolution $dr$. Such a disordered pattern is shown in panels (b) and (c) of
Figure 6.

ii) Coherent Structure: In physical space, these patterns emerge after the
infectious population spreads across the entire spatial domain. The infections
are space-filling and no lattice remains free from infections within the
lattice. Thus, the $I(x,y,t)/N_{o}$ in space is primarily dominated by high
values interspersed with infrequent but relatively smaller infectious spots of
low (but finite) $I(x,y,t)/N_{o}$. In spectral space, a cluster of several
contiguous dominant energetic modes with characteristic scales smaller than
$L$ emerge. However, the energy content $\sigma kE^{*}(k)$ of the large
wavenumbers remains high and quasi-random when compared to wavenumbers
commensurate with $k_{o}$. This feature is shown in (d) and (e) of Figure 6.

iii) Pseudo-diffusive Structure: In physical space, the transitions between
infected and non-infected zones are almost obscured corresponding to the panel
(i) in Figure 4. The $I(x,y,t)/N_{o}$ in physical space spans numerous values.
In spectral space, the pattern is associated with $k_{e}$ approaching $k_{o}$,
meaning the spectral energy reaches its maximum at the lowest wavenumber
(space-filling), and intermediate to high wavenumber modes are not energetic.
The domain size clearly impacts $k_{o}$ here (as expected) - but the
qualitative aspect of the pattern remains less affected by the domain size
(i.e. as long as the pattern is almost space-filling irrespective of the
size). Since it is a transient structure and can be impacted by domain size,
it is not the main focus here. For completeness, its spectral properties are
discussed in the Appendix for three domain sizes to illustrate the qualitative
aspects of such pattern type and its sensitivity to domain size.

In all calculations, the dimensionless energy $E^{*}(k)$ at dimensionless time
$\gamma 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=\sqrt{k_{x}^{2}+k_{y}^{2}}$, $k=1/r_{o}$ is the wavelength and
$dr<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.

###  Impact of Kernel Size on Pattern Formation

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
($\gamma 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
$\sigma/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}$.

Figure 7: Spatial map for the pattern transition at $t_{p}$. The occurrence of
different patterns is marked with different patches. The blue line
(coefficient of determination near unity) is the interface separating the
disordered pattern (cyan, phase I) from the coherent pattern (white, phase
II). The pattern is identified as coherent when the energetic mode becomes one
of the last 10 nodes in the spectral energy diagram. The simulation is
conducted with $N=512$, $\phi=0.01$ and $R_{o}=4.5$.

As shown in Figure 7, kernel size can influence the type of spatial pattern
setting at $t_{p}$. When $R_{a}/dr<5$, disordered spatial patterning always
forms initially (Phase I) for a wide range of $R_{a}/\sigma$. For
$R_{a}/\sigma>3$ values, the separation between the disordered (Phase I) and
coherent (Phase II) patterns follows an approximate linear relation with
$\sigma/dr\approx 5$. Hence, the transition between disordered and coherent
pattern at $t_{p}$ depends on the shape (i.e. $\sigma/dr$, normalized standard
deviation) of the kernel for such kernel-based IDE [Fuentes et al., 2003].

##  Discussion

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: Left column: Variation of the spreading area $A_{s}$ in
dimensionless time $\gamma t$ normalized by the domain size $L^{2}$. Middle
column: maximal expansion speed variations with , $R_{o}$, and kernel
spreading property . Right column: Normalized time to peak infections as a
function of , $R_{o}$, and . Randomized initial conditions are used. The other
input model parameters are identical to Figure 2.

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_{max}$ as before. Unsurprisingly, the time to
peak infection also decreases accordingly with increasing $c_{max}$ and
increasing . The coefficient $c_{k}$ is higher than those featured in Figure
2, indicating dependency on initial conditions.

While is small (but finite), the developing infectious area occupies the
entire domain. Hence, infectious individuals acting through may be termed as
geographical spreaders because only their mobility results in disease
expansion that impacts the entire spread in the domain area. Last, variations
in $\sigma/dr$ lead to a linear regime for intermediate values of
$\sigma/dr>5$ again consistent with predictions from equation 19. Around
$\sigma/dr\approx 5$, $c_{max}$ in both Figure 2 and Figure 8 suggest a
departure from linearity. This threshold $\sigma/dr$ is consistent with the
pattern transition results in Figure 7. Hence, while $\langle I\rangle$
dynamics are not sensitive to $\sigma/dr$, $c_{max}$ does sense the spatial
pattern in SIR dynamics. To test the stability of this transition at
$\sigma/dr\approx 5$, systems with different domain sizes are also
investigated. The simulation results with a domain size of $N=$ 1024 and 2048
are shown in Figure 9. It is noted here the randomly initialized conditions
and $R_{o}=4.5$ are the same as Figure 7. Figure 9 shows that the transition
between the disordered pattern and the coherent structure remains persistent
at $\sigma/dr\approx 5$, regardless of the domain size enlargement from N=512
to 1024 and 2048.

Figure 9: Evolution of patterns at $t_{p}$ when the domain size is enlarged.
It is noted that the kernel sizes in the three rows are $R_{a}/dr=256$, $64$
and $256$, respectively. The patterns shown in the first two columns have
clear disordered structures where as the patterns observed in the last two
columns are coherent.

##  Summary and Conclusion

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 speed and retarding
the fast expansion of the epidemics. Similarly, controlling the variance of
the spatial kernel can act to prevent the formation of coherent patterns
(simultaneous breakouts) in the entire domain, which can significantly
leverage the burden of the hospital system across communities.

## Acknowledgments

G. Katul acknowledges support from the U.S. National Science Foundation (NSF-
AGS-1644382, NSF-AGS-2028633, and NSF-IOS-1754893). S. Li acknowledges the
financial support from the Nicholas School of the Environment at Duke
University.

## Appendix

As shown in Figure 4, the pseudo-diffusive pattern appears at $t_{i}$ when
$d\langle I\rangle/dt$ 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
pseudo-diffusive 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.

Figure A1: The spectral density $E^{*}(k)$ and cumulative energy spectra
$R_{a}kE^{*}(k)$ of the pseudo-diffusive pattern at $t_{i}$. The simulation
parameters are identical to plot (i) of Figure 4.

Figure A2: Spatial map for the pattern transition at time $t_{i}$. The
occurrence of different patterns is marked with with different patches. The
blue line (coefficient of determination approximately unity) is the interface
separating the disordered pattern (green) from the pseudo-diffusive pattern
(white). The simulation is conducted with $N=512$, $\phi=0.01$ and
$R_{o}=4.5$.

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.

Figure A3: Comparison of the pseudo-diffusive patterns when the domain size
doubles. Plots (a) and (b) show two different pseudo-diffusive patterns at
$t_{i}$ when $N=512$. The two patterns in plots (c) and (d) correspond to a
disordered (panel c) and pseudo-diffusive (panel d) when the domain size
doubles (N=1024).

It is shown in Figure A3 that when the domain size changes, the maximal
energetic modes at $t_{i}$ 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.

## References

  * Belik et al. [2011] V. Belik, T. Geisel, and D. Brockmann.  Natural human mobility patterns and spatial spread of infectious diseases.  _Physical Review X_ , 1(1):011001, 2011.
  * Boerlijst and Van Ballegooijen [2010] M. Boerlijst and W. Van Ballegooijen.  Spatial pattern switching enables cyclic evolution in spatial epidemics.  _PLoS Computational Biology_ , 6(12):e1001030, 2010.
  * Brockmann and Helbing [2013] D. Brockmann and D. Helbing.  The hidden geometry of complex, network-driven contagion phenomena.  _Science_ , 342(6164):1337–1342, 2013.
  * Cadoni and Gaeta [2020] M. Cadoni and G. Gaeta.  Size and timescale of epidemics in the SIR framework.  _Physica D: Nonlinear Phenomena_ , 411:132626, 2020.
  * Colizza and Vespignani [2007] V. Colizza and A. Vespignani.  Invasion threshold in heterogeneous metapopulation networks.  _Physical Review Letters_ , 99(14):148701, 2007\.
  * Da Cunha et al. [2009] J. Da Cunha, A. Penna, M. Vainstein, R. Morgado, and F. Oliveira.  Self-organization analysis for a nonlocal convective Fisher equation.  _Physics Letters A_ , 373(6):661–667, 2009.
  * Dietz and Heesterbeek [2002] K. Dietz and J. Heesterbeek.  Daniel Bernoulli’s epidemiological model revisited.  _Mathematical Biosciences_ , 180(1-2):1–21, 2002\.
  * Eisinger and Thulke [2008] D. Eisinger and H.-H. Thulke.  Spatial pattern formation facilitates eradication of infectious diseases.  _The Journal of Applied Ecology_ , 45(2):415, 2008\.
  * Fisher [1937] R. Fisher.  The wave of advance of advantageous genes.  _Annals of Eugenics_ , 7(4):355–369, 1937.
  * Fuentes et al. [2003] M. Fuentes, M. Kuperman, and V. Kenkre.  Nonlocal interaction effects on pattern formation in population dynamics.  _Physical Review Letters_ , 91(15):158104, 2003\.
  * Fuentes et al. [2004] M. Fuentes, M. Kuperman, and V. Kenkre.  Analytical considerations in the study of spatial patterns arising from nonlocal interaction effects.  _The Journal of Physical Chemistry B_ , 108(29):10505–10508, 2004.
  * Galvani and May [2005] A. P. Galvani and R. M. May.  Dimensions of superspreading.  _Nature_ , 438(7066):293–295, 2005.
  * Geng et al. [2021] X. Geng, G. Katul, F. Gerges, E. Bou-Zeid, H. Nassif, and M. Boufadel.  A kernel-modulated sir model for COVID-19 contagious spread from county to continent.  _Proceedings of the National Academy of Sciences_ , 118(21), 2021.  10.1073/pnas.2023321118.
  * Girvan and Newman [2002] M. Girvan and M. E. Newman.  Community structure in social and biological networks.  _Proceedings of the National Academy of Sciences_ , 99(12):7821–7826, 2002.
  * Gleeson et al. [2016] J. P. Gleeson, K. P. O’Sullivan, R. A. Baños, and Y. Moreno.  Effects of network structure, competition and memory time on social spreading phenomena.  _Physical Review X_ , 6(2):021019, 2016.
  * Grassly and Fraser [2008] N. C. Grassly and C. Fraser.  Mathematical models of infectious disease transmission.  _Nature Reviews Microbiology_ , 6(6):477–487, 2008.
  * Guo et al. [2020] Z.-G. Guo, G.-Q. Sun, Z. Wang, Z. Jin, L. Li, and C. Li.  Spatial dynamics of an epidemic model with nonlocal infection.  _Applied Mathematics and Computation_ , 377:125158, 2020\.
  * Hethcote [2000] H. Hethcote.  The mathematics of infectious diseases.  _SIAM Review_ , 42(4):599–653, 2000.
  * Hooshyar et al. [2020] M. Hooshyar, C. Wagner, R. Baker, C. Metcalf, B. Grenfell, and A. Porporato.  Cyclic epidemics and extreme outbreaks induced by hydro-climatic variability and memory.  _Journal of the Royal Society Interface_ , 17(171):20200521, 2020.
  * Jo et al. [2014] H.-H. Jo, J. I. Perotti, K. Kaski, and J. Kertész.  Analytically solvable model of spreading dynamics with non-Poissonian processes.  _Physical Review X_ , 4(1):011041, 2014.
  * Katul et al. [2020] G. Katul, A. Mrad, S. Bonetti, G. Manoli, and A. Parolari.  Global convergence of COVID-19 basic reproduction number and estimation from early-time sir dynamics.  _PLOS One_ , 15:e0239800, 2020.  https://doi.org/10.1371/journal.pone.0239800.
  * Keeling and Rohani [2011] M. Keeling and P. Rohani.  _Modeling Infectious Diseases in Humans and Animals_.  Princeton University Press, 2011.
  * Keeling et al. [2009] M. Keeling, L. Danon, et al.  Mathematical modelling of infectious diseases.  _British Medical Bulletin_ , 92(1):33–42, 2009\.
  * Keeling et al. [2001] M. J. Keeling, P. Rohani, and B. T. Grenfell.  Seasonally forced disease dynamics explored as switching between attractors.  _Physica D: Nonlinear Phenomena_ , 148(3-4):317–335, 2001.
  * Kefi et al. [2008] S. Kefi, M. Rietkerk, and G. G. Katul.  Vegetation pattern shift as a result of rising atmospheric CO2 in arid ecosystems.  _Theoretical Population Biology_ , 74(4):332–344, 2008.
  * Kermack and McKendrick [1927] W. Kermack and A. McKendrick.  A contribution to the mathematical theory of epidemics.  _Proceedings of the Royal Society of London. Series A-Containing papers of a mathematical and physical character_ , 115(772):700–721, 1927.
  * Koher et al. [2019] A. Koher, H. H. Lentz, J. P. Gleeson, and P. Hövel.  Contact-based model for epidemic spreading on temporal networks.  _Physical Review X_ , 9(3):031017, 2019.
  * Kot et al. [1996] M. Kot, M. Lewis, and P. van den Driessche.  Dispersal data and the spread of invading organisms.  _Ecology_ , 77(7):2027–2042, 1996.
  * Kucharski et al. [2020] A. J. Kucharski, T. W. Russell, C. Diamond, Y. Liu, J. Edmunds, S. Funk, R. M. Eggo, F. Sun, M. Jit, J. D. Munday, et al.  Early dynamics of transmission and control of COVID-19: a mathematical modelling study.  _The Lancet Infectious Diseases_ , 20(5):553–558, 2020.
  * Kuebart and Stabler [2020] A. Kuebart and M. Stabler.  Infectious diseases as socio-spatial processes: The COVID-19 outbreak in Germany.  _Tijdschrift voor Economische en Sociale Geografie_ , 111(3):482–496, 2020.
  * Lang et al. [2018] J. C. Lang, H. De Sterck, J. L. Kaiser, and J. C. Miller.  Analytic models for SIR disease spread on random spatial networks.  _Journal of Complex Networks_ , 6(6):948–970, 2018.
  * Lefever and Lejeune [1997] R. Lefever and O. Lejeune.  On the origin of Tiger Bush.  _Bulletin of Mathematical Biology_ , 59(2):263–294, 1997.
  * Meron [2018] E. Meron.  From patterns to function in living systems: Dryland ecosystems as a case study.  _Annual Reviews: Condensed Matter Physics_ , 9:79––103, 2018.  https://doi.org/10.1146/annurev-conmatphys-033117-053959.
  * Murray [2001] J. Murray.  _Mathematical Biology II: Spatial Models and Biomedical Applications_.  Springer New York, 2001.
  * Murray [2007] J. Murray.  _Mathematical Biology: I. An Introduction_ , volume 17.  Springer Science & Business Media, 2007.
  * Pastor-Satorras et al. [2015] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani.  Epidemic processes in complex networks.  _Reviews of Modern Physics_ , 87(3):925, 2015\.
  * Riley [2007] S. Riley.  Large-scale spatial-transmission models of infectious disease.  _Science_ , 316(5829):1298–1301, 2007.
  * Riley et al. [2015] S. Riley, K. Eames, V. Isham, D. Mollison, and P. Trapman.  Five challenges for spatial epidemic models.  _Epidemics_ , 10:68–71, 2015.
  * Sanz et al. [2014] J. Sanz, C.-Y. Xia, S. Meloni, and Y. Moreno.  Dynamics of interacting diseases.  _Physical Review X_ , 4(4):041005, 2014.
  * Segal et al. [2013] B. Segal, V. Volpert, and A. Bayliss.  Pattern formation in a model of competing populations with nonlocal interactions.  _Physica D: Nonlinear Phenomena_ , 253:12–22, 2013.
  * Sherratt and Smith [2008] J. A. Sherratt and M. J. Smith.  Periodic travelling waves in cyclic populations: field studies and reaction–diffusion models.  _Journal of the Royal Society Interface_ , 5(22):483–505, 2008.
  * Small et al. [2006] M. Small, C. K. Tse, and D. M. Walker.  Super-spreaders and the rate of transmission of the SARS virus.  _Physica D: Nonlinear Phenomena_ , 215(2):146–158, 2006.
  * Soriano-Paños et al. [2018] D. Soriano-Paños, L. Lotero, A. Arenas, and J. Gómez-Gardeñes.  Spreading processes in multiplex metapopulations containing different mobility networks.  _Physical Review X_ , 8(3):031039, 2018.
  * Sun et al. [2016] G.-Q. Sun, M. Jusup, Z. Jin, Y. Wang, and Z. Wang.  Pattern transitions in spatial epidemics: Mechanisms and emergent properties.  _Physics of Life Reviews_ , 19:43–73, 2016.
  * Thompson and Katul [2008] S. Thompson and G. Katul.  Plant propagation fronts and wind dispersal: an analytical model to upscale from seconds to decades using superstatistics.  _The American Naturalist_ , 171(4):468–479, 2008\.
  * Thompson and Katul [2009] S. Thompson and G. Katul.  Secondary seed dispersal and its role in landscape organization.  _Geophysical Research Letters_ , 36(2):L02402, 2009.
  * Thompson et al. [2008] S. Thompson, G. Katul, and S. M. McMahon.  Role of biomass spread in vegetation pattern formation within arid ecosystems.  _Water Resources Research_ , 44(10):W10421, 2008\.
  * Thompson et al. [2009] S. Thompson, G. Katul, J. Terborgh, and P. Alvarez-Loayza.  Spatial organization of vegetation arising from non-local excitation with local inhibition in tropical rainforests.  _Physica D: Nonlinear Phenomena_ , 238(13):1061–1067, 2009.
  * Tian and Yuan [2017] B. Tian and R. Yuan.  Traveling waves for a diffusive seir epidemic model with non-local reaction.  _Applied Mathematical Modelling_ , 50:432–449, 2017.
  * Turkyilmazoglu [2021] M. Turkyilmazoglu.  Explicit formulae for the peak time of an epidemic from the sir model.  _Physica D: Nonlinear Phenomena_ , 422:132902, 2021.
  * Valdano et al. [2015] E. Valdano, L. Ferreri, C. Poletto, and V. Colizza.  Analytical computation of the epidemic threshold on temporal networks.  _Physical Review X_ , 5(2):021005, 2015.
  * Volpert and Petrovskii [2009] V. Volpert and S. Petrovskii.  Reaction–diffusion waves in biology.  _Physics of Life Reviews_ , 6(4):267–310, 2009\.