# The propagation of infection fronts in spatially distributed compartment
models in epidemiology

Joseph Rudnick1, David Jasnow2, and Jorge Viñals3 1 Department of Physics and
Astronomy, University of California Los Angeles, CA 90095, 2 Department of
Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260, (3)
School of Physics and Astronomy, and Institute for Health Informatics,
University of Minnesota, Minneapolis, MN 55455, USA

###### Abstract

Spatio-temporal extensions of familiar compartment models for disease
transmission incorporating diffusive behavior, or interactions between
individuals at separate locations, are explored. The models considered have
the character of reaction-diffusion systems, which allow familiar techniques
to be applied. The focus is largely on the appearance of soliton-like moving
fronts that spread infection to previously uninfected regions. Near threshold
dynamical critical behavior and a degree of universality are revealed.
Extending two of the models to include a simple nonlinearity in the strength
of the binary interaction between a susceptible individual and an infected
one, we find the possibility of static coexistence between spatial regions
with different levels of infection and an analogy with first-order transitions
in thermodynamics.

##  I Introduction

The so called compartmental models in epidemiology follow the foundational
research of Kermack and McKendrick Kermack and McKendrick (1927), who
formulated kinetic equations to describe the temporal evolution of an epidemic
in terms of various subpopulations and their interactions. For example,
subpopulations can refer to individuals susceptible to a disease, those
exposed to it, infected, or recovered Anderson and May (1979); Murray (2001);
Li (2018). Initially space-independent, the solutions of these models provided
key insights into epidemic progression mechanisms, and into the existence of
steady states that were either disease free, or associated with persistent
infection. The models also focused on the characteristic rates of the various
stages of disease progression, and led to predictions of the characteristic
time scales and infection levels that followed. The emergence of SARS-CoV-2
has spurred new research into compartment models of disease spread, specially
given the wealth of fine grained, time and space dependent, empirical data.
Generalizations of compartment models in this case include, for example, the
consideration of fine grained compartments Parolini  _et al._ (2021); Yan  _et
al._ (2021), a cell based network that also incorporates subpopulation fluxes
to mimic population mobility Le Treut  _et al._ (2021), or the addition of
delays relating the newly infected and the newly recovered Huang  _et al._
(2021)

More generally, a key development in the study of population dynamics is the
work of Skellam Skellam (1951), which was inspired by earlier work of Fisher
Fisher (1937) on evolutionary genetics, and who introduced diffusive processes
(Brownian motion) to model species dispersal in Europe. The concept of
diffusive dispersal has been widely applied in Biology Belgacem (2001); Fife
(2013), including in the study of disease spreading and a theory of epidemics
Hoppensteadt and Peskin (2012); Postnikov and Sokolov (2007); Abdelheq  _et
al._ (2019). Beginning with the work of Noble Noble (1974), attention was paid
to the details of the geographic spread of epidemics by allowing and
explicitly modeling the spatio-temporal evolution of the subpopulations
involved. Therefore spatially dependent generalizations of compartment models
play an important role in describing the spread of disease as they intend to
capture the effect of mobile disease carriers across geographical regions.
Communication occurs through migration Noble (1974), or through contacts
between infected and susceptible individuals, and therefore predominantly
between adjacent geographical areas Naether  _et al._ (2008). This research
highlighted the possible existence of traveling fronts separating regions with
high prevalence of the disease from regions with healthier populations.
Results were obtained on the conditions for fronts to develop and propagate,
as well as their speed. Noteworthy is the analysis of Källén et al. Källén
(1984) of a model of the spread of rabies in Europe. The model studied is a
variant of the so called SI model (definitions of the various models are given
below), and it was shown to lead to a propagating front of infection, with a
velocity which is substantially equivalent to the results which we present
below. We mention here the general limitation that the spatial coupling or
correlation among compartments in these models is rather simple, as they are
generally restricted to a small spatial neighborhood of any given segment of
the population, and the dynamics follows rather idealized assumptions (e.g.,
simple diffusion).

More recently, compartment models have been extended to include complex
networks that incorporate a large number of distinct nodes (generalized
compartments). The nodes form an interconnected network with some specified
topology, and are allowed to interact in more general ways than simple
geographical proximity Rock  _et al._ (2014); Verdasca  _et al._ (2005); Eksin
_et al._ (2021). Recent examples of application of the methodology include the
study of waves in the spread of influenza Viboud  _et al._ (2006), an
empirical determination of a network of social contacts, and its consequences
on disease incidence Mossong  _et al._ (2008), or the spread of SARS-CoV-2
associated with a number of connectedness and mobility indices Holtz  _et al._
(2020); Paltiel  _et al._ (2020); Le Treut  _et al._ (2021). A modified SIR(S)
model has been introduced to study the effect of changing climate on basic
infectious disease reproduction number, and modeled heterogeneity by
introducing a set of coupled compartments into the model Saad-Roy  _et al._
(2020). Many current SARS-CoV-2 studies involve contact networks with a very
large number of nodes, on the order of one million (the population of a
typical metropolitan area) Nande  _et al._ (2021). The role of explicit delay
terms in the differential equations defining a compartment model has also been
considered to describe SARS-CoV-2 outbreaks in Illinois in Ref. Wong  _et al._
(2020). These methods are considered predictive in that model parameters are
estimated by statistical fits to past empirical data, and then models can be
propagated into the future in order to predict the evolution of an epidemic.
Finally, we mention results concerning front propagation in more general
models that consider nonlocal interaction kernels and diffusivity in SEIR
compartment models. The conditions for the existence of traveling solutions in
the model have been established Wu  _et al._ (2018), and later extended to
include the effects of temporal delays in the interaction kernels Wu and Ma
(2021).

The models that we analyze below are written in terms of coarse grained
variables describing the number of individuals belonging to a certain
subpopulation. For example, and following standard convention, we denote by
$S(x,t)$ the number density of susceptible individuals in an element of volume
around $x$ at time $t$, $E(x,t)$ is the number of individuals who have been
exposed to the infectious disease, leading to a number, $I(x,t)$, of infected
and $R(x,t)$ recovered individuals. The various models are described in their
corresponding section, but they all involve rate functions describing
transitions among the subpopulations, and either spatial diffusion of
individuals belonging to the various subpopulations, or space-dependent
disease transmission among them depending on the model. As a consequence of
their definition, all the models studied here belong to the class of reaction-
diffusion equations, and we therefore use the tools of bifurcation theory and
nonlinear dynamics to analyze solutions of interest. We restrict our analysis
to the simpler case of one dimensional spatial variations, although the models
are readily generalized to higher dimensions. We find in all cases that the
evolution of some of the dynamical variables satisfies the Fisher equation
(Appendix A) Fisher (1937); Kolmogorov  _et al._ (1991) _near onset_ of
infection propagation. This connection allows us to find analytic expressions
for the shape of the propagating front and its velocity in some cases, whereas
in others we proceed numerically to confirm the existence of fronts and obtain
their velocity.

In Section II we begin by considering a simple SI model that includes only two
populations, susceptible and infected, as first studied by Noble Noble (1974).
This is a model that addresses the evolution of a disease that spreads through
a population but that, at asymptotically long times, there are no infected
individuals remaining. We show in this case that a wave of infection exists
originating from an initial condition in which a few individuals are assumed
to be infected in some small region. If the infection rate is larger than the
recovery rate, a wave of infection sweeps across the entire system. The wave
envelope has constant shape as it propagates as a front of constant velocity,
which we compute numerically. The same method of analysis is used next to
describe the SIR model with vitality, with similar conclusions when the rate
of infection is larger that the rate of recovery plus the rate of population
renewal. The same results obtain for the more complex SEIR model, which we
show reduces, near threshold, to the Fisher-Kolmogorov model for each of the
subpopulations. We also look at a model which ignores spread through diffusion
of individual populations but takes into account non-locality of the infection
process. Finally we examine extended SIR and SEIR models in which the
infectivity rate coefficient is assumed to depend on the number of infected
individuals. The added nonlinearity allows bistability between two fixed
points, one featuring persistent infection, the other disease free. In this
case we find a unique front solution separating both fixed points, and compute
the front velocity for a range of parameter models.

##  II Noble’s SI model

In a seminal paper Noble (1974), J. V. Noble investigated a simple
epidemiological model for the propagation of an infection, which spreads by
means of a moving front. The model divides individuals into two populations,
consisting of those who are disease-free and susceptible to the infection
($S$) and those who have been infected ($I$). The ultimate fate of the members
of the $I$ population is to be removed from the system, either through death
or by virtue of acquired immunity. This model is appropriate to an epidemic
such as the bubonic plague in the pre-antibiotic era which rapidly sweeps
through a population, causing most of the deaths that occur in the time-frame
of its duration. The one-dimensional equations for this system are taken to be

$\displaystyle\frac{\partial S(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\frac{\partial^{2}S(x,t)}{\partial x^{2}}-\beta S(x,t)I(x,t)$
(1) $\displaystyle\frac{\partial I(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\frac{\partial^{2}I(x,t)}{\partial x^{2}}+\beta
S(x,t)I(x,t)-\mu I(x,t)$ (2)

The population of susceptible individuals, $S(x,t)$, changes through diffusive
motion, quantified by the coefficient $D$, and attrition due to infection, the
rate being determined by the extent of contact between the two populations and
an infectivity coefficient, $\beta$. The infected population, $I(x,t)$ changes
as a result of the process of infection, diffusion and the rate at which
individuals are removed through death or the acquisition of immunity, governed
by the coefficient $\mu$.

We start with the case of a uniform (spatially mixed) set of populations, in
which the model can be analytically solved with the use of quadratures, and
which guides the exploration of the transition between the conditions under
which the infection dies out locally and those in which it spreads and
develops into an epidemic.

###  II.1 The equations assuming uniform populations

In this case the spatial dependence in Eqs. (1) and (2) can be ignored. They
then reduce to

$\displaystyle\frac{dS(t)}{dt}$ $\displaystyle=$ $\displaystyle-\beta
S(t)I(t)$ (3) $\displaystyle\frac{dI(t)}{dt}$ $\displaystyle=$
$\displaystyle\beta S(t)I(t)-\mu I(t)$ (4)

Dividing (4) by (3), we end up with the single differential equation

$\frac{dI}{dS}=-1+\frac{\mu}{\beta S}$ (5)

Integrating this equation up, we find

$I=I_{0}-(S-S_{0})+\frac{\mu}{\beta}\ln(S/S_{0})$ (6)

where $I_{0}$ and $S_{0}$ are the initial values of the infected and
susceptible populations. For further discussion of the analysis of this model
see Brauer  _et al._ (2008).

Note the properties of the equations above. First, given (3), $S(t)$ is
monotonically decreasing. In light of this fact, the change with $t$ of $I(t)$
for small values of $t$ is controlled by the sign of the right hand side of
(5); if it is negative $I(t)$ increases with $t$, and if it is positive $I(t)$
decreases with $t$. Given the form of (5) and the monotonic behavior of $S(t)$
we know that eventually $I(t)$ will decrease with time. Whether it ever
increases is determined by the values of $\beta$, $\mu$ and $S_{0}$. When
$\beta>S_{0}\mu$, $I(t)$ decreases monotonically with $t$, and when
$\beta<S_{0}\mu$, there is an initial increase with time before the inevitable
decline of $I(t)$. As background for the analysis to follow, we set the
initial conditions to $I_{0}=0$ and $S_{0}=1$. This corresponds to an
uninfected system. However, we will use it as the limit of a system with an
extremely small infected population. The parametric solution (6) then becomes

$I=1-S+\frac{\mu}{\beta}\ln(S)$ (7)

At $S=1$, the slope of the curve—see Fig. 1—is

$\left.\frac{dI}{dS}\right|_{S=1}=\frac{\mu}{\beta}-1$ (8)

If $\beta>\mu$, the $I$ vs. $S$ curve has a negative slope. Given the time
dependence of $S(t)$ this means $I(t)$ increases from its initial value,
(near) zero, with time. If $\beta<\mu$ there is no such growth, and the
infection never spreads.

Figure 1: A plot of $I(t)$ versus $S(t)$ for various values of $\mu/\beta$.
The initial value of $S$ is 1 and the initial value of $I$ is 0 The time flow
for all plots is to the left, in the direction of decreasing $S(t)$. The plots
terminate when $I(t)=0$, at which point there is no change in either $S(t)$ or
$I(t)$. The points at which a curve intersects the horizontal axis are both
steady states. The state at the right is unstable and the state at the left is
stable.

As for the late-time behavior of $S(t)$ and $I(t)$, the last term on the right
hand side of (6) guarantees that there will be a value of $t$ for which $I(t)$
passes through zero. At that time, the absence of an infected population
causes the time derivative of $S(t)$ to vanish, while Eq. (5) guarantees that
$I(t)$ is also stationary, A dynamically stable steady state has thus been
reached. The value of $S$ at this point is determined by the ratio
$\mu/\beta$. There is an explicit formula for this value, obtained by solving
(7) with the left hand side set equal to 0. One finds for that value

$S_{f}(\mu/\beta)=-\frac{\mu}{\beta}W\left(-\frac{\beta}{\mu}e^{-\frac{\beta}{\mu}}\right)$
(9)

where $W(x)$ is the Lambert, or product log, function.Corless  _et al._ (1996)
Figure 2 shows how the final value of $S$ depends on the ratio $\mu/\beta$.

Figure 2: Plot of $S_{f}(\mu/\beta)$.

Finally, the full solutions of Eqs.(1) and (2) can be numerically completed by
quadratures. Substituting the right hand side of (7) for $I$ in (3), an
implicit solution for $S(t)$ is obtained by straightforward integration. Then,
inserting that solution into the right hand side of (4) we recover $I(t)$ by
numerical integration.

###  II.2 The spread of the infection as a propagating front: determination of
the front’s velocity

According to numerical evidence, after transients have died off a spreading
infection propagates as a front which retains its shape as it translates; i.e.
it is soliton-like. Figure 3 displays the results of the numerical solution of
Eqs. (1) and (2).

Figure 3: The solutions of the Noble equations Eqs. (1) and (2) for $I(x,t)$
and $S(x,t)$ plotted for three equally spaced times. The coefficients in those
equations are $D=300$, $\beta=1$ and $\mu=0.2$. The fronts propagate from the
left to the right. See also corresponding plots in Noble Noble (1974). As the
advancing front sweeps through the system the initial state of $S=1,I=0$
transforms to a final steady state in which $I=0$ and $S$ is reduced to a
smaller, but non-zero, value.

If the propagation of the fronts is indeed soliton-like, then we can write

$\displaystyle S(x,t)$ $\displaystyle=$ $\displaystyle S_{s}(x-vt)$ (10)
$\displaystyle I(x,t)$ $\displaystyle=$ $\displaystyle I_{s}(x-vt)$ (11)

where the subscript $s$ indicates soliton-like structure. If we define

$\displaystyle O_{S}$ $\displaystyle=$ $\displaystyle
v\frac{d}{dx}+D\frac{d^{2}}{dx^{2}}$ (12) $\displaystyle O_{I}$
$\displaystyle=$ $\displaystyle v\frac{d}{dx}+D\frac{d^{2}}{dx^{2}}+\mu$ (13)

the equations for $S_{s}$ and and $I_{s}$ that follow from (1) and (2) can be
reduced to

$O_{S}S_{s}(x)=-O_{I}I_{s}(x)=\beta S_{s}(x)I_{s}(x)$ (14)

This allows us to utilize the numerical solution of equation set (1) and (2)
to extract the value of the velocity, $v$ of propagation of the contagion. An
example of the determination of the propagation velocity of the fronts is
shown in Fig. 4. The three functions coincide for only one of the three
assumed velocities ($v=30.7$).

Figure 4: The three functions in (14), plotted for three assumed front
velocities, $v$. The three sets of curves are offset along the $x$ axis for
clarity of presentation. For only one of the three velocities—the actual
propagation velocity of the front—do the three functions coincide. The
coefficients are the same as in Fig. 3.

###  II.3 The transition between a dissipating and a propagating infection

As in the case of the uniform system considered in Section II.1, the evolution
of the system described by the spatiotemporal equations (1) and (2) changes
dramatically as the ratio $\beta/\mu$ passes through its threshold value. When
$\beta/\mu$ is larger than the value of the threshold, the infection
propagates through the system as a soliton-like front as noted above, while if
the ratio is less than the threshold the infection remains localized and
eventually dies off. This fundamental change in the nature of the contagion is
a dynamical phase transition, and like many phase transitions, it displays
distinctive properties, including diverging space and time scales governed by
critical exponents.

###  II.4 The uniform population just above threshold

Referring back to Section II.1, particularly Eq. (6), we can expand about the
threshold for increasing infected population. We set $I_{0}=0{+}$, $S_{0}=1$,
$S=1-\Delta S$, with $\Delta S\ll 1$ and $\mu/\beta=1-\epsilon$ with positive
$\epsilon\ll 1$ measuring the deviation from threshold. (Note that $I_{0}$, a
positive infinitesimal, and $S_{0}$ will refer to a time in the distant past.)
In this case, expanding to second order in $\Delta S$ we find

$I=\epsilon\Delta S-\frac{\Delta S^{2}}{2}$ (15)

We now turn to (3), the equation for the time dependence of $S(t)$. Using the
expansion above, we find

$\frac{d\Delta S}{dt}=\beta(\epsilon\Delta S-\Delta S^{2}/2)$ (16)

which is readily integrated up. We find

$\Delta S(t)=\frac{2\epsilon e^{\beta\epsilon t}}{1+e^{\beta\epsilon t}}$ (17)

If we plug this form into Eq. (15), we find

$I(t)=\frac{2\epsilon^{2}e^{\beta\epsilon t}}{(1+e^{\beta\epsilon t})^{2}}$
(18)

The net change in the susceptible population as a result of the propagating
infection vanishes as $\epsilon\rightarrow 0$ while the maximum value of the
infected population over the course of the infection scales as $\epsilon^{2}$.
Furthermore, the length of time over which the changes arising from the
infection take place diverges as $\epsilon^{-1}$.

###  II.5 A representation of the SI system immediately above threshold

We now know three things about the solution to the SI equations in the uniform
system near onset of propagating fronts. First, the difference between $S(t)$
and its uninfected value is of order $\epsilon$. Second, we know that the
magnitude of the infected population goes as $\epsilon^{2}$. Third, we know
that the time dependence slows down as $\epsilon\rightarrow 0$, the time scale
going as $1/\epsilon$.

We now set the initial susceptible population infinitesimally smaller than
one. We then set the infectivity coefficient $\beta$ equal to one and the
death rate $\mu$ equal to $1-\epsilon$. Then, we re-parameterize the
populations as follows

$\displaystyle S(x,t)$ $\displaystyle=$ $\displaystyle
1-\epsilon+\epsilon\sigma(\epsilon t,\sqrt{\epsilon}x)$ (19) $\displaystyle
I(x,t)$ $\displaystyle=$ $\displaystyle\epsilon^{2}\iota(\epsilon
t,\sqrt{\epsilon}x)$ (20)

The particular form of the presumed $x$-dependence is motivated by the desire
to assure that the diffusion terms behave in the same way as the time
derivative contributions as $\epsilon\rightarrow 0$. Substituting into Eqs.
(1) and (2) and extracting the leading order terms—which means we can neglect
the difference between $S(x,t)$ and 1 on the right hand side of the scaled
version of Eq. (1)—we end up with the following two equations describing the
behavior of the two functions $\sigma$ and $\iota$

$\displaystyle\frac{\partial\sigma(\tau,\chi)}{\partial\tau}$ $\displaystyle=$
$\displaystyle
D\frac{\partial^{2}\sigma(\tau,\chi)}{\partial\chi^{2}}-\iota(\tau,\chi)$ (21)
$\displaystyle\frac{\partial\iota(\tau,\chi)}{\partial\tau}$ $\displaystyle=$
$\displaystyle
D\frac{\partial^{2}\iota(\tau,\chi)}{\partial\chi^{2}}+\iota(\tau,\chi)\sigma(\tau,\chi)$
(22)

Here,

$\displaystyle\tau$ $\displaystyle=$ $\displaystyle\epsilon t$ (23)
$\displaystyle\chi$ $\displaystyle=$ $\displaystyle\sqrt{\epsilon}x$ (24)

The equation set (21) and (22) can be readily solved numerically.

Figure 5: The solutions to Eqs. (21) and (22) for three equally-spaced values
of the rescaled time $\tau$. The diffusion constant $D$ has been set equal to
300.

Figure 5 shows the result of implementing those equations for the case
$D=300$. This is clear evidence for soliton-like behavior. The behavior of the
curves for $\sigma(\tau,\chi)$ at the far left reflects the relaxation of
vestiges of initial conditions to the stable steady state at very large values
of $\tau$. We can check for soliton-like behavior by assuming that
$\sigma(\tau,\chi)$ has the form $\sigma_{s}(\chi-v^{\prime}\tau)$, and
similarly for $\iota$. , i.e.,

$\displaystyle\sigma(\chi,\tau)$ $\displaystyle=$
$\displaystyle\sigma_{s}(\chi-v^{\prime}\tau)$ (25)
$\displaystyle\iota(\chi,\tau)$ $\displaystyle=$
$\displaystyle\iota_{s}(\chi-v^{\prime}\tau)$ (26)

Then, (21) and (22) take the form

$\displaystyle
v^{\prime}\sigma_{s}^{\prime}(\chi-v^{\prime}\tau)+D\sigma_{s}^{\prime\prime}(\chi-v^{\prime}\tau)$
$\displaystyle=$ $\displaystyle\iota_{s}(\chi-v^{\prime}\tau)$ (27)
$\displaystyle
v^{\prime}\iota_{s}^{\prime}(\chi-v^{\prime}\tau)+D\iota_{s}^{\prime\prime}(\chi-v^{\prime}\tau)$
$\displaystyle=$
$\displaystyle-\iota_{s}(\chi-v^{\prime}\tau)\sigma_{s}(\chi-v^{\prime}\tau)$
(28)

Both equations will be satisfied with the correct choice of the selected
velocity $v^{\prime}$. The velocity can be determined inserting numerical
solutions of (21) and (22) for a fixed value of $t$ into the right and left
hand sides of (27) and (28) and then determining the value of $v^{\prime}$ for
which those equations are satisfied. Numerics show that the proper choice for
the velocity, $v^{\prime}$, is

$v^{\prime}\simeq 2\sqrt{D}$ (29)

Given the relationships (23) and (24) this corresponds to a velocity in the
original coordinates scaling as

$\displaystyle v$ $\displaystyle=$ $\displaystyle\epsilon^{1/2}v^{\prime}$
(30) $\displaystyle\simeq$ $\displaystyle 2\sqrt{\epsilon D}$

.

##  III SIR model with vitality

It is possible to amend the Noble SI equation set (1), (2) by adding an
equation for the dynamics of the population of “removed” individuals,
$R(x,t)$. This new equation is

$\frac{\partial R(x,t)}{\partial t}=D\frac{\partial^{2}R(x,t)}{dx^{2}}+\mu
I(x,t)$ (31)

The resulting equation set is an SIR model for the evolution of a contagion.
No new behavior results from the expansion of the model to incorporate the
population of individuals who have contracted the infection and either died or
survived with protection against reinfection. However, there is an interesting
variant of the model consisting of Eqs. (1), (2) and (31). This is the SIR
model with vitality, consisting of the following equations,

$\displaystyle\frac{\partial S(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\frac{\partial^{2}S(x,t)}{\partial x^{2}}+\mu N(x,t)-\mu
S(x,t)-\beta\frac{S(x,t)I(x,t)}{N(x,t)}$ (32) $\displaystyle\frac{\partial
I(x,t)}{\partial t}$ $\displaystyle=$ $\displaystyle
D\frac{\partial^{2}I(x,t)}{\partial
x^{2}}+\beta\frac{S(x,t)I(x,t)}{N(x,t)}-\gamma I(x,t)-\mu I(x,t)$ (33)
$\displaystyle\frac{\partial R(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\frac{\partial^{2}R(x,t)}{\partial x^{2}}+\gamma I(x,t)-\mu
R(x,t)$ (34)

where

$S(x,t)+I(x,t)+R(x,t)=N(x,t)$ (35)

$N(x,t)$ being the total population (density) at location $x$ and time $t$.
This SIR model describes a collection of individuals with a birth rate, $\mu$,
equal to the death rate, which is the same for all three populations, $S$, $I$
and $R$; every individual in this system is born susceptible. All other terms
refer to diffusion of the populations and the transition of individuals from
one population to another. For instance the term $\gamma I(x,t)$ describes the
rate at which infected individuals recover and are hence removed.

Adding (32)–(34) we find

$\frac{\partial N(x,t)}{\partial t}=D\frac{\partial^{2}N(x,t)}{\partial
x^{2}}$ (36)

This tells us that the overall population tends to an $x$ and $t$ independent
distribution. As a consequence we can assume a constant value for $N(x,t)=N$.
This simplifies matters, and we can replace $S(x,t)$ with $Ns(x,t)$ and
similarly for $I(x,t)$ and $R(x,t)$. We then have the relationship

$s(x,t)+i(x,t)+r(x,t)=1$ (37)

###  III.1 Uniform version of the SIR equations with vitality

The three equations in this case are

$\displaystyle\frac{ds(t)}{dt}$ $\displaystyle=$
$\displaystyle\mu(1-s(t))-\beta s(t)i(t)$ (38) $\displaystyle\frac{di(t)}{dt}$
$\displaystyle=$ $\displaystyle\beta s(t)i(t)-(\gamma+\mu)i(t)$ (39)
$\displaystyle\frac{dr(t)}{dt}$ $\displaystyle=$ $\displaystyle\gamma i(t)-\mu
r(t)$ (40)

Equation (34) is actually redundant, as we can also use Eq. (37) to determine
$r(t)$. Furthermore, the variable $r(t)$ appears nowhere in Eqs. (38) and
(39). If we assume as an initial condition $s(t)\approx 1$, with $r(0)=0$ and
$i(0)=1-s(0)-r(0)$ very small, we see from (39) that the requirement for a
growing infection is

$\beta>\gamma+\mu$ (41)

###  III.2 Propagation of infection as a moving front: determination of the
front’s velocity

As in the case of the SI equations of Section II, when $\beta$ exceeds a
threshold, in this case the inequality (41) being satisfied, and after
transients have subsided, an initially localized region of infection will
spread through the system as a soliton-like front. Figure 6 is a snapshot of
the front, which propagates by translation from left to right.

Figure 6: A snapshot of the propagating infection front, in particular the
values of $S(x,t)$, $I(x,t)$ and $R(x,t)$. In Eqs. (32)–(34). $N(x,t)$ has
been set equal to 1, and the coefficients in the equation set are $D=300$,
$\gamma=0.3$, $\mu=0.1$ and $\beta=1$. All functions translate at the same
fixed speed from left to right.

The infected population, $I(x,t)$, differs in its behavior from the the
corresponding population in the Noble SI model, in that it does not relax to
zero in the aftermath of the spreading front. Rather, as a result of the
sustaining birth rate, the system relaxes to an “endemic” state in which a
non-zero fraction of all individuals are infected. The remainder is either
uninfected and susceptible or previously infected and now recovered. The
values of the coefficients $\gamma$, $\mu$ and $\beta$ control the
distribution of the individuals among those three populations. Solving for a
dynamically stable steady state solution to the constitutive uniform system
equations, i.e., (32)–(34) with left hand sides set equal to zero and $N(x,t)$
set equal to 1, we find

$\displaystyle S_{SS}$ $\displaystyle=$
$\displaystyle\frac{\gamma+\mu}{\beta}$ (42) $\displaystyle I_{SS}$
$\displaystyle=$
$\displaystyle\frac{\mu(\beta-\gamma-\mu)}{\beta(\gamma+\mu)}$ (43)
$\displaystyle R_{SS}$ $\displaystyle=$
$\displaystyle\frac{\gamma(\beta-\gamma-\mu)}{\beta(\gamma+\mu)}$ (44)

where the subscript $SS$ stands for “steady state.”

To find the velocity of the moving front, we follow the steps in Section II.2.
We re-express $S(x,t)$ and $I(x,t)$ as in Eqs. (10) and (11). Then, defining

$\displaystyle O_{S}$ $\displaystyle=$ $\displaystyle
v\frac{d}{dx}+D\frac{d^{2}}{dx^{2}}+\mu$ (45) $\displaystyle O_{I}$
$\displaystyle=$ $\displaystyle
v\frac{d}{dx}+D\frac{d^{2}}{dx^{2}}+\mu+\gamma$ (46)

the following equations will be satisfied for the proper choice of the front
velocity, $v$,

$O_{S}S_{s}(x)-\mu=-O_{I}I_{s}(x)=\beta S_{s}(x)I_{s}(x)$ (47)

The expression for $R_{s}(x-vt)$ follows from the results above and

$S_{s}(x-vt)+I_{s}(x-vt)+R_{s}(x-vt)=N$ (48)

###  III.3 Behavior of the propagating infection front just above onset:
uniform system analysis

In the vicinity of onset, we write

$\beta=\gamma+\mu+\epsilon$ (49)

with $\epsilon\ll 1$. Anticipating the behavior in the vicinity of onset, we
make use of the re-parameterization in Section II.5, writing

$\displaystyle s(t)$ $\displaystyle=$ $\displaystyle 1+\epsilon\sigma(\tau)$
(50) $\displaystyle i(t)$ $\displaystyle=$ $\displaystyle\epsilon\iota(\tau)$
(51) $\displaystyle r(t)$ $\displaystyle=$ $\displaystyle\epsilon\rho(\tau)$
(52)

where $\tau=\epsilon t$. In light of the sum rule (37), the following
relationship now holds

$\sigma(\tau)+\iota(\tau)+\rho(\tau)=0$ (53)

(Note that, as defined, $\sigma<0$.) We start by focusing on the equation for
the evolution of the infected population, which can now be written

$\epsilon^{2}\frac{d\iota(\tau)}{d\tau}=\epsilon(\beta-(\gamma+\mu))\iota(\tau)+\epsilon^{2}\beta\sigma(\tau)\iota(\tau)$
(54)

Using (49) we transform (54) into

$\frac{d\iota(\tau)}{d\tau}=\iota(\tau)+\beta\iota(\tau)\sigma(\tau)$ (55)

We now look at Eq. (38) under the substitutions (50)–(52). Keeping leading-
order-in-$\epsilon$ contributions to the right hand side of that equation we
find

$\displaystyle\epsilon\frac{d\sigma(\tau)}{d\tau}$ $\displaystyle=$
$\displaystyle-\mu\sigma(\tau)-\beta\iota(\tau)$ (56)
$\displaystyle\rightarrow$
$\displaystyle-\mu\sigma(\tau)-(\gamma+\mu)\iota(\tau)$

Given the factor of $\epsilon$ in front of the derivative above, the leading-
order-in-$\epsilon$ solution of this equation for $\sigma(\tau)$ is

$\sigma(\tau)=-\frac{\gamma+\mu}{\mu}\iota(\tau)$ (57)

By a similar process we find, to leading order,

$\rho(\tau)=\frac{\gamma}{\mu}\iota(\tau)$ (58)

These last two relations are entirely consistent with the sum rule (53).
Making use of (57), Eq. (55) becomes, to leading order in $\epsilon$,

$\frac{d\iota(\tau)}{d\tau}=\epsilon\iota(\tau)\left(1-\frac{(\gamma+\mu)^{2}}{\mu}\iota(\tau)\right)$
(59)

If we introduce the spatial dependence, add diffusion and re-parameterize as
in Section II.5, we end up with the well-studied Fisher-KPP equation (see Eq.
(163) and accompanying discussion). This version of the equation is

$\frac{\partial\iota(\chi,\tau)}{\partial\tau}=D\frac{\partial\iota(\chi,\tau)}{\partial\chi^{2}}+\iota(\chi,\tau)\left(1-\frac{(\gamma+\mu)^{2}}{\mu}\iota(\chi,\tau)\right)$
(60)

Once again, the velocity of propagation of the front can be established
numerically by assuming an $\iota(\chi,\tau)$ of the form
$\iota_{s}(\chi-v^{\prime}\tau)$, so that the time derivative in (60) is
replaced by $-v^{\prime}\partial\iota/\partial\chi$, and then determining the
value of $v^{\prime}$ for which the equation is satisfied. The end result of
this analysis is

$v^{\prime}=c_{2}\sqrt{D}$ (61)

where $c_{2}\simeq 2$ in accord with well known results Fisher (1937);
Kolmogorov  _et al._ (1991); Fife and McLeod (1977) for the FKPP equation.
Returning to the original variables $x$ and $t$, we find for the front’s
velocity in the vicinity of onset

$\displaystyle v$ $\displaystyle=$ $\displaystyle\epsilon^{1/2}v^{\prime}$
$\displaystyle=$ $\displaystyle 2\sqrt{\epsilon D}$

which exhibits the characteristic slowing down of the front velocity in the
vicinity of onset.

##  IV A model incorporating non-locality of the infection process

Naether et. al. Naether  _et al._ (2008) proposed an epidemiological model
that takes into account the mobility of the various populations and which
results in a spatially extended infection process. Assuming that the the
disease can be transmitted to a susceptible individual from members of the
infected population in the neighborhood of that individual, and making use of
a gradient expansion, the following continuum model is proposed,

$\displaystyle\frac{\partial S(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle-\beta I(x,t)S(x,t)-BS(x,t)\frac{\partial^{2}I(x,t)}{\partial
x^{2}}$ (62) $\displaystyle\frac{\partial I(x,t)}{\partial t}$
$\displaystyle=$ $\displaystyle\beta
I(x,t)S(x,t)+BS(x,t)\frac{\partial^{2}I(x,t)}{\partial x^{2}}-\mu I(x,t)$ (63)
$\displaystyle\frac{\partial R(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle\mu I(x,t)$ (64)

For the uniform system, i.e., in the absence of spatial derivatives on the
right hand sides, the system above is just Noble’s SI model for uniform
populations (See Section II.1) with an added equation governing the removed
individuals—the SIR equations without vitality. Numerical solution of Eqs.
(62)–(64) yields propagating fronts when, as in the case of the SI model,
$\beta>\mu$. Such a propagating front is shown in Fig. 7

Figure 7: A snapshot of the propagating front as generated by Eqs. (62)–(64).
The coefficients are as follows: $\beta=1$, $\mu=0.5$, $B=300$. The fronts
move from left to right. The higher rate of attrition of the infected
population ($\mu=0.5$ here as opposed to $\mu=0.2$ in Fig. 3) leads to a
larger surviving susceptible population.

###  IV.1 Onset analysis

The uniform equations as noted, are the just the SIR model without vitality.
To investigate onset, we write $\mu=\beta(1-\epsilon)$ and make use of the
parameterization in (19) and (20). Under the assumption of a scaling solution
the “reduced” versions of the constitutive equations corresponding to (27) and
(28) are

$\displaystyle\frac{\partial\sigma(\chi,\tau)}{\partial\tau}$ $\displaystyle=$
$\displaystyle-\beta\iota(\chi,\tau)$ (65)
$\displaystyle\frac{\partial\iota(\chi,\tau)}{\partial\tau}$ $\displaystyle=$
$\displaystyle\beta\sigma(\chi,\tau)\iota(\chi,\tau)+B\frac{\partial^{2}\iota(\chi,\tau)}{\partial\chi^{2}}$
(66)

We then assume soliton-like behavior, as embodied in

$\displaystyle-v^{\prime}\sigma_{s}^{\prime}(\chi-v^{\prime}\tau)$
$\displaystyle=$ $\displaystyle-\beta\iota_{s}(\chi-v^{\prime}\tau)$ (67)
$\displaystyle-v^{\prime}\iota_{s}^{\prime}(\chi-v^{\prime}\tau)$
$\displaystyle=$
$\displaystyle\beta\sigma_{s}(\chi-v^{\prime}\tau)\iota_{s}(\chi-v^{\prime}\tau)+B\iota_{s}^{\prime\prime}(\chi-v^{\prime}\tau)$
(68)

Finding, numerically, the value of $v^{\prime}$ that satisfies Eqs. (67) and
(68) we discover that the velocity, $v^{\prime}$, in the vicinity of onset is
given by

$v=c\sqrt{\epsilon\beta B}$ (69)

where $c\simeq 2.0$.

##  V The SEIR model

An extension of the SIR model with vitality adds a population of exposed
individuals ($E$) and an associated, additional time scale. The members of
this population are assumed to have taken up the pathogen through the
infection process but are not symptomatic or able to pass it on, which is to
say they are neither sick nor carriers. However, over the course of time, some
fraction of the population does become truly infected. The rest are subject to
the same other hazards as the susceptibles so have the same (background) death
rate as the individuals in every other population. The equations governing
this model are

$\displaystyle\frac{\partial S(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\frac{\partial^{2}S(x,t)}{dx^{2}}+\mu N(x,t)-\mu
S(x,t)-\frac{\beta S(x,t)I(x,t)}{N}$ (70) $\displaystyle\frac{\partial
E(x,t)}{\partial t}$ $\displaystyle=$ $\displaystyle
D\frac{\partial^{2}E(x,t)}{dx^{2}}+\frac{\beta S(x,t)I(x,t)}{N}-(\mu+a)E(x,t)$
(71) $\displaystyle\frac{\partial I(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\frac{\partial^{2}I(x,t)}{dx^{2}}+aE(x,t)-(\mu+\gamma)I(x,t)$
(72) $\displaystyle\frac{\partial R(x,t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\frac{\partial^{2}R(x,t)}{dx^{2}}+\gamma E(x,t)-\mu R(x,t)$
(73)

with

$N(x,t)=S(x,t)+E(x,t)+I(x,t)+R(x,t)$ (74)

The coefficient $a$ in Eqs. (71) and (72) quantifies the rate of transition of
exposed individuals to the state of being infected. It can readily be
established that if $N(x,t)$ as given by (74) is spatially constant, then it
remains unchanged under the action of Eqs. (70)–(73). This means that $R(x,t)$
obeys the relationship

$R(x,t)=N-S(x,t)-E(x,t)-I(x,t)$ (75)

with the total population density, $N$, fixed in space and time. From here on
we will set $N=1$.

###  V.1 Onset: the ODEs

The analysis of onset begins with the consideration of the results of the
fully mixed model, i.e., the first order differential equations that ignore
spatial dependence. Those equations are (see, e.g., (70)–(73))

$\displaystyle\frac{dS(t)}{dt}$ $\displaystyle=$
$\displaystyle\mu(1-S(t))-\beta S(t)I(t)$ (76) $\displaystyle\frac{dE(t)}{dt}$
$\displaystyle=$ $\displaystyle\beta S(t)I(t)-(\mu+a)E(t)$ (77)
$\displaystyle\frac{dI(t)}{dt}$ $\displaystyle=$ $\displaystyle
aE(t)-(\mu+\gamma)I(t)$ (78) $\displaystyle\frac{dR(t)}{dt}$ $\displaystyle=$
$\displaystyle\gamma I(t)-\mu R(t)$ (79)

We can ignore the last equation for two reasons. First, given the comments
above, we can simply use the relation $S(t)+E(t)+I(t)+R(t)=1$ to solve for
$R(t)$ given the values of the other populations. Second, the population
$R(t)$ does not enter into the first three equations, which can be treated as
a self-contained system.

As a first step, we solve for fixed point solutions to (77)–(79).
Straightforward algebra leads to two solutions; the first fixed point
corresponds to a population free of infection: $S=1,E=0,I=0,R=0$. The second
set of solutions is

$\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{(a+\mu)(\gamma+\mu)}{a\beta}$ (80) $\displaystyle E$
$\displaystyle=$
$\displaystyle\frac{\mu(a\beta-(\mu+a)(\gamma+\mu))}{a\beta(a+\mu)}$ (81)
$\displaystyle I$ $\displaystyle=$
$\displaystyle\frac{\mu(a\beta-(\mu+a)(\gamma+\mu))}{\beta(a+\mu)(\gamma+\mu)}$
(82) $\displaystyle R$ $\displaystyle=$
$\displaystyle\frac{\gamma(a\beta-(\mu+a)(\gamma+\mu))}{\beta(a+\mu)(\gamma+\mu)}$
(83)

The solutions above are consistent with realizable population values – in that
$E$, $I$ and $R$ are positive – only if $\beta>(\mu+a)(\mu+\gamma)/a$. When
this inequality is satisfied, the relations listed above describe the
dynamically stable steady state of the equations (80)–(83).

Now, consider the implications when $\beta=(1+\epsilon)(\mu+a)(\mu+\gamma)/a$,
with $\epsilon\ll 1$. To first order in $\epsilon$, the fixed point solutions
are

$\displaystyle S$ $\displaystyle=$ $\displaystyle 1-\epsilon$ (84)
$\displaystyle E$ $\displaystyle=$ $\displaystyle\frac{\mu\epsilon}{a+\mu}$
(85) $\displaystyle I$ $\displaystyle=$
$\displaystyle\frac{a\mu\epsilon}{(a+\mu)(\gamma+\mu)}$ (86) $\displaystyle R$
$\displaystyle=$ $\displaystyle\frac{a\gamma\epsilon}{(a+\mu)(\gamma+\mu)}$
(87)

Next, we reparameterize the populations. Writing

$\displaystyle S(t)$ $\displaystyle=$ $\displaystyle 1+\epsilon s(t)$ (88)
$\displaystyle E(t)$ $\displaystyle=$
$\displaystyle\frac{\mu\epsilon}{\mu+a}e(t)$ (89) $\displaystyle I(t)$
$\displaystyle=$ $\displaystyle\frac{a\mu\epsilon}{(\mu+a)(\mu+\gamma)}i(t)$
(90)

with a similar—but unneeded—reparameterization for $R(t)$, we insert the right
hand sides (88)–(90) along with $\beta=(1+\epsilon)(\mu+a)(\mu+\gamma)/a$ into
Eqs. (76)–(78). Expanding to second order in $\epsilon$, we find for the three
relevant equations of motion,

$\displaystyle\epsilon\frac{ds(t)}{dt}$ $\displaystyle=$
$\displaystyle-\epsilon\mu(s(t)+i(t))-\epsilon^{2}\mu i(t)(1+s(t))$ (91)
$\displaystyle\epsilon\frac{\mu}{\mu+\alpha}\frac{de(t)}{dt}$ $\displaystyle=$
$\displaystyle\epsilon\mu(i(t)-e(t))+\epsilon^{2}\mu i(t)(1+s(t))$ (92)
$\displaystyle\epsilon\frac{a\mu}{(\mu+a)(\mu+\gamma)}\frac{di(t)}{dt}$
$\displaystyle=$
$\displaystyle\epsilon\frac{a\mu}{(\mu+a)(\mu+\gamma)}(e(t)-i(t))$ (93)

Note that there is no term going as $\epsilon^{2}$ on the right hand side of
(93). Under most circumstances we would be justified in ignoring all terms
going as $\epsilon^{2}$. However, in this case, with the addition of an
additional subpopulation (E), the higher-order in $\epsilon$ contributions to
the SEIR equations near onset prove to be quite consequential, and the
analysis of onset is more involved.

We start by looking at the the leading order terms in (91)–(93). Dividing by
$\epsilon$, the equations are

$\displaystyle\frac{ds(t)}{dt}$ $\displaystyle=$
$\displaystyle-\mu(s(t)+i(t))$ (94)
$\displaystyle\frac{\mu}{\mu+\alpha}\frac{de(t)}{dt}$ $\displaystyle=$
$\displaystyle\mu(i(t)-e(t))$ (95)
$\displaystyle\frac{a\mu}{(\mu+a)(\mu+\gamma)}\frac{di(t)}{dt}$
$\displaystyle=$ $\displaystyle\frac{a\mu}{(\mu+a)}(e(t)-i(t))$ (96)

Recasting the variables $s(t)$, $e(t)$ and $r(t)$ as the components of a
column vector and simplifying notation, the relations above can be rewritten
as

$\frac{d}{dt}\left(\begin{array}[]{c}s(t)\\\ e(t)\\\
i(t)\end{array}\right)=\left(\begin{array}[]{ccc}-A&0&-A\\\ 0&-B&B\\\
0&C&-C\\\ \end{array}\right)\left(\begin{array}[]{c}s(t)\\\ e(t)\\\
i(t)\end{array}\right)$ (97)

where

$\displaystyle A$ $\displaystyle=$ $\displaystyle\mu$ (98) $\displaystyle B$
$\displaystyle=$ $\displaystyle\mu+a$ (99) $\displaystyle C$ $\displaystyle=$
$\displaystyle\mu+\gamma$ (100)

To obtain the solution to this equation, we express the matrix on the right
hand side as a “spectral representation” in terms of its eigenvalues and
associated projection operators. The three eigenvalues are 0, $-A$ and $-B-C$.
The projection operators satisfy ${\bf P}_{i}\cdot{\bf P}_{j}=\delta_{i,j}{\bf
P}_{i}$, and are, in the order that the eigenvalues are listed,

$\displaystyle{\bf P}_{1}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccc}0&-\frac{C}{B+C}&-\frac{B}{B+C}\\\
0&\frac{C}{B+C}&\frac{B}{B+C}\\\ 0&\frac{C}{B+C}&\frac{B}{B+C}\\\
\end{array}\right)$ (104) $\displaystyle{\bf P}_{2}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccc}1&\frac{C}{-A+B+C}&-\frac{A-B}{-A+B+C}\\\
0&0&0\\\ 0&0&0\\\ \end{array}\right)$ (108) $\displaystyle{\bf P}_{3}$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccc}0&\frac{AC}{(B+C)(A-B-C)}&-\frac{AC}{(B+C)(A-B-C)}\\\
0&\frac{B}{B+C}&-\frac{B}{B+C}\\\ 0&-\frac{C}{B+C}&\frac{C}{B+C}\\\
\end{array}\right)$ (112)

The solution to (97) can then be written in the form

$\left(\begin{array}[]{c}s(t)\\\ e(t)\\\ i(t)\end{array}\right)=\left({\bf
P}_{1}+e^{-At}{\bf P_{2}}+e^{-(B+C)t}{\bf
P}_{3}\right)\cdot\left(\begin{array}[]{c}s(0)\\\ e(0)\\\
i(0)\end{array}\right)$ (113)

Note that at large values of $t$, the solution (113) becomes

${\bf P}_{1}\cdot\left(\begin{array}[]{c}s(0)\\\ e(0)\\\
i(0)\end{array}\right)=\left(\begin{array}[]{c}-\frac{Ce(0)}{B+C}-\frac{Bi(0)}{B+C}\\\
\frac{Ce(0)}{B+C}+\frac{Bi(0)}{B+C}\\\
\frac{Ce(0)}{B+C}+\frac{Bi(0)}{B+C}\end{array}\right)$ (114)

The solutions converge to $s(t)=-e(t)$ and $e(t)=i(t)$ but to limits
determined by initial values and not to the dynamical fixed point,
$s(t)=-1,e(t)=i(t)=1$. Approach to those fixed points is driven by the next-
to-leading-order terms.

Restoring the next-to-leading order terms, in generality we now have the set
of equations

$\frac{d}{dt}\left(\begin{array}[]{c}s(t)\\\ e(t)\\\
i(t)\end{array}\right)=\left(\begin{array}[]{ccc}-A&0&-A\\\ 0&-B&B\\\
0&C&-C\\\ \end{array}\right)\left(\begin{array}[]{c}s(t)\\\ e(t)\\\
i(t)\end{array}\right)+\epsilon\left(\begin{array}[]{c}-\mu i(t)(1+s(t))\\\
\mu i(t)(1+s(t))\\\ 0\end{array}\right)$ (115)

We then write

$\left(\begin{array}[]{c}s(t)\\\ e(t)\\\ i(t)\end{array}\right)=\left({\bf
P}_{1}+e^{-At}{\bf P_{2}}+e^{-(B+C)t}{\bf
P}_{3}\right)\cdot\left(\begin{array}[]{c}s^{\prime}(t)\\\ e^{\prime}(t)\\\
i^{\prime}(t)\end{array}\right)$ (116)

This allows us to cancel the leading order term on the right hand side of
(115). The end result is

$\frac{d}{dt}\left(\begin{array}[]{c}s^{\prime}(t)\\\ e^{\prime}(t)\\\
i^{\prime}(t)\end{array}\right)=\epsilon\left({\bf P}_{1}+e^{At}{\bf
P_{2}}+e^{(A+B)t}{\bf P}_{3}\right)\cdot\left(\begin{array}[]{c}-\mu
i(t)(1+s(t))\\\ \mu i(t)(1+s(t))\\\ 0\end{array}\right)$ (117)

Integrating up and making use of (116) we find

$\displaystyle\left(\begin{array}[]{c}s(t)\\\ e(t)\\\ i(t)\end{array}\right)$
(121) $\displaystyle=$ $\displaystyle\epsilon\int_{0}^{t}{\bf
P}_{1}\cdot\left(\begin{array}[]{c}-\mu i(t^{\prime})(1+s(t^{\prime}))\\\ \mu
i(t^{\prime})(1+s(t^{\prime}))\\\ 0\end{array}\right)\,dt^{\prime}$ (133)
$\displaystyle+\epsilon\int_{0}^{t}e^{-A(t-t^{\prime})}{\bf
P}_{2}\cdot\left(\begin{array}[]{c}-\mu i(t^{\prime})(1+s(t^{\prime}))\\\ \mu
i(t^{\prime})(1+s(t^{\prime}))\\\ 0\end{array}\right)\,dt^{\prime}$
$\displaystyle+\epsilon\int_{0}^{t}e^{-(B+C)(t-t^{\prime})}{\bf
P}_{3}\cdot\left(\begin{array}[]{c}-\mu i(t^{\prime})(1+s(t^{\prime}))\\\ \mu
i(t^{\prime})(1+s(t^{\prime}))\\\ 0\end{array}\right)\,dt^{\prime}$

Given that we expect the populations to evolve slowly, with a time constant of
order $1/\epsilon$, we can find the leading order contributions to the second
and third terms on the right hand side of (133). They yield

$(\frac{\epsilon}{A}{\bf P}_{2}+\frac{\epsilon}{B+C}{\bf
P}_{3})\cdot\left(\begin{array}[]{c}-\mu i(t)(1+s(t))\\\ \mu i(t)(1+s(t))\\\
0\end{array}\right)$ (134)

Finally, we take the time derivative of the resulting equation. In light of
(134) and the expectation of slow evolution of the densities, the time
derivatives of the second and third terms in (133) will introduce a further
factor of $\epsilon$ and can thus be discarded. This leaves us with

$\displaystyle\frac{d}{dt}\left(\begin{array}[]{c}s(t)\\\ e(t)\\\
i(t)\end{array}\right)$ $\displaystyle=$ $\displaystyle\epsilon{\bf
P}_{1}\cdot\left(\begin{array}[]{c}-\mu i(t)(1+s(t))\\\ \mu i(t)(1+s(t))\\\
0\end{array}\right)$ (141)

Given Eq. (104), (141) becomes

$\frac{d}{dt}\left(\begin{array}[]{c}s(t)\\\ e(t)\\\
i(t)\end{array}\right)=\epsilon\frac{\mu(\mu+\gamma)}{2\mu+a+\gamma}\left(\begin{array}[]{c}-i(t)(1+s(t))\\\
i(t)(1+s(t))\\\ i(t)(1+s(t))\end{array}\right)$ (142)

where we have used (98)–(100). Given that the equations are consistent with
$e(t)=i(t)=-s(t)$, we can recast the three resulting equations as follows

$\displaystyle\frac{ds(t)}{dt}$ $\displaystyle=$
$\displaystyle\epsilon\frac{\mu(\mu+\gamma)}{2\mu+a+\gamma}s(t)(1+s(t))$ (143)
$\displaystyle\frac{de(t)}{dt}$ $\displaystyle=$
$\displaystyle\epsilon\frac{\mu(\mu+\gamma)}{2\mu+a+\gamma}e(t)(1-e(t))$ (144)
$\displaystyle\frac{di(t)}{dt}$ $\displaystyle=$
$\displaystyle\epsilon\frac{\mu(\mu+\gamma)}{2\mu+a+\gamma}i(t)(1-i(t))$ (145)

These are all consistent with the Fisher-KPP equation as a description of the
scaled densities in the vicinity of onset, as in the SIR equation with
vitality – recall that $s(t)$ is negative. See Section III for the full
narrative leading to that conclusion.

##  VI Nonlinear infectivity constant and coexistence of two stable states

All the models discussed so far assume that the likelihood of transmission
goes as the product of the two population densities $S(x,t)$ and $I(x,t)$.
This amounts to ignoring possible “threshold” effects, in which the rate of
spread of the infection increases nonlinearly with the level of infection.
Such effects can be encoded in an $I$-dependence of the transmission
coefficient $\beta$, so that $\beta\rightarrow\beta(I(x,t))$. As we shall see,
a simple $I$ dependence leads in the modified SIR and SEIR models to the
possibility of two dynamically stable steady states, one in which there is an
endemic level of infection while the other is free of the disease. We start by
considering a spatially uniform system.

###  VI.1 The spatially uniform SIR system with an $I$-dependent $\beta$

We return to Eqs. (38)–(40) in which, now, $\beta$ takes the form

$\beta(I)=\beta_{0}+\alpha\frac{I}{N}$ (146)

We then use the conservation of $s(t)+i(t)+r(t)$ to replace the sum by the
constant value 1. Next, we choose a set of values for the coefficients in the
resulting equations. Here, they are

$\displaystyle\gamma=0.3,\mu=0.2,\beta_{0}=0.375,\alpha=5$ (147)

One important criterion for the choice of coefficients is that the value of
$\beta_{0}$, the rate of infection at vanishingly small $I$, is below the
threshold set by (41). in this case that threshold is 0.5.

We then seek the stationary solutions of the first two of Eqs. (38)–(40) (the
solution of the fourth one following from $s(t)+e(t)+i(t)=1$) with $\beta$
replaced by $\beta(I)$ given by (146). In this case one finds three sets of
stationary solutions

  * •

Solution 1:

$\displaystyle s=1$ (148)

  * •

Solution 2:

$\displaystyle
s=\frac{\alpha\mu+\sqrt{\left(\alpha\mu-\beta_{0}\beta_{\text{th}}\right){}^{2}-4\alpha\mu\beta_{\text{th}}\left(\beta_{\text{th}}-\beta_{0}\right)}+\beta_{0}\beta_{\text{th}}}{2\alpha\mu}$
(149)

  * •

Solution 3:

$\displaystyle s$ $\displaystyle=$
$\displaystyle\frac{\alpha\mu-\sqrt{\left(\alpha\mu-\beta_{0}\beta_{\text{th}}\right){}^{2}-4\alpha\mu\beta_{\text{th}}\left(\beta_{\text{th}}-\beta_{0}\right)}+\beta_{0}\beta_{\text{th}}}{2\alpha\mu}$
(150)

In the expressions above, $\beta_{\rm th}$ is the threshold value for the
coefficient $\beta$, in the case of an $I$-independent $\beta$ given by

$\beta_{\rm th}=\gamma+\mu$ (151)

(See Eq. (41) and accompanying text.) In all three cases above, the stationary
solution for $i$ and $r$ are given by

$\displaystyle\ i$ $\displaystyle=$ $\displaystyle\mu\frac{1-s}{\beta_{\rm
th}}$ (152) $\displaystyle r$ $\displaystyle=$
$\displaystyle\gamma\frac{1-s}{\beta_{\rm th}}$ (153)

Finally, we recover the steady state solution for $r$ from the sum rule above.

One absolute requirement that must be satisfied by the steady state solutions
for the sub-populations is that they are real and positive. In the case of the
second and third solutions listed above, reality of the population values
requires a positive argument in the square roots in (149) and (150). Coupled
with the demand for positive sub-populations we are led to leads to the
following inequality,

$\alpha>\frac{2\sqrt{\mu^{2}\beta_{\text{th}}^{3}\left(\beta_{\text{th}}-\beta_{0}\right)}-\mu\left(\beta_{0}-2\beta_{\text{th}}\right)\beta_{\text{th}}}{\mu^{2}}$
(154)

Next, we determine the stability of those solutions. At this point, we replace
coefficients by the representative values listed above. Expanding about the
three solutions, we find for the eigenvalues of the linear stability matrix of
the solutions,

  * •

Solution 1: $\lambda_{1}=-0.2,\lambda_{2}=-0.125$

  * •

Solution 2: $\lambda_{1}=-0.182396,\lambda_{2}=0.120832$

  * •

Solution 3: $\lambda_{1}=-0.166874+0.397813i,\lambda_{2}=-0.166874-0.397813i$

We see that the eigenvalues of the linear stability matrix for Solution 1,
which corresponds to an uninfected population, are all negative; the
uninfected state is stable with respect to the introduction of an infection,
assuming the initial rate of infection is sufficiently low. Additionally,
Solution 3, in which there is an endemic infection rate, is also stable with
respect to small perturbations, in that the real part of the eigenvalue is
negative. On the other hand in the case of Solution 2—in which the rate of
infection lies between 0 and the steady state value in Solution 3—one of the
eigenvalues of the linear stability matrix is positive, so it is dynamically
unstable.

We can explore the basins of attraction of the two stable stationary states by
solving the full set of spatially uniform equations, (76)–(79). Figures 8 and
9 show what we find for coefficients $\gamma$, $\mu$ and $\beta_{0}$ fixed
according to (147), and two different values of $\alpha$. Figure 8 displays
the basins of attraction when $\alpha=5$.

Figure 8: The basins of attraction for the SIR model with modified $\beta$ as
given by (146), coefficients $\gamma$, $\mu$ and $\beta_{0}$ fixed according
to (147) and $\alpha=5$. The axes are the initial values of $s(t)$ and $i(t)$,
which both range from 0 to 1 subject to the restriction that $s(0)+i(0)\leq
1$. The value of $r(0)$ is given by $r(0)=1-s(0)-i(0)$. In the blue region the
fixed point to which the system converges is endemically infected,
corresponding to solution 3 for $s$, with $i$ and $r$ given by (152) and
(153). In the red region the limiting values correspond to the uninfected
state: $s=1$, $i=r=0$.

Figure 9 displays the basins of attraction when $\alpha=3$, close to the
minimum value of 2.815 allowed by the requirement of reality and positivity of
the dynamically stable fixed points, as expressed in (154). As $\alpha$
approaches that value the blue basin of attraction continuously shrinks to
zero.

Figure 9: The basins of attraction with all coefficients the same as in Fig.
8, with the exception of $\alpha$, which is now 3. Again, the red region is
corresponds to an uninfected limiting state and the blue region to the limit
of an endemic infection.

###  VI.2 Two regions with different levels of infection in the SIR model with
an $I$-dependent $\beta$

The results in Section VI.1 point to the prospect of competition—and possibly
coexistence—between two regions with different levels of infection,
specifically an uninfected region and a region sustaining an endemic
infection.

To see how this can come about, we turn our focus to the full set of spatio-
temporal eqquations, (32)–(34) with the $I(x,t)$-dependent $\beta$ given by
(146), $I$ on the right hand side being replaced by $I(x,t)$. We can replace
$N(x,t)$ by a constant value, which we set equal to 1, and we have the sum
rule

$S(x,t)+I(x,t)+R(x,t)=1$ (155)

which means we need only concern ourselves with the first two equations.

Numerical simulations using these equations with parameters $\gamma$, $\mu$
and $\beta_{0}$ set according to (147) and $\alpha$ in the immediate vicinity
of 3 yield two regions separated by a compact front, as shown in Fig. 10.

Figure 10: Snapshot of the sub-populations in two regions. On the left the
population is infection free and on the right, there is infection at an
endemic level. The coefficient $\alpha$ in (146) has been set equal to 2.9,
and the front between the two regions propagates to the right, with a velocity
0.66776, so the system is evolving towards a uniform uninfected state.

The front can propagate in either direction, depending on the choice for
$\alpha$. For $\alpha<\alpha_{0}\approx 2.9267$ the front propagates to the
right, corresponding to an evolution to a uniform uninfected state, while if
$\alpha$ is greater than that transition value the propagation is to the left,
and the system ends up converting to a uniformly endemic state. The front
remains compact as $\alpha$ is assigned values in the immediate vicinity of
$\alpha_{0}$. When $\alpha=\alpha_{0}$, the front is stationary, and we have
perfect coexistence of the regions. When the fronts propagate, they do so in a
solitonic manner. This can be verified, and the dependence of the velocity on
the value of $\alpha$ determined, by testing the assumption that
$S(x)=S_{s}(x-vt)$ and similarly for $E(x)$, $I(x)$ and $R(x)$. Then the SIR
equations reduce to

$\displaystyle-v\frac{dS_{s}(x)}{dx}$ $\displaystyle=$ $\displaystyle
D\frac{d^{2}S_{s}(x)}{dx^{2}}+\mu-\mu
S_{s}(x)-\beta(I_{s}(x))S_{s}(x)I_{s}(x)$ (156)
$\displaystyle-v\frac{dI_{s}(x)}{dx}$ $\displaystyle=$ $\displaystyle
D\frac{d^{2}I_{s}(x)}{dx^{2}}+aE_{s}(x)-(\mu+\gamma)I_{s}(x)$ (157)
$\displaystyle-v\frac{dR_{s}(x)}{dx}$ $\displaystyle=$ $\displaystyle
D\frac{d^{2}R_{s}(x)}{dx^{2}}+\gamma E_{s}(x)-\mu R_{s}(x)$ (158)

where we have set the total population, $N_{s}(x)=1$. We then solve the
original set of spatio-temporal equations (32)–(34) numerically, take a
snapshot of the system, such as the one shown in Fig. 10, and check whether
there is a value of $v$ for which Eqs. (156)–(158) are all satisfied. To
satisfactory numerical tolerance (agreement to five significant figures in the
velocities) we find this to be the case for coefficients as given by (147)
with $\alpha$ varying from 2.85 to 3. The dependence of $v$ on $\alpha$ is
shown in Fig. 11.

Figure 11: The velocity $v(\alpha)$ of the propagation of the fronts, as shown
in Fig. 10, as a function of the coefficient $\alpha$ in (146). The value of
$\alpha$ for which the fronts do not propagate either to the left, as they do
for higher values of $\alpha$, or to the right, as they do for lower values,
is $\alpha_{0}\approx 2.927$.

###  VI.3 The case of the SEIR model with an $I$-dependent $\beta$

The analysis of an SEIR model in which the coefficient $\beta$ has the $I$
dependence (146) is similar, and in some cases identical, to the discussion
above. The homogeneous system admits three fixed points. The values of $s$ at
those fixed points are given by (148)–(150), the one change being $\beta_{\rm
th}$, used which in this case is

$\beta_{\rm th}=\frac{(\mu+a)(\mu+\gamma)}{a}$ (159)

Given the fixed point value of $s$, the populations $e$, $i$ and $r$ are

$\displaystyle e$ $\displaystyle=$ $\displaystyle\frac{\mu(1-s)}{\mu+a}$ (160)
$\displaystyle i$ $\displaystyle=$ $\displaystyle\frac{\mu(1-s)}{\beta_{\rm
th}}$ (161) $\displaystyle r$ $\displaystyle=$
$\displaystyle\frac{\gamma(1-s)}{\beta_{\rm th}}$ (162)

The stability properties of the three fixed points are the same as in the case
of the SIR model. The basins of attraction can be determined, with the
complication that they are now three dimensional manifolds. In order for the
system to admit of an uninfected stable fixed point, it is necessary to adjust
$\beta_{0}$ to be less than $\beta_{\rm th}$. With $\gamma=0.3$, $\mu=0.2$ and
$a=0.4$, $\beta_{\rm th}=0.75$. Choosing $\beta_{0}=0.675$ we find that for
$\alpha\approx 5.0533$ a state with endemic infection coexists with an
uninfected state, while for $\alpha$ slightly different than that value the
front between those states propagates and one of the states grows at the
expense of the other.

There are a number of interesting aspects to the possibility of a stationary
kink, since there may be real world examples in which two neighboring regions
may ”coexist” with disparate infection levels. In fact, the existence of a
stationary kink, which requires precise tuning of parameters ($\alpha$, for
fixed $\beta_{0},\gamma,\mu,a$), and an advancing front with even slight de-
tuning, is analogous to a first-order transition in thermodynamics. For
example, when the pressure $P$ is on the vapor pressure curve $P_{\sigma}(T)$
with $T<T_{c},$ the critical temperature, a liquid can coexist with its vapor
with any volume fractions for the two phases. In principle, with
$P-P_{\sigma}(T)=0\pm$, all of the material will convert to one phase or the
other in equilibrium occupying the system. If the pressure is adjusted
rapidly, the interface separating the liquid and its vapor (a kink solution)
will move one way or the other Cross and Hohenberg (1993) until one phase
fully occupies the system. This dynamical behavior can be seen in SEIR
compartmental models as discussed here, with $\alpha_{0}-\alpha_{c}$ playing
the role of the appropriate chemical potential difference (or magnetic field)
the context of condensed matter. Appendix B outlines the analysis of such a
model.

##  VII Summary

We have examined several extensions to compartment models that incorporate
spatial dependence of the relevant variables so as to model either individual
mobility or disease spreading through contact between individuals at different
geographical locations. All models considered here can be written as simple
reaction-diffusion systems, and we find that stable propagating front
solutions generically exist in certain ranges of model parameters that depend
on the specifics of each model. In that case, we also find the velocity of the
propagating front. These results potentially provide an additional avenue for
obtaining model parameters from empirical infection data, independently from
estimates of the various rates from single compartment data.

The various models studied have to be considered mean-field approximations in
which local statistical fluctuations in the subpopulations have been
neglected. This is only realistic in cases in which all the subpopulations are
large so that the various rates of interaction introduced in the models can be
meaningfully defined. While this may well be the case for some diseases, it is
expected to fail for the SARS-CoV-2 pandemic. The rates of infection are quite
low relative to the size of the population, and transmission and spread is
expected to be dominated by large fluctuations. As a consequence, one would
not expect, and indeed does not observe in the empirical data, compact fronts
separating regions at different stages of infection The New York Times (2022);
Li  _et al._ (2022). Current extensions of compartment models therefore
involve networks with a very large number of nodes to analyze the observed
effects of subpopulation heterogeneity. Whether such heterogeneity can be
modeled via stochastic extensions of the compartment models is currently under
study.

## References

  * Kermack and McKendrick (1927) W. O. Kermack and A. G. McKendrick, Proceedings of the Royal Society of London. Series A 115, 700 (1927).
  * Anderson and May (1979) R. M. Anderson and R. M. May, Nature 280, 361 (1979).
  * Murray (2001) J. Murray,  _Mathematical biology II: spatial models and biomedical applications_ , Vol. 3 (Springer-Verlag, 2001).
  * Li (2018) M. Y. Li,  _An introduction to mathematical modeling of infectious diseases_ , Vol. 2 (Springer, 2018).
  * Parolini  _et al._ (2021) N. Parolini, L. Dede’, P. F. Antonietti, G. Ardenghi, A. Manzoni, E. Miglio, A. Pugliese, M. Verani, and A. Quarteroni, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477, 20210027 (2021).
  * Yan  _et al._ (2021) H. Yan, Y. Zhu, J. Gu, Y. Huang, H. Sun, X. Zhang, Y. Wang, Y. Qiu, and S. Xi Chen, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477, 20200440 (2021), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2020.0440 .
  * Le Treut  _et al._ (2021) G. Le Treut, G. Huber, M. Kamb, K. Kawagoe, A. McGeever, J. Miller, R. Pnini, B. Veytsman, and D. Yllanes, medRxiv (2021).
  * Huang  _et al._ (2021) N. E. Huang, F. Qiao, Q. Wang, H. Qian, and K.-K. Tung, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477, 20210551 (2021), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2021.0551 .
  * Skellam (1951) J. G. Skellam, Biometrika 38, 196 (1951).
  * Fisher (1937) R. Fisher, Annals of Eugenics 7, 353 (1937).
  * Belgacem (2001) F. Belgacem, Int. J. Appl. Math. 5, 85 (2001).
  * Fife (2013) P. C. Fife,  _Mathematical aspects of reacting and diffusing systems_ , Vol. 28 (Springer Science & Business Media, 2013).
  * Hoppensteadt and Peskin (2012) F. C. Hoppensteadt and C. S. Peskin,  _Modeling and simulation in medicine and the life sciences_ , Vol. 10 (Springer Science & Business Media, 2012).
  * Postnikov and Sokolov (2007) E. B. Postnikov and I. M. Sokolov, Mathematical biosciences 208, 205 (2007).
  * Abdelheq  _et al._ (2019) M. Abdelheq, O. Belhamiti, L. Bouzid, D. Y. Trejos, and J. C. Valverde, Journal of Theoretical Biology 480, 192 (2019).
  * Noble (1974) J. V. Noble, Nature 250, 726 (1974).
  * Naether  _et al._ (2008) U. Naether, E. Postnikov, and I. Sokolov, The European Physical Journal B 65, 353 (2008).
  * Källén (1984) A. Källén, Nonlinear Analysis-Theory Methods & Applications 8, 851 (1984).
  * Rock  _et al._ (2014) K. Rock, S. Brand, J. Moir, and M. J. Keeling, Reports on Progress in Physics 77, 026602 (2014).
  * Verdasca  _et al._ (2005) J. Verdasca, M. Telo da Gama, A. Nunes, N. Bernardino, J. Pacheco, and M. Gomes, Journal of Theoretical Biology 233, 553 (2005).
  * Eksin  _et al._ (2021) C. Eksin, M. Ndeffo-Mbah, and J. S. Weitz, Journal of Theoretical Biology 520, 110632 (2021).
  * Viboud  _et al._ (2006) C. Viboud, O. N. Bjørnstad, D. L. Smith, L. Simonsen, M. A. Miller, and B. T. Grenfell, Science 312, 447 (2006).
  * Mossong  _et al._ (2008) J. Mossong, N. Hens, M. Jit, P. Beutels, K. Auranen, R. Mikolajczyk, M. Massari, S. Salmaso, G. S. Tomba, J. Wallinga, J. Heijne, M. Sadkowska-Todys, M. Rosinska, and W. J. Edmunds, PLOS Medicine 5, 1 (2008).
  * Holtz  _et al._ (2020) D. Holtz, M. Zhao, S. G. Benzell, C. Y. Cao, M. A. Rahimian, J. Yang, J. Allen, A. Collis, A. Moehring, T. Sowrirajan, _et al._ , Proceedings of the National Academy of Sciences 117, 19837 (2020).
  * Paltiel  _et al._ (2020) A. D. Paltiel, A. Zheng, and R. P. Walensky, JAMA Network Open 3, e2016818 (2020).
  * Saad-Roy  _et al._ (2020) C. M. Saad-Roy, C. E. Wagner, R. E. Baker, S. E. Morris, J. Farrar, A. L. Graham, S. A. Levin, M. J. Mina, C. J. E. Metcalf, and B. T. Grenfell, Science 370, 811 (2020).
  * Nande  _et al._ (2021) A. Nande, B. Adlam, J. Sheen, M. Z. Levy, and A. L. Hill, PLoS computational biology 17, e1008684 (2021).
  * Wong  _et al._ (2020) G. N. Wong, Z. J. Weiner, A. V. Tkachenko, A. Elbanna, S. Maslov, and N. Goldenfeld, Phys. Rev. X 10, 041033 (2020).
  * Wu  _et al._ (2018) X. Wu, B. Tian, and R. Yuan, Mathematical Methods in the Applied Sciences 41, 8136 (2018), https://onlinelibrary.wiley.com/doi/pdf/10.1002/mma.5273 .
  * Wu and Ma (2021) X. Wu and Z. Ma, Bound. Value Probl. 87 (2021), https://doi.org/10.1186/s13661-021-01564-z.
  * Kolmogorov  _et al._ (1991) A. Kolmogorov, I. Petrovskii, and N. Piskunov, in  _Selected Works of A. N. Kolmogorov I_ , edited by V. M. Tikhomirov (Kluwer, 1991) p. 248.
  * Brauer  _et al._ (2008) F. Brauer, P. Van den Driessche, J. Wu, and L. J. S. Allen,  _Mathematical epidemiology_, Mathematical biosciences subseries (Springer, Berlin, 2008) pp. xviii, 408 p.
  * Corless  _et al._ (1996) R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, Advances in Computational Mathematics 5, 329 (1996).
  * Fife and McLeod (1977) P. C. Fife and J. B. McLeod, Archive for Rational Mechanics and Analysis 65, 335 (1977).
  * Cross and Hohenberg (1993) M. C. Cross and P. C. Hohenberg, Reviews of Modern Physics 65, 851 (1993).
  * The New York Times (2022) The New York Times, (2022), Covid 19 data project. County level data times series at: https://github.com/nytimes/covid-19-data/blob/master/us-counties.csv.
  * Li  _et al._ (2022) X. Li, Q. Lin, and M. Kolak, “The US COVID-19 Atlas,” https://www.uscovidatlas.org (2022).
  * Landau (2008) L. Landau, Ukrainian Journal of Physics 53, 25 (2008).
  * Ginzburg (2004) V. L. Ginzburg, Physics-Uspekhi 47, 1155 (2004).
  * Coleman (1988) S. Coleman,  _Aspects of symmetry : selected Erice lectures of Sidney Coleman_ , 1st ed. (Cambridge University Press, Cambridge, England ; New York, 1988) pp. xiv, 402.
  * Jasnow  _et al._ (1979) D. Jasnow, T. Ohta, and J. Rudnick, Physical Review B 20, 2774 (1979).

##  Appendix A The FKPP model

The one dimensional Fisher-Kolmogorov–Petrovsky–Piskunov equation Fisher
(1937); Kolmogorov  _et al._ (1991) determines the evolution of a scalar field
$u(x,t)$ which obeys the equation,

$\frac{\partial u}{\partial t}=D\frac{\partial^{2}u}{\partial x^{2}}+ru(1-u),$
(163)

where both $D>0$ and $r>0$ are model parameters. The equation was originally
introduced to describe the spatial propagation of an advantageous gene, but
given the generality of the reaction term (the last term in the right hand
side of Eq. (163)), this equation has been used to describe front propagation
in many different contexts in Physics, Chemistry, and Biology, including
Ecology, Physiology, combustion, and solidification. The equation as written
has two fixed points $u=0$, and $u=1$. The fixed point $u=0$ is linearly
unstable, whereas the fixed point $u=1$ is linearly stable to perturbations of
the form $u(x,t)=1+\epsilon\cos(qx)$, $\epsilon\ll 1$ while
$\sigma=-r+Dq^{2}<0$. The equation has front solutions $u(x,t)=f(x-vt)$ that
satisfy the boundary conditions $u(x\rightarrow-\infty)=0$,
$u(x\rightarrow\infty)=1$ so that the front advances at constant speed
$v=2\sqrt{Dr}$.

##  Appendix B Motion of a one-dimensional interface separating two different
phases

We start with an equation that closely resembles the Fisher-
Kolmogorov–Petrovsky–Piskunov (FKPP) Fisher (1937); Kolmogorov  _et al._
(1991) equation in (163).

$\frac{\partial u}{\partial t}=D\frac{\partial^{2}u}{\partial
x^{2}}+ru(1-u^{2})$ (164)

This equation, unlike the FKPP equation, supports a stationary front. To see
that this is so, set the left hand side of (164) equal to zero, and insert the
solution

$u_{s}(x)=\tanh(B(x-x_{0}))$ (165)

into the new equation. Taking the second derivative of $u_{s}(x)$, we find

$\displaystyle D\frac{d^{2}u_{s}(x)}{dx^{2}}$ $\displaystyle=$
$\displaystyle-2DB^{2}\tanh(B(x-x_{0})){\mathop{\rm sech}}^{2}(B(x-x_{0}))$
(166) $\displaystyle=$ $\displaystyle-2DB^{2}u_{s}(x)(1-u_{s}(x)^{2})$

The stationary solution of (164)—i.e. $u$ is independent of $t$—is thus
satisfied if

$B=\sqrt{\frac{r}{2D}}$ (167)

There are two instructive physical interpretations of the stationary version
of Eq. (164). First, thinking of $u$ as a position variable and $x$ as a time
variable, the equation

$D\frac{d^{2}u(x)}{dx^{2}}=-ru(x)(1-u(x)^{2})$ (168)

is the equation of motion of a non-relativistic point particle with a
potential energy of the form

$V(u)=r\left(\frac{u^{2}}{2}-\frac{u^{4}}{4}\right)$ (169)

Second, thinking $u(x)$ as a scalar property, like the local magnetization
along an easy axis, of a system with spatial variation in one dimension, Eq.
(168) is the minimization equation for the free energy of a system, if that
free energy takes the form

$\displaystyle\int_{-\infty}^{\infty}\left[\frac{1}{2}\left(\frac{du(x)}{dx}\right)^{2}-r\left(\frac{u(x)^{2}}{2}-\frac{u(x)^{4}}{4}\right)\right]dx$
(170) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\left[\frac{1}{2}\left(\frac{du(x)}{dx}\right)^{2}-V(u(x))\right]dx$

This is an example of a Ginzburg-Landau-like free energy Landau (2008);
Ginzburg (2004). A graph of $V(u)$ is shown in Fig. 12.

Figure 12: The potential function $V(u)$ as given by (169), with $r=1$.

From the dynamical viewpoint, the solution (165), with (167), describes the
motion of a point mass that starts out at rest on the leftmost peak at $u=-1$
and moves through the valley, eventually coming to rest on the rightmost peak
at $u=1$. From the statistical mechanical viewpoint, the solution describes
the profile of a domain wall separating a region in which the system is
uniform with a negative magnetization corresponding to a minimum of $-V(u)$
and a region of positive magnetization corresponding to the other minimum of
that function.

Before proceeding to the case of a propagating front solution to Eq. (164),
let’s note the arbitrariness of the term $x_{0}$ in the solution (165), which
reflects the translational invariance of the solution. If we were to change
$x_{0}$ to $x_{0}+\delta x$. The solution, to first order in $\delta x$, is

$u_{s}(x)-u^{\prime}_{s}(x)\delta x$ (171)

This new function is still a solution to the stationary equation, which means
substituting (171) into (168) we find no term going as $\delta x$. The
function $u_{s}^{\prime}(x)$ is an example of what is known as a translation
mode Coleman (1988); Jasnow  _et al._ (1979).

Now, suppose we change the term $V(u)$ from its original form to $V(u)+\Delta
V(u)$. Anticipating that this change induces motion in the stationary front
and that to lowest order in $\Delta V(u)$ the moving front has the same
general form as the stationary one, we replace $u(x)$ with $u_{s}(x-vt)$ with
$v$ small. Substituting this in Eq. (164), expanding to first order in $\Delta
V(u)$ and $u$ we obtain

$-vu_{s}^{\prime}(x-vt)=-\Delta V^{\prime}(u_{s}(x-vt))-\delta V^{\prime}$
(172)

The term $\delta V^{\prime}$ consists of first order changes in
$V^{\prime}(u)$ in $u(x-vt)-u_{s}(x-vt)$. This term can be shown to be
mathematically orthogonal to $u^{\prime}_{s}(x-vt)$. Multiplying both sides of
this equation by $u_{s}^{\prime}(x-vt)$ and integrating we find

$\displaystyle v$ $\displaystyle=$
$\displaystyle\frac{\int_{-\infty}^{\infty}\Delta
V^{\prime}(u_{s}(x-vt))u_{s}^{\prime}(x-vt)dx}{\int_{-\infty}^{\infty}\phi_{s}^{\prime\,2}(x-vt)dx}$
(173) $\displaystyle=$ $\displaystyle\frac{\Delta V(x=\infty)-\Delta
V(x=-\infty)}{\int_{-\infty}^{\infty}\phi_{s}^{\prime\,2}(x-vt)dx}$

The velocity of propagation of the front depends on the difference between the
two limiting values of the change in the local free energy term.