# The contact process with an asymptomatic state

Nicolas Lanchier and Max Mercer Nicolas Lanchier was partially supported by
NSF grant CNS-2000792.

###### Abstract

In order to understand the role of a potentially high infectiousness of
symptomatic individuals or, on the contrary, social distancing, quarantine,
etc. in the course of an infectious disease, this paper introduces a natural
variant of the popular contact process that distinguishes between asymptomatic
and symptomatic individuals. Infected individuals all recover at rate one but
infect nearby individuals at a rate that depends on whether they show the
symptoms of the disease or not. Newly infected individuals are always
asymptomatic and may or may not show the symptoms before they recover. The
analysis of the corresponding mean-field model reveals that, in the absence of
local interactions, regardless of the rate at which asymptomatic individuals
become symptomatic, there is an epidemic whenever at least one of the
infection rates is sufficiently large. In contrast, our analysis of the
interacting particle system shows that, when the rate at which asymptomatic
individuals become symptomatic and the rate at which asymptomatic individuals
infect nearby individuals are small, there cannot be any epidemic even when
symptomatic individuals are highly infectious.

††footnotetext: AMS 2000 subject classifications: Primary 60K35 Keywords and
phrases: Contact process, coupling, block construction, COVID-19, influenza,
HIV.

##  1 Introduction

Infectious diseases commonly spread through the transfer of bacteria from one
person to another either directly through talking, coughing, sneezing, etc. or
indirectly through touching and subsequently infecting objects which are then
touched by another person. This process is referred to as viral shedding, in
which diseased cells create copies of themselves that leave the body and
infect other individuals. It is also common that infected individuals do not
immediately exhibit the symptoms of the disease, which leads to two categories
of infected individuals that we shall call asymptomatic individuals (the ones
who do not exhibit the symptoms yet) and symptomatic individuals (the ones who
exhibit the symptoms of the disease). Typically, asymptomatic individuals are
less contagious but unaware that they have contracted the disease whereas
symptomatic individuals are more contagious but can take precautions (use of a
mask, social distancing, quarantine, use of a condom, etc.) to limit the
spread of the disease. In particular, the characteristics of the disease and
the social behavior of the population induce a variability in the rate at
which the disease spreads from asymptomatic versus symptomatic individuals.
COVID-19. A recent example of an infectious disease is COVID-19. Upon
infection, individuals experience a period during which they do not show any
symptoms of the virus, and it is estimated that the mean incubation period is
8.29 days [19]. It is even possible for individuals to spend their entire
disease cycle in the asymptomatic stage in which case they might never know
that they have had contracted the disease at some point. Similarly,
individuals may never know about the time at which they got infected. However,
just because some individuals are asymptomatic does not mean that they are not
contagious, and some studies report that most COVID-19 cases have been spread
from asymptomatic people: according to the study [12], 59% of all
transmissions came from asymptomatic individuals, including 35% from
presymptomatic individuals and 24% from individuals who never developed any
symptoms while, according to the study [11], the estimated proportion of
presymptomatic transmissions was 44%. In contrast, due to the enforcement of
quarantines, symptomatic individuals cannot easily spread the disease, showing
that the presymptomatic spread is far from being a negligible factor in the
case of COVID-19.
Influenza. Another example of an infectious disease is the influenza virus
(flu). In humans, it is known that viral shedding is greater in symptomatic
individuals than in asymptomatic individuals according to [18] which cites two
studies [2, 3] that established this result. In this case, the rate of
symptomatic spread is higher than the rate of asymptomatic spread. This is in
sharp contrast with the case of COVID-19 discussed above in which the
estimated mean incubation period and the smaller estimated proportion of
presymptomatic transmissions along with the result that most adults have
replication-competent SARS-CoV-2 for 10–14 days after symptom onset [5]
suggest that there is a higher rate of asymptomatic spread than symptomatic
spread.
HIV. Another important example is the human immunodeficiency virus (HIV) which
has been studied in the context of interacting particle systems before for its
interesting dynamics. HIV is a sexually-transmitted disease that can spread
during both the asymptomatic and the symptomatic stages of the infection.
Individuals with symptomatic HIV infection, especially those with advanced
disease or opportunistic infections, may have higher levels of the virus in
their body fluids, which can increase the risk of transmission [17]. However,
as people become aware of their infection, they may take steps to avoid
spreading the infection to others they interact with. In any case, the rate of
asymptomatic spread may differ significantly from the rate of symptomatic
spread.

To understand the effects of a potentially high infectiousness of symptomatic
individuals due to an increased viral shedding or on the contrary the benefit
of social distancing, quarantine, etc., we study a natural variant of Harris’
contact process [10] in which infected individuals are divided into two
categories: asymptomatic individuals and symptomatic individuals. Like in the
contact process, individuals are located on the $d$-dimensional integer
lattice, and infected individuals recover at rate one and infect their healthy
neighbors by contact. However, the rate at which infected individuals infect
their neighbors now depends on whether they are asymptomatic or symptomatic,
with newly infected individuals being always asymptomatic before possibly
transitioning to the symptomatic state. More precisely, each site of the
integer lattice is in state

$0=\hbox{healthy}\qquad 1=\hbox{infected/asymptomatic}\qquad\hbox{or}\qquad
2=\hbox{infected/symptomatic}$

so the state of the enire system at time $t$ is a spatial configuration

$\xi_{t}:\mathbb{Z}^{d}\longrightarrow\\{0,1,2\\}\quad\hbox{where}\quad\xi_{t}(x)=\hbox{state
at site\leavevmode\nobreak\ $x$ at time\leavevmode\nobreak\ $t$}.$

We assume for simplicity that individuals can only interact with their nearest
neighbors but most of our results easily extend to larger ranges of
interactions. To describe the dynamics (local interactions) in the nearest
neighbor case, we let

$N_{x}=\\{y\in\mathbb{Z}^{d}:|\\!|x-y|\\!|=1\\}\quad\hbox{for all}\quad
x\in\mathbb{Z}^{d},$

where $|\\!|\cdot|\\!|$ refers to the Euclidean norm, be the set of nearest
neighbors of site $x$. Then, letting for each type $i$ of infected
individuals, each site $x$ and each configuration $\xi$,

$\begin{array}[]{l}f_{i}(x,\xi)=(1/2d)\,\sum_{y\in
N_{x}}\mathbf{1}\\{\xi(y)=i\\}\end{array}$

be the fraction of neighbors of site $x$ that are in state $i$, the contact
process with an asymptomatic state evolves according to the local transition
rates

$\begin{array}[]{rclcrcl}0\to 1&\hbox{at
rate}&\beta_{1}f_{1}(x,\xi)+\beta_{2}f_{2}(x,\xi)&&1\to 2&\hbox{at
rate}&\gamma\vspace*{4pt}\\\ 1\to 0&\hbox{at rate}&1&&2\to 0&\hbox{at
rate}&1.\end{array}$

The first transition rate indicates that healthy individuals get infected from
nearby asymptomatic individuals at rate $\beta_{1}$ and from nearby
symptomatic individuals at rate $\beta_{2}$. Though the infection rate
$\beta_{2}$ is expected to be larger than $\beta_{1}$, factors such as social
distancing and quarantine may contribute to make $\beta_{2}$ smaller than the
infection rate $\beta_{1}$. The two transition rates at the bottom indicate
that infected individuals recover at rate one regardless of whether they are
symptomatic or not. Finally, the second transition rate at the top indicates
that asymptomatic individuals may become symptomatic. For small $\gamma$, most
of the infected individuals will recover before showing any symptoms and may
be unaware that they were ever infected whereas, for large $\gamma$, most of
the infected individuals will quickly show the symptoms of the disease.
Assuming that the population is homogeneously mixing, and letting $u_{1}$ be
the density of infected individuals who are asymptomatic and $u_{2}$ be the
density of infected individuals who are symptomatic, the system is described
by the mean-field model

$u_{1}^{\prime}=(\beta_{1}u_{1}+\beta_{2}u_{2})(1-u_{1}-u_{2})-(1+\gamma)u_{1}\quad\hbox{and}\quad
u_{2}^{\prime}=\gamma u_{1}-u_{2}.$

The analysis of this system shows that, in addition to the disease-free
equilibrium $p_{0}=(0,0)$, there is a nontrivial fixed point $p_{12}$ if and
only if $\beta_{1}+\gamma\beta_{2}>1+\gamma$, in which case this fixed point
is globally stable. In particular, for all $0<\gamma<\infty$, even if one of
the two infection rates is equal to zero, there is an epidemic whenever the
other infection rate is sufficiently large. Our analysis of the spatial model
again shows that, for all $0<\gamma<\infty$ and all $\beta_{2}$, there is an
epidemic provided the infection rate $\beta_{1}$ is sufficiently large.
However, when both $\gamma$ and $\beta_{1}$ are small, the population always
recovers regardless of the infection rate $\beta_{2}$, which strongly
contrasts with the behavior of the mean-field model. Starting with a single
infection at 0, we say that the process

$\begin{array}[]{rcl}\hbox{survives strongly}&\hbox{when}&P(\xi_{t}(0)=1\
\hbox{infinitely often})>0\vspace*{4pt}\\\
\hbox{survives}&\hbox{when}&P(\xi_{t}\not\equiv 0\ \hbox{for all times}\
t)>0\end{array}$

and goes extinct when it does not survive. Intuitively, the larger the
asymptomatic infection rate $\beta_{1}$ and/or the larger the symptomatic
infection rate $\beta_{2}$, the larger the set of infected sites:

$\begin{array}[]{l}\hbox{(a) \ \ the set of infected sites is nondecreasing
in\leavevmode\nobreak\ $\beta_{1}$}\vspace*{4pt}\\\ \hbox{(b) \ \ the set of
infected sites is nondecreasing in\leavevmode\nobreak\
$\beta_{2}$}.\end{array}$ (1)

In addition, because increasing $\gamma$ decreases the number of asymptomatic
individuals but increases the number of symptomatic, we also expect that

$\begin{array}[]{rcl}\hbox{(a)}\ \
\beta_{1}<\beta_{2}&\Longrightarrow&\hbox{the set of infected sites is
nondecreasing in\leavevmode\nobreak\ $\gamma$}\vspace*{4pt}\\\ \hbox{(b)}\ \
\beta_{1}>\beta_{2}&\Longrightarrow&\hbox{the set of infected sites is
nonincreasing in\leavevmode\nobreak\ $\gamma$}\end{array}$ (2)

and we refer to the simulation pictures in Figure 1 for an illustration of
this result.

$\beta_{1}=8$, $\beta_{2}=0$, $\gamma=1$$\beta_{1}=8$, $\beta_{2}=0$,
$\gamma=2$$\beta_{1}=0$, $\beta_{2}=4$, $\gamma=1$$\beta_{1}=0$,
$\beta_{2}=4$, $\gamma=2$

Figure 1:  Snapshots of the process on the $300\times 300$ torus at time 200.
In all the pictures, the black sites represent the individuals who are
infected but asymptomatic, the gray sites the individuals who are infected and
symptomatic, and the white sites the individuals who are healthy. The
simulations show the reversed monotonicity of the parameter $\gamma$ depending
on the order of the two infection rates $\beta_{1}$ and $\beta_{2}$.

Using a graphical representation to couple processes with different
parameters, we will prove that the two monotonicity results (1) and (2) indeed
hold when $\beta_{1}<\beta_{2}$. However, though we believe (and conjecture)
that the probability of survival of the infection is monotone with respect to
each of the parameters, the same coupling argument fails to prove this result
when $\beta_{1}>\beta_{2}$. Monotonicity implies that, two of the parameters
being fixed, there exists at most one phase transition (between extinction and
survival) in the direction of the third parameter. To prove the existence of
phase transitions, we first compare the process with the basic contact process
[10], the interacting particle system in which sites can be 0 = healthy or 1 =
infected with transition rates

$\begin{array}[]{rclcrcl}0\to 1&\hbox{at rate}&\beta
f_{1}(x,\xi)&\quad\hbox{and}&1\to 0&\hbox{at rate}&1.\end{array}$

The main result about the contact process states that there exists a
nondegenerate (positive and finite) critical value $\beta_{c}$ that depends on
the spatial dimension such that

$\begin{array}[]{rcl}\hbox{the process survives
strongly}&\Longleftrightarrow&\beta>\beta_{c}.\end{array}$ (3)

Bezuidenhout and Grimmett [1] also proved that the contact process dies out at
$\beta=\beta_{c}$. For a detailed review of the process, we refer to Liggett
[15, 16]. Note that the contact process with an asymptomatic state can be
viewed as a mixture of the contact process with parameter $\beta_{1}$ and the
contact process with parameter $\beta_{2}$, with the relative weight of the
second contact process being measured by the rate $\gamma$ at which
asymptomatic individuals become symptomatic. Observing that the rate at which
infected individuals infect their neighbors is always between

$\min(\beta_{1},\beta_{2})\quad\hbox{and}\quad\max(\beta_{1},\beta_{2})$

and that the rate at which they recover is equal to one regardless of their
state (asymptomatic or symptomatic), it follows from (3) and a standard
coupling argument that

$\begin{array}[]{rcl}\hbox{(a)}\ \
\beta_{1}\leq\beta_{c}\quad\hbox{and}\quad\beta_{2}\leq\beta_{c}&\Longrightarrow&\hbox{the
infection dies out}\vspace*{4pt}\\\ \hbox{(b)}\ \
\beta_{1}>\beta_{c}\quad\hbox{and}\quad\beta_{2}>\beta_{c}&\Longrightarrow&\hbox{the
infection survives strongly}.\end{array}$ (4)

We now look at the parameter regions where one of the two infection rates is
subcritical and the other one supercritical in which case the behavior of the
process is much less obvious.

Asymptomatic = subcritical and symptomatic = supercritical
To begin with, note that taking $\gamma=0$ blocks the access to the
symptomatic state, so the process reduces to the basic contact process with
parameter $\beta_{1}$. On the contrary, taking $\gamma=\infty$ makes the
transition $1\to 2$ instantaneous, thus removing the asymptomatic state,
therefore the process reduces to the basic contact process with parameter
$\beta_{2}$. Using these observations together with block constructions and
perturbation techniques, we obtain the following theorem.

###### Theorem 1

Starting from a translation invariant distribution with infinitely many 1s and
2s,

$\begin{array}[]{rcl}\beta_{1}<\beta_{c}\quad\hbox{and}\quad\beta_{2}>\beta_{c}&\Longrightarrow&\left\\{\begin{array}[]{l}\hbox{(a)
strong survival for all\leavevmode\nobreak\ $\gamma<\infty$
large}\vspace*{4pt}\\\ \hbox{(b) extinction for all\leavevmode\nobreak\
$\gamma>0$ small}.\end{array}\right.\end{array}$ (5)

To study the process for intermediate values of the parameter $\gamma$, we now
look at the $\beta_{1}=0$ case in which asymptomatic individuals are not
contagious. In this extreme case, the transition rates of the contact process
with an asymptomatic state reduce to

$\begin{array}[]{rclcrcl}0\to 1&\hbox{at rate}&\beta_{2}f_{2}(x,\xi)&&1\to
2&\hbox{at rate}&\gamma\vspace*{4pt}\\\ 1\to 0&\hbox{at rate}&1&&2\to
0&\hbox{at rate}&1.\end{array}$ (6)

This process is a particular case of Krone’s two-stage contact process [14].
In the demographic context of the two-stage contact process, the states are
interpreted as

$0=\hbox{empty}\qquad 1=\hbox{occupied by a child}\qquad\hbox{and}\qquad
2=\hbox{occupied by an adult}.$

The transition rates (6) indicate that only adults can give birth, and when
they do, they give birth to a child. In addition, children become adults at a
maturation rate $\gamma$. In the original model, children die at a possibly
higher rate than the adults. Using an idea inspired from Foxall [7] to compare
the process with a subcritical branching process, we will prove that

###### Theorem 2

Starting from a translation invariant distribution with infinitely many 1s and
2s,

$\begin{array}[]{rcl}\gamma<1/(2d-1)\quad\hbox{and}\quad\beta_{1}=0&\Longrightarrow&\hbox{extinction
for all \,$\beta_{2}\geq 0$}.\end{array}$

Letting $\beta_{1}=0$ and $\beta_{+}>\beta_{c}$ be fixed, it follows from
Theorem 1 and the monotonicity of the set of infected sites with respect to
the parameters that there exists $\gamma_{+}=\gamma_{+}(\beta_{+})<\infty$
such that there is an epidemic whenever $\gamma>\gamma_{+}$ and
$\beta_{2}\geq\beta_{+}$. In contrast, it follows from Theorem 2 that there
exists $\gamma_{-}>0$ such that the infection dies out when
$\gamma<\gamma_{-}$ and $\beta_{1}=0$ regardless of the value of the
symptomatic infection rate. Using again the monotonicity of the set of
infected sites with respect to the parameter $\gamma$ when
$\beta_{1}<\beta_{2}$, we deduce the existence of a critical value

$0<\gamma_{-}\leq\gamma_{c}\leq\gamma_{+}=\gamma_{+}(\beta_{+})<\infty$

such that the following two implications hold:

$\begin{array}[]{rcl}\hbox{(a)}\ \
\beta_{1}=0\quad\hbox{and}\quad\gamma>\gamma_{c}&\Longrightarrow&\hbox{survival
for\leavevmode\nobreak\ $\beta_{2}<\infty$ large}\vspace*{4pt}\\\ \hbox{(b)}\
\
\beta_{1}=0\quad\hbox{and}\quad\gamma<\gamma_{c}&\Longrightarrow&\hbox{extinction
even when\leavevmode\nobreak\ $\beta_{2}=\infty$}.\end{array}$ (7)

This is in sharp contrast with the behavior of the mean-field model in which,
provided $\gamma>0$, there is an epidemic for $\beta_{2}$ sufficiently large
even when $\beta_{1}=0$.

Asymptomatic = supercritical and symptomatic = subcritical
Recall that in this case the coupling arguments used when
$\beta_{1}<\beta_{2}$ fail to prove monotonicity of the set of infected sites
with respect to each of the three parameters. However, using some obvious
symmetry and exchanging the role of the two infection rates also accounting
for the reverse effect of the parameter $\gamma$, the proof of Theorem 1
easily extends to obtain

###### Theorem 3

Starting with infinitely many infected sites,

$\begin{array}[]{rcl}\beta_{1}>\beta_{c}\quad\hbox{and}\quad\beta_{2}<\beta_{c}&\Longrightarrow&\left\\{\begin{array}[]{l}\hbox{(a)
strong survival for all\leavevmode\nobreak\ $\gamma>0$ small}\vspace*{4pt}\\\
\hbox{(b) extinction for all\leavevmode\nobreak\ $\gamma<\infty$
large}.\end{array}\right.\end{array}$ (8)

To study the behavior of the process for intermediate values of $\gamma$, we
now assume $\beta_{2}=0$, meaning that the symptomatic individuals cannot
spread the disease, which happens when all the symptomatic individuals are
placed in quarantine. Assuming in addition that an infected individual who is
asymptomatic cannot recover but turns instead into a symptomatic individual
results in the three-state interacting particle system with local transition
rates

$\begin{array}[]{rclcrclcrcl}0\to 1&\hbox{at rate}&\beta_{1}f_{1}(x,\xi)&&1\to
2&\hbox{at rate}&1+\gamma&&2\to 0&\hbox{at rate}&1.\end{array}$ (9)

This stochastic process is a time change of the standard forest fire model in
which time is sped up by the factor $1+\gamma$. In this context, the states
are interpreted as

$0=\hbox{live tree}\qquad 1=\hbox{burning tree}\qquad\hbox{and}\qquad
2=\hbox{burnt tree}.$

The transition rates (9) indicate that the fire spread from burning trees to
nearby live trees, trees keep burning for an exponential amount of time, and
burnt trees are eventually replaced by a new live tree. Durrett and Neuhauser
[6] proved that, in $d=2$, the fire keeps burning when the rate at which
burning trees emit sparks exceeds some critical value $\alpha_{c}$. In
addition, because the transition $1\to 0$ occurring at rate one in our process
is replaced by the transition $1\to 2$ occurring at the same rate in the
forest fire model, we expect that our epidemic model has more 1s and less 2s
than the forest fire model which, accounting for the time change and using
[6], suggests that there is an epidemic when
$\beta_{1}>(1+\gamma)\,\alpha_{c}$. Though we believe (and conjecture) that
this result is true, standard coupling arguments again fail to compare both
processes. However, using a construction from Cox and Durrett [4], we can
compare the set of sites that ever get infected with the set of wet sites in a
certain site percolation process to prove the following.

###### Theorem 4

Assuming that $d>1$ and starting the process with a single 1 at the origin in
an otherwise healthy population,

$\begin{array}[]{rcl}\gamma<\infty\quad\hbox{and}\quad\beta_{2}=0&\Longrightarrow&\hbox{the
infection survives for all \,$\beta_{1}<\infty$ large}.\end{array}$ (10)

Contrary to Theorem 2, this result is in agreement with the behavior of the
mean-field model in which there is an epidemic whenever one infection rate is
sufficiently large even if the other infection rate is equal to zero.

$\beta_{1}$$\beta_{2}$epidemic(4.b)$\infty$(10)$\gamma>\gamma_{c}$(7.a)11(4.a)$\gamma=\infty$$\beta_{c}$$\beta_{c}$(10)$\gamma<\gamma_{c}$(7.b)

$\gamma=0$

$\displaystyle\frac{1}{\gamma}+1$$1+\gamma$recovery(5.a)(8.a)(5.b)(8.b)

Figure 2:  Phase structure of the contact process with an asymptomatic state
and its mean-field model. The thick black lines show the phase transitions of
the interacting particle system for various values of $\gamma$. The solid
vertical and horizontal lines show the phase transitions in the limiting cases
$\gamma=0$ and $\gamma=\infty$. The dashed lines are the straight lines with
equation $\beta_{1}+\gamma\beta_{2}=1+\gamma$ showing the corresponding phase
transitions for the mean-field model.

This shows that, while the asymptomatic and the symptomatic infection rates
play similar roles in the mean-field model, they have significantly different
effects in the spatial model. This is due to the presence of local
interactions. In the case $\beta_{1}=0$, meaning that the 1s are not
contagious, the 2s are quickly surrounded and blocked by the nearby 1s they
create, and if in addition these 1s take too much time to turn into 2s, the
infection dies out. In contrast, when $\beta_{2}=0$, meaning that the 2s
cannot spread the disease, because the 1s produce 1s, even when the 1s quickly
turn into 2s, if these 1s spread fast enough they can still escape from the 2s
and form traveling waves drifting off to infinity. Figure 2 shows a picture of
the phase structure of the process and its mean-field model summarizing these
results.
The rest of the paper is devoted to the proofs. Section 2 focuses on the mean-
field model, using in particular Bendixson–Dulac theorem to give a complete
description of the global stability of the fixed points. Section 3 explains
how to construct the process from a graphical representation and couple
processes with different parameters to prove monotonicity when
$\beta_{1}<\beta_{2}$. The last four sections give the proofs of the four
theorems. Sections 4 and 5 use block constructions and perturbation techniques
to compare the process when $\gamma$ is either very small or very large with a
supercritical (Section 4) and a subcritical (Section 5) contact process and
deduce Theorems 1 and 3. Section 6 relies on a comparison of the process with
a subcritical Galton–Watson branching process to prove extinction, as stated
in Theorem 2. Finally, the proof of Theorem 4 is carried out in Section 7 by
showing that the set of sites that ever get infected dominates stochastically
a certain site percolation process on the integer lattice in dimension $d>1$.

##  2 Mean-field analysis

This section is devoted to the analysis of the mean-field model (existence,
local stability and global stability of the fixed points) described above.
Recall that, assuming that the population is homogeneously mixing, and letting
$u_{i}$ be the density of individuals in state $i$, the mean-field model
consists of the system of coupled differential equations

$\begin{array}[]{rcl}u_{1}^{\prime}=F_{1}(u_{1},u_{2})&=&(\beta_{1}u_{1}+\beta_{2}u_{2})(1-u_{1}-u_{2})-(1+\gamma)u_{1}\vspace*{4pt}\\\
u_{2}^{\prime}=F_{2}(u_{1},u_{2})&=&\gamma u_{1}-u_{2}.\end{array}$ (11)

Setting the right-hand side of both equations equal to zero, some basic
algebra shows that, in addition to the disease-free equilibrium $p_{0}=(0,0)$
which is always a fixed point, there is also a nontrivial fixed point
$p_{12}=(\bar{u}_{1},\bar{u}_{2})$ whose coordinates are given by

$\bar{u}_{1}=\frac{1}{1+\gamma}-\frac{1}{\beta_{1}+\gamma\beta_{2}}\quad\hbox{and}\quad\bar{u}_{2}=\gamma\bar{u}_{1}=\frac{\gamma}{1+\gamma}-\frac{\gamma}{\beta_{1}+\gamma\beta_{2}}.$
(12)

This fixed point is biologically relevant in the sense that it belongs to the
two-dimensional simplex if and only if $\beta_{1}+\gamma\beta_{2}>1+\gamma$.
To study the local stability of the fixed points, notice first that the
Jacobian matrix associated to the system (11) is given by

$J(u_{1},u_{2})=\left(\begin{array}[]{ccc}\beta_{1}(1-2u_{1})-(\beta_{1}+\beta_{2})u_{2}-(1+\gamma)&&\beta_{2}(1-2u_{2})-(\beta_{1}+\beta_{2})u_{1}\vspace*{2pt}\\\
\gamma&&-1\end{array}\right).$

Taking $u_{1}=u_{2}=0$, the disease-free equilibrium is locally stable if and
only if

$\operatorname{Tr}\,(J(0,0))=\beta_{1}-\gamma-2<0\quad\hbox{and}\quad\operatorname{Det}\,(J(0,0))=1+\gamma-\beta_{1}-\gamma\beta_{2}>0,$

which is equivalent to $\beta_{1}+\gamma\beta_{2}<1+\gamma$. To study the
local stability of the nontrivial fixed point, we let $D_{1}$ and $D_{2}$ be
the two denominators in the expression (12) of the fixed point:

$D_{1}=1+\gamma\quad\hbox{and}\quad D_{2}=\beta_{1}+\gamma\beta_{2}.$

Multiplying the first coefficient of $J(\bar{u}_{1},\bar{u}_{2})$ by
$D_{1}D_{2}$, we get

$\begin{array}[]{l}D_{1}D_{2}(\beta_{1}(1-2\bar{u}_{1})-(\beta_{1}+\beta_{2})\bar{u}_{2}-(1+\gamma))\vspace*{4pt}\\\
\hskip
40.0pt=D_{1}D_{2}(\beta_{1}-(1+\gamma)-((\beta_{1}+\gamma\beta_{2})+(1+\gamma)\beta_{1})\bar{u}_{1})\vspace*{4pt}\\\
\hskip
40.0pt=D_{1}D_{2}(\beta_{1}-D_{1}-(D_{2}+D_{1}\beta_{1})(1/D_{1}-1/D_{2}))\vspace*{4pt}\\\
\hskip
40.0pt=D_{1}D_{2}(\beta_{1}-D_{1})-(D_{2}+D_{1}\beta_{1})(D_{2}-D_{1}),\end{array}$
(13)

while multiplying the second coefficient of the matrix by $D_{1}D_{2}$, we get

$\begin{array}[]{l}D_{1}D_{2}(\beta_{2}(1-2\bar{u}_{2})-(\beta_{1}+\beta_{2})\bar{u}_{1})\vspace*{4pt}\\\
\hskip
40.0pt=D_{1}D_{2}(\beta_{2}-((\beta_{1}+\gamma\beta_{2})+(1+\gamma)\beta_{2})\bar{u}_{1})\vspace*{4pt}\\\
\hskip
40.0pt=D_{1}D_{2}(\beta_{2}-(D_{2}+D_{1}\beta_{2})(1/D_{1}-1/D_{2}))\vspace*{4pt}\\\
\hskip
40.0pt=D_{1}D_{2}\beta_{2}-(D_{2}+D_{1}\beta_{2})(D_{2}-D_{1}).\end{array}$
(14)

Using (13), we deduce that the trace is negative:

$\begin{array}[]{rcl}D_{1}D_{2}\,\operatorname{Tr}\,(J(\bar{u}_{1},\bar{u}_{2}))&=&D_{1}D_{2}(\beta_{1}-D_{1})-(D_{2}+D_{1}\beta_{1})(D_{2}-D_{1})-D_{1}D_{2}\vspace*{4pt}\\\
&=&-D_{1}^{2}D_{2}-D_{2}^{2}+D_{1}^{2}\beta_{1}=-D_{2}^{2}-D_{1}^{2}\gamma\beta_{2}<0,\end{array}$

which shows that the Jacobian matrix has at least one negative eigenvalue.
Similarly, using both expressions (13) and (14), we deduce that the
determinant is positive:

$\begin{array}[]{l}D_{1}D_{2}\,\operatorname{Det}\,(J(\bar{u}_{1},\bar{u}_{2}))\vspace*{4pt}\\\
\hskip
20.0pt=-D_{1}D_{2}(\beta_{1}-D_{1})+(D_{2}+D_{1}\beta_{1})(D_{2}-D_{1})-\gamma
D_{1}D_{2}\beta_{2}+\gamma(D_{2}+D_{1}\beta_{2})(D_{2}-D_{1})\vspace*{4pt}\\\
\hskip
20.0pt=-D_{1}D_{2}(\beta_{1}+\gamma\beta_{2}-D_{1})+((1+\gamma)D_{2}+D_{1}(\beta_{1}+\gamma\beta_{2}))(D_{2}-D_{1})\vspace*{4pt}\\\
\hskip
20.0pt=-D_{1}D_{2}(D_{2}-D_{1})+(D_{1}D_{2}+D_{1}D_{2})(D_{2}-D_{1})=D_{1}D_{2}(D_{2}-D_{1})>0\end{array}$

when $D_{2}=\beta_{1}+\gamma\beta_{2}>1+\gamma=D_{1}$, which is exactly the
condition for the nontrivial fixed point to belong to the two-dimensional
simplex. This shows that both eigenvalues have the same sign therefore that
both eigenvalues are negative and the interior fixed point $p_{12}$ is locally
stable.

$\beta_{1}=2,\ \beta_{2}=3,\ \gamma=2$$\beta_{1}=2,\ \beta_{2}=3,\ \gamma=0.5$

Figure 3:  Solution curves of the mean-field model for different values of
$\gamma$.

Finally, to prove that the local stability of the fixed points under the
appropriate conditions implies their global stability, the last step is to
apply Bendixson–Dulac theorem to exclude the existence of periodic orbits.
Using the Dulac function $\phi=1/(u_{1}u_{2})$, we obtain

$\begin{array}[]{rcl}\nabla\cdot(\phi F_{1},\phi
F_{2})&=&\displaystyle\frac{\partial}{\partial
u_{1}}\bigg{(}\bigg{(}\frac{\beta_{1}}{u_{2}}+\frac{\beta_{2}}{u_{1}}\bigg{)}(1-u_{1}-u_{2})-\frac{1+\gamma}{u_{2}}\bigg{)}+\displaystyle\frac{\partial}{\partial
u_{2}}\bigg{(}\frac{\gamma}{u_{2}}-\frac{1}{u_{1}}\bigg{)}\vspace*{8pt}\\\
&=&\displaystyle-\bigg{(}\frac{\beta_{2}(1-u_{1}-u_{2})}{u_{1}^{2}}\bigg{)}-\bigg{(}\frac{\beta_{1}}{u_{2}}+\frac{\beta_{2}}{u_{1}}\bigg{)}-\frac{\gamma}{u_{2}^{2}}<0\end{array}$

for all $(u_{1},u_{2})$ in the interior of the simplex. In conclusion,

$\begin{array}[]{rcl}\hbox{there is an
epidemic}&\Longleftrightarrow&\beta_{1}+\gamma\beta_{2}>1+\gamma\end{array}$

where epidemic in the context of the mean-field model means that, starting
from any initial condition in the simplex other than $p_{0}$, the system
converges to the interior fixed point $p_{12}$. In contrast, when the
inequality is flipped, the trivial fixed point $p_{0}$ is the only fixed point
in the simplex and, starting from any initial condition in the simplex, the
system converges to $p_{0}$. Figure 3 shows a picture of the solution curves
when there is an epidemic while Figure 2 shows a visualization of the
parameter region $\beta_{1}+\gamma\beta_{2}>1+\gamma$ (above the oblique
dashed line when $0<\gamma<\infty$).

##  3 Graphical representation and monotonicity

The first step to study the spatial model is to use a key idea due to Harris
[9] that consists in constructing the process graphically from a collection of
independent Poisson processes/exponential clocks. The concept of graphical
representation not only shows that the process on the infinite lattice is
well-defined (even though the time to the first/next update does not exist)
but also will be used later to couple processes with different parameters in
order to prove monotonicity results, and apply block
constructions/perturbation techniques to compare the process properly rescaled
in space and time with oriented site percolation. In the case of the contact
process with an asymptomatic state, the most natural approach to construct the
process is to use four collections of exponential clocks corresponding to
infection from an asymptomatic individual, infection from a symptomatic
individual, transition from the asymptomatic state to the symptomatic state,
and recovery.

12211112112243210$-1$$-2$$-3$$-4$

Figure 4:  Graphical representation of the contact process with an
asymptomatic state. The black lines represent the infected asymptomatic
individual and the gray lines the infected symptomatic individuals.

More precisely, we think of the integer lattice as the vertex set of a
directed graph in which there is a directed edge $\vec{xy}$ if and only if
vertices $x$ and $y$ are nearest neighbors, and attach the following
exponential clocks to the directed edges and vertices as follows.

  * •

Infection from asymptomatic. Place an exponential clock with rate
$\beta_{1}/2d$ along each directed edge $\vec{xy}$. Each time the clock rings,
say at time $t$, draw an arrow $(x,t)\overset{\hbox{\tiny
1}}{\longrightarrow}(y,t)$ to indicate that if the tail of the arrow is
occupied by a 1 and the head of the arrow is healthy then the head of the
arrow becomes of type 1.

  * •

Infection from symptomatic. Place an exponential clock with rate
$\beta_{2}/2d$ along each directed edge $\vec{xy}$. Each time the clock rings,
say at time $t$, draw an arrow $(x,t)\overset{\hbox{\tiny
2}}{\longrightarrow}(y,t)$ to indicate that if the tail of the arrow is
occupied by a 2 and the head of the arrow is healthy then the head of the
arrow becomes of type 1.

  * •

Transition to symptomatic. Place an exponential clock with rate $\gamma$ at
each vertex $x$. Each time the clock rings, say at time $t$, put a dot
$\bullet$ at point $(x,t)$ to indicate that if site $x$ is asymptomatic then
it becomes symptomatic.

  * •

Recovery. Place an exponential clock with rate one at each vertex $x$. Each
time the clock rings, say at time $t$, put a cross $\times$ at point $(x,t)$
to indicate that if site $x$ is infected (asymptomatic or symptomatic) then it
becomes healthy.

Harris [9] proved for a large class of interacting particle systems including
our contact process with an asymptomatic state that the process is well-
defined and can indeed be constructed using the collection of exponential
clocks and algorithm above, and we refer to Figure 4 for a picture of the
graphical representation and an illustration of how to construct the process.
We now use different variations of the graphical representation together with
the superposition property for Poisson processes to couple/construct jointly
processes with different parameters and prove monotonicity with respect to
each of the three parameters when $\beta_{1}<\beta_{2}$. Let
$\Xi(\beta_{1},\beta_{2},\gamma)$ be the process with infection rate from
asymptomatic individuals $\beta_{1}$, infection rate from symptomatic
individuals $\beta_{2}$, and transition rate to the symptomatic state
$\gamma$.

Monotonicity in $\beta_{1}$. Fix $\beta_{1}<\beta_{1}^{\prime}<\beta_{2}$, let

$\xi_{t}=\Xi(\beta_{1},\beta_{2},\gamma)\quad\hbox{and}\quad\xi_{t}^{\prime}=\Xi(\beta_{1}^{\prime},\beta_{2},\gamma),$

and construct the two processes as follows.

  * •

Place arrows $x\longrightarrow y$ along each directed edge $\vec{xy}$ at rate
$\beta_{1}$. These arrows have the same effect as the 1-arrows and the
2-arrows defined above on both processes.

  * •

Place arrows $x\overset{\hbox{\tiny 1'}}{\longrightarrow}y$ along each
directed edge $\vec{xy}$ at rate $\beta_{1}^{\prime}-\beta_{1}$. These arrows
have the same effect as the 1-arrows defined above on the process
$\xi_{t}^{\prime}$ only and the same effect as the 2-arrows defined above on
both processes.

  * •

Place arrows $x\overset{\hbox{\tiny 2}}{\longrightarrow}y$ along each directed
edge $\vec{xy}$ at rate $\beta_{2}-\beta_{1}^{\prime}$. These arrows have the
same effect as the 2-arrows defined above on both processes.

  * •

Place black dots $\bullet$ at each site $x$ at rate $\gamma$. These black dots
have the same effect as the dots defined above on both processes.

  * •

Place crosses $\times$ at each site $x$ at rate one. These crosses have the
same effect as the crosses defined above on both processes.

Starting both processes from the same configuration, it can be proved that all
the possible pairs of states generated by the coupling
$(\xi_{t},\xi_{t}^{\prime})$ are $S=\\{00,01,02,11,12,22\\}$, i.e., the set
$S$ is closed under the dynamics.

|  $\longrightarrow$ |  00 |  01 |  02 |  11 |  12 |  22 |  $\overset{\hbox{\tiny 1'}}{\longrightarrow}$ |  00 |  01 |  02 |  11 |  12 |  22 |  $\overset{\hbox{\tiny 2}}{\longrightarrow}$ |  00 |  01 |  02 |  11 |  12 |  22 | $\bullet$ | $\times$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
00 |  |  00 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 | 00 | 00
01 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 | 02 | 00
02 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 | 02 | 00
11 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 | 22 | 00
12 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 | 22 | 00
22 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  11 |  11 |  12 |  11 |  12 |  22 | 22 | 00
Table 1: Coupling showing the monotonicity with respect to $\beta_{1}$ when
$\beta_{1}<\beta_{2}$

See Table 1 where the pairs of states in the first column are the possible
pairs of states at site $x$ before a potential update, the pairs of states in
the first row are the possible pairs of states at $y$ before a potential
update, and the pairs of states in the table are the pairs of states resulting
from the dynamics. This implies that the process $\xi_{t}^{\prime}$ has more
infected sites than the process $\xi_{t}$ therefore the set of infected sites
is nondecreasing with respect to $\beta_{1}$.

Monotonicity in $\beta_{2}$. Fix $\beta_{1}<\beta_{2}<\beta_{2}^{\prime}$, let

$\xi_{t}=\Xi(\beta_{1},\beta_{2},\gamma)\quad\hbox{and}\quad\xi_{t}^{\prime}=\Xi(\beta_{1},\beta_{2}^{\prime},\gamma),$

and construct the two processes as follows.

  * •

Place arrows $x\longrightarrow y$ along each directed edge $\vec{xy}$ at rate
$\beta_{1}$. These arrows have the same effect as the 1-arrows and the
2-arrows defined above on both processes.

  * •

Place arrows $x\overset{\hbox{\tiny 2}}{\longrightarrow}y$ along each directed
edge $\vec{xy}$ at rate $\beta_{2}-\beta_{1}$. These arrows have the same
effect as the 2-arrows defined above on both processes.

  * •

Place arrows $x\overset{\hbox{\tiny 2'}}{\longrightarrow}y$ along each
directed edge $\vec{xy}$ at rate $\beta_{2}^{\prime}-\beta_{2}$. These arrows
have the same effect as the 2-arrows defined above on the process
$\xi_{t}^{\prime}$ only.

  * •

Place black dots $\bullet$ at each site $x$ at rate $\gamma$. These black dots
have the same effect as the dots defined above on both processes.

  * •

Place crosses $\times$ at each site $x$ at rate one. These crosses have the
same effect as the crosses defined above on both processes.

|  $\longrightarrow$ |  00 |  01 |  02 |  11 |  12 |  22 |  $\overset{\hbox{\tiny 2}}{\longrightarrow}$ |  00 |  01 |  02 |  11 |  12 |  22 |  $\overset{\hbox{\tiny 2'}}{\longrightarrow}$ |  00 |  01 |  02 |  11 |  12 |  22 | $\bullet$ | $\times$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
00 |  |  00 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 | 00 | 00
01 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 | 02 | 00
02 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 | 02 | 00
11 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 | 22 | 00
12 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 | 22 | 00
22 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  11 |  11 |  12 |  11 |  12 |  22 | 22 | 00
Table 2: Coupling showing the monotonicity with respect to $\beta_{2}$ when
$\beta_{1}<\beta_{2}$

As previously, Table 2 shows that the set $S$ is closed under the dynamics of
the coupling $(\xi_{t},\xi_{t}^{\prime})$, which implies monotonicity with
respect to $\beta_{2}$.

Monotonicity in $\gamma$. Fix $\beta_{1}<\beta_{2}$ and
$\gamma<\gamma^{\prime}$, let

$\xi_{t}=\Xi(\beta_{1},\beta_{2},\gamma)\quad\hbox{and}\quad\xi_{t}^{\prime}=\Xi(\beta_{1},\beta_{2},\gamma^{\prime}),$

and construct the two processes as follows.

  * •

Place arrows $x\longrightarrow y$ along each directed edge $\vec{xy}$ at rate
$\beta_{1}$. These arrows have the same effect as the 1-arrows and the
2-arrows defined above on both processes.

  * •

Place arrows $x\overset{\hbox{\tiny 2}}{\longrightarrow}y$ along each directed
edge $\vec{xy}$ at rate $\beta_{2}-\beta_{1}$. These arrows have the same
effect as the 2-arrows defined above on both processes.

  * •

Place black dots $\bullet$ at each site $x$ at rate $\gamma$. These black dots
have the same effect as the dots defined above on both processes.

  * •

Place white dots $\circ$ at each site $x$ at rate $\gamma^{\prime}-\gamma$.
These white dots have the same effect as the dots defined above on the process
$\xi_{t}^{\prime}$ only.

  * •

Place crosses $\times$ at each site $x$ at rate one. These crosses have the
same effect as the crosses defined above on both processes.

|  $\longrightarrow$ |  00 |  01 |  02 |  11 |  12 |  22 |  $\overset{\hbox{\tiny 2}}{\longrightarrow}$ |  00 |  01 |  02 |  11 |  12 |  22 | $\bullet$ | $\circ$ | $\times$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
00 |  |  00 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 | 00 | 00 | 00
01 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 | 02 | 02 | 00
02 |  |  01 |  01 |  02 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 | 02 | 02 | 00
11 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  00 |  01 |  02 |  11 |  12 |  22 | 22 | 12 | 00
12 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  01 |  01 |  02 |  11 |  12 |  22 | 22 | 12 | 00
22 |  |  11 |  11 |  12 |  11 |  12 |  22 |  |  11 |  11 |  12 |  11 |  12 |  22 | 22 | 22 | 00
Table 3: Coupling showing the monotonicity with respect to $\gamma$ when
$\beta_{1}<\beta_{2}$

Table 3 confirms that the set $S$ is again closed under the dynamics of the
coupling $(\xi_{t},\xi_{t}^{\prime})$ therefore the set of infected sites is
nondecreasing with respect to $\gamma$.
Using similar joint graphical representations in the case
$\beta_{1}>\beta_{2}$ and starting both processes from the same configuration
first create the same pairs of states as before, including 12. However,
because the asymptomatic individuals are more infectious, the graphical
representation now contains type 1 arrows through which only the 1s can infect
nearby healthy individuals. In particular, the presence of a type 1 arrow
going from a site in state 12 to a site in state 00 creates a site in state
10, thus showing that a standard coupling fails to prove monotonicity when
$\beta_{1}>\beta_{2}$.

##  4 Proof of Theorems 1.a and 3.a (strong survival)

This section is devoted to the survival parts of Theorems 1 and 3. The
intuition behind these results is simple. When $\gamma=0$, the epidemic model
reduces to a contact process with parameter $\beta_{1}$ so there is an
epidemic when $\beta_{1}>\beta_{c}$. Similarly, when $\gamma=\infty$, the
epidemic model reduces to a contact process with parameter $\beta_{2}$ so
there is an epidemic when $\beta_{2}>\beta_{c}$. To extend these results to
$\gamma$ small/large, the idea is to apply a standard perturbation argument to
the block construction introduced by Bezuidenhout and Grimmett [1] to study
the supercritical contact process.

Proof of Theorem 3.a. Let $\beta_{1}>\beta_{c}$ and think of the contact
process

$\zeta_{t}^{0}:\mathbb{Z}^{d}\longrightarrow\\{0,1\\}\quad\hbox{where}\quad
0=\hbox{healthy}\quad\hbox{and}\quad 1=\hbox{infected}$

with parameter $\beta_{1}$ as being constructed from the graphical
representation depicted in Figure 4 by using the type 1 infection arrows and
the recovery crosses (but ignoring the type 2 arrows and the dots). Note that,
denoting by $P_{\gamma}$ the probability associated to the epidemic model with
parameter $\gamma$, the probability associated to the contact process is
simply $P_{0}$. Now let

$\mathscr{L}_{1}=\\{(m,n)\in\mathbb{Z}\times\mathbb{Z}_{+}:m+n\ \hbox{is
even}\\}$

which we turn into a directed graph $\vec{\mathscr{L}}_{1}$ by placing arrows

$(m,n)\to(m^{\prime},n^{\prime})\quad\hbox{if and only
if}\quad|m-m^{\prime}|=1\ \ \hbox{and}\ \ n^{\prime}=n+1.$

The idea of the block construction is to compare the contact process properly
rescaled in space and time with oriented site percolation on this directed
graph. Having four integers $J,K,L,T$ playing the role of space and time
scales, we let $D=[\\![-J,J]\\!]^{d}\times\\{0\\}$ and

$\begin{array}[]{rcl}A_{m,n}&=&(2mLe_{1},nKT)+[-2L,2L]^{d}\times[0,(K+1)T]\vspace*{4pt}\\\
B_{m,n}&=&(2mLe_{1},(n+1)KT)+[-2L,2L]^{d}\times[0,T]\end{array}$

for all $(m,n)\in\mathscr{L}_{1}$ where $e_{1}$ is the first unit vector in
$\mathbb{Z}^{d}$. Let also

$E_{m,n}=\hbox{there is a translate of\leavevmode\nobreak\ $D$ contained
in\leavevmode\nobreak\ $B_{m,n}$ that is fully infected}$

and declare site $(m,n)$ to be a good site if $E_{m,n}$ occurs. Bezuidenhout
and Grimmett proved that, for all infection rates $\beta_{1}>\beta_{c}$ and
all $\epsilon>0$ small, the scale parameters can be chosen large enough that
the set of good sites dominates stochastically the set of wet sites in an
oriented site percolation process on $\vec{\mathscr{L}}_{1}$ with parameter
$1-\epsilon$. More precisely, there exists a collection of good events
$G_{m,n}$ that only depend on the graphical representation in $A_{m,n}$ such
that

$P_{0}(G_{m,n})\geq 1-\epsilon\quad\hbox{and}\quad E_{m,n}\cap G_{m,n}\subset
E_{m-1,n+1}\cap E_{m+1,n+1}.$ (15)

In other words, a translate of $D$ fully infected in $B_{m,n}$ will produce
with probability at least $1-\epsilon$ two translates of $D$ fully infected in
the two blocks immediately above. Choosing $\epsilon>0$ small so that the
percolation process with parameter $1-\epsilon$ is supercritical, this
construction shows that the contact process with parameter
$\beta_{1}>\beta_{c}$ survives strongly. Returning to the epidemic model, we
denote by $G_{m,n}^{0}$ the event that there is no dots in the block
$A_{m,n}$. On this event, the epidemic model and the contact process share the
same graphical representation in the block so

$E_{m,n}\cap(G_{m,n}\cap G_{m,n}^{0})\subset E_{m-1,n+1}\cap E_{m+1,n+1}.$
(16)

In addition, the scale parameters being fixed, letting $S=(4L+1)^{d}(K+1)T$ be
the size (number of sites times length of the time interval) of $A_{m,n}$,
there exists $\gamma_{-}>0$ small such that

$P_{\gamma}(G_{m,n}^{0})=P(\operatorname{Poisson}(\gamma S)=0)=\exp(-\gamma
S)\geq 1-\epsilon$ (17)

for all $\gamma<\gamma_{-}$. In particular, for all $\gamma<\gamma_{-}$,

$P_{\gamma}(G_{m,n}\cap G_{m,n}^{0})\geq
1-P_{0}(G_{m,n}^{c})-P_{\gamma}((G_{m,n}^{0})^{c})\geq 1-2\epsilon.$ (18)

Combining (16) and (18) shows that (15) also holds (with $2\epsilon$ instead
of $\epsilon$ and different good events) for the epidemic model with
$\gamma<\gamma_{-}$. Because $\epsilon$ can be chosen arbitrarily small, this
shows strong survival (there is an epidemic) and proves the survival part of
Theorem 3.

Proof of Theorem 1.a. Let $\beta_{2}>\beta_{c}$ and think of the contact
process

$\zeta_{t}^{\infty}:\mathbb{Z}^{d}\longrightarrow\\{0,2\\}\quad\hbox{where}\quad
0=\hbox{healthy}\quad\hbox{and}\quad 2=\hbox{infected}$

with parameter $\beta_{2}$ as being constructed from the graphical
representation depicted in Figure 4 by only using the type 2 infection arrows
and the recovery crosses. Note that the probability associated to this contact
process is now $P_{\infty}$. In addition, because the process is
supercritical, the two statements in (15) still hold but for the probability
$P_{\infty}$ instead of $P_{0}$, i.e.,

$P_{\infty}(G_{m,n})\geq 1-\epsilon\quad\hbox{and}\quad E_{m,n}\cap
G_{m,n}\subset E_{m-1,n+1}\cap E_{m+1,n+1}.$ (19)

To prove strong survival of the epidemic model when $\gamma$ is large, note
that the process behaves like its coupled contact process if each time a 2 in
the contact process tries to infect a neighbor then the same happens in the
epidemic model, and if none of the 1s in the epidemic model tries to infect a
neighbor. From the point of view of the graphical representation, this occurs
in the block $A_{m,n}$ if each time there is an arrow pointing at $(x,t)$ in
the block, there is a dot at site $x$ before the next arrow starting at $x$.
In particular, denoting this event by $G_{m,n}^{\infty}$, the analog of (16)
holds:

$E_{m,n}\cap(G_{m,n}\cap G_{m,n}^{\infty})\subset E_{m-1,n+1}\cap
E_{m+1,n+1}.$ (20)

Following the same approach as before, the last step to prove the theorem is
to show that the probability of $G_{m,n}^{\infty}$ is close to one when
$\gamma$ large. To do this, let

$H_{m,n}=\hbox{there are less than\leavevmode\nobreak\
$2H=2(\beta_{1}+\beta_{2})S$ arrows pointing at\leavevmode\nobreak\
$A_{m,n}$}$

where $S$ is the size of $A_{m,n}$. Because the arrows occur at rate
$\beta_{1}+\beta_{2}$ at each site,

$P_{\gamma}(H_{m,n}^{c})=P(\operatorname{Poisson}((\beta_{1}+\beta_{2})S)\geq
2(\beta_{1}+\beta_{2})S)\leq\epsilon/2$ (21)

for all $K,L,T$ large. Now, let $X$ and $Y$ be exponentially distributed with
rate $\beta_{1}+\beta_{2}$ and $\gamma$, respectively. The scale parameters
being fixed so that (19) and (21) hold, because the dots appear at rate
$\gamma$ at each site, it follows from the superposition property that

$P_{\gamma}((G_{m,n}^{\infty})^{c}\,|\,H_{m,n})\leq
2HP(X<Y)=2H\bigg{(}\frac{\beta_{1}+\beta_{2}}{\beta_{1}+\beta_{2}+\gamma}\bigg{)}\leq\epsilon/2$
(22)

for all $\gamma$ larger than some $\gamma_{+}<\infty$. Combining (21) and
(22), we deduce that

$\begin{array}[]{rcl}P_{\gamma}((G_{m,n}^{\infty})^{c})&=&P_{\gamma}((G_{m,n}^{\infty})^{c}\,|\,H_{m,n})P_{\gamma}(H_{m,n})+P_{\gamma}((G_{m,n}^{\infty})^{c}\,|\,H_{m,n}^{c})P_{\gamma}(H_{m,n}^{c})\vspace*{4pt}\\\
&\leq&P_{\gamma}((G_{m,n}^{\infty})^{c}\,|\,H_{m,n})+P_{\gamma}(H_{m,n}^{c})\leq\epsilon\end{array}$

which, together with (19), implies that, for all $\gamma>\gamma_{+}$,

$P_{\gamma}(G_{m,n}\cap G_{m,n}^{\infty})\geq
1-P_{\infty}(G_{m,n}^{c})-P_{\gamma}((G_{m,n}^{\infty})^{c})\geq 1-2\epsilon.$
(23)

As previously, combining (20) and (23) shows strong survival (there is an
epidemic) for all $\gamma>\gamma_{+}$, which completes the proof of the
survival part of Theorem 1.

##  5 Proof of Theorems 1.b and 3.b (extinction)

This section is devoted to the extinction parts of Theorems 1 and 3. The proof
relies on the same three ingredients as for the survival parts: coupling of
the epidemic model with the contact processes $\zeta_{t}^{0}$ and
$\zeta_{t}^{\infty}$, block construction, and perturbation arguments at
$\gamma=0$ and at $\gamma=\infty$. The main difference with the proof of
survival is the block construction.

Proof of Theorem 1.b. As previously, let $\zeta_{t}^{0}$ be the basic contact
process with parameter $\beta_{1}$ coupled with the epidemic model, and assume
that $\beta_{1}<\beta_{c}$. To compare the process properly rescaled in space
and time with oriented site percolation, we now let
$\mathscr{L}_{2}=\mathbb{Z}^{d}\times\mathbb{Z}_{+}$ which we turn into a
directed graph $\vec{\mathscr{L}}_{2}$ by placing arrows

$\begin{array}[]{l}(m,n)\to(m^{\prime},n^{\prime})\quad\hbox{if and only
if}\vspace*{4pt}\\\ \hskip
20.0pt|m_{1}-m_{1}^{\prime}|+\cdots+|m_{d}-m_{d}^{\prime}|+|n-n^{\prime}|=1\ \
\hbox{and}\ \ n\leq n^{\prime}.\end{array}$

One can think of the horizontal arrows as expansion in space and of the
vertical arrows as persistence in time. Having an integer $L$ playing the role
of the space and time scales, we let

$\begin{array}[]{rcl}A_{m,n}&=&(2mL,nL)+[-2L,2L]^{d}\times[0,2L]\vspace*{4pt}\\\
B_{m,n}&=&(2mL,nL)+[-L,L]^{d}\times[L,2L]\end{array}$

for all $(m,n)\in\mathscr{L}_{2}$. Then, we define the events

$E_{m,n}=\hbox{the space-time block\leavevmode\nobreak\ $B_{m,n}$ is fully
healthy}$

and declare site $(m,n)$ to be a good site if $E_{m,n}$ occurs. To prove that
the set of good sites dominates the set of wet sites in an oriented site
percolation process on $\vec{\mathscr{L}}_{2}$ with parameter arbitrarily
close to one, the idea is to show that, with high probability when $L$ is
large, none of the infection paths starting from the bottom or the periphery
of $A_{m,n}$ reaches the block $B_{m,n}$. For $s<t$, we say that there is an
infection path from $(y,s)$ to $(x,t)$ if we can go from the first space-time
point to the second one going up in the graphical representation, crossing the
infection arrows from tail to head but avoiding the recovery crosses.
Bezuidenhout and Grimmett [Bezuidenhout_Grimmett_1991] proved the following
exponential decay for the length of the infection paths in the subcritical
phase: starting with a single infection at site $x$ and letting
$|\zeta_{s}^{0}|$ be the number of infected sites at time $s$,

$P\bigg{(}\int_{0}^{\infty}|\zeta_{s}^{0}|\,ds>t\bigg{)}\leq e^{-\alpha
t}\quad\hbox{for some}\quad\alpha>0.$ (24)

Let $\epsilon>0$ be small, and let

$S_{b}=(4L+1)^{2}\quad\hbox{and}\quad S_{p}=2L((4L+1)^{d}-(4L-1)^{d})$

be the size (number of sites) of the bottom of the block $A_{m,n}$ and the
size (number of sites times length of the time interval) of its periphery.
Consider also the events

$H_{m,n}=\hbox{there are less than\leavevmode\nobreak\ $2\beta_{1}S_{p}$
arrows starting from the pheriphery of\leavevmode\nobreak\ $A_{m,n}$}.$

Using that an infection path starting from the bottom or the periphery of the
block $A_{m,n}$ must have length at least $L$ to reach $B_{m,n}$, it follows
from (24) that

$\begin{array}[]{l}P(\hbox{path from bottom or periphery
of\leavevmode\nobreak\ $A_{m,n}$ reaching\leavevmode\nobreak\
$B_{m,n}$})\vspace*{4pt}\\\ \hskip 30.0pt\leq P(\hbox{path from
bottom})+P(\hbox{path from periphery}\,|\,H_{m,n})P(H_{m,n})\vspace*{4pt}\\\
\hskip 60.0pt+\ P(\hbox{path from
periphery}\,|\,H_{m,n}^{c})P(H_{m,n}^{c})\vspace*{4pt}\\\ \hskip 30.0pt\leq
P(\hbox{path from bottom})+P(\hbox{path from
periphery}\,|\,H_{m,n})+P(H_{m,n}^{c})\vspace*{4pt}\\\ \hskip
30.0pt\leq(S_{b}+2\beta_{1}S_{p})\,e^{-\alpha
L}+P(\operatorname{Poisson}(\beta_{1}S_{p})\geq
2\beta_{1}S_{p})\leq\epsilon\end{array}$ (25)

for all $L$ sufficiently large. In particular, there exists a collection of
good events $G_{m,n}$ that only depend on the graphical representation in
$A_{m,n}$ such that

$P_{0}(G_{m,n})\geq 1-\epsilon\quad\hbox{and}\quad G_{m,n}\subset E_{m,n}.$
(26)

In other words, for all infection rates $\beta_{1}<\beta_{c}$ and all
$\epsilon>0$ small, the scale parameter $L$ can be chosen large enough that
the set of good sites dominates stochastically the set of wet sites in an
oriented site percolation process on $\vec{\mathscr{L}}_{2}$ with parameter
$1-\epsilon$. For $\epsilon>0$ small enough, not only the set of open sites
percolates but also the probability of a path of closed (not open) sites with
length at least $n$ starting at $(0,0)$ decays exponentially with $n$. Because
the infection cannot appear spontaneously, this shows extinction. To extend
this result to the epidemic model with parameter $\gamma>0$ small, the idea is
again to use a perturbation argument similar to (17).

Proof of Theorem 3.b. As in the previous section, let $\zeta_{t}^{\infty}$ be
the (subcritical) contact process with parameter $\beta_{2}<\beta_{c}$ coupled
with the epidemic model. Then, the two properties in (26) again hold for this
process and can be extended to the epidemic model with $\gamma<\infty$ large
(depending on the scale parameter $L$) using the estimates (21) and (22).

##  6 Proof of Theorem 2

Throughout this section, we assume that $\beta_{1}=0$, meaning that the 1s are
not contagious. Because only the 2s can infect nearby healthy individuals, in
order to prove extinction of the disease when $\gamma$ is small, the idea is
to show that the set of 2s is dominated by a subcritical Galton–Watson
branching process. More precisely, the first objective is to prove that,
starting with a single 2 in an otherwise healthy population, the infection
paths that originate from this 2 cannot expand too far or live too long. The
next lemma shows how to control the radius of the infection.

###### Lemma 5

Assume that $\beta_{1}=0$ and $\gamma<1/(2d-1)$. Then, for all $K>0$, starting
with a single 2 at the origin in an otherwise healthy population

$P(\xi_{t}(x)\neq 0\ \hbox{for some}\
(x,t)\notin[-K,K]^{d}\times\mathbb{R}_{+})\leq(2d\gamma/(1+\gamma))^{K}.$

###### Proof.

By monotonicity of the set of infected sites with respect to the infection
rates when $\beta_{1}<\beta_{2}$, it suffices to prove the result in the
limiting case $\beta_{2}=\infty$. In this case, the process is dominated
stochastically by a subcritical Galton–Watson branching process. Each time a 2
at site $x$ gives birth to a 1 at site $y\in N_{x}$ that turns into a 2,

  * •

we think of the 2 at site $x$ as the parent of the 2 at site $y$ and

  * •

we think of the 2 at site $y$ as the child of the 2 at site $x$,

which induces a partition of the set of 2s into generations, starting with the
initial single 2 belonging to generation zero. Note that, in the limiting case
$\beta_{2}=\infty$, each 2 gives birth instantaneously to 1s onto each of the
healthy neighbors, whose number is bounded by $2d$. In addition, by the
superposition property for the exponential distribution, each of these 1s
turns into a 2 independently with probability $\gamma/(1+\gamma)$. In
particular, the number of 2s at generation $n$ is dominated by the
Galton–Watson branching process $(X_{n})$ that has offspring distribution

$\operatorname{Binomial}(2d,\gamma/(1+\gamma))\quad\hbox{with
mean}\quad\mu=2d\gamma/(1+\gamma)<1.$

Using that the process $(\mu^{-n}X_{n})$ is a martingale, the probability that
the Galton–Watson branching process lives for at least $n$ generations is

$P(X_{n}\neq 0)\leq
E(X_{n})=\mu^{n}E(X_{0})=\mu^{n}=(2d\gamma/(1+\gamma))^{n}.$ (27)

To prove that the infection generated by the initial 2 cannot spread too far,
note that, because it takes at least $K$ generations to bring a 2 outside the
spatial box $[-K,K]^{d}$,

$P(\xi_{t}(x)=2\ \hbox{for some}\
(x,t)\notin[-K,K]^{d}\times\mathbb{R}_{+})\leq P(X_{K}\neq
0)\leq(2d\gamma/(1+\gamma))^{K}$

according to (27), which completes the proof. ∎

Using the previous branching process, we now control the duration of the
infection.

###### Lemma 6

Assume that $\beta_{1}=0$ and $\gamma<1/(2d-1)$. Then, for all $K>0$, starting
with a single 2 at the origin in an otherwise healthy population

$E(T_{0})\leq\frac{2d}{1-2d\gamma/(1+\gamma)}<\infty\quad\hbox{where}\quad
T_{0}=\inf\\{t:\xi_{t}\equiv 0\\}.$

###### Proof.

Let $N$ be the total number of individuals in the Galton–Watson branching
process introduced in the previous proof. By the monotone convergence theorem
and (27),

$E(N)=E\bigg{(}\sum_{n=0}^{\infty}\,X_{n}\bigg{)}=\sum_{n=0}^{\infty}\,E(X_{n})=\sum_{n=0}^{\infty}\,\mu^{n}=\frac{1}{1-2d\gamma/(1+\gamma)}<\infty.$
(28)

In addition, because the number of 1s created by each 2 is bounded by $2d$,
and each infected site recovers at rate one, the time to extinction of the
infection $T_{0}$ is bounded by the sum of $2dN$ independent exponential
random variables $Y_{i}$ with parameter one:

$T_{0}\leq Y_{1}+Y_{2}+\cdots+Y_{2dN}.$ (29)

Combining (28) and (29), we deduce that

$\begin{array}[]{rcl}E(T_{0})&=&E(E(T_{0}\,|\,N))\leq
E(E(Y_{1}+\cdots+Y_{2dN}\,|\,N))\vspace*{4pt}\\\
&=&E(2dNE(Y_{1}))=2dE(N)E(Y_{1})=2d/(1-2d\gamma/(1+\gamma))<\infty,\end{array}$

which proves the lemma. ∎

Combining Lemmas 5 and 6, we can now deduce Theorem 2 which states that the
population recovers. The basic idea is to use the exponential decay of the
radius of the infection to show that, regardless of the initial configuration,
the probability of an invasion path starting outside a large box centered at
the origin and reaching the origin can be made arbitrarily small. Because this
box is bounded, using that the mean duration of the infection is finite, we
also have that the infection paths starting from the box die out collectively
in a finite time. In particular, the probability that the origin is infected
at a large time is small.

###### Lemma 7

Let $\beta_{1}=0$ and $\gamma<1/(2d-1)$. Then,

$\begin{array}[]{l}\lim_{t\to\infty}P(\xi_{t}(0)\neq 0)=0\quad\hbox{for
all}\quad\beta_{2}\geq 0.\end{array}$

###### Proof.

By monotonicity, it suffices to prove the result for $\beta_{2}=\infty$ and
starting from the all 2 configuration. In this case, the infected region
generated by the 2 initially at site $y$ is dominated by its counterpart when
starting with a single 2 at site $y$ in an otherwise healthy population. In
particular, letting $\epsilon>0$ be small, it follows from Lemma 5 that the
probability of an infection path starting from outside the cube
$\Lambda_{K}=[-K,K]^{d}$ and ending at the origin is bounded by

$\begin{array}[]{l}\displaystyle\sum_{i=0}^{\infty}\
\Big{(}\operatorname{Vol}\,(\Lambda_{K+i+1})-\operatorname{Vol}\,(\Lambda_{K+i})\Big{)}\bigg{(}\frac{2d\gamma}{1+\gamma}\bigg{)}^{K}\vspace*{0pt}\\\
\hskip 100.0pt\leq\displaystyle\sum_{i=0}^{\infty}\
2d(2K+2i+3)^{d-1}\bigg{(}\frac{2d\gamma}{1+\gamma}\bigg{)}^{K}\leq\epsilon/2\end{array}$
(30)

for some $K<\infty$ large fixed from now on. Let $T_{y}$ be the extinction
time of the infected cluster starting at site $y$. The integer $K$ being
fixed, it follows from Lemma 6 that the time $\tau_{K}$ at which the infected
region generated by all the 2s initially in $\Lambda_{K}$ satisfies

$\begin{array}[]{l}E(\tau_{K})\leq
E(\max_{y\in\Lambda_{K}}T_{y})\leq\displaystyle\sum_{y\in\Lambda_{K}}E(T_{y})\leq\frac{2d(2K+1)^{d}}{1-2d\gamma/(1+\gamma)}<\infty.\end{array}$

This implies that time $\tau_{K}$ is almost surely finite therefore

$P(\tau_{K}>T)\leq\epsilon/2\quad\hbox{for some}\quad T<\infty.$ (31)

Combining (30) which controls the infection paths starting from outside
$\Lambda_{K}$ and (31) which controls the infection paths starting from inside
$\Lambda_{k}$, we deduce that the probability of the origin being infected at
some time larger than $T$ is bounded by $\epsilon$, which proves the lemma. ∎

Because the evolution rules of the process are translation invariant, the
previous lemma also applies to each site of the lattice, which implies that,
regardless of the initial configuration, the density of infected sites
converges to zero, and shows the theorem.

##  7 Proof of Theorem 4

Recall that, under the assumptions of Theorem 4, the rate $\beta_{2}=0$, the
rate $\gamma$ is fixed and finite, and the process starts with a single 1 at
the origin in an otherwise healthy population. To prove survival of the
infection, the idea is to use a construction from Cox and Durrett [4] to show
that the set of sites that ever get infected dominates a certain independent
site percolation process that is supercritical provided $\beta_{1}$ is
sufficiently large. To each site $x$ and each oriented edge $\vec{xy}$ where
site $y$ is a nearest neighbor of site $x$, we attach independent random
variables

$\tau(x,y)=\operatorname{Exponential}(\beta_{1}/2d)\quad\hbox{and}\quad\sigma(x)=\operatorname{Exponential}(1+\gamma).$

Then, we say that directed edge $\vec{xy}$ is

$\hbox{open when}\ \tau(x,y)<\sigma(x)\quad\hbox{and}\quad\hbox{closed when}\
\tau(x,y)\geq\sigma(x)$

and define the percolation cluster

$C_{0}=\\{x\in\mathbb{Z}^{d}:\hbox{there is a directed
path\leavevmode\nobreak\ $0\to x$ of open edges}\\}.$

Thinking of $\tau(x,y)$ as the time it takes for a 1 at site $x$ to spread
germs to $y$ and thinking of $\sigma(x)$ as the time it take for a 1 at site
$x$ to either recover, which occurs at rate one, or become symptomatic, which
occurs at rate $\gamma$, the presence of a directed open path $0\to x$ in the
percolation process implies that the infection spreads up to site $x$ in the
interacting particle system. Using also the memoryless property of the
exponential distribution, the contact process and the directed bond
percolation process can be coupled in such a way that the sites that either
get infected contains the percolation cluster $C_{0}$ from which it follows
that the probability of survival is

$P_{0}(\xi_{t}\not\equiv 0\ \hbox{for all}\ t)\geq P(|C_{0}|=\infty).$ (32)

Therefore, to prove the theorem, it suffices to show that percolation occurs
with positive probability when the infection rate $\beta_{1}$ is large enough.
Note that whether the directed edges starting from the same site $x$ are open
or closed are not independent events because they all depend on the
realization of the random variable $\sigma(x)$. However, declaring site $x$ to
be open if

$\tau(x,y)<\sigma(x)\quad\hbox{for all}\quad y\in N_{x}$

and closed otherwise, whether different sites are open or closed are now
independent events because they depend on disjoint collections of independent
random variables. In addition, because a directed edge $\vec{xy}$ is open
whenever site $x$ is open, we have

$\bar{C}_{0}\subset
C_{0}\quad\hbox{where}\quad\bar{C}_{0}=\\{x\in\mathbb{Z}^{d}:\hbox{there is a
directed path\leavevmode\nobreak\ $0\to x$ of open sites}\\}.$ (33)

In view of (32)–(33), to prove survival of the infection, it suffices to show
that, for $\beta_{1}$ large enough, sites are (independently) open with a
probability that exceeds the critical value $p_{c}(\mathbb{Z}^{d})$ of site
percolation, which is known to be strictly less than one when $d\geq 2$. Let

$\tau_{1},\tau_{2},\ldots,\tau_{2d}=\operatorname{Exponential}(\beta_{1}/2d)\quad\hbox{and}\quad\sigma=\operatorname{Exponential}(1+\gamma)$

be independent, and let $\tau=\max(\tau_{1},\tau_{2},\ldots,\tau_{2d})$. Then,
conditioning on the value of $\sigma$, we obtain that the probability
$p(\beta_{1},\gamma)$ of a given site being open is equal to

$\begin{array}[]{rcl}\displaystyle
P(\tau<\sigma)&=&\displaystyle\int_{0}^{\infty}P(\tau<t)\,(1+\gamma)\,e^{-(1+\gamma)t}\,dt\vspace*{6pt}\\\
&=&\displaystyle\int_{0}^{\infty}(1-e^{-\beta_{1}t/2d})^{2d}\,(1+\gamma)\,e^{-(1+\gamma)t}\,dt\vspace*{2pt}\\\
&=&\displaystyle\int_{0}^{\infty}\sum_{k=0}^{2d}{2d\choose
k}(-1)^{k}(1+\gamma)\,e^{-(1+\gamma+k\beta_{1}/2d)t}\,dt=\displaystyle\sum_{k=0}^{2d}{2d\choose
k}\frac{(-1)^{k}}{1+\frac{k\beta_{1}}{2d(1+\gamma)}}.\end{array}$

For each $\gamma$, the function $p(\,\cdot\,,\gamma)$ is continuous increasing
on $\mathbb{R}_{+}$ and

$\begin{array}[]{l}p(0,\gamma)=0\qquad\hbox{and}\qquad\lim_{\beta_{1}\to\infty}p(\beta_{1},\gamma)=1.\end{array}$

In particular, there exists a unique $\bar{\beta}$ such that
$p(\bar{\beta},\gamma)=p_{c}(\mathbb{Z}^{d})$ and this $\bar{\beta}$ is
finite. In addition, it follows from (32)–(33) and strict monotonicity that

$P_{0}(\xi_{t}\not\equiv 0\ \hbox{for all}\ t)\geq P(|C_{0}|=\infty)\geq
P(|\bar{C}_{0}|=\infty)>0\quad\hbox{for all}\quad\beta_{1}>\bar{\beta}.$

This complete the proof of Theorem 4. For instance, in the two-dimensional
case, Kesten [13] proved that the critical value of undirected bond
percolation is equal to one-half. This, together with a general result due to
Grimmett and Stacey [8], implies that

$p_{c}(\mathbb{Z}^{2})\leq 1-\bigg{(}1-\frac{1}{2}\bigg{)}^{3}=\frac{7}{8}$

and one can easily check that, for all $r>15.27$,

$\sum_{k=0}^{4}{4\choose
k}\frac{(-1)^{k}}{1+kr}=1-\frac{4}{1+r}+\frac{6}{1+2r}-\frac{4}{1+3r}+\frac{1}{1+4r}>\frac{7}{8}.$

This shows that the infection survives when $\beta_{1}>4\cdot
15.27\cdot(1+\gamma)=61.08\cdot(1+\gamma)$.

## References

  * [1] C. Bezuidenhout and G. Grimmett.  The critical contact process dies out.  Ann. Probab., 18(4):1462–1482, 1990.
  * [2] A. B. Bjornson, M. A. Mellencamp, and G. M. Schiff.  Complement is activated in the upper respiratory tract during influenza virus infection.  American Review of Respiratory Disease, 143(5 Pt 1):1062–1066, 1991\.
  * [3] R. B. Couch, R. G. Jr. Douglas, D. S. Fedson, and J. A. Kasel.  Correlated studies of a recombinant influenza-virus vaccine. 3. protection against experimental influenza in man.  Journal of Infectious Diseases, 124:473–480, 1971.
  * [4] J. T. Cox and R. Durrett.  Limit theorems for the spread of epidemics and forest fires.  Stochastic Process. Appl., 30(2):171–191, 1988.
  * [5] P. K. Drain et al.  Duration of viral infectiousness and correlation with symptoms and diagnostic testing in non-hospitalized adults during acute sars-cov-2 infection: A longitudinal cohort study.  Journal of Clinical Virology, 161:105420, April 2023.
  * [6] R. Durrett and C. Neuhauser.  Epidemics with recovery in $D=2$.  Ann. Appl. Probab., 1(2):189–206, 1991.
  * [7] E. Foxall.  New results for the two-stage contact process.  J. Appl. Probab., 52(1):258–268, 2015.
  * [8] G. Grimmett and A. Stacey.  Critical probabilities for site and bond percolation models.  Ann. Probab., 26(4):1788–1812, 1998.
  * [9] T. E. Harris.  Nearest-neighbor Markov interaction processes on multidimensional lattices.  Advances in Math., 9:66–89, 1972.
  * [10] T. E. Harris.  Contact interactions on a lattice.  Ann. Probability, 2:969–988, 1974.
  * [11] X. He, E. H. Y. Lau, P. Wu, and et al.  Temporal dynamics in viral shedding and transmissibility of COVID-19.  Nat Med, 26(5):672–675, 2020.
  * [12] M. A. Johansson et al.  Sars-cov-2 transmission from people without covid-19 symptoms.  JAMA network open, 4(1):e2035057–e2035057, 2021.
  * [13] H. Kesten.  The critical probability of bond percolation on the square lattice equals ${1\over 2}$.  Comm. Math. Phys., 74(1):41–59, 1980.
  * [14] S. M. Krone.  The two-stage contact process.  Ann. Appl. Probab., 9(2):331–351, 1999.
  * [15] T. M. Liggett.  Interacting particle systems, volume 276 of Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, New York, 1985.
  * [16] T. M. Liggett.  Stochastic interacting systems: contact, voter and exclusion processes, volume 324 of Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, Berlin, 1999.
  * [17] W. C. Miller, N. E. Rosenberg, S. E. Rutstein, and K. A. Powers.  Role of acute and early hiv infection in the sexual transmission of hiv.  Current Opinion in HIV and AIDS, 5(4):277–282, 2010.
  * [18] E. Patrozou and L. A. Mermel.  Does influenza transmission occur from asymptomatic infection or prior to symptom onset?  Public Health Reports, 124(2):193–196, Mar-Apr 2009.
  * [19] J. Qin, C. You, Q. Lin, T. Hu, S. Yu, and X.-H. Zhou.  Estimation of incubation period distribution of covid-19 using disease onset forward time: a novel cross-sectional and forward follow-up study.  Science advances, 6(33):eabc1202, 2020.


# A flexible and fast PyTorch toolkit for simulating training and inference on
analog crossbar arrays111submitted to AICAS 2021

Malte J. Rasch  IBM Research  Diego Moreda  IBM Research  Tayfun Gokmen  IBM
Research  Manuel Le Gallo  IBM Research  Fabio Carta  IBM Research  Cindy
Goldberg  IBM Research  Kaoutar El Maghraoui  IBM Research  Abu Sebastian  IBM
Research  Vijay Narayanan  IBM Research

###### Abstract

We introduce the IBM Analog Hardware Acceleration Kit, a new and first of a
kind open source toolkit to simulate analog crossbar arrays in a convenient
fashion from within PyTorch (freely available at
https://github.com/IBM/aihwkit). The toolkit is under active development and
is centered around the concept of an “analog tile” which captures the
computations performed on a crossbar array. Analog tiles are building blocks
that can be used to extend existing network modules with analog components and
compose arbitrary artificial neural networks (ANNs) using the flexibility of
the PyTorch framework. Analog tiles can be conveniently configured to emulate
a plethora of different analog hardware characteristics and their non-
idealities, such as device-to-device and cycle-to-cycle variations, resistive
device response curves, and weight and output noise. Additionally, the toolkit
makes it possible to design custom unit cell configurations and to use
advanced analog optimization algorithms such as Tiki-Taka. Moreover, the
backward and update behavior can be set to “ideal" to enable hardware-aware
training features for chips that target inference acceleration only. To
evaluate the inference accuracy of such chips over time, we provide
statistical programming noise and drift models calibrated on phase-change
memory hardware. Our new toolkit is fully GPU accelerated and can be used to
conveniently estimate the impact of material properties and non-idealities of
future analog technology on the accuracy for arbitrary ANNs.

##  1 Introduction

To cope with the ever increasing demand for computing resources in artificial
intelligence applications, dedicated hardware accelerators for artificial
neural network (ANN) training and inference have been proposed recently. One
promising future technology is the use of memristive crossbar arrays for
accelerating the ubiquitous matrix-vector multiply and rank-update operations
in ANNs by employing in-memory computation of matrices stored as analog
quantities in tunable resistive elements [1, 7, 8, 13]. While considerable
run-time and performance gain over today’s technology is projected for analog
hardware in principle [5], noisy and non-linear device characteristics as well
as design restrictions such as stationary weight matrices, demand for a new
breed of ANN topologies and algorithmic advances adapted to analog technology
to unleash its full potential. An ideal outcome of such “analog AI”
algorithmic endeavour would be reminiscent of the increase in popularity of
convolutional networks once it became clear that their compute is ideally
suited to be performed on graphic processing units (GPUs).

To be viable, these new analog algorithms and optimized network architectures
need to be tested in a simulator that takes analog hardware constraints and
realistic material properties into account. Moreover, these simulations need
to be applied on realistically sized deep neural networks (DNNs) and datasets
to make informed decisions about hardware design and the expected accuracy.

Figure 1: Repository and code structure of the IBM Analog Hardware
Acceleration Kit.

While there exists a number of current software packages that are dedicated to
simulating AI workloads on non-volatile memory elements, notably MLP-NeuroSim
[2] and RxNN [9], they are usually not tied to modern AI frameworks, such as
PyTorch [11]. Such frameworks have become indispensable tools for AI
researchers to build new ANN architectures and implement AI workloads. For
instance, while MLP-NeuroSim [2] has the ability to handle a great number of
realistic hardware models, it is built around a small (by today’s standard)
MNIST network example, hence lacking the tools to flexibly implement new ANN
topologies or simulate larger DNNs. RxNN [9] is based on the rather outdated
Caffe framework and only caters to analog chips dedicated to inference,
lacking more advanced algorithms and pulse update schemes that are needed for
training-enabled chip designs. A recent PyTorch re-implementation of a subset
of the MLP-NeuroSim package using PyTorch, called DNN+NeuroSim [12], is
closest to our framework. However, it does not support fast custom CUDA
kernels for in-situ parallel analog update, and also lacks many of the aspects
of modern code development, such as modularity, object oriented programming,
proper Python binding of fast C++ elements, unit testing, easy Python package
installation and so on. To address these limitations, we decided to open
source IBM’s internal analog AI simulation core with a convenient PyTorch
interface, hence providing adequate software tools for the advancement of
algorithmic analog AI.

##  2 Toolkit overview

We provide a professionally maintained software kit that is easy to install
(as simple as `pip install aihwkit`), has an extensive
documentation222https://aihwkit.readthedocs.io/en/latest/ with a growing
number of examples and tutorials, and has a well integrated PyTorch experience
providing a number of new `Analog` layers, such as linear (fully-connected)
layer and convolutions (see Fig. 1). Under the hood, these layers perform the
analog hardware simulations using a highly optimized C++/CUDA core library
(RpuCuda) that is part of the package and made available under an Apache-2.0
Open Source license. The core library is bound to the Python layer so that no
re-compilation is needed for defining models and parameters. Instead, all C++
functions and analog hardware defining parameters are conveniently accessible
from Python. Figure 1 illustrates the architecture and code structure of the
toolkit.

Among others, we provide an `AnalogLinear` layer, which uses a single analog
tile to compute the forward pass of a fully-connected layer within a DNN.
Defining a simple analog training simulations (with default device parameters)
is similar to how PyTorch to define DNN layers (see example in Fig. 2).

⬇

# Define crossbar (RPU) config

rpu_config = SingleRPUConfig(device=ReRamESPresetDevice())

# Define a single-layer network

model = AnalogLinear(4, 2, bias=True, rpu_config=config)

# Define analog-aware optimizer

opt = AnalogSGD(model.parameters(), lr=0.1)

opt.regroup_param_groups(model)

# Run the training

for epoch in range(100):

pred = model(x) # forward pass

loss = mse_loss(pred, y)

loss.backward() # backward pass

opt.step() # (analog pulsed) update

Figure 2: Simple network definition and training with parallel pulsed update.
Note that imports and data preparation are omitted for brevity (see example 01
in our GitHub repository).

##  3 Analog tile model simulation

The toolkit is centered around an “analog tile” that corresponds to a 2D
weight matrix that is stored on a non-volatile resistive crossbar array and
includes most pre- and post-processing (such as dynamic input scaling and
analog-digital converters) that are based on (and extends) the resistive
processing unit (RPU) definitions in [5]. The basic operation of this tile is
a matrix-vector product $y_{i}=\sum_{j}w_{ij}x_{j}$, where the matrix elements
$w_{ij}$ are assumed to be stored in the crossbar (optionally together with
the biases). The analog matrix-vector product is, however, corrupted by noise
and non-idealities that are defined by the user through parameter
configurations of the tile during creation (called `rpu_config`).

For instance, input/output/weight noise processes, analog-to-digital
conversion resolutions, or other non-idealities can be added so that the
effective computation simulated is e.g.333Consult the online documentation for
a complete overview at https://aihwkit.readthedocs.io.  given by

$y_{i}=f_{\text{adc}}\big{(}\sum_{j}(w_{ij}+\sigma_{\text{w}}\xi_{ij})(f_{\text{dac}}(x_{j})+\sigma_{\text{inp}}\xi_{j})+\sigma_{\text{out}}\xi_{i}\big{)}$
(1)

where $f_{\text{adc}}$ and $f_{\text{dac}}$ are discretization functions
(together with e.g. dynamic scaling and clipping) of the digital-analog
conversion process (as given by the hardware design), and $\xi$ are Gaussian
noise processes. In our tool, we generally opt for a more abstract functional
representation of the underlying hardware (i.e. largely working with
normalized units), since in this manner the cause of accuracy loss of
different hardware designs become more apparent and specifications can be
derived based on parameters relevant to the DNN.

The RPU C++ library is the core of the high-performing implementation of Eq.
(1), which extensively uses dedicated and fused GPU kernels. It provides the
above basic model, however, internal pre- and post processing (e.g. analog to
digital discretization) can be simply turned off or adapted by the user using
PyTorch’s functionality.

Within a DNN that is considered to be accelerated with analog crossbar arrays,
we generally assume that digital operations can be separated from analog
operations. For instance, by default, activation functions are performed with
digital computing using PyTorch’s library of floating point functions and
assuming that analog result values of a matrix vector product are digitized.
Moreover, we currently support only parallel read-out schemes and linear
analog-digital conversion. However, iterative read-out schemes and shift-and-
add for more complicated analog-digital converters could be straightforwardly
added in a similar manner using PyTorch.

While the forward pass includes cycle-to-cycle variations of the resistive
materials and circuit noise sources (see Eq. 1), for chips dedicated to ANN
training, backward and (potentially) update passes need to be similarly
simulated to be performed on analog as well. Thus, in our framework, backward
pass (the transposed version of the forward pass) can be configured to contain
noise and non-linearity similar to the forward pass. In general, these
parameters do not need to be the same, to reflect hardware designs that differ
in the noisiness in both directions.

A particular strength of our new software framework is the highly-optimized
pulse update with tune-able material pulse response properties for chips that
are dedicated to perform ANN training. For that, we consider many algorithmic
advances [5, 3] and use stochastic (or deterministic) pulse trains to
incrementally perform the theoretical rank update

$w_{ij}^{(t+1)}=w_{ij}^{(t)}+\lambda d_{i}x_{j}$ (2)

where $d_{i}$ are the back-propagated error signals. One possible way to
implement the update pass is to use stochastic pulse sequences of variable or
fixed length, where the probability of having a pulse is proportional to the
size of $x_{j}$ and $d_{i}$. In case of pulse coincidences, the weight element
at crosspoint $ij$ is updated by a finite amount $\Delta w_{ij}$ (see [5]).

Importantly, due to physical constraints, the amount $\Delta w_{ij}$ can vary
in a systematic way from crosspoint to crosspoint (device-to-device
variations) or might have a systematic bias in one or the other direction (up
versus down pulses). Moreover, the step $\Delta w_{ij}$ might depend in a non-
linear way on the actual conductance of the crosspoint value, depending on the
physical properties of the material used as resistive element.

We thus provide a flexible way to choose the correct update response and fit
the update behavior to the material properties under investigation, such as
the amount of device-to-device variations in various parameters and additional
pulse-to-pulse variations. We also provide a number of presets calibrated on
hardware data (e.g. for ReRAM [6] characteristics, see Fig. 3). Example
response curves are shown in Fig. 3.

Note that our framework differs significantly from the update approach
implemented by DNN+NeuroSim. For instance, DNN+NeuroSim uses native PyTorch
operators to compute the accumulated gradient, which in case of batch size
larger than one or generally for convolutions will accumulate $p$ outer
products in digital, and then only discretize the accumulated result to
calculate the pulse number and add pulse variations. This means that the
outer-product is actually done in digital and not analog. This will result in
a dramatic under-estimation of the update noise properties for a system that
performs the rank-1 update in place in parallel in analog (as suggested by
[5]). Thus, we re-implemented convolution and linear operators in our C++ core
to ensure that gradient accumulation happens using parallel update in analog
memory even for batch sizes larger than 1 or in case of convolutions.

Moreover, our framework draws the stochastic pulse trains and applies the
update pulses (in case of a coincidence) sequentially on the devices, which
are then updated according to the device model pulse response curve. These
pulsed updates are highly optimized for GPU compute, where we provide hundreds
of custom CUDA kernel templates that are automatically tuned for each tile,
thus selecting the fastest template on the fly. Consequently, the full analog
training simulation often takes only 2-5x longer than the conventional
floating-point training444Depending on DNN and device setting, e.g. 60s per
epoch for VGG-8/CIFAR10 with parallel pulsed update versus 15s with floating
point on a V100 GPU.

Figure 3: A) Schematic illustration of ReRAM and PCM devices. B) Pulse
response of the experimental [6] and simulated ReRAM device (with device-to-
device variations, write noise and, cycle-to-cycle variations). C)
Experimentally obtained temporal evolution of PCM conductance [10] compared to
that simulated by the statistical noise model of our inference tile. Note that
we assume currently that all weights are programmed at the same time in the
simulation, whereas in the experiment devices converge at different iterations
of programming and there will be jitter in times at which drift starts in
individual devices. This effect leads to the slight difference noticeable for
small times (less than 1000 seconds) between model and experiment.

##  4 Advanced device models and optimizers

All the tile parameters are chosen by setting the resistive device
configuration. The default is having a single device per crosspoint. However,
in many hardware designs, it might be advantageous to have multiple conductive
elements within a “unit cell”. Therefore, we provide a mechanism to define
unit cells having an arbitrary number of devices (each potentially having
different material properties), which can be updated one-by-one or all
together. This allows defining two uni-directional devices per cell, a multi-
device cell or a bi-directional single device cell with reference devices, by
configuring the parameters of our analog tile accordingly.

Moreover, some analog optimization schemes have been proposed, where each
analog tile of SGD is replaced with two analog tiles to create a coupled
system by exchanging information between the two [3] tiles. We provide such
constructs that can be configured to test different and new analog inspired
optimization algorithms (see Fig. 4).

Lastly, a number of temporal noise processes are available, such as
conductance decay, reset, and diffusion, which are all implemented with
systematic device-to-device variability. These temporal processes are applied
once per mini-batch and are thus assumed to be negligibly small within the
physical time a one mini-batch computation takes.

##  5 Training for analog inference chips

It has become increasingly apparent that accuracy on chips that are designed
for inference benefit greatly from “hardware-aware” training, where the DNN
workload is trained (or fine-tuned) in software with some analog noise and
non-idealities in the forward pass, but with otherwise perfect backward and
update pass [4, 14, 8, 10]. Our framework supports hardware-aware training for
inference. In particular, backward and update of a tile can simply be
configured as being “perfect”. Moreover, we provide weight noise functions for
inference training, that reversibly add noise onto weights during a mini-batch
(in forward and backward, but not during update pass) to train for noise
resiliency.

Finally, we also provide a specific inference “tile” setting, which supports
adding carefully calibrated, conductance-dependent programming noise, weight
read noise and conductance drift onto a trained networks’ analog weights
during inference time and provides automatic global drift compensation (see
eg. [10]). The noise strength and temporal profile is calibrated on a 1M
phase-change memory (PCM) device array [10] (see Fig. 3).

Thus, any network can be tested for realistic inference accuracy and how the
accuracy degrades over time.

⬇

# Define more complicated crossbar (RPU) config

rpu_config = UnitCellRPUConfig(

device=TransferCompound(

# Devices that compose the Tiki-taka compound.

unit_cell_devices=[

ReRamSBPresetDevice(dw_min_dtod=0.1),

ReramSBPresetDevice(dw_min_std=0.2)],

# Some adjustments to how to perform Tiki-Taka

units_in_mbatch=True,

transfer_every=2))

Figure 4: A more complex device configuration. This example implements the
Tiki-taka modified SGD rule[3] for analog training. Once the rpu_config is
defined, the DNN training is identical to Fig. 2.

##  6 Extendability and open source contributions

The reason of using PyTorch as an AI framework is to have a system to
flexibility define arbitrary DNN architectures. Our simulator inherits this
flexibility in defining networks with both analog and digital components,
where for the latter all features of PyTorch can be used. Moreover, because of
the modular structure of our code base, it is easy to extend and incorporate
new features. For instance, the inference noise models use Python classes and
could be easily derived to implement other drift characteristics.

Device update behavior, on the other hand, is internally computed in C++. This
core C++ RpuCuda library is designed for speed but retains a modular structure
to make it also extendable. While a number of abstract device models are
already provided, including a number of parameter presets, new device models
can be implemented relatively easily by inheriting most properties of existing
devices. In particular, if only the update response curve of a new material
needs to be changed, many optimized kernel routines can be re-used, as only
the update function itself needs to be written and added to the code (in C++
for CPU and GPU).

In general, we follow an open source development model using GitHub standards,
where we encourage contributions and discussions via an issue tracker, while
using modern industry best-practices for software development (continuous
integration, automated unit testing, code reviews) and packaging (semantic
versioning, documentation, PEP adherence).

##  7 Conclusion

We present a new flexible open source software tool for simulating the
functional aspects of analog AI with a clear features roadmap. It is solely
focused on the algorithmic development and functional verification for ANN
training and inference on emerging analog chips. This toolkit is _not_
intended to estimate run time, latency, or power performance of a hardware
chip. While this is very important in its own right, performance aspects are
(largely) irrelevant to the functional characterization and need dedicated and
complimentary tools. We hope that this toolkit will make it easier for
researchers from AI hardware development and deep learning communities, to
develop new algorithms and design new ANN topologies that are optimized for
analog AI, hence experiencing the exciting new area of analog AI.

## References

  * [1] G. W. Burr et al.  Neuromorphic computing using non-volatile memory.  Advances in Physics: X, 2(1):89–124, 2017.
  * [2] P.-Y. Chen et al.  NeuroSim+: An integrated device-to-algorithm framework for benchmarking synaptic devices and array architectures.  In 2017 IEEE International Electron Devices Meeting (IEDM), pages 6–1, 2017.
  * [3] T. Gokmen and W. Haensch.  Algorithm for training neural networks on resistive device arrays.  Frontiers in Neuroscience, 14, 2020.
  * [4] T. Gokmen, M. J. Rasch, and W. Haensch.  The marriage of training and inference for scaled deep learning analog hardware.  In 2019 IEEE International Electron Devices Meeting (IEDM), pages 22–3. IEEE, 2019.
  * [5] T. Gokmen and Y. Vlasov.  Acceleration of deep neural network training with resistive cross-point devices: Design considerations.  Frontiers in neuroscience, 10:333, 2016.
  * [6] N. Gong and other.  Signal and noise extraction from analog memory elements for neuromorphic computing.  Nature communications, 9(2102), 2018.
  * [7] W. Haensch, T. Gokmen, and R. Puri.  The next generation of deep learning hardware: Analog computing.  Proceedings of the IEEE, 107(1):108–122, 2018.
  * [8] S. Jain et al.  Neural network accelerator design with resistive crossbars: Opportunities and challenges.  IBM Journal of Research and Development, 63(6):10–1, 2019.
  * [9] S. Jain, A. Sengupta, K. Roy, and A. Raghunathan.  RxNN: A framework for evaluating deep neural networks on resistive crossbars.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 2020.
  * [10] V. Joshi et al.  Accurate deep neural network inference using computational phase-change memory.  Nature Communications, 11(2473), 2020.
  * [11] A. Paszke et al.  Pytorch: An imperative style, high-performance deep learning library.
  * [12] X. Peng, S. Huang, H. Jiang, A. Lu, and S. Yu.  DNN+NeuroSim V2.0: An end-to-end benchmarking framework for compute-in-memory accelerators for on-chip training.  arXiv:2003.06471, 2020.
  * [13] A. Sebastian, M. Le Gallo, R. Khaddam-Aljameh, and E. Eleftheriou.  Memory devices and applications for in-memory computing.  Nature Nanotechnology, 15:529–544, 2020.
  * [14] H. Tsai et al.  Inference of long-short term memory networks at software-equivalent accuracy using 2.5 M analog phase change memory devices.  In 2019 Symposium on VLSI Technology, pages T82–T83. IEEE, 2019\.