# When ideas go viral - complex bifurcations
in a two-stage transmission model

J. HEIDECKE and M. V. BARBAROSSA  Frankfurt Institute for Advanced Studies,
Ruth-Moufang-Straße 1
60438 Frankfurt, Germany
<EMAIL_ADDRESS>
###### Abstract

We consider the qualitative behavior of a mathematical model for transmission
dynamics with two nonlinear stages of contagion. The proposed model is
inspired by phenomena occurring in epidemiology (spread of infectious
diseases) or social dynamics (spread of opinions, behaviors, ideas), and
described by a compartmental approach. Upon contact with a promoter
(contagious individual), a naive (susceptible) person can either become
promoter himself or become weakened, hence more vulnerable. Weakened
individuals become contagious when they experience a second contact with
members of the promoter group. After a certain time in the contagious
compartment, individuals become inactive (are insusceptible and cannot spread)
and are removed from the chain of transmission. We combine this two-stage
contagion process with renewal of the naive population, modeled by means of
transitions from the weakened or the inactive status to the susceptible
compartment. This leads to rich dynamics, showing for instance coexistence and
bistability of equilibria and periodic orbits. Properties of (nontrivial)
equilibria are studied analytically. In addition, a numerical investigation of
the parameter space reveals numerous bifurcations, showing that the dynamics
of such a system can be more complex than those of classical epidemiological
ODE models.

##  1 Introduction

Social contagion is the spread of behaviors or attitudes through (physical or
virtual) groups of people [1]. From a mathematical point of view, modeling
social contagion in large communities is very similar to modeling the
transmission of an infectious disease in a population. Hence, it seems natural
that methods from the field of mathematical epidemiology, such as
compartmental models [2, 3, 4, 5, 6, 7, 8], are used to model social contagion
phenomena. In certain cases social contagion and disease spread even have to
be considered coupled to one another, as when a group in a social network
criticizes vaccination [9]. Despite the analogies, social contagion differs
from biological contagion in various aspects. For example, intellectual
epidemics could be advantageous [3], ideas do not require interpersonal
contact to spread [10], or people might be asked to choose between opposite
opinions [7]. Thus, to mathematically describe social contagion processes,
models from theoretical epidemiology might require adaptation to the specific
context.

The classical SEIR (susceptible-exposed-infective-recovered) model in
mathematical epidemiology describes the transmission of an infectious disease
in a population [11]. When a susceptible individual comes in contact with an
infective one, there is a certain probability that contagion occurs and the
susceptible moves to the exposed compartment. After a latent period exposed
individuals become infectious themselves and can infect others. Once the
infectious period is over, the individual recovers, cannot transmit the
disease to others anymore and becomes immune. Waning of immunity, i.e.,
transitions from the recovered to the susceptible compartment, is possible for
certain diseases [11, 12].

In the context of social contagion, the spreading of a specific behavior or
opinion in a population could be described as the transition of individuals
from the ”naive” (susceptible) status to the ”promoter” (infectious) one. This
transition might require several steps and depend on repeated exposure to
promoters [8, 13]. Therefore, we classify individuals as:

  * •

naive/susceptible ($S$), those who have not yet been exposed to the considered
behavior/opinion,

  * •

weakened ($W$), those who came in contact with the considered
behavior/opinion, but are not yet spreading to others

  * •

promoters/infectious ($I$), those who have embraced the considered
behavior/opinion and are able to transmit it to others

  * •

inactive/resistant ($R$), those who have been sharing the considered
behavior/opinion earlier but are now neither transmitting to others nor can be
re-exposed.

In contrast to the classical SEIR approach, transition from the
exposed/weakened stage to the promoter/infected stage depends on contacts with
infectives/promoters. Promoters and susceptibles make contact sufficient to
transmit the opinion/behavior at rate $\beta$. We assume that upon such a
contact the susceptibles have a certain probability $\rho\in(0,1)$ to enter
the $I$-compartment directly, becoming a promoter themselves (perfect
contact). With probability $1-\rho$ susceptible individuals enter the
$W$-compartment (imperfect contact). When a $W$-individual comes in contact
with promoters, the weakened enters the $I$-compartment at rate
$\beta_{2}\geq\beta$. In this sense, an individual in the $W$-compartment is
more vulnerable to the opinion/behavior than an individual in the
$S$-compartment. Over time promoters might reject the considered
opinion/behavior. Thus, we assume that promoters leave the promoting class at
rate $\alpha$ and become inactive/resistant. As we assume that the
inactivity/resistance of a $R$-individual towards the considered
opinion/behavior might wane over time, we allow transition from $R$ back to
$S$ at rate $\eta$. Moreover, we assume that the infection can fade away in
weakened individuals, whereby transitions from $W$ to $S$ are occuring at rate
$\gamma$. We assume that $\gamma\geq\alpha$, suggesting that the average time
an individual is able to promote a certain opinion is not shorter than his
exposure time.

Figure 1: Flowchart of the two-stage contagion model (1) with waning of
immunity and fading of infection in weakened individuals.

A model sketch is given in Fig. 1 and the corresponding differential equations
system is

$\displaystyle S^{\prime}$ $\displaystyle=-\beta IS+\gamma W+\eta R$ (1)
$\displaystyle W^{\prime}$ $\displaystyle=(1-\rho)\beta IS-\beta_{2}IW-\gamma
W$ $\displaystyle I^{\prime}$ $\displaystyle=\beta_{2}IW+\rho\beta IS-\alpha
I$ $\displaystyle R^{\prime}$ $\displaystyle=\alpha I-\eta R.$

Limit cases of this system lead, on the one hand, to the standard SIRS model
(cf. [11]) if $\rho=1$ and the $W$-compartment is empty at the beginning of
observations. On the other hand, if the transition from $W$ to $I$ would be a
linear one and $\gamma=\rho=0$, the model would be equivalent to the classical
SEIRS model (cf. [11]).

Observing that the total population $N=S+W+I+R$ does not vary over time, we
set $N\equiv 1$. Clearly, system (1) has a unique non-negative global solution
for every choice of non-negative initial values, and the set
$\\{(S,I,W,R)\in\mathbb{R}^{4}_{\geq 0}\,|\,S+W+I+R=1\\}$ is forward
invariant. For analytical simplicity, we restrict to the case
$\beta_{2}=\beta$. Using the conservation relation $W=1-S-I-R$, we consider
the reduced system

$\displaystyle S^{\prime}$ $\displaystyle=-\beta IS+\gamma(1-S-I-R)+\eta R$
(2) $\displaystyle I^{\prime}$ $\displaystyle=\rho\beta IS-\alpha I+\beta
I(1-S-I-R)$ $\displaystyle R^{\prime}$ $\displaystyle=\alpha I-\eta R.$

A similar compartmental model for two-stage contagion was previously proposed
by Guy Katriel [8]. System (2) differs from Katriel’s work in two aspects.
First we include waning immunity and fading of infection, that is transitions
from $R$, respectively $W$, to the susceptible compartment. Second, we do not
consider population demography (births/deaths). Such differences lead to major
analytical challenges with respect to Katriel’s study [8]. In the rest of this
work we study the qualitative properties of system (2) by means of analytical
and numerical methods.

##  2 Existence and local stability of equilibria

To understand the long-term behavior of the system we first investigate its
equilibria. The criteria on the stability of the disease-free equilibrium,
where the $I$ compartment is empty, are commonly related to the so-called
basic reproduction number, $\mathcal{R}_{0}$. We remark that the following
results are derived assuming that $\gamma\geq\alpha$ holds.

######  Theorem 2.1.

  1. a)

System (2) has a unique disease-free (DFE) equilibrium
$\mathcal{E}_{0}:=(1,0,0)$, which exists for any choice of
$\beta,\,\gamma,\,\eta,\,\alpha>0$, $\rho\in(0,1)$.

  2. b)

The DFE $\mathcal{E}_{0}$ is locally asymptotically stable (shortly, LAS) if
$\mathcal{R}_{0}:=\frac{\rho\beta}{\alpha}<1$, and unstable if
$\mathcal{R}_{0}>1$.

###### Proof.

a): Equilibria of system (2) are determined setting the right-hand side of the
system equal to zero,

$\displaystyle 0=$ $\displaystyle-\beta
I^{*}S^{*}+\gamma(1-S^{*}-I^{*}-R^{*})+\eta R^{*}$ (3) $\displaystyle 0=$
$\displaystyle I^{*}\left(\rho\beta
S^{*}-\alpha+\beta(1-S^{*}-I^{*}-R^{*})\right)$ $\displaystyle 0=$
$\displaystyle\alpha I^{*}-\eta R^{*}.$

If $I^{*}=0$ (DFE condition), then the third relation in (3) implies that
$R^{*}=0$ and the first relation in (3) implies that $S^{*}=1$. Thus,
$\mathcal{E}_{0}$ is the unique DFE of the system. b): For studying the local
stability of the DFE we linearize about $\mathcal{E}_{0}$. The Jacobian matrix
of (2) is

$J=\left({\begin{array}[]{ccc}-\beta I-\gamma&-\beta S-\gamma&\eta-\gamma\\\
\left(\rho-1\right)\beta
I&\beta\left(\left(\rho-1\right)S+1-2I-R\right)-\alpha&-\beta I\\\
0&\alpha&-\eta\\\ \end{array}}\right).$ (4)

Evaluation of $J$ at $\mathcal{E}_{0}$ yields the eigenvalues
$\lambda_{1}=-\gamma<0$, $\lambda_{2}=\rho\beta-\alpha$ and
$\lambda_{3}=-\eta<0$. The condition $\lambda_{2}<0\iff\mathcal{R}_{0}<1$
guarantees local stability of the DFE. ∎

The above expression for $\mathcal{R}_{0}$ can be interpreted as the number of
secondary infections that a single promoter introduced into a completely
susceptible population produced through perfect contacts ($\rho\beta$) over
the duration of its promoting period ($1/\alpha$). The secondary infections
produced through imperfect contacts require two nonlinear transitions over the
$W$-compartment, which is empty at $\mathcal{E}_{0}$, hence do not contribute
to determining the local stability of the DFE.

Theorem 2.1 provides conditions for the local stability of the DFE. What
exactly happens at the bifurcation value $\mathcal{R}_{0}=1$ can be determined
with the help of Theorem A.1, first introduced by Castillo-Chavez and Song
[14]. For convenience of notation we define
$\kappa:=\left(1+\alpha/\eta\right)$.

######  Theorem 2.2.

Let $\rho^{*}:=\frac{-\alpha+\sqrt{\alpha\gamma\kappa}}{\gamma\kappa-\alpha}$.
If $\rho>\rho^{*}$, then a transcritical bifurcation of forward type occurs at
$\mathcal{R}_{0}=1$. If $\rho<\rho^{*}$, then a transcritical bifurcation of
backward type occurs at $\mathcal{R}_{0}=1$.

###### Proof.

We use Theorem A.1 and take $\beta$ as our bifurcation parameter, with
$\beta^{*}:=\frac{\alpha}{\rho}$ corresponding to $\mathcal{R}_{0}=1$. The
Jacobian (4) of system (2) evaluated at $\mathcal{E}_{0}$ for
$\beta=\beta^{*}$, $\mathcal{A}:=J|_{\mathcal{E}_{0},\beta^{*}}$, has a simple
zero eigenvalue and two eigenvalues with negative real part. To compute a
right eigenvector $w$ of $\mathcal{A}$ corresponding to the zero eigenvalue,
solve $\mathcal{A}w=0$, which is an underdetermined system. We fix $w_{2}=1$
and obtain $w_{3}=\frac{\alpha}{\eta}w_{2}=\frac{\alpha}{\eta}$ and
$w_{1}=-\frac{\alpha}{\rho\gamma}-1+\frac{\alpha}{\gamma}-\frac{\alpha}{\eta}$.
Analogously, to find a left eigenvector, we solve the system $v\mathcal{A}=0$,
obtaining $v=(0,1,0)$.

Next, we denote the vector field of system (2) by $f=(f_{1},f_{2},f_{3})$ and
we calculate second order partial derivatives to verify the conditions of
Theorem A.1. As $v_{1}=0=v_{3}$ the derivatives of $f_{1}$ and $f_{3}$ are not
needed. All second order derivatives of $f_{2}$ are zero, except for

$\displaystyle\frac{\partial^{2}f_{2}}{\partial^{2}I}$
$\displaystyle=-2\beta^{*},$ $\displaystyle\frac{\partial^{2}f_{2}}{\partial
I\partial R}$ $\displaystyle=-\beta^{*},$
$\displaystyle\frac{\partial^{2}f_{2}}{\partial S\partial I}$
$\displaystyle=\rho\beta^{*}-\beta^{*},$
$\displaystyle\frac{\partial^{2}f_{2}}{\partial I\partial\beta^{*}}$
$\displaystyle=\rho.$

Therefore, the quantities $a$ and $b$ in Theorem A.1 are given by

$b=\sum_{i=1}^{3}v_{2}w_{i}\frac{\partial^{2}f_{2}}{\partial
x_{i}\partial\beta}(0,0)=\rho v_{2}w_{2}=\rho>0,$

and

$\displaystyle a=$
$\displaystyle\sum_{i,j=1}^{3}v_{2}w_{i}w_{j}\frac{\partial^{2}f_{2}}{\partial
x_{i}\partial x_{j}}(0,0)$ $\displaystyle=$
$\displaystyle-2\frac{\alpha}{\rho}v_{2}w_{2}^{2}-2\frac{\alpha}{\rho}v_{2}w_{2}w_{3}+2\left(\alpha-\frac{\alpha}{\rho}\right)v_{2}w_{1}w_{2}$
$\displaystyle=$
$\displaystyle-2\frac{\alpha}{\rho}-2\frac{\alpha^{2}}{\rho\eta}+2\left(\alpha-\frac{\alpha}{\rho}\right)\left(-\frac{\alpha}{\rho\gamma}-1+\frac{\alpha}{\gamma}-\frac{\alpha}{\eta}\right),$

where $x=(x_{1},x_{2},x_{3})=(S,I,R)$. The condition $a>0$ is equivalent to

$\displaystyle\psi(\rho):=\rho^{2}\left(\gamma\kappa-\alpha\right)+2\rho\alpha-\alpha<0.$
(5)

Because of the assumption $\gamma\geq\alpha$, we have
$\gamma\kappa=\gamma(1+\frac{\alpha}{\eta})>\alpha$. The open up parabola
$\psi(\rho)$ has negative intercept and vertex on the left half-plane. Thus,
condition (5) is fulfilled only if $0<\rho<\rho^{*}$, where

$\rho^{*}=\frac{-\alpha+\sqrt{\alpha\gamma\kappa}}{\gamma\kappa-\alpha}$

is the positive zero of $\psi(\rho)$. ∎

Endemic equilibria of (2) are determined by solving (3) for $I^{*}\neq 0$.
From the third and the first relation in (3) we get

$R^{*}=\frac{\alpha}{\eta}I^{*}\;\mbox{ and }\;S^{*}=\frac{\gamma(1-\kappa
I^{*})+\alpha I^{*}}{\beta I^{*}+\gamma}.$

Note here that $1-I^{*}-\frac{\alpha}{\eta}I^{*}=1-I^{*}-R^{*}\geq 0$, hence
$S^{*}\geq 0$ if $I^{*}\geq 0$. The equilibria condition thus reduces to the
quadratic equation

$\displaystyle 0=a{I^{*}}^{2}+bI^{*}+c=:\Phi(I^{*}),$ (6)

with $a:=\beta\kappa>0$, $b:=\alpha\left(2-\rho\right)+\rho\gamma\kappa-\beta$
and $c:=\gamma\left(\frac{\alpha}{\beta}-\rho\right)$. Solutions to (6) are
given by

$I^{*}_{\pm}=\frac{1}{2\beta\kappa}\left[\beta-\rho\gamma\kappa-\alpha\left(2-\rho\right)\pm\sqrt{\Delta}\right],$

where $\Delta$ is the discriminant of (6). Hence, system (2) has at most two
endemic equilibria
$\mathcal{E}_{\pm}:=\left(S_{\pm}^{*},I_{\pm}^{*},R_{\pm}^{*}\right)$, whereby
$I^{*}_{\pm}$ are the positive roots of the open up parabola $\Phi(I^{*})$. We
write the coefficients $b$ and $c$ of $\Phi(I^{*})$ as functions of $\beta$,

$\displaystyle b(\beta):=$
$\displaystyle\alpha(2-\rho)+\rho\gamma\kappa-\beta$ $\displaystyle
c(\beta):=$ $\displaystyle\gamma\left(\frac{\alpha}{\beta}-\rho\right).$

Note that $b(\beta)$ is a strictly decreasing linear function with zero at

$\displaystyle\beta_{b}:=\alpha(2-\rho)+\rho\gamma\kappa$

and $c(\beta)$ is a strictly decreasing curve with zero at

$\displaystyle\beta_{c}:=\frac{\alpha}{\rho}.$

Observe further that

$c(\beta)\left\\{\begin{matrix}>0\iff\mathcal{R}_{0}<1,\\\
=0\iff\mathcal{R}_{0}=1,\\\ <0\iff\mathcal{R}_{0}>1.\end{matrix}\right.$

Theorem 2.2 suggest to consider the cases $\rho^{*}<\rho$ and $\rho^{*}>\rho$
separately.

######  Theorem 2.3.

If $\rho^{*}<\rho<1$, then:

  * •

For $\mathcal{R}_{0}\leq 1$ there are no endemic equilibria.

  * •

For $\mathcal{R}_{0}>1$ there is a unique endemic equilibrium
$\mathcal{E}_{+}$.

At $\mathcal{R}_{0}=1$ a transcritical bifurcation of forward type occurs, and
a branch of stable endemic equilibria $\mathcal{E}_{+}$ emerges from
$\mathcal{E}_{0}$.

###### Proof.

Observe that

$\beta_{b}>\beta_{c}\,\iff\,2\rho\alpha-\alpha\rho^{2}+\rho^{2}\gamma\kappa>\alpha\;\iff\,\psi(\rho)>0\,\underbrace{\iff}_{\mbox{Thm.~{}\ref{theo:bifdir}}}\,\rho^{*}<\rho.$
(7)

Thus, if $\rho^{*}<\rho$ we have that the zero of $b(\beta)$ lies on the right
of the zero of $c(\beta)$. In other words, if $c(\beta)\geq 0$ (that is
$\mathcal{R}_{0}\leq 1$) then necessarily $b(\beta)>0$. This means that if the
open up parabola $\Phi(I^{*})$ has a positive y-intercept, then it also has
vertex on the left half-plane, hence no positive root $I^{*}$. If $c(\beta)<0$
(that is $\mathcal{R}_{0}>1$), $\Phi(I^{*})$ has a negative y-intercept, thus
only its larger root $I^{*}_{+}$ is positive. The rest follows from Theorem
2.2. ∎

Next, define

$\mathcal{R}_{C}:=\rho(2-\rho)-\rho^{2}\frac{\gamma\kappa}{\alpha}+2\frac{\rho(1-\rho)}{\alpha}\sqrt{\alpha\gamma\kappa}$
(8)

and consider the case $\rho^{*}>\rho$.

######  Theorem 2.4.

If $0<\rho<\rho^{*}$, then the following results hold:

  1. a)

If $\mathcal{R}_{0}<\mathcal{R}_{C}$, there are no endemic equilibria.

  2. b)

If $\mathcal{R}_{C}\leq\mathcal{R}_{0}<1$, there are two endemic equilibria
$\mathcal{E}_{\pm}$, which coincide if $\mathcal{R}_{C}=\mathcal{R}_{0}$.

  3. c)

If $1\leq\mathcal{R}_{0}$, there is a unique endemic equilibrium
$\mathcal{E}_{+}$.

At $\mathcal{R}_{0}=1$ a transcritical bifurcation of backward type occurs,
and a branch of unstable endemic equilibria $\mathcal{E}_{-}$ emerges from
$\mathcal{E}_{0}$.

###### Proof.

From $\rho<\rho^{*}$ it follows from (7) that $\beta_{b}<\beta_{c}$. Hence if
$c(\beta)\leq 0$ (that is $\mathcal{R}_{0}\geq 1$) then necessarily
$b(\beta)<0$. This means that if the open up parabola $\Phi(I^{*})$ has a non-
positive y-intercept, it has vertex on the right half-plane as well, and hence
a unique positive root $I_{+}^{*}$. This proofs statement c). For two positive
roots $I_{\pm}^{*}$ the conditions $c(\beta)>0$, $b(\beta)<0$ and
$\Delta(\beta)\geq 0$, where

$\displaystyle\Delta(\beta)=$
$\displaystyle\beta^{2}+\beta{\left(2\rho\alpha+2\rho\gamma\kappa-4\alpha\right)}$
$\displaystyle\;+\underbrace{4\alpha^{2}\left(1-\rho\right)+\rho^{2}\alpha^{2}+\rho^{2}\gamma^{2}\kappa^{2}+4\rho\alpha\gamma\kappa-2\rho^{2}\alpha\gamma\kappa-4\alpha\gamma\kappa}_{:=c_{\Delta}},$

are necessary. For $\Delta(\beta)=0$ the two roots coincide. The discriminant
$\Delta(\beta)$ as a function of $\beta$ is an open up parabola itself. Short
computation shows that

$c_{\Delta}<0\;\iff\;\rho<2\rho^{*},$

hence $c_{\Delta}$ is always negative under the assumptions of the theorem.
Thus $\Delta(\beta)\geq 0$ for $\beta$ larger than the positive zero of the
discriminant-parabola, that is

$\displaystyle\beta$
$\displaystyle\geq\alpha(2-\rho)-\rho\gamma\kappa+2(1-\rho)\sqrt{\alpha\gamma\kappa}=:\beta_{\Delta},$

and $\Delta(\beta)=0\iff\beta=\beta_{\Delta}$. Observe that

$\displaystyle\beta_{b}<\beta_{\Delta}$
$\displaystyle\iff\;\rho^{2}\gamma^{2}\kappa^{2}<(1-\rho)^{2}\alpha\gamma\kappa$
(9) $\displaystyle\iff\;\rho^{2}(\gamma\kappa-\alpha)+\rho 2\alpha-\alpha<0$
$\displaystyle\iff\;\rho<\rho^{*}.$

Further

$\displaystyle\beta_{\Delta}<\beta_{c}$
$\displaystyle\iff\alpha(2-\rho)-\rho\gamma\kappa+2(1-\rho)\sqrt{\alpha\gamma\kappa}<\frac{\alpha}{\rho}$
$\displaystyle\iff\underbrace{\rho^{2}\left(-\alpha-\gamma\kappa-2\sqrt{\alpha\gamma\kappa}\right)+\rho
2\left(\alpha+\sqrt{\alpha\gamma\kappa}\right)-\alpha}_{=:\chi(\rho)}<0.$

Note that $\chi(\rho)$ is an open down parabola with vertex on the right half-
plane, negative y-intercept, and zero discriminant. Therefore, $\chi(\rho)$
has a unique zero at

$\displaystyle\rho=\frac{(\alpha+\gamma\kappa+2\sqrt{\alpha\gamma\kappa})(-\alpha+\sqrt{\alpha\gamma\kappa})}{(\gamma\kappa-\alpha)(\alpha+\gamma\kappa+2\sqrt{\alpha\gamma\kappa})}=\rho^{*}.$

Hence,

$\rho<\rho^{*}\;\implies\;\chi(\rho)<0\;\iff\;\beta_{\Delta}<\beta_{c}.$

With (9) this means that for $\rho<\rho^{*}$ it holds that
$\beta_{b}<\beta_{\Delta}<\beta_{c}$. As the positive root $\beta_{\Delta}$ of
$\Delta(\beta)$ lies on the right of the zero of $b(\beta)$, then for
$\beta\geq\beta_{\Delta}$ we have $\Delta(\beta)\geq 0$ and necessarily
$b(\beta)<0$. As $c(\beta)>0$ (that is $\mathcal{R}_{0}<1$) for
$\beta<\beta_{c}$, the condition on $\beta$ to have two positive roots
$I_{\pm}^{*}$ of $\Phi(I^{*})$ is $\beta_{\Delta}\leq\beta<\beta_{c}$.
Dividing by $\beta_{c}$ and using the definition (8) of $\mathcal{R}_{C}$, we
obtain the condition on $\mathcal{R}_{0}$ as in statement b) of the theorem.
The rest follows from Theorem 2.2. ∎

Figure 2: Visualization of statements in Theorem 2.3 and Theorem 2.4: Endemic
equilibria of system (2) in dependence on $\rho$.

The results of Theorem 2.3 and Theorem 2.4 are summarized in Fig. 2. Local
stability of the endemic equilibria can be determined evaluating the Jacobian
matrix (4) of the system at $\mathcal{E}_{\pm}$,

$J|_{\mathcal{E}_{\pm}}=\left({\begin{array}[]{ccc}-\beta
I^{*}_{\pm}-\gamma&-\beta S_{\pm}^{*}-\gamma&\eta-\gamma\\\
\left(\rho-1\right)\beta
I^{*}_{\pm}&\beta\left(\left(\rho-1\right)S_{\pm}^{*}+1-\left(1+\kappa\right)I^{*}_{\pm}\right)-\alpha&-\beta
I^{*}_{\pm}\\\ 0&\alpha&-\eta\\\ \end{array}}\right).$

The characteristic polynomial of $J|_{\mathcal{E}_{\pm}}$ is given by

$\Pi(\lambda)=\lambda^{3}+\lambda^{2}a_{2}+\lambda a_{1}+a_{0},$ (10)

where

$\displaystyle a_{2}=$ $\displaystyle\beta
I^{*}_{\pm}\left(2+\kappa\right)+\gamma+\alpha+2\beta
S_{\pm}^{*}\left(1-\rho\right)-\beta+\eta$ $\displaystyle a_{1}=$
$\displaystyle
2\beta^{2}I^{*}_{\pm}S_{\pm}^{*}\left(1-\rho\right)+\beta^{2}{I^{*}_{\pm}}^{2}\left(1+\kappa\right)+\gamma\alpha+\gamma\beta
S_{\pm}^{*}\left(1-\rho\right)$ $\displaystyle\;+{}\gamma\beta
I^{*}_{\pm}\left(3-\rho\right)+\gamma\beta\frac{\alpha}{\eta}I^{*}_{\pm}+\eta\alpha+\eta\beta
S^{*}\left(1-\rho\right)+3\eta\beta I^{*}_{\pm}$
$\displaystyle\;+{}\alpha\beta
I^{*}_{\pm}+\eta\gamma-\beta^{2}I^{*}_{\pm}-\eta\beta-\gamma\beta$
$\displaystyle a_{0}=$
$\displaystyle\eta\beta^{2}I^{*}_{\pm}\left(2\left(S_{\pm}^{*}\left(1-\rho\right)+I^{*}_{\pm}\right)-1\right)+\gamma\eta\beta(S_{\pm}^{*}(1-\rho)$
$\displaystyle\;+{}I^{*}_{\pm}\left(3-\rho\right)-1)+(\eta-\gamma)\alpha\beta
I^{*}_{\pm}\left(2-\rho\right)+\gamma\eta\alpha.$

With the Routh–Hurwitz criteria [15] it follows that:

######  Theorem 2.5.

Let $\mathcal{E}^{*}$ be an endemic equilibrium of system (2) and let
$J|_{\mathcal{E}^{*}}$ the coefficient matrix of the linearization of the
system (2) about $\mathcal{E}^{*}$. Then $\mathcal{E}^{*}$ is LAS if the
coefficients of the characteristic polynomial (10) satisfy $a_{2}>0$,
$a_{0}>0$, and $a_{1}a_{2}>a_{0}$.

We enrich the analytical results obtained so far by means of numerical
investigations.

##  3 Numerical bifurcation analysis

From the analysis in Sect. 2 it becomes clear that $\beta$ and $\rho$ are
critical parameters affecting the qualitative behavior of the system. In this
section we extensively investigate the $(\beta,\rho)$-parameter plane and
numerically identify bifurcations of codimension 1 and 2. The remaining
parameters are fixed $\eta=0.02$, $\alpha=0.2$, $\gamma=0.3$, such that
$\rho^{*}\approx 0.1976$. This choice is not motivated by a specific
application or data, but highlights the rich dynamics the model can produce.
All parametric portraits and simulations shown in what follows are produced
using the numerical bifurcation software MATCONT [16]. Background on
bifurcation theory can be found e.g. in the seminal books by Yuri Kuznetsov
[17] or Stephen Wiggins [18].

At first we investigate regions of the $(\beta,\rho)$-plane where $\rho$ is
constant and look for bifurcations in $\beta$ only. Figure 3 shows the
$(\beta,I)$-plane, where curves of (the $I$-component of) endemic equilibria
in dependence of $\beta$ are plotted for six different values of $\rho$,
ordered by ascending $\rho$. Branches in blue (respectively, yellow) represent
locally asymptotically stable (respectively, unstable) equilibria. Red
asterisks mark either a transcritical bifurcation point (BP), a fold
bifurcation point (LP), a supercritical Hopf bifurcation point
($\text{H}^{\text{sup}}$), a subcritical Hopf bifurcation point
($\text{H}^{\text{sub}}$) or a neutral saddle equilibrium (NS, no bifurcation
point). Figure 3 shows that as $\rho$ increases there is a continuous change
from a backward bifurcation (Fig. 3(a–d)) to a forward bifurcation (Fig.
3(e–f)), as the fold bifurcation point LP crosses the $\beta$-axis to the
lower half-plane. This is in line with our analytical results from the
previous section (cf. Theorems 2.2–2.4). Further, Hopf bifurcations occur on
the equilibria branch corresponding to $\mathcal{E}_{+}$. We observe a change
from a subcritical Hopf bifurcation (Fig. 3(a)) to a supercritical Hopf
bifurcation (Fig. 3(b–e)). The neutral saddle equilibrium (Fig. 3(a–c)) gets
closer to the fold bifurcation point, finally collides with it and turns into
a second supercritical Hopf bifurcation point on the branch of
$\mathcal{E}_{+}$ (Fig. 3(d–e)). Increasing $\rho$, the two supercritical Hopf
bifurcation points approach each other (Fig. 3(d–e)), collapse into one point
and eventually vanish (Fig. 3(f)).

Figure 3: Curves of (the $I$-component of) endemic equilibria in dependence of
$\beta$ are plotted for different values of $\rho$. Branches in blue
(respectively, yellow) represent LAS (respectively, unstable) equilibria.

MATCONT provides additional information on the eigenvalues of the
linearization about the endemic equilibria. In case of a backward bifurcation
with a single Hopf bifurcation point (Fig. 3(a–c)), both endemic equilibria
are born unstable at the fold bifurcation. The unstable manifold of
$\mathcal{E}_{+}$ is two-dimensional and the one of $\mathcal{E}_{-}$ is one-
dimensional. The two eigenvalues of the linearization about $\mathcal{E}_{+}$
with positive real part cross the imaginary axis at the Hopf bifurcation point
so that $\mathcal{E}_{+}$ becomes LAS. In the situation where we observe two
(supercritical, Fig. 3(d–e)) Hopf bifurcation points on the branch of
$\mathcal{E_{+}}$, this equilibrium was born LAS at the fold bifurcation. Two
eigenvalues cross the imaginary axis to the right half-plane at the first Hopf
bifurcation point and cross it again back to the left half-plane at the second
Hopf bifurcation point.

Figure 4: (a) Continuation of fold and transcritical bifurcation points with
respect to $\beta$ and $\rho$; (b) Continuation of Hopf bifurcation points
with respect to $\beta$ and $\rho$ starting from the Bogdanov–Takens point BT.

The above observations help us to understand the parametric portrait in
$\beta$ and $\rho$, and the arising codimension 2 bifurcations. Continuation
of fold and transcritical bifurcation points with respect to $\beta$ and
$\rho$ leads to Fig. 4(a). The curve of fold bifurcations (straight line) and
the curve of transcritical bifurcations (curved line) meet at the point
labeled CP. This point marks a cusp bifurcation where the normal form
coefficient of the fold bifurcation vanishes. It divides the curve of
transcritical bifurcations into the branch of forward bifurcations (upper
part) and backward bifurcations (lower part). In addition, a Bogdanov–Takens
(BT) bifurcation occurs on the curve of fold bifurcation points, where an
additional eigenvalue approaches the imaginary axis. At the BT point the curve
of Hopf bifurcation points meets the curve of fold bifurcations tangentially.
Continuation of the Hopf bifurcation curve from BT leads to the red curve in
Fig. 4(b) (neutral saddle equilibria are not shown here). The Bogdanov–Takens
point BT marks the point where the neutral saddle equilibrium turns into a
supercritical Hopf bifurcation point or vice versa. The semi-elliptic shape of
the Hopf curve leads to the occurrence of two Hopf bifurcations, in the
continuation with respect to $\beta$, for values of $\rho$ above the
Bogdanov–Takens point (cf. Fig. 3(d–e)). Another codimension 2 bifurcation
arises on the curve of Hopf bifurcations. A generalized Hopf bifurcation,
where the first Lyapunov coefficient is zero and changes its sign, labeled GH
marks the point where the Hopf bifurcation changes from subcritical (lower
part) to supercritical (upper part). In what follows, we study the branches of
limit cycles emerging from these Hopf bifurcation points.

Figure 5: (a) Continuation of (the I-component of) limit cycles with respect
to $\beta$ starting from subcritical Hopf bifurcation point for $\rho=0.02$;
(b) Zoom of (a); (c) Continuation of (the I-component of) limit cycles with
respect to $\beta$ starting from supercritical Hopf bifurcation point for
$\rho=0.03$. The branch of stable limit cycles reaches a limit point cycle,
folds back and becomes unstable; (d) Zoom of (c).

In Fig. 5(a–b) we see that, for sufficiently small $\rho$, as $\beta$
increases, the amplitudes of the unstable limit cycles (with respect to the
I-component) born at the subcritical Hopf bifurcation increase as well. The
limit cycles get closer and closer to $\mathcal{E}_{-}$. This increases the
period of the closed orbits, and a point moving on such orbits spends more and
more time in the proximity of the equilibrium. The branch of limit cycles
disappears by colliding with the one-dimensional unstable manifold of
$\mathcal{E_{-}}$ leading to a saddle homoclinic bifurcation (more details are
provided in the master thesis [19] of the first author).

Figure 5(c–d) shows the same plot for a slightly larger value of $\rho$ such
that the involved Hopf bifurcation is supercritical. The branch of stable
limit cycles reaches a limit point cycle (LPC) then folds back and becomes
unstable. Thus, there is a small interval of values of $\beta$ in which a
stable limit cycle and an unstable one coexist. The unstable branch of limit
cycles again vanishes through a saddle homoclinic bifurcation. Increasing
$\rho$ further, the interval of values of $\beta$ where a stable limit cycle
exists becomes larger, while the interval where an unstable limit cycle exists
gets smaller and eventually vanishes (i.e. there is no LPC involved anymore).
In this case it is the branch of stable limit cycles that vanishes through a
saddle homoclinic bifurcation (not shown here, cf. [19]).

Figure 6: Continuation of (the I-component of) limit cycles with respect to
$\beta$ starting from two supercritical Hopf bifurcation points for (a)
$\rho=0.2$ and (b) $\rho=0.24$. The two branches of stable limit cycles meet
and form an endemic bubble.

Increasing $\rho$ further, two Hopf Bifurcation points, hence two branches of
stable limit cycles exist, both vanishing through a saddle homoclinic
bifurcation (not shown here, cf. [19]). For even larger values of $\rho$ there
exists a bubble of limit cycles, starting and ending at the Hopf Bifurcation
points (cf. Fig. 6(a)). Increasing $\rho$ the bubble becomes smaller (cf. Fig.
6(b)), until the two Hopf points collide and the limit cycles disappear (cf.
Fig. 3(f)).

Figure 7: Parametric portrait with respect to $\beta$ and $\rho$. The
different bifurcation curves subdivide the $(\beta,\rho)$-plane into regions.
The green and yellow balloons show zoomed regions around the BT and GH point,
respectively.

Continuation of the curve of saddle homoclinic bifurcations (light green) and
fold bifurcations of limit cycles (magenta) yield the complete parametric
portrait in Fig. 7. The different bifurcation curves subdivide the
$(\beta,\rho)$-plane into regions where the system shows different qualitative
behavior. We conclude this numerical study by a discussion of typical solution
trajectories and orbits for each region.

Figure 8: Dynamics in region I: Convergence to $\mathcal{E}_{0}$ shown in (a)
$(t,I)$-plane and (b) $(S,I,R)$-phase space. For these simulations we set
parameter values: $\beta=1,\,\rho=0.15,\,\alpha=0.2,\,\gamma=0.3$, and
$\eta=0.02$; and initial conditions: $R_{0}=0,\,I_{0}=0.5$ and
$I_{0}=0.1,\,S_{0}=1-I_{0}$.

Region I: The disease-free equilibrium $\mathcal{E}_{0}$ is the only
equilibrium and it is LAS (Fig. 8). Crossing the boundary to region IX leads
to a transcritical bifurcation through which $\mathcal{E}_{0}$ loses stability
and the LAS $\mathcal{E}_{+}$ arises. Crossing the boundary to region II
results in a fold bifurcation, with the LAS equilibrium $\mathcal{E_{+}}$ and
the unstable equilibrium $\mathcal{E}_{-}$ (with one-dimensional unstable
manifold). Passing from region I to region III also leads to a fold
bifurcation, through which both $\mathcal{E}_{+}$ (two-dimensional unstable
manifold) and $\mathcal{E}_{-}$ (one-dimensional unstable manifold) arise
unstable.

Figure 9: Dynamics in region II: Bistability of $\mathcal{E}_{0}$ and
$\mathcal{E}_{+}$. For these simulations we set parameter values:
$\beta=1.1,\,\rho=0.181,\,\alpha=0.2,\,\gamma=0.3$, and $\eta=0.02$, and
initial conditions: (a) $R_{0}=0.1,\,S_{0}=0.85,\,I_{0}=0.015$; (b)
$R_{0}=0,\,S_{0}=0.9,\,I_{0}=0.1$.

Region II: The disease-free equilibrium $\mathcal{E}_{0}$ is LAS. Both endemic
equilibria $\mathcal{E}_{\pm}$ exist, $\mathcal{E}_{+}$ being LAS and
$\mathcal{E}_{-}$ being unstable with one-dimensional unstable manifold.
Figure 9 shows the bistability of $\mathcal{E}_{0}$ and $\mathcal{E}_{+}$.
Crossing the boundary to region IX leads to a transcritical backward
bifurcation through which $\mathcal{E}_{0}$ loses stability and
$\mathcal{E}_{-}$ moves out of the first quadrant. Thus, we observe a
discontinous jump of orbits from the DFE $\mathcal{E}_{0}$ to the already
quite large endemic equilibrium $\mathcal{E}_{+}$, which is known as
hysteresis. Passing from region II to region IV leads to a Hopf bifurcation
through which $\mathcal{E}_{+}$ becomes unstable with two-dimensional unstable
manifold, and a stable periodic orbit emerges. Crossing the boundary between
region II and region VI (intersection of Hopf and transcritical bifurcation
curve) the above bifurcations happen simultaneously.

Figure 10: Dynamics in region III: Convergence to the DFE $\mathcal{E}_{0}$.
Both endemic equilibria $\mathcal{E}_{\pm}$ are unstable. For these
simulations we set parameter values:
$\beta=2.5,\,\rho=0.05,\,\alpha=0.2,\,\gamma=0.3$, and $\eta=0.02$; and
initial conditions: (a) $R_{0}=0.65,\,S_{0}=0.207,\,I_{0}=0.065$; (b)
$R_{0}=0.05,\,S_{0}=0.9,\,I_{0}=0.006$.

Region III: The system has three equilibria (cf. Fig. 10). The DFE
$\mathcal{E}_{0}$ is LAS and both endemic equilibria $\mathcal{E}_{\pm}$ are
unstable, the unstable manifold of $\mathcal{E}_{+}$ being two-dimensional,
that of $\mathcal{E}_{-}$ one-dimensional. We observe the phenomenon of
excitability: Solutions starting close to $\mathcal{E}_{-}$ follow its
unstable manifold, which leads to a large epidemic outbreak before converging
to $\mathcal{E}_{0}$ (cf. Fig. 10(b)). Crossing the boundary to region IV
results into a saddle homoclinic bifurcation, which produces a stable periodic
orbit. Passing the boundary to region V, a fold bifurcation of limit cycles
takes place and a stable and an unstable periodic orbit are born. Moving from
region III to region VII leads to a subcritical Hopf bifurcation,
$\mathcal{E}_{+}$ becomes LAS and an unstable limit cycle arises.

Figure 11: Dynamics in region IV: Bistability of $\mathcal{E}_{0}$ and a limit
cycle. For these simulations we set parameter values:
$\beta=1.14,\,\rho=0.174,\,\alpha=0.2,\,\gamma=0.3$ and $\eta=0.02$; and
initial conditions: (a) $R_{0}=0.15,\,S_{0}=0.82,\,I_{0}=0.015$; (b)
$R_{0}=0,\,S_{0}=0.9,\,I_{0}=0.1$.

Region IV: The disease-free equilibrium $\mathcal{E}_{0}$ is LAS and both
endemic equilibria $\mathcal{E}_{\pm}$ are unstable, the unstable manifold of
$\mathcal{E}_{+}$ being two-dimensional, that of $\mathcal{E}_{-}$ one-
dimensional. Additionally, a stable limit cycle exists. This leads to a
bistability of the disease-free equilibrium and a periodic solution (cf. Fig.
11). Crossing the boundary to region V results into a saddle homoclinic
bifurcation, which produces a second (unstable) periodic orbit. Passing the
boundary to region VI, a transcritical bifurcation occurs, through which
$\mathcal{E}_{0}$ loses stability and $\mathcal{E}_{-}$ moves out of the first
quadrant (hysteresis, cf. Region II). Moving from region IV to region VIII
leads to a Hopf bifurcation, $\mathcal{E}_{+}$ becomes LAS and the stable
periodic orbit vanishes. Crossing the boundary to region IX (intersection of
Hopf and transcritical bifurcation curve) the last two bifurcations happen
simultaneously.

Figure 12: Dynamics in region V: Bistability of $\mathcal{E}_{0}$ and a limit
cycle. For these simulations we set parameter values:
$\beta=2.952,\,\rho=0.046,\,\alpha=0.2,\,\gamma=0.3,\,\eta=0.02$ and initial
conditions: (a) $R_{0}=0.03,\,S_{0}=0.94,\,I_{0}=0.003$; (b)
$R_{0}=0.3,\,S_{0}=0.3,\,I_{0}=0.2$.

Region V: The dynamical behavior in this region is similar to that in Region
IV (cf. Fig. 12). Additionally to the stable limit cycle, an unstable limit
cycle exists, separating the basin of attraction of $\mathcal{E}_{0}$ and of
the stable limit cycle. Like in Region III, we observe the phenomenon of
excitability for orbits starting close to $\mathcal{E}_{-}$ (cf. Fig. 12(a)).
Moving to region VII leads to a supercritical Hopf bifurcation. The stable
limit cycle disappears and $\mathcal{E}_{+}$ becomes LAS. Crossing the
boundary to region VIII (intersection of Hopf and saddle homoclinic
bifurcation curve) the supercritical Hopf bifurcation happens simultaneously
with a saddle homoclinic bifurcation through which also the unstable limit
cycle vanishes.

Figure 13: Dynamics in region VI: Convergence to a stable limit cycle shown in
(a) $(t,I)$-plane and (b) $(S,I,R)$-phase space. For these simulations we set
parameter values: $\beta=1.4,\,\rho=0.2,\,\alpha=0.2,\,\gamma=0.3,\,\eta=0.02$
and initial conditions: $R_{0}=0,\,S_{0}=0.9,\,I_{0}=0.1$.

Region VI: The disease-free equilibrium $\mathcal{E}_{0}$ is unstable, and so
is also the only endemic equilibrium $\mathcal{E}_{+}$. There exists a stable
limit cycle to which solutions converge (cf. Fig. 13). Moving from region VI
to region IX leads to a supercritical Hopf bifurcation, the stable limit cycle
vanishes and $\mathcal{E}_{+}$ becomes LAS.

Figure 14: Dynamics in region VII: Bistability of $\mathcal{E}_{0}$ and
$\mathcal{E}_{+}$. For these simulations we set parameter values:
$\beta=3.27,\,\rho=0.01,\,\alpha=0.2,\,\gamma=0.3,\,\eta=0.02$ and initial
conditions: (a) $R_{0}=0.5,\,S_{0}=0.3,\,I_{0}=0.2$; (b)
$R_{0}=0.06,\,S_{0}=0.88,\,I_{0}=0.006$.

Region VII: The equilibria $\mathcal{E}_{0}$ and $\mathcal{E}_{+}$ are both
LAS (cf. Fig. 14), whereas the third equilibrium $\mathcal{E}_{-}$ is unstable
with one-dimensional unstable manifold. Furthermore, there exists an unstable
periodic orbit that separates the basin of attraction of $\mathcal{E}_{0}$ and
$\mathcal{E}_{+}$. Like in Region III, we observe the phenomenon of
excitability for orbits starting close to $\mathcal{E}_{-}$ (cf. Fig. 14(b)).
Crossing the boundary to region VIII leads to a saddle homoclinic bifurcation
and the unstable periodic orbit vanishes.

Region VIII: In this region, we have bistability of $\mathcal{E}_{0}$ and
$\mathcal{E}_{+}$ (cf. Fig. 15(a)). The second endemic equilibrium
$\mathcal{E}_{-}$ is unstable with one-dimensional unstable manifold
connecting to $\mathcal{E}_{+}$. Crossing the boundary to region IX, a
transcritical bifurcation takes place, so that $\mathcal{E}_{0}$ loses
stability and $\mathcal{E}_{-}$ leaves the positive quadrant (hysteresis, cf.
Region II).

Figure 15: (a) Dynamics in region VIII: Bistability of $\mathcal{E}_{0}$ and
$\mathcal{E}_{+}$. For these simulations we set parameter values:
$\beta=3.4,\rho=0.01,\,\alpha=0.2,\,\gamma=0.3,\,\eta=0.02$ and initial
values: $R_{0}=0,\,S_{0}=0.9,\,I_{0}=0.001$ and $I_{0}=0.01$. (b) Dynamics in
region IX: Convergence to the unique endemic equilibrium $\mathcal{E}_{+}$.
For these simulations we set parameter values:
$\beta=2.5,\,\rho=0.25,\,\alpha=0.2,\,\gamma=0.3,\,\eta=0.02$ and initial
conditions: $R_{0}=0,\,I_{0}=0.1$ and $I_{0}=0.3,\,S_{0}=1-I_{0}$.

Region IX: The DFE $\mathcal{E}_{0}$ is unstable, and the unique endemic
equilibrium $\mathcal{E}_{+}$ is LAS (cf. Fig. 15(b)).

So far we have investigated the $(\beta,\rho)$-parameter plane fixing the
values of $\alpha,\,\gamma,\,\eta$. One might ask how the qualitative behavior
of system (2) is affected by the particular choice of these three parameters.
For instance, we constructed the parametric portrait from Fig. 7 for different
values of $\eta$ (same was done for $\alpha$ and $\gamma$, though not shown
here). The result is shown in Fig. 16. Reducing $\eta$ (Fig. 16(a)) region VI
gets significantly larger. Slowing down the transition from $R$ to $S$ makes
the model closer to a SIRS system with delay (cf. [20, 21]) enhancing the
occurrence of oscillations. Moreover, the generalized Hopf point in Fig. 16(b)
moves to the fourth quadrant. In contrast, increasing $\eta$ (Fig. 16(c))
region VI gets significantly smaller, whereas the points GH, CP, BT move up
along the bifurcation curves.

Figure 16: Parametric portrait in $\beta$ and $\rho$, for $\alpha=0.2$,
$\gamma=0.3$, and (a) $\eta=0.016$; (b) $\eta=0.02$; (c)=0.024.

##  4 Discussion

In this work we have presented a mathematical model for transmission dynamics
with two nonlinear stages of contagion, applying analytical and numerical
methods to analyze its qualitative behavior. The two-stage contagion is
combined with a renewal of the susceptible compartment, modeled by (i) the
waning of immunity in inactive/resistant individuals and (ii) the fading of
infection in those individuals who are weakened after the first contact with
the promoting/infected community. This leads to rich dynamics, including
bistability of equilibria or bistability of an equilibrium and a periodic
solution, discontinuous regime shifts through hysteresis effects, and
excitability. Thus, the multi-stage nature of social contagion processes might
explain some of the complex phenomena observed in social dynamics (see e.g.
the irregular [22] or periodic [23] outcomes in political elections, or the
emergence of new trends in the usage of social media [24]). In a previous
study by Guy Katriel [8] similar properties were determined for a two-stage
contagion model with demographic turnover. This lets us conjecture that the
rich dynamics observed in our work and in [8] is due to the coupling of a two-
stage contagion process with any (demographic or ”immunological”) source of
renewal of the susceptible population. Despite of the analogies with Katriel’s
work, the combination of waning of immunity and fading of infection in our
model leads to additional analytical complexity. The analytical advantage of
Katriel’s model was possibly due to the choice of the birth rate in the
$S$-compartment matching with the death rates of all compartments. Katriel’s
model shows also both backward and forward bifurcation, and in the latter case
it behaves like a one-stage contagion model for any choice of $\beta$. In
contrast, our model shows Hopf bifurcations, hence periodic solutions, also in
case of a forward bifurcation (cf. Fig. 7(e)). Single phenomena which can be
observed in our model have been previously found also for variations of the
classical one-stage contagion models. For example SIRS models with delayed
loss of immunity [20] naturally show stable periodic solutions, whereas
bistability of equilibria was found e. g. in models with exogenous reinfection
[14] or imperfect vaccination [25].

In this work we have focused on the investigation of the qualitative
properties of a simple two-stage contagion model. We have numerically
investigated the parameter space, focusing in particular on the effects of the
transmission rate ($\beta$) and the probability of a perfect contact $(\rho)$.
Of course the study could be repeated deriving the parametric portrait of the
system (2) with respect to the other model parameters as well. Moreover, we
see three possible generalizations of our model: (i) the assumption that
$\beta=\beta_{2}$ could be relaxed, and e.g. $W$-individuals might be assumed
to have a higher susceptibility than $S$-individuals; (ii) individuals in the
$W$-compartment could also be contagious (cf. also [26]); (iii) the model
could also include births and deaths, in addition to waning/fading processes.
All these variations would make the analytical investigations more
challenging; however, a similar numerical investigation as presented in this
work could be performed. Thinking of applications and comparison with data,
the major limitation of our work is due to the deterministic approach.
Dividing a population into a few homogeneous compartments, without taking into
account interpersonal variability, is indeed a major simplification of
reality. Refining our approach, agent-based modeling [27] and complex networks
[28] could be used. In certain cases, previous works based on these methods
also included two-stage contagion [29, 30].

## Acknowledgement

The authors are supported by the LOEWE focus CMMS.

## Appendix A

In the proof of Theorem 2.2 we referred to the following result by Castillo-
Chavez and Song [14] proved using center-manifold theory.

######  Theorem A.1.

Let $f\in C^{2}(\mathbb{R}^{n}\times\mathbb{R},\mathbb{R}^{n})$. Consider the
following general system of ODEs with a parameter $\beta$

$\displaystyle\frac{dx}{dt}=f(x,\beta).$ (11)

Without loss of generality assume that $x_{0}=0$ is an equilibrium point of
the system, that is, $f(0,\beta)=0$ for all $\beta\in\mathbb{R}$. Assume the
following:

  * •

The linearization of the system (11)
$\mathcal{A}:=D_{x}f(0,0)=\left(\frac{\partial f_{i}}{\partial
x_{j}}(0,0)\right)$ has zero as a simple eigenvalue and all other eigenvalues
of $\mathcal{A}$ have negative real parts.

  * •

The matrix $\mathcal{A}$ has a non-negative right eigenvector $w$ and a left
eigenvector $v$ each corresponding to the zero eigenvalue.

Let $f_{k}$ be the $k$-th component of $f$ and

$\displaystyle a$
$\displaystyle=\sum_{k,i,j=1}^{n}v_{k}w_{i}w_{j}\frac{\partial^{2}f_{k}}{\partial
x_{i}\partial x_{j}}(0,0)$ $\displaystyle b$
$\displaystyle=\sum_{k,i=1}^{n}v_{k}w_{i}\frac{\partial^{2}f_{k}}{\partial
x_{i}\partial\beta}(0,0).$

Then, the local dynamics of the system around $0$ is completely determined by
the signs of $a$ and $b$:

  1. 1.

$a>0,\,b>0$. When $\beta<0$ with $|\beta|\ll 1$, $0$ is locally asymptotically
stable and there exists a positive unstable equilibrium; when $0<\beta\ll 1$,
$0$ is unstable and there exists a negative and locally asymptotically stable
equilibrium.

  2. 2.

$a<0,\,b<0$. When $\beta<0$ with $|\beta|\ll 1$, $0$ is unstable; when
$0<\beta\ll 1$, $0$ is locally asymptotically stable and there exists a
positive unstable equilibrium.

  3. 3.

$a>0,\,b<0$. When $\beta<0$ with $|\beta|\ll 1$, $0$ is unstable and there
exists a locally asymptotically stable negative equilibrium; when $0<\beta\ll
1$, $0$ is stable and a positive unstable equilibrium appears.

  4. 4.

$a<0,\,b>0$. When $\beta$ changes from negative to positive, $0$ changes its
stability from stable to unstable. Correspondingly, a negative unstable
equilibrium becomes positive and locally asymptotically stable.

In particular, if $a>0$ and $b>0$, a backward bifurcation occurs at $\beta=0$,
and if $a<0$ and $b>0$, a forward bifurcation occurs at $\beta=0$.

######  Remark.

Remark 1 in [14] suggests that if the equilibrium of interest in Theorem A.1
is a non-negative equilibrium $x_{0}$, then the requirement that w is non-
negative is not necessary. When some components in w are negative, one can
still apply Theorem A.1 provided that the $j$-th component of $w$ is positive
whenever the $j$-th component of $x_{0}$ is zero. If the $j$-th component of
$x_{0}$ is positive, then the $j$-th component of $w$ need not be positive.

## References

  * [1] American Psychological Association. Dictionary of Psychology, https://dictionary.apa.org/social-contagion (last accessed 10 Sept, 2020).
  * [2] J. Sooknanan, and D. M. Comissiong, When behaviour turns contagious: the use of deterministic epidemiological models in modeling social contagion phenomena, Int. J. Dyn. Control 5(4), 1046 (2017).
  * [3] W. Goffman, and V. A. Newill, Generalization of epidemic theory: An application to the transmission of ideas, Nature 204(4955), 225 (1964).
  * [4] L. M. Bettencourt, A. Cintrón-Arias, D. I. Kaiser and C. Castillo-Chávez, The power of a good idea: Quantitative modeling of the spread of ideas from epidemiological models, Physica A 364, 513 (2006).
  * [5] F. Jin, E. Dougherty, P. Saraf, Y. Cao and N. Ramakrishnan, Epidemiological modeling of news and rumors on twitter, In SNAKKD’13: Proceedings of The 7th Workshop on Social Network Mining and Analysis (2013).
  * [6] J. Woo and H. Chen, Epidemic model for information diffusion in web forums: experiments in marketing exchange and political dialog, Springer (2016).
  * [7] S. Wang, L. Rong and J. Wu, Bistability and multistability in opinion dynamics models, Appl. Math. Comput. 289, 388 (2016).
  * [8] G. Katriel, The dynamics of two-stage contagion, Chaos Solitons Fract. 2, 100010 (2019).
  * [9] C. T. Bauch and A. P. Galvani, Social factors in epidemiology, Science 342(6154), 47 (2013).
  * [10] A. L. Hill, D. G. Rand, M. A. Nowak and N. A. Christakis, Infectious disease modeling of social contagion in networks, PLoS Comp. Biol. 6(11), e1000968 (2010).
  * [11] M. Martcheva, An introduction to mathematical epidemiology, Springer (2015).
  * [12] M.V. Barbarossa and G. Röst, Immuno-epidemiology of a population structured by immune status: a mathematical study of waning immunity and immune system boosting, J. Math. Biol. 71(6), 1737 (2015).
  * [13] E. M. Rogers, Diffusion of Innovations, Free Press (2010).
  * [14] C. Castillo-Chavez and B. Song, Dynamical models of tuberculosis and their applications, Math. Bioscie. Eng. 1(2), 361 (2004).
  * [15] J. D. Murray, Mathematical biology: I. An introduction, Springer (2007).
  * [16] A. Dhooge, W. Govaerts and Y. A. Kuznetsov, MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs, ACM Trans. Math. Software 29(2), 141 (2003).
  * [17] Y. A. Kuznetsov, Elements of applied bifurcation theory, Springer (2013).
  * [18] S. Wiggins, Introduction to applied nonlinear dynamical systems and chaos, Springer (2003).
  * [19] J. Heidecke, Mathematical modelling of transmission processes with two-stage contagion, Master Thesis, Heidelberg University, Heidelberg (2020).
  * [20] M. L. Taylor and T. W. Carr, An SIR epidemic model with partial temporary immunity modeled with delay, J. Math. Biol. 59(6), 841 (2009).
  * [21] M. V. Barbarossa, M. Polner and G. Röst, Stability switches induced by immune system boosting in an SIRS model with discrete and distributed delays, SIAM J. Appl. Math. 77(3), 905 (2017).
  * [22] Politico. Italy – National parliament voting intention https://www.politico.eu/europe-poll-of-polls/italy/ (last accessed Oct 14 2020)
  * [23] Wikipedia. United States presidential election https://en.wikipedia.org/wiki/United_States_presidential_election (last accessed Oct 14 2020)
  * [24] E. Ortiz-Ospina. The rise of social media. Our world in data, Sept 28 2019 https://ourworldindata.org/rise-of-social-media (last accessed Oct 14 2020)
  * [25] F. Brauer, Backward bifurcations in simple vaccination models, J. Math. Anal. Appl. 298(2), 418 (2004).
  * [26] N. M. Crisosto, C. M. Kribs-Zaleta, C. Castillo-Chávez and S. Wirkus, Community resilience in collaborative learning, Discr. Cont. Dyn. Sys. B 14(1), 17 (2010).
  * [27] G. Deffuant, D. Neau, F. Amblard and G. Weisbuch, Mixing beliefs among interacting agents, Adv. Complex Sys. 3(01n04), 87 (2000).
  * [28] G. Albi, L. Pareschi and M. Zanella, Opinion dynamics over complex networks: kinetic modeling and numerical methods, Kin. Rel. Models 10(1) (2017).
  * [29] T. Hasegawa and K. Nemoto, Discontinuous transition of a multistage independent cascade model on networks, J. Stat. Mech. 2014(11), P11024 (2014).
  * [30] S. Melnik, J. A. Ward, J. P. Gleeson and M. A. Porter, Multi-stage complex contagions, Chaos 23(1), 013124 (2013).