# Global stability of SIR model with heterogeneous transmission rate modeled
by the Preisach operator

Ruofei Guan, Jana Kopfová, Dmitrii Rachinskii

###### Abstract

In recent years, classical epidemic models, which assume stationary behavior
of individuals, have been extended to include an adaptive heterogeneous
response of the population to the current state of the epidemic. However, it
is widely accepted that human behavior can exhibit history-dependence as a
consequence of learned experiences. This history-dependence is similar to
hysteresis effects that have been well-studied in control theory. To
illustrate the importance of history-dependence for epidemic theory, we study
dynamics of a variant of the SIRS model where individuals exhibit lazy-switch
responses to prevalence dynamics. The resulting model, which includes the
Preisach hysteresis operator, possesses a continuum of endemic equilibrium
states characterized by different proportions of susceptible, infected and
recovered populations. We discuss stability properties of the endemic
equilibrium set and relate them to the degree of heterogeneity of the adaptive
response. Our results support the argument that public health responses during
the emergence of a new disease can have long-term consequences for subsequent
management efforts. The main mathematical contribution of this work is a
method of global stability analysis, which uses a family of Lyapunov functions
corresponding to different branches of the hysteresis operator.

Keywords: SIR model, switching transmission rate, endemic equilibrium,
connected set of equilibrium states, Preisach hysteresis operator, Lyapunov
function

MSC: 34D23, 92D30, 92D25

##  1 Introduction

The standard SIR model assumes constant transmission rate (Kermack and
McKendrick 1927). However, during epidemics, the transmission rate can change
due to interventions of the health authorities and an adaptive response of the
population to dynamics of the infectious disease. The role of the modulation
of the transmission rate due to a human response to dynamics of epidemics has
been widely studied (Capasso and Serio 1978; d’Onofrio and Manfredi 2009;
d’Onofrio and Manfredi 2020). Preventive measures that decrease the
transmission rate can include media campaigns raising awareness in the
population about the current severity of the epidemic, massive vaccination
campaigns, access to effective and affordable tests and medicines provided by
the authorities, face covering in public places, quarantine and social
distancing measures, transition to online teaching at schools and colleges,
gathering and travel limitations, business restrictions and stay-at-home
orders (Gostin and Wiley 2020; Marquioni and de Aguiar 2020; Center of Disease
Control and Prevention 2021). Public participation in these measures can be
partly mandatory and partly voluntary.

One major objective of preventive measures is flattening the curve, i.e.
slowing down the spread of infection in order to keep the number of active
disease cases at a manageable level dictated by the capacity of the health
care system (Matrajt and Leung 2020). However, restrictive preventive measures
can be typically introduced only for a limited period of time due to the
economic and social cost they incur (Fairlie 2020). Because of these
constraints, many prevention policies (such as the intervention protocols
which have been recently implemented by the health authorities in order to
contain the Covid-19 epidemic (Wilder-Smith et al. 2020) can be thought of, at
least simplistically, as threshold based. As such, an intervention starts when
a certain variable such as the number of daily new infections, the percentage
of the occupied hospital beds or the effective reproduction number $R$ reaches
a critical threshold value set by the health or government authority
(Department of Health and Human Services, Nebraska 2020; DeBenedetto and Ruiz
2021); the intervention is revoked when this variable drops below the level
deemed safer. If the thresholds at which the intervention begins and ends are
different, then hysteresis effect is present.

In their earliest form, hysteresis effects, which are commonly called a
‘backward bifurcation’ (a subcritical steady state bifurcation) in the
epidemic theory literature, were clearly identified by Hadeler and Chavez
(1995) in their study of a core group. Hysteresis in adopting a ‘healthy
behavior’ can also result from social contagion processes such as social
interaction, reinforcement and imitation (Su et al. 2017; Wang et al. 2017).
For instance, in the presence of imperfect vaccine, hysteresis loops of
vaccination rate arise with respect to dynamic changes in the perceived cost
of vaccination (Chen and Fu 2019) as individuals revisit their vaccination
decision in response to dynamics of the epidemic and through a social learning
processes under peer influence. These hysteresis loops form a complex
structure with minor loops nested within the maximal loop.

Hysteresis is a type of dependence of the state of a system on its past states
(the term was first coined to describe magnetic hysteresis as lagging of the
magnetization of a ferromagnetic material behind variations of the magnetizing
field (Ewing 1881)). Characteristically, this lagging is determined by the
past values of the phase variables and the ordering of the past events rather
than the timing of those events. In particular, hysteresis is not associated
with any explicit time scale and, as such, differs from the delayed response
of systems with deviating argument (delayed and functional differential
equations). In fact, hysteresis is defined as a rate-independent operator in
(Visintin 1994) and is characterized by a durable effect of temporary stimuli
in biological systems and population processes (Pimenov et al. 2012). Thus,
the rate-independence of hysteresis and the prescribed time scale of systems
with a deviating argument (and closely associated convolution operators) can
be viewed as two idealizations which present a modeler with a choice of
complementary mathematical tools for modeling a particular process.

The effectiveness of the prevention measures depends on the response of the
population to the intervention policies and, in particular, the degree of
homogeneity of this response. The ability and willingness of an individual to
receive immunization or to follow community isolation policies depends on the
level of trust to the authorities in the community, personal beliefs,
perceived risk of contracting the disease, risk of possible complications and
other factors, which vary with age, health condition, living circumstances,
profession and social group (Guidry et al. 2021; Lazarus et al. 2021). All
these variations contribute to the heterogeneity of the response of the
society to the advent of an epidemic including behavioral changes. Effects of
switching self-protecting behavior on the transmission of a disease were
analyzed, for example, by Geoffard and Philipson (1996) and Chen (2004) who
modeled a population of agents with heterogeneous preferences for taking
protective actions using rational choice theory. Oscillatory vaccination
dynamics in heterogeneous populations were identified by Reluga, Bauch and
Galvani (2006). An important finding by Chen and Fu (2019) who studied a
relaxation oscillator limit of these dynamics is that hysteresis in adopting
healthy behavior becomes more pronounced with increasing heterogeneity of the
population.

In this paper, we attempt to model the effect of hysteresis and heterogeneity
in the response of the population to preventive measures on the epidemic
trajectory. Our focus is on the global stability of the set of equilibrium
states. We use a variant of the SIR model where the transmission rate depends
on dynamics of the infected population. As a starting point, we adopt the
approach used by Chladná et al. (2020) to modeling the homogeneous switched
response of the population to the varying number of infected individuals by a
two-threshold two-state relay operator. In this model, it is assumed that the
health authorities implement a two-threshold intervention policy whereby the
intervention starts when the number of infected individuals exceeds a critical
threshold value and stops whenever this number drops below a different (lower)
threshold. The objective of the two-threshold policy is to navigate the system
to the endemic equilibrium simultaneously keeping the number of infected
individuals in check and not committing to continuous intervention. If the
response is ideally homogeneous, then prevention measures are assumed to
translate immediately to reduced values of the transmission coefficient and
the effective reproduction number during the intervention. This response is
characterized by the simplest rectangular hysteresis loop in the relationship
between the density of the infected population and the transmission
coefficient. The rectangular hysteresis loop can be considered as a model of
the bi-stable response associated with the backward bifurcation.

Next, to reflect the heterogeneity of the response, the population is divided
into multiple subpopulations, each characterized by a different pair of
switching thresholds. In order to keep the model relatively simple, we apply
averaging under additional simplifying assumptions as in Kopfová et al.
(2021); Rouf and Rachinskii (2021). This leads to a differential model with
just two variables, $S$ and $I$, but with a complex operator relationship
between the transmission rate and the density of the infected population $I$.
As such, this operator relationship, known as the Preisach hysteresis operator
(see, for example, Krasnosel’skii and Pokrovskii 1989; Mayergoyz 2003)
accounts for hysteresis in the heterogeneous response. On the other hand,
applications of differential models with a population of two-state relay
operators to population dynamics date back to the pioneering work by Jäger and
Hoppensteadt (1980).

We show that due to hysteresis, the heterogeneous model has a connected
continuum of endemic equilibrium states characterized by different proportions
of the susceptible, infected and recovered populations. We consider the
convergence of the epidemic trajectory to this continuum. In mechanical
systems, hysteresis incurs energy losses thus promoting convergence to an
equilibrium state. This convergence can be shown using the energy potential as
a Lyapunov function (for the construction of hysteretic thermodynamic
potentials and examples, see Brokate and Sprekels (1996); Krejčí et al.
(2013); Krejčí et al. (2019); related methods of the construction of the
Lyapunov function for the problem of absolute stability of Lurie systems with
hysteresis in the feedback loop can be found in Gelig and Yakubovich (1978)
and the review Leonov et al.(2017). However, for SIR systems, hysteresis is a
source of instability, which can lead to periodic oscillations corresponding
to recurrent waves of infection (Chladná et al. 2020; Kopfová et al. 2021).
The objective of this work is to show the convergence to the equilibrium set
if hysteresis has a limited magnitude (as quantified by an appropriate
measure) and provide estimates for the magnitude of the hysteresis effect,
which guarantee this stable scenario globally. A numerical case study of the
dynamical scenarios exhibited by the heterogeneous SIR model considered below
has been provided by Rachinskii and Rouf (2021).

To analyze stability, we interpret the SIR model with the Preisach operator as
a switched system (di Bernardo et al. 2008) associated with an infinite family
of nonlinear planar vector fields. Between the switching moments, a trajectory
of the system is an integral curve of a particular globally stable vector
field with the associated global Lyapunov function. The Preisach operator
imposes non-trivial rules for switching from one vector field to another
corresponding to switching from one to another branch of the hysteresis
nonlinearity. We prove stability by carefully controlling (a) the increments
of the Lyapunov functions along the trajectory between the switching points;
and, (b) the difference of the Lyapunov functions corresponding to different
vector fields (associated with different hysteresis branches) at the switching
points. This methodology has been previously applied to a simpler special
situation where the switching surface was the same for all the hysteresis
branches (Kopfová et al. 2021[27]). The main mathematical contribution of this
paper is the extension of this method of global stability analysis to the more
typical generic case when the switching surface changes dynamically depending
on the history of the state.

The paper is organized as follows. In the next section, we discuss an SIR
model with a hysteretic heterogeneous response. Stability results are
presented in Section 3 and discussed in Section 4. The Appendix contains the
proofs.

##  2 Model

2.1. The systems considered below adapt the standard dimensionless SIR model

$\displaystyle\dot{I}$ $\displaystyle=\beta SI-(\gamma+\mu)I,$ (1)
$\displaystyle\dot{S}$ $\displaystyle=-\beta SI-\mu S+\mu$

with simple immigration demography, where $I$ and $S$ are the densities of the
infected and susceptible populations, respectively; $\beta$ is the
transmission coefficient; $\gamma$ is the recovery rate; and $\mu$, is the
departure rate due to the disease unrelated death and emigration. We assume a
constant total population scaled to unity, hence the density of the recovered
population $R=1-I-S$ can be removed from the system. The domain $I,S\geq 0$,
$I+S\leq 1$ is positively invariant for system (1).

Re-scaling the time, one obtains the normalized system

$\displaystyle\dot{I}$ $\displaystyle=R_{0}SI-I,$ (2) $\displaystyle\dot{S}$
$\displaystyle=-R_{0}SI-\rho S+\rho,$

with the dimensionless parameters

$R_{0}=\frac{\beta}{\gamma+\mu},\qquad\rho=\frac{\mu}{\gamma+\mu}<1.$

Here $R_{0}$ is the basic reproduction number and $R=R_{0}S$ is the effective
reproduction number for any given density $S$ of the susceptible population.
If $R_{0}<1$, then the infection free equilibrium $(I^{*},S^{*})=(0,1)$ is
globally stable in the closed positive quadrant. On the other hand, if
$R_{0}>1$, then the infection free equilibrium is a saddle, and the positive
endemic equilibrium $(I_{*},S_{*})$ defined by

$I_{*}=\left(1-\frac{1}{R_{0}}\right)\rho,\qquad S_{*}=\frac{1}{R_{0}}$ (3)

is globally stable in the open positive quadrant. If

$\rho(R_{0})^{2}<4(R_{0}-1),$ (4)

then the endemic equilibrium is of focus type.

2.2. A switched system with two flows of the form (1), which have different
values of the basic reproduction number, $R_{0}=R_{0}^{nat}$ and
$R_{0}=R_{0}^{int}$, was considered in Chladná et al., 2020; Rouf and
Rachinskii, 2021. Following this work, we assume that the parameter $\rho$ is
the same for both flows and a switch from one flow to the other occurs when
the density of the infected population reaches certain thresholds,
$I=\alpha_{1}$ and $I=\alpha_{2}$, satisfying $0<\alpha_{1}<\alpha_{2}<1.$
More precisely, it is postulated that the basic reproduction number
instantaneously switches from the value $R_{0}^{nat}$ to the value
$R_{0}^{int}$ when the variable $I=I(t)$ reaches the upper threshold value
$\alpha_{2}$. On the other hand, $R_{0}$ switches back from the value
$R_{0}^{int}$ to the value $R_{0}^{nat}$ as $I(t)$ reaches the lower threshold
value $\alpha_{1}$, see Figure 1.

$R_{0}$$I$$R_{0}^{nat}$$R_{0}^{int}$$\alpha_{1}$$\alpha_{2}$ Figure 1:
Switching rule (10) for the basic reproduction number. The vertical segments
correspond to instantaneous transitions of the basic reproduction number from
the value $R_{0}^{nat}$ to $R_{0}^{int}$ and back.

This two-threshold switching rule can be formalized using the standard non-
ideal relay operator, which is also known as a rectangular hysteresis loop or
a lazy switch (see e.g. Visintin 1994). The relay is characterized by the
threshold parameters $\alpha_{1},\alpha_{2}$, and we use the notation
$\alpha=(\alpha_{1},\alpha_{2})$. The function
$I(t):\mathbb{R}_{+}\to\mathbb{R}$ is called the input of the relay. The state
of the relay, denoted $r_{\alpha}(t)$, equals either $0$ or $1$ at any moment
$t\in\mathbb{R}_{+}$. Specifically, given any continuous input
$I(t):\mathbb{R}_{+}\to\mathbb{R}$ and an initial value of the state,
$r_{\alpha}(0)=r_{\alpha}^{0}$, which satisfies the constraints

$r_{\alpha}^{0}\in\\{0,1\\}\quad\text{if}\quad\alpha_{1}<I(0)<\alpha_{2};$ (5)
$r_{\alpha}^{0}=0\quad\text{if}\quad I(0)\leq\alpha_{1};\qquad
r_{\alpha}^{0}=1\quad\text{if}\quad I(0)\geq\alpha_{2},$ (6)

the state of the relay at the future moments $t>0$ is defined by the relations

$r_{\alpha}(t)=\left\\{\begin{array}[]{cl}0&\text{if there is $t_{1}\in[0,t]$
such that $I(t_{1})\leq\alpha_{1}$}\\\ &\text{and $I(\tau)<\alpha_{2}$ for all
$\tau\in(t_{1},t]$;}\\\ 1&\text{if there is $t_{1}\in[0,t]$ such that
$I(t_{1})\geq\alpha_{2}$}\\\ &\text{and $I(\tau)>\alpha_{1}$ for all
$\tau\in(t_{1},t]$};\\\ r_{\alpha}(0)&\text{if $\alpha_{1}<I(\tau)<\alpha_{2}$
for all $\tau\in[0,t]$.}\end{array}\right.$ (7)

Function (7) will be denoted by

$r_{\alpha}(t)=({\mathcal{R}}_{\alpha}[r_{\alpha}^{0}]I)(t),\qquad t\geq 0.$
(8)

This standard short notation adopted from Brokate and Sprekels, 1996 stresses
that the relay’s state is dependent on the history of $I(t)$, i.e. (7) defines
an operator which takes the function $I(t):\mathbb{R}_{+}\to\mathbb{R}$ rather
than an instantaneous value of $I$ as an input. Moreover, the function
$r_{\alpha}(t):\mathbb{R}_{+}\to\\{0,1\\}$ depends both on the input
$I(t):\mathbb{R}_{+}\to\mathbb{R}$ and the initial state
$r_{\alpha}(0)=r_{\alpha}^{0}$ of the relay. By definition (7) of the input-
to-state map (8), the state satisfies the constraints

$r_{\alpha}(t)=1\quad\text{whenever}\quad I(t)\geq\alpha_{2};\quad
r_{\alpha}(t)=0\quad\text{whenever}\quad I(t)\leq\alpha_{1}$ (9)

at all times. Further, the function (7) has at most a finite number of jumps
between the values $0$ and $1$ on any finite time interval $t_{0}\leq t\leq
t_{1}$.

Using formula (7) and notation (8), the switching rule for the basic
reproduction number depicted in Figure 1 can be expressed as

$R_{0}(t)=R_{0}^{nat}-(R_{0}^{nat}-R_{0}^{int})\cdot(\mathcal{R}_{\alpha}[r_{\alpha}^{0}]I)(t),\qquad
t\geq 0,$ (10)

Combining equations (2) with formula (10) for $R_{0}=R_{0}(t)$, one obtains a
switched system. According to the interpretation discussed in the
Introduction, the switched system (2), (10) models the dynamics of the
epidemic under the assumption that the health authorities implement the two-
threshold intervention policy — an intervention begins when the density of
infection (the number of active cases) exceeds the threshold value
$\alpha_{2}$ and is revoked when the density of infection drops below the
lower threshold value $\alpha_{1}$. It is assumed that the intervention
quickly translates into the reduction of the basic reproduction number; once
the intervention stops, $R_{0}$ returns to the larger value.

As shown in Chladná et al. (2020); Rouf and Rachinskii (2021), each trajectory
of switched system (2), (10) from the invariant domain $I,S\geq 0$, $I+S\leq
1$ converges to either an equilibrium state or a periodic orbit. In
particular, depending on the positioning of the thresholds
$\alpha_{1},\alpha_{2}$, the system can have either one attractor (an
equilibrium state or a periodic orbit), or two attractors (two co-existing
stable equilibrium states or a stable equilibrium state co-existing with a
stable periodic orbit), or three attractors (two equilibrium states and one
periodic orbit). For a periodic orbit with $I(t)$ oscillating between the
values $\alpha_{1}$ and $\alpha_{2}$, the point $(I(t),R_{0}(t))$ moves
clockwise along the rectangular hysteresis loop shown in Figure 1.

2.3. Model (10) corresponds to an ideally homogeneous response of the
population when all the individuals change their behavior simultaneously
following the recommendation of the health authority. Now, we consider a
model, in which several responses of the form (10), with different thresholds
$\alpha$, are combined because different groups of individuals respond
differently to the advent and dynamics of an epidemic and to the interventions
of the health authorities.

In particular, the ability and willingness to follow the recommendations of
the health authority can vary significantly from one to another group of
individuals for the same level of threat of contracting the disease. Multiple
factors are at play such as the occupation, age, living environment and health
conditions of an individual, to mention a few. In order to account for the
heterogeneity of the individual response, let us divide the susceptible
population into non-intersecting subpopulations $S_{\alpha}$ parameterized by
points $\alpha$ of a subset
$\Pi\subset\\{\alpha=(\alpha_{1},\alpha_{2}):\alpha_{1}<\alpha_{2}\\}$ of the
$\alpha$-plane, assuming a homogeneous response within each group
$S_{\alpha}$. As a simplification, let us assume that the response of the
infected population is homogeneous, i.e. all infected individuals act the same
and don’t change their behaviors over time, hence the heterogeneity of the
transmission coefficient is due to the behavior of the susceptible individuals
only. Further, assume that the basic reproduction number for the subpopulation
$S_{\alpha}$ is given by (10) and the contact structure is separable. Then,
the average basic reproduction number for the entire population at time $t$
equals

$R_{0}(t)=R_{0}^{nat}-(R_{0}^{nat}-R_{0}^{int})\iint_{\Pi}({\mathcal{R}}_{\alpha}[r_{\alpha}^{0}]I)(t)\,dF(\alpha),\qquad
t\in\mathbb{R}_{+},$ (11)

where the cumulative distribution function $F(\alpha)$ describes the
distribution of the susceptible population $S$ over the index set $\Pi$ (the
set of threshold pairs $\alpha=(\alpha_{1},\alpha_{2})$). The effective
reproduction number at time $t$ is then $R(t)=R_{0}(t)S(t)$. Finally, we
assume for simplicity that the distribution function $F(\alpha)$ is
independent of time, i.e. $F(\alpha)$ does not change with variations of $I$.
In this case, the mapping of the space of continuous inputs
$I(t):\mathbb{R}_{+}\to\mathbb{R}$ to the space of outputs
$R_{0}(t):\mathbb{R}_{+}\to\mathbb{R}$ defined by (11) is known as the
Preisach operator (see e.g. Krasnosel’skii and Pokrovskii 1989).

Under the above assumptions, the dynamics of the epidemic is modeled by system
(2) where $R_{0}(t)$ is related to $I(t)$ by the Preisach operator (11), which
accounts for the heterogeneity of the transmission coefficient among different
subgroups of the susceptible population.

A similar system results from the assumption that the health authorities have
multiple intervention policies (numbered $n=1,\ldots,N$) in place, each
decreasing the transmission coefficient by a certain amount $\Delta\beta_{n}$
while the intervention is implemented (Rouf and Rachinskii 2021). The
authorities aim to provide an adaptive response, which is adequate to the
severity of the epidemic. Let us suppose that each intervention policy is
guided by the two-threshold start/stop rule, such as in (10), associated with
a particular pair of thresholds $\alpha^{n}=(\alpha_{1}^{n},\alpha_{2}^{n})$.
The response of the population is assumed to be homogeneous. Under these
assumptions, the basic reproduction number of system (2) is given by

$R_{0}(t)=R_{0}^{nat}-\sum_{n=0}^{N-1}\hat{q}_{n}\cdot({\mathcal{R}}_{\alpha^{n}}[r_{\alpha^{n}}^{0}]I)(t),$
(12)

where $\hat{q}_{n}=\Delta\beta_{n}/(\gamma+\mu)$; hence, $R_{0}$ is set to
change at multiple thresholds $\alpha_{1}^{n},\alpha_{2}^{n}$. Operator (12),
known as the discrete Preisach model, is a particular case of (11) with an
atomic measure $F$.

Below we consider absolutely continuous distributions $F(\alpha)$ (see
Appendix 5.1 for details). The corresponding operator (11), which is called
the continuous Preisach model, can be written in the equivalent form

$R_{0}(t)=R_{0}^{nat}-(R_{0}^{nat}-R_{0}^{int})\iint_{\Pi}q(\alpha)\,\bigl{(}{\mathcal{R}}_{\alpha}[r_{\alpha}^{0}]I\bigr{)}(t)\,d\alpha_{1}d\alpha_{2},\qquad
t\geq 0,$ (13)

where $0<R_{0}^{int}<R_{0}^{nat}$ and $q=q(\alpha):\Pi\to\mathbb{R}_{+}$ is a
strictly positive measurable and bounded probability density function, i.e.

$\iint_{\Pi}q(\alpha)\,d\alpha_{1}d\alpha_{2}=1.$ (14)

This operator can be approximated by discrete operators (12). We will assume
that $q$ is bounded. Also, in what follows,

$\Pi=\\{\alpha=(\alpha_{1},\alpha_{2}):0\leq\alpha_{1}<\alpha_{2}\leq 1\\}$
(15)

and

$R_{0}^{nat}>R_{0}^{int}>1.$ (16)

The input-output diagram of the Preisach operator in Figure 2 illustrates the
relation between $I(t)$ and $R_{0}(t)$ defined by equation (13). One can see
that the point $(I,R_{0})$ switches from one branch $R_{0}(I)$ to another at
every turning point of the input $I(t)$. The presence of multiple branches and
hysteresis loops is a manifestation of the history-dependence, i.e. the output
value $R_{0}(t)$ at a time $t$ is determined not only by the concurrent value
$I(t)$ of the input but also by some of the past input values $I(\tau)$
achieved on the time interval $\tau\in[0,t]$. Fundamental properties of the
history-dependence, which are typical of hysteresis operators, are briefly
discussed in Appendix 5.1.

Figure 2: An example of the input-output diagram for the Preisach operator
(13) on panel (a) for the input $I(t)$ shown on panel (b). The turning points
on panel (a) correspond to local extrema of the input $I(t)$ on panel (b).

##  3 Results

3.1. Trajectories of system (2) coupled with the operator relationship (13)
lie in the infinite-dimensional phase space $\mathfrak{U}$ of triplets
$(I,S,r_{\alpha})$, where the measurable function
$r_{\alpha}=r(\alpha):\Pi\to\\{0,1\\}$ of the variable
$\alpha=(\alpha_{1},\alpha_{2})$ describes the states of the relays
${\mathcal{R}}_{\alpha}$ at a given moment (see Appendix 5.2 for details).
Slightly abusing the notation, we will also refer to the two-dimensional curve
$(I(t),S(t))$ as a trajectory, omitting the component (8) in the state space
of the Preisach operator.

Let us first consider equilibrium states of system (2), (13). The components
$I$, $S$, $r_{\alpha}$ of the solution and the basic reproduction number
$R_{0}$ at an equilibrium state are constant, and $R_{0}$ is related to the
function $r_{\alpha}=r(\alpha):\Pi\to\\{0,1\\}$ by the equation

$R_{0}=R_{0}^{nat}-(R_{0}^{nat}-R_{0}^{int})\iint_{\Pi}q(\alpha)r_{\alpha}\,d\alpha_{1}d\alpha_{2}.$
(17)

Due to assumption (15) and the compatibility constraint (9), the inclusion
$\alpha\in\Pi$ implies that all the relays are in state $r_{\alpha}=0$ when
$I=0$. Therefore, system (2), (13) has a unique infection free equilibrium
state

$E^{*}=(I^{*},S^{*})=(0,1),$ (18)

in which the function $r_{\alpha}:\Pi\to\\{0,1\\}$ is the identical zero and
the basic reproduction number equals $R_{0}=R_{0}^{nat}$ according to (17). On
the other hand, it is easy to see that due to assumption (16), system (2),
(13) also has a connected continuum of endemic equilibrium states

$E^{\theta}_{*}=(I_{*}^{\theta},S_{*}^{\theta})=\left(\left(1-\frac{1}{R_{0}^{\theta}}\right)\rho,\frac{1}{R_{0}^{\theta}}\right),\qquad
0\leq\theta\leq 1,$ (19)

which are characterized by different proportions of the infected, susceptible
and recovered populations and different values $R_{0}^{\theta}$ of the basic
reproduction number (Rachinskii and Rouf 2021).

3.2. The phase space $\mathfrak{U}$ of system (2), (13) is naturally embedded
into the space $\mathbb{R}^{2}\times L_{1}(\Pi;\mathbb{R})$ and is endowed
with a metric by this embedding (see Appendix 5.2 for details). As such, local
and global stability properties of trajectories of system (2), (13) are
defined in the standard way. We consider the positively invariant set

$\mathfrak{U}_{0}=\\{(I,S,r_{\alpha})\in\mathfrak{U}:I>0,S>0,I+S\leq 1\\}.$

###### Theorem 1

For any given $\rho>0$, $R_{0}^{nat}>1$ and a strictly positive probability
density function $q=q(\alpha):\Pi\to\mathbb{R}_{+}$ of the Preisach operator,
there is an $\varepsilon\in(0,R_{0}^{nat}-1)$ such that if $1+\varepsilon\leq
R_{0}^{int}<R_{0}^{nat}$, then the set of endemic equilibrium states of system
(2), (13) is globally asymptotically stable in $\mathfrak{U}_{0}$. Moreover,
each trajectory from $\mathfrak{U}_{0}$ converges to an endemic equilibrium
state. In particular, the components $I(t),S(t)$ of any trajectory from
$\mathfrak{U}_{0}$ satisfy $(I(t),S(t))\to E_{*}^{\theta}$ for some
$\theta\in[0,1]$.

The proof presented in the Appendix provides an explicit estimate for the
interval of values of $R_{0}^{int}\in[1+\varepsilon,R_{0}^{nat})$, for which
the set of endemic equilibrium states is globally stable. In particular, the
global stability is guaranteed if the quantities $\epsilon_{0}$ and $\kappa$
defined by (27), (64) using (23), (62) are both positive.

##  4 Conclusion

We considered an SIR model where the transmission coefficient changes in
response to dynamics of the epidemic. We aimed to emphasize that the adaptive
response of the population to the state of the epidemic is typically
heterogeneous and that this response can feature history-dependence in a form,
which is similar to a hysteresis effect. Specifically, we adapted a two-state
two-threshold hysteretic switch as a model of the adaptive response of an
individual to the varying number of active cases. In an ideally homogeneous
population, the two switching thresholds of the transmission coefficient can
be imposed by the health authority which starts the intervention when the
number of active cases exceeds a threshold $\alpha_{1}$ and ends the
intervention when the number of active cases drops below another threshold
$\alpha_{2}$. In order to account for the possibility of a heterogeneous
response among the susceptible individuals, we allowed a distribution of
switching thresholds and modeled the aggregate response of the susceptible
population by the Preisach operator. The variance of the distribution,
$\sigma^{2}$, measured the degree of heterogeneity of the public response to
the interventions.

The resulting heterogeneous model is shown to have a continuum of endemic
equilibrium states differing by the proportions of susceptible, infected and
recovered populations. Numerical simulations of the model provided earlier by
Rouf and Rachinskii, 2021 suggested that increasing the heterogeneity of the
response by widening the spread of thresholds of different population groups
(a larger variance $\sigma^{2}$) leads to a significant increase of the peak
of infection during the epidemic. Moreover, a higher degree of heterogeneity
of the public response also tends to steer the epidemic trajectory to an
endemic equilibrium state with a higher prevalence and a lower proportion of
the susceptible population. In other words, a more homogeneous response of the
public to transmission prevention measures helps “flattening the curve” and
can lead to lower prevalence and transmission rate when an endemic equilibrium
state is reached after the epidemic. These results suggest that intervention
programs are more effective when accompanied by education campaigns, which
convey the importance of the intervention measures to the public thus ensuring
a more homogeneous response from different population groups.

However, when the history-dependent two-threshold response in our model is
ideally homogeneous, i.e. the transmission coefficient behaves as a lazy-
switch, some threshold pairs lead to the convergence of the epidemic
trajectory to a periodic orbit predicting recurrent outbreaks of the epidemic
(Chladná et al. 2020). Emergence of sustained periodic oscillations of the
prevalence was also shown numerically in a similar SIR model with vaccination
when the heterogeneity of the response (measured by the variance $\sigma^{2}$
of the probability density function of the Preisach operator) decreases
(Kopfová et al. 2021). This scenario is in line with earlier findings by
Reluga, Bauch and Galvani (2006) who identified stable oscillatory vaccination
dynamics in heterogeneous populations and concluded that the more homogeneous
the response of a population, the less likely vaccine uptake is to be stable.
The stability theorem presented in this paper agrees with all these results
suggesting that the heterogeneity of the adaptive response promotes the
convergence of the epidemic trajectory to an equilibrium state. Indeed, the
conditions of the theorem, which ensure the global stability of the endemic
equilibrium set, require a sufficient degree of heterogeneity of the response.
On the other hand, for a sufficiently homogeneous response, these conditions
are violated leaving a possibility of undesirable non-equilibrium dynamics,
i.e. a too homogeneous response can potentially trigger epidemic instability.

Thus, our results support the argument that a more homogeneous public
response, which helps “flattening the curve” and tends to lead to a lower
prevalence at the endemic state after the epidemic, comes at the price of a
slower transitioning to the endemic state and a higher risk of epidemic
instability. A higher degree of heterogeneity of the public response can
stabilize and accelerate transitioning to an endemic state but after a higher
infection peak; it also potentially results in a higher prevalence at the
endemic equilibrium state after the epidemic.

The mathematical contribution of this paper is a method of global stability
analysis, which uses a family of Lyapunov functions corresponding to different
branches of the hysteresis operator. We applied this method for the first time
to an epidemic model with an adaptive transmission rate. The method is
potentially useful for analysis of other systems with hysteresis operators,
which have a continuum of equilibrium states.

##  5 Appendix

5.1. Continuous Preisach model. Let us briefly recall a rigorous definition of
the continuous Preisach operator (13) (Krasnosel’skii and Pokrovskii 1989). It
involves a collection of non-ideal relays ${\mathcal{R}}_{\alpha}$, which
respond to the same continuous input $I(t):\mathbb{R}_{+}\to\mathbb{R}$
independently according to formula (7). The relays contributing to the system
have different pairs of thresholds $\alpha=(\alpha_{1},\alpha_{2})\in\Pi$,
where the subset $\Pi$ of the half-plane
$\\{\alpha=(\alpha_{1},\alpha_{2}):\alpha_{1}<\alpha_{2}\\}$ is assumed to be
measurable and bounded; the $\alpha$-plane is called the Preisach plane. The
output of the continuous Preisach model is the scalar-valued function
$R_{0}(t):\mathbb{R}_{+}\to\mathbb{R}$ defined by (13), where
$q(\alpha):\Pi\to\mathbb{R}_{+}$ is a positive bounded measurable function
(probability density function) representing the weights of the relays; and,
$r_{\alpha}^{0}$ is the initial state of the relay ${\mathcal{R}}_{\alpha}$
for any given $\alpha\in\Pi$. The function
$r_{\alpha}^{0}=r^{0}(\alpha):\Pi\to\\{0,1\\}$ of the variable
$\alpha=(\alpha_{1},\alpha_{2})$ is referred to as the initial state function
of the Preisach operator. It is assumed to be measurable and satisfy the
constraints (5), (6), in which case the initial state-input pair is called
compatible. These requirements ensure that the integral in (13) is well-
defined for each $t\geq 0$ and, furthermore, the output
$R_{0}(t):\mathbb{R}_{+}\to\mathbb{R}$ of the Preisach model is a continuous
function of time.

The function (7) with a fixed $t\geq 0$ and varying $\alpha\in\Pi$ is
interpreted as the state function of the Preisach model at the moment $t$
because it describes the states of all the relays at this moment; this state
function $r_{\alpha}(t)=r(t;\alpha):\Pi\to\\{0,1\\}$ is an element of the
space $L_{1}(\Pi;\mathbb{R})$ for each $t\geq 0$. In this sense, formula (7)
defines the evolution of the state function in the space
$L_{1}(\Pi;\mathbb{R})$ for a given input provided that the initial state-
input pair is compatible. Moreover, the state function
$r_{\alpha}(t)=r(t;\alpha):\Pi\to\\{0,1\\}$ and the concurrent input value
$I(t)$ are compatible for each $t\geq 0$ and the semigroup property holds:

$r_{\alpha}(\tau)=({\mathcal{R}}_{\alpha}[r^{0}_{\alpha}]I)(\tau)\ \ \
\Rightarrow\ \ \ r_{\alpha}(t)=({\mathcal{R}}_{\alpha}[r_{\alpha}(\tau)]I)(t)\
\ \ \text{for any}\ \ \ \alpha\in\Pi,\ t\geq\tau\geq 0.$ (20)

For brevity, let us denote the input-to-output operator of the Preisach model
defined by (13) as

$R_{0}(t)=({\mathcal{P}}[r^{0}_{\alpha}]I)(t),\qquad t\geq 0,$ (21)

where both the input $I(t):\mathbb{R}_{+}\to\mathbb{R}$ and the initial state
function $r^{0}_{\alpha}=r^{0}({\alpha}):\Pi\to\\{0,1\\}$ (which is compatible
with the initial value of the input) are the arguments; the value of this
operator is the output $R_{0}(t):\mathbb{R}_{+}\to\mathbb{R}$. An important
property of the Preisach operator (13) is that it is Lipschitz continuous if
$q(\alpha):\Pi\to\mathbb{R}_{+}$ is bounded (Krasnosel’skii et al. 1989). More
precisely, the relations

$R_{0}^{k}(t)=({\mathcal{P}}[r^{0,k}_{\alpha}]I^{k})(t),\qquad t\geq 0,{\qquad
k=1,2,}$

and $0\leq I^{1}(t),I^{2}(t)\leq 1$ ($t\geq 0$) imply

$\|R_{0}^{1}(t)-R_{0}^{2}(t)\|_{C([0,\tau];\mathbb{R})}\leq
q_{0}\Big{(}\|r^{0,1}_{\alpha}-r^{0,2}_{\alpha}\|_{L_{1}(\Pi;\mathbb{R})}+\|I^{1}(t)-I^{2}(t)\|_{C([0,\tau];\mathbb{R})}\Big{)}$
(22)

for any $\tau\geq 0$ with the Lipschitz constant

$q_{0}:=(R_{0}^{nat}-R_{0}^{int})\,\sup_{\alpha\in\Pi}q(\alpha).$ (23)

5.2. Phase space of the system with the Preisach operator. Let $\mathfrak{U}$
be the set of all triplets $(I_{0},S_{0},r_{\alpha}^{0})$, where
$(I_{0},S_{0})\in\mathfrak{D}=\\{(I,S):0\leq I,S;I+S\leq 1\\}$ and the state
function $r^{0}_{\alpha}=r^{0}(\alpha):\Pi\to\\{0,1\\}$ of the Preisach
operator is compatible with $I_{0}$. The set $\mathfrak{U}$ is a natural phase
space for system (2), (13). In particular, the global Lipschitz estimate (22)
ensures (for example, using the Picard-Lindelöf type of argument) that for a
given $(I_{0},S_{0},r^{0}_{\alpha})\in\mathfrak{U}$, system (2), (13) has a
unique local solution with the initial data $I(0)=I_{0},S(0)=S_{0}$ and the
initial state function $r^{0}_{\alpha}$ of the Preisach opeartor (see e.g. the
survey Leonov et al. 2017). This solution satisfies
$(I(t),S(t),r_{\alpha}(t))\in\mathfrak{U}$ for each $t$ in its domain.
Further, the positive invariance of $\mathfrak{D}$ implies that each solution
is extendable to the whole semi-axis $t\in\mathbb{R}_{+}$. Due to the
semigroup property (20), these solutions induce a continuous semi-flow in the
phase space $\mathfrak{U}$, which is endowed with a metric by the natural
embedding into the space $\mathbb{R}^{2}\times L_{1}(\Pi;\mathbb{R})$. This
construction leads to the standard definition of local and global stability
including stability of equilibrium states and periodic solutions. In
particular, an equilibrium is a triplet
$(I_{0},S_{0},r^{0}_{\alpha})\in\mathfrak{U}$ and a periodic solution is a
periodic function $(I(t),S(t),r_{\alpha}(t)):\mathbb{R}_{+}\to\mathfrak{U}$
where the last component, viewed as a function
$r_{\alpha}(t)=r(t;\alpha):\mathbb{R}_{+}\times\Pi\to\\{0,1\\}$ of two
variables $t\in\mathbb{R}_{+}$ and $\alpha\in\Pi$, is given by (8). The basic
reproduction number (13) at an equilibrium is constant, while for a periodic
solution the basic reproduction number $R_{0}(t)$ and the state function of
the Preisach operator change periodically with the period of $I(t)$ and
$S(t)$.

5.3. Switched system.  A fundamental property of hysteresis operators, which
is often used as their definition (see, e.g. Visisntin 2006; Kluwer
Encyclopedia of Mathematics 1997), is rate-independence. It means that for any
increasing continuous function $\theta(t):\mathbb{R}_{+}\to\mathbb{R}_{+}$
(transformation of time) inputs $I(t)$ and $I(\theta(t))$ produce identical
input-output diagrams such as in Figure LABEL:fig_loops. Due to the rate-
independence property of the Preisach operator, given an initial state
function $r^{0}=r^{0}_{\alpha}:\Pi\to\\{0,1\\}$, there is a continuous
function $R_{r^{0}}:\mathbb{R}\to\mathbb{R}$ such that every monotone input
$I(t):\mathbb{R}_{+}\to\mathbb{R}$ (which is compatible with the state
function at the initial moment) and the corresponding output (21) of the
Preisach operator are related by

$R_{0}(t)=R_{r^{0}}(I(t)),\qquad t\geq 0.$

The function $R_{r^{0}}(I)$ will be referred to as a branch of the Preisach
operator. Hence, the state function $r^{0}=r^{0}_{\alpha}$ parameterizes the
infinite set of all possible branches. Due to the assumption that the density
$q=q(\alpha)$ of the Preisach operator is positive, every branch
$R_{r^{0}}(I)$ is a decreasing function on the interval $0\leq I\leq 1$. Along
with $R_{r^{0}}(I)$, the functions

$f_{r^{0}}(I)=IR_{r^{0}}(I)$ (24)

will play an important role.

With the above notation, if the first component $I(t)$ of a solution
$(I(t),S(t),r(t)):\mathbb{R}_{+}\to\mathfrak{U}$ of system (2), (13) is
monotone on a time interval $[t_{k},t_{k+1}]$, then the pair $(I(t),S(t))$ is
a solution of system

$\displaystyle\dot{I}$ $\displaystyle=R(I)SI-I,$ (25) $\displaystyle\dot{S}$
$\displaystyle=-R(I)SI-\rho S+\rho$

with

$R(I)=R_{r^{k}}(I),$ (26)

where $r^{k}=r(t_{k})$, on the same time interval. This ordinary differential
system will be called a branch of the operator-differential system (2), (13);
again, the state function $r^{k}=r^{k}_{\alpha}$ of the Preisach operator
parameterizes the set of ordinary differential systems (25), (26).

We see that dynamics of system (2), (13) can be interpreted as switched
dynamics of planar systems (25), (26), where the switching moments are the
turning points of the component $I$ of a trajectory. As such, a switch occurs
when a trajectory $(I(t),S(t))$ of system (25), (26) (with a fixed $r^{k}$)
reaches the nullcline $\dot{I}=0$, i.e. at the moment when
$R_{r^{k}}(I(t))S(t)=1$ (notice that the nullcline depends on $r^{k}$).

In what follows, the pair $(I(t),S(t))$ along with the triplet
$(I(t),S(t),r(t)):\mathbb{R}_{+}\to\mathfrak{U}$ will be referred to as a
solution of system (2), (13) whenever this doesn’t cause confusion.

5.4. Proofs.  Throughout the proofs we assume that $R_{0}^{nat}-R_{0}^{int}>0$
is sufficiently small to ensure that the quantities $\epsilon_{0},\kappa$
defined below by (27), (64) using (23), (62) satisfy $\epsilon_{0},\kappa>0$.

We first recall well-known stability properties of the planar ordinary
differential system (25), which are ensured by the following lemmas.

###### Lemma 1

If

$\epsilon_{0}:=R_{0}^{int}-q_{0}>0$ (27)

with $q_{0}$ given by (23), then for any branch $R(I)$ of the Preisach
operator, the derivative of the function $f(I)=IR(I)$ satisfies

$f^{\prime}(I)\geq\epsilon_{0}.$

Proof.  From (22), it follows that $-q_{0}\leq R^{\prime}(I)\leq 0$. Combining
this with $R(I)\geq R_{0}^{int}$ and $I\leq 1$, we obtain
$f^{\prime}(I)=R(I)+IR^{\prime}(I)\geq\epsilon_{0}$. $\Box$

###### Lemma 2

System (25), where $R(I)$ is any branch of the Preisach operator, has a unique
infection-free equilibrium $(I^{*},S^{*})=(0,1)$ and a unique positive
(endemic) equilibrium $(I_{*},S_{*})$, where $I_{*}$ is a unique solution of

$\frac{1}{R(I)}=1-\frac{I}{\rho}$ (28)

and $S_{*}=1/R(I_{*})$ (cf. (18), (19)). The infection-free equilibrium is
unstable, while the positive equilibrium $(I_{*},S_{*})$ is globally stable in
the positive quadrant $I,S>0$.

Proof.  By assumption, $1/R(0)=1/R_{nat}<1$ and $1/R(\rho)>0$, hence equation
(28) has a solution $I_{*}$ in the interval $[0,\rho]\subset[0,1]$. The
measure density $q(\alpha)$ of the Preisach operator (13) is assumed to be
positive, which implies that $R(I)$ is a decreasing function. Since the
function $1/R(I)$ increases, while the function $1-I/\rho$ decreases, this
solution is unique.

The Jacobi matrix of system (25) at the infection-free equilibrium $(0,1)$ has
eigenvalues $\lambda_{1}=R_{nat}-1,\lambda_{2}=-R_{nat}-\rho$, hence this
equilibrium is a saddle due to the assumption $R_{nat}>1.$

At the endemic equilibrium $(I_{*},S_{*})$, the characteristic polynomial of
the Jacobi matrix of system (25) equals

$P(\lambda)=\lambda^{2}+\lambda\bigl{(}I_{*}R(I_{*})+\rho-
I_{*}R^{\prime}(I_{*})S_{*}\bigr{)}+I_{*}R(I_{*})-\rho
I_{*}R^{\prime}(I_{*})S_{*}.$

Since $R(I)$ decreases, all coefficients of this polynomial are positive,
hence all eigenvalues have negative real parts, and the endemic equilibrium is
asymptotically stable.

We claim that

$V(I,S)=\int_{I_{*}}^{I}\left(1-\frac{f(I_{*})}{f(i)}\right)\,di+\int_{S_{*}}^{S}\left(1-\frac{S_{*}}{s}\right)\,ds$
(29)
$=I-I_{*}+S-S_{*}-f(I_{*})\int_{I_{*}}^{I}\frac{di}{f(i)}-S_{*}\ln{\frac{S}{S_{*}}}.$

is a Lyapunov function for system (25) in the positive quadrant.

Indeed, $V(I_{*},S_{*})=0$, and Lemma 1 implies that $V(I,S)>0$ for $I,S>0$
except at the point $(I_{*},S_{*})$. Further, $V\to\infty$ if either
$I+S\to\infty$ or $S\to 0$ or $I\to 0$ (the latter follows from $f(0)=0$). On
the other hand, using the notation

$R_{*}=R(I_{*}),\qquad f_{*}=f(I_{*})=I_{*}R(I_{*})$

and $R=R(I)$, $f=f(I)$, the time derivative of $V$ along a trajectory of
system (25) equals

$\begin{array}[]{lll}\dot{V}&=&\displaystyle\left(1-\frac{S_{*}}{S}\right)\rho(1-S)-I+S_{*}f+\frac{f_{*}I}{f}-Sf_{*}\\\
&=&\displaystyle\frac{\rho}{S}(S-S_{*})(1-S)-f_{*}(S-S_{*})+S_{*}(f-f_{*})-I\left(1-\frac{f_{*}}{f}\right)\\\
&=&\displaystyle\frac{\rho}{S}(S-S_{*})(1-S)-(S-S_{*})\frac{\rho(1-S_{*})}{S_{*}}+(S_{*}f-I)\left(1-\frac{f_{*}}{f}\right),\end{array}$

where we used that $I_{*}=\rho(1-S_{*})=\rho(1-I_{*}/f_{*})$ due to Lemma 2.
Equivalently,

$\begin{array}[]{lll}\dot{V}&=&\displaystyle\rho(S-S_{*})\left(\frac{1}{S}-\frac{1}{S_{*}}\right)+\left(\frac{I_{*}f}{f_{*}}-I\right)\left(1-\frac{f_{*}}{f}\right)\\\
&=&\displaystyle-\frac{\rho}{SS_{*}}(S-S_{*})^{2}+\frac{I_{*}}{f}(f-f_{*})\cdot\frac{I}{f_{*}}\left(\frac{f}{I}-\frac{f_{*}}{I_{*}}\right)\\\
&=&\displaystyle-\frac{\rho}{SS_{*}}(S-S_{*})^{2}+\frac{1}{RR_{*}}(f-f_{*})(R-R_{*}).\end{array}$

Since the function $R=R(I)$ decreases and the function $f=f(I)$ increases
(Lemma 1), we conclude that

$\dot{V}\leq-\frac{\rho}{SS_{*}}(S-S_{*})^{2}<0,\qquad(I,S)\neq(I_{*},S_{*}).$
(30)

Hence, the equilibrium $(I_{*},S_{*})$ is globally stable in the positive
quadrant. $\Box$

It turns out that the function $V$ is convex.

###### Lemma 3

All the level sets $\\{(I,S):V(I,S)<c\\}$, $c>0$, of the function $V$ are
convex.

Proof.  Denoting the partial derivatives of the function $V$ by $V_{I},V_{S}$,
the curvature of the level line $\\{(I,S):V(I,S)=c\\}$ of $V$ for any $c>0$ is
given by

$\kappa=-\frac{V_{SS}(V_{I})^{2}+V_{II}(V_{S})^{2}-2V_{IS}V_{I}V_{S}}{(V_{S}^{2}+V_{I}^{2})^{3/2}}=-\frac{V_{SS}(V_{I})^{2}+V_{II}(V_{S})^{2}}{(V_{S}^{2}+V_{I}^{2})^{3/2}},$

because $V_{IS}=0$. Using that

$V_{I}=1-\frac{f(I_{*})}{f(I)},\quad V_{S}=1-\frac{S_{*}}{S},\quad
V_{II}=\frac{f^{\prime}(I)f(I_{*})}{f^{2}(I)},\quad
V_{SS}=\frac{S_{*}}{S^{2}}$

and that $f^{\prime}(I)>0$ due to Lemma 1, we obtain $\kappa<0$, which
completes the proof. $\Box$

Let us consider the trajectory
$(I(\cdot),S(\cdot),r(\cdot)):\mathbb{R}_{+}\to\mathfrak{U}$ of system (2),
(13) and the corresponding planar curve
$(I(\cdot),S(\cdot)):\mathbb{R}_{+}\to\mathbb{R}^{2}$. Consider the sequence
$\mathbb{T}=\\{t_{k}\\}$ of successive switching moments
$0=t_{0}<t_{1}<t_{2}<\cdots$ such that

  * •

The restriction of $(I(\cdot),S(\cdot)):\mathbb{R}_{+}\to\mathbb{R}^{2}$ to
each interval $(t_{k},t_{k+1})$, $t_{k},t_{k+1}\in\mathbb{T}$ is a solution of
system (25) with $R(I)=R_{r(t_{k})}(I)$ satisfying

$R_{r(t_{k})}(I(t))S(t)\neq 1,\quad t\in(t_{k},t_{k+1});\qquad
R(I(t_{k+1}))S(t_{k+1})=1.$ (31)

  * •

In addition, if the set $\mathbb{T}=\\{t_{0},\ldots,t_{N}\\}$ is finite, then
the restriction of $(I(\cdot),S(\cdot)):\mathbb{R}_{+}\to\mathbb{R}^{2}$ to
the infinite interval $(t_{N},\infty)$ is a solution of system (25) with
$R(I)=R_{r(t_{N})}(I)$ satisfying

$R_{r(t_{N})}(I(t))S(t)\neq 1,\quad t>t_{N}.$

If the set $\mathbb{T}$ is finite, the trajectory of system (2), (13) follows
the flow (25), (26) on the interval $[t_{N},\infty)$ and due to Lemma 2
converges to an equilibrium. In this scenario, this equilibrium of (25), (26)
is a node.

We continue the proof assuming that $\mathbb{T}=\\{t_{k}\\}$ is infinite, i.e.
there is an infinite sequence of switches along the trajectory. Let

$I_{k}=I(t_{k}),\quad S_{k}=S(t_{k}).$

Denote by $V_{k}(I,S)$ the Lyapunov function (36) for system (25), (26) (where
$r^{k}=r(t_{k})$) and by $(I_{*}^{k},S_{*}^{k})$ the positive (endemic)
equilibrium of this system. We will also use the notation
$R_{k}(I)=R_{r(t_{k})}(I)$. The next lemma provides an estimate for the
(negative) increment

$\Delta V_{k}:=V_{k}(I_{k+1},S_{k+1})-V_{k}(I_{k},S_{k})$ (32)

of this function along the trajectory on the time interval $[t_{k},t_{k+1}]$.

###### Lemma 4

If $R_{r(t_{k})}(I(t))S(t)>1$ on $(t_{k},t_{k+1})$ (cf. (31)), then the
increment (32) of the Lyapunov function $V_{k}$ along the trajectory satisfies

$\Delta
V_{k}\leq-\frac{\rho}{4}(I_{k+1}-I_{*}^{k})\sqrt{S_{M}\int_{I_{*}^{k}}^{I_{k+1}}\left(1-\frac{f(I_{*}^{k})}{f(i)}\right)di},$
(33)

where $M=(I_{M},S_{M})$ is the highest point of the level line
$V_{k}(I,S)=V_{k}(I_{k+1},S_{k+1})$ of the Lyapunov function $V_{k}$, see
Figure 3.

Figure 3: An arc of a trajectory of system (2), (13) satisfying
$R_{k}(I(t))S(t)\geq 1$.

Proof. At the point $M=(I_{M},S_{M})$,

$0=\left.\frac{\partial V_{k}}{\partial
I}\right|_{M}=1-\frac{f(I_{*}^{k})}{f(I_{M})},$

which is equivalent to $f(I_{*}^{k})=f(I_{M})$. Since $f(I)$ increases, it
follows that $I_{M}=I_{*}^{k}$, hence $M=(I_{*}^{k},S_{M})$ and

$V_{k}(I_{*}^{k},S_{M})=S_{M}-S_{*}^{k}-S_{*}^{k}\ln{\frac{S_{M}}{S_{*}^{k}}}.$

Consider a point $(I_{b},S_{b})$ on the arc of the level line of $V_{k}$
connecting the points $M$ and $(I_{k+1},S_{k+1})$. Due to Lemma 3,

$\frac{S_{b}-S_{k+1}}{I_{k+1}-I_{b}}\geq\frac{S_{M}-S_{k+1}}{I_{k+1}-I_{*}^{k}}$

(see Figure 3). Equivalently,

$I_{k+1}-I_{b}\leq\frac{S_{b}-S_{k+1}}{S_{M}-S_{k+1}}(I_{k+1}-I_{*}^{k}).$

This implies the estimate

$I_{b}-I_{*}^{k}\geq
I_{k+1}-I_{*}^{k}-\frac{S_{b}-S_{k+1}}{S_{M}-S_{k+1}}(I_{k+1}-I_{*}^{k})=\frac{S_{M}-S_{b}}{S_{M}-S_{k+1}}(I_{k+1}-I_{*}^{k}).$
(34)

In the domain $I\geq I_{*}^{k},S\geq{1}/{R_{k}(I)},$ we have
${1}/{R_{k}(I)}\geq{1}/{R_{k}(I_{*}^{k})}=S_{*}^{k},$ which implies $0\leq
S-{1}/{R_{k}(I)}\leq S-S_{*}^{k}$. Hence, relation (30) implies that along the
trajectory,

$\displaystyle\frac{dV_{k}}{dI}=\frac{\dot{V}_{k}}{\dot{I}}\leq-\frac{\rho}{SS_{*}^{k}IR_{k}(I)}\cdot\frac{(S-S_{*}^{k})^{2}}{(S-1/R_{k}(I))}\leq-\frac{\rho}{SS_{*}^{k}IR_{k}(I_{*}^{k})}(S-S_{*}^{k})$
$\displaystyle=-\frac{\rho}{SI}(S-S_{*}^{k})\leq-\frac{\rho}{I}(S-S_{*}^{k}).$

Since $V_{k}$ decreases along the trajectory, the segment of the trajectory
with $t\in(t_{k},t_{k+1})$ lies above the level line of $V_{k}$, hence for
this segment $I\in[I_{*}^{k},I_{b}]$ implies $S\geq S_{b}$ and

$\frac{dV_{k}}{dI}\leq-\frac{\rho}{I}(S_{b}-S_{*}^{k}).$ (35)

The relation $\dot{I}(t)>0$, $t\in(t_{k},t_{k+1})$, implies that the arc of
the trajectory for this time interval admits a parameterization
$(I,\hat{S}(I))$, $I\in[I_{k},I_{k+1}]$. Therefore, integrating (35) over the
interval $[I_{k}^{*},I_{b}]$, we obtain

$V_{k}(I_{b},\hat{S}(I_{b}))-V_{k}(I_{*}^{k},\hat{S}(I_{*}^{k}))\leq-\rho(S_{b}-S_{*}^{k})\ln{\frac{I_{b}}{I_{*}^{k}}}.$

Further, using the identity $\ln(1+x)\geq{x}/(1+x),$ $x\geq-1$, we arrive at

$V_{k}(I_{b},\hat{S}(I_{b}))-V_{k}(I_{*}^{k},\hat{S}(I_{*}^{k}))\leq-\rho(S_{b}-S_{*}^{k})\frac{I_{b}-I_{*}^{k}}{I_{b}}\leq-\rho(S_{b}-S_{*}^{k})(I_{b}-I_{*}^{k})$

and due to (34),

$V_{k}(I_{b},\hat{S}(I_{b}))-V_{k}(I_{*}^{k},\hat{S}(I_{*}^{k}))\leq-\rho\frac{(S_{b}-S_{*}^{k})(S_{M}-S_{b})(I_{k+1}-I_{*}^{k})}{S_{M}-S_{k+1}}.$
(36)

The maximum of the function
$\frac{(S_{b}-S_{*}^{k})(S_{M}-S_{b})}{S_{M}-S_{k+1}}$ on the interval
$[S_{k+1},S_{M}]$ is achieved at $S_{k+1}=(S_{M}+S_{*}^{k})/{2},$ hence

$\max_{S_{b}\in[S_{k+1},S_{M}]}\frac{(S_{b}-S_{*}^{k})(S_{M}-S_{b})}{S_{M}-S_{k+1}}\geq\frac{(S_{M}-S_{*}^{k})^{2}}{4(S_{M}-S_{k+1})}\geq\frac{(S_{M}-S_{*}^{k})^{2}}{4(S_{M}-S_{*}^{k})}=\frac{S_{M}-S_{*}^{k}}{4}.$

Therefore, from (36) it follows that

$V_{k}(I_{b},\hat{S}(I_{b}))-V_{k}(I_{*}^{k},\hat{S}(I_{*}^{k}))\leq-\frac{\rho}{4}(I_{k+1}-I_{*}^{k})(S_{M}-S_{*}^{k}).$

Since $V_{k}$ decreases along the trajectory, this estimate implies

$\Delta V_{k}\leq-\frac{\rho}{4}(I_{k+1}-I_{*}^{k})(S_{M}-S_{*}^{k}).$ (37)

Finally, let us estimate $S_{M}-S_{*}^{k}.$ Since $V_{k}$ has the same value
at the points $(I_{*}^{k},S_{M})$ and $(I_{k+1},S_{k+1})$, we have

$S_{M}-S_{*}^{k}-S_{*}^{k}\ln{\frac{S_{M}}{S_{*}^{k}}}=I_{k+1}-I_{*}^{k}-f(I_{*}^{k})\int_{I_{*}^{k}}^{I_{k+1}}\frac{di}{f(i)}+S_{k+1}-S_{*}^{k}-S_{*}^{k}\ln{\frac{S_{k+1}}{S_{*}^{k}}},$

which implies

$S_{M}-S_{k+1}-S_{*}^{k}\ln{\frac{S_{M}}{S_{k+1}}}=\int_{I_{*}^{k}}^{I_{k+1}}\left(1-\frac{f(I_{*}^{k})}{f(i)}\right)di.$

We estimate the left hand side of this equation using
$\ln{(1+x)}\geq{x}/{1+x}$, $x\geq-1$:

$S_{M}-S_{k+1}-S_{*}^{k}\ln{\frac{S_{M}}{S_{k+1}}}\leq
S_{M}-S_{k+1}-\frac{S_{*}^{k}(S_{M}-S_{k+1})}{S_{M}}=\left(1-\frac{S_{*}^{k}}{S_{M}}\right)\left(S_{M}-S_{k+1}\right),$

hence

$\frac{1}{S_{M}}(S_{M}-S_{*}^{k})(S_{M}-S_{k+1})\geq\int_{I_{*}^{k}}^{I_{k+1}}\left(1-\frac{f(I_{*}^{k})}{f(i)}\right)di.$

Since $S_{*}^{k}<S_{k+1}$, this implies

$(S_{M}-S_{*}^{k})^{2}\geq
S_{M}\int_{I_{*}^{k}}^{I_{k+1}}\left(1-\frac{f(I_{*}^{k})}{f(i)}\right)di.$
(38)

Combining this relation with (37), we obtain (33). $\Box$

Next, we estimate the difference

$\tilde{\Delta}V_{k}:=V_{k+1}(I_{k+1},S_{k+1})-V_{k}(I_{k+1},S_{k+1})$ (39)

between the values of Lyapunov functions $V_{k+1}$ and $V_{k}$ at the point
$(I_{k+1},S_{k+1})$. We begin with the following statement.

###### Lemma 5

Under the assumptions of Lemma 4,

$0\leq I_{*}^{k}-I_{*}^{k+1}=\rho(S_{*}^{k+1}-S_{*}^{k})\leq\frac{\rho
q_{0}}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})}\,(I_{k+1}-I_{*}^{k}).$ (40)

Proof.  Consider the branch $R_{k}(I):=R_{r(t_{k})}(I)$ of the Preisach
operator (see Section 5.3), the corresponding function $f_{k}(I):=IR_{k}(I)$
and the endemic equilibrium $(I_{*}^{k},S_{*}^{k})$ of system (25), (26) for
this branch:

$S_{*}^{k}=1-\frac{I_{*}^{k}}{\rho}=\frac{1}{R_{k}(I_{*}^{k})}.$

These equations imply

$S_{*}^{k+1}-S_{*}^{k}=\frac{1}{\rho}(I_{*}^{k}-I_{*}^{k+1})$

and

$\displaystyle\frac{I_{*}^{k}-I_{*}^{k+1}}{\rho}=\frac{1}{R_{k+1}(I_{*}^{k+1})}-\frac{1}{R_{k}(I_{*}^{k})}=\frac{R_{k}(I_{*}^{k})-R_{k+1}(I_{*}^{k})+R_{k+1}(I_{*}^{k})-R_{k+1}(I_{*}^{k+1})}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})}$
$\displaystyle=\frac{1}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})}\bigr{(}R_{k}(I_{*}^{k})-R_{k+1}(I_{*}^{k})+R_{k+1}^{\prime}(\hat{I})(I_{*}^{k}-I_{*}^{k+1})\bigr{)}$

for an intermediate value of $\hat{I}$ between $I_{*}^{k}$ and $I_{*}^{k+1}$.

If we assume that $I_{*}^{k+1}>I_{*}^{k}$, then we arrive at a contradiction:

$0>\frac{1}{\rho}(I_{*}^{k}-I_{*}^{k+1})=\frac{1}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})}\bigr{(}R_{k}(I_{*}^{k})-R_{k+1}(I_{*}^{k})+R_{k+1}^{\prime}(\hat{I})(I_{*}^{k}-I_{*}^{k+1})\bigr{)}\geq
0$

because $R_{k+1}(I)$ decreases and (22) implies that

$0\leq R_{k}(I_{*}^{k})-R_{k+1}(I_{*}^{k})\leq q_{0}(I_{k+1}-I_{*}^{k}).$ (41)

Hence, $I_{*}^{k+1}\leq I_{*}^{k}$, and consequently $S_{*}^{k+1}\geq
S_{*}^{k}$ as well as

$0\leq
S_{*}^{k+1}-S_{*}^{k}=\frac{1}{\rho}(I_{*}^{k}-I_{*}^{k+1})\leq\frac{R_{k}(I_{*}^{k})-R_{k+1}(I_{*}^{k})}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})}.$
(42)

Combining (41) and (42), we obtain (40). $\Box$

We use Lemma 5 to prove the following estimate for the difference (39).

###### Lemma 6

Under the assumptions of Lemma 4, the difference (39) satisfies

$\tilde{\Delta}V_{k}\leq\frac{q_{0}(I_{k+1}-I_{*}^{k})^{2}}{R_{k}(I_{*}^{k})}\left(\frac{\rho^{2}q_{0}}{R_{k+1}^{2}(I_{*}^{k+1})f_{k}(I_{*}^{k})}+\frac{\rho}{f_{k}(I_{*}^{k})}+1\right).$
(43)

Proof. By the definition of the Lyapunov function,
$\tilde{\Delta}V_{k}=\Delta_{S}+\Delta_{I}$, where

$\displaystyle\Delta_{S}$ $\displaystyle=$
$\displaystyle\int_{S_{*}^{k+1}}^{S_{k+1}}\left(1-\frac{S_{*}^{k+1}}{s}\right)\,ds-\int_{S_{*}^{k}}^{S_{k+1}}\left(1-\frac{S_{*}^{k}}{s}\right)\,ds$
$\displaystyle=$
$\displaystyle-\int_{S_{*}^{k}}^{S_{*}^{k+1}}\left(1-\frac{S_{*}^{k}}{s}\right)\,ds-\int_{S_{*}^{k+1}}^{S_{k+1}}\left(\frac{S_{*}^{k+1}-S_{*}^{k}}{s}\right)\,ds,$
$\displaystyle\Delta_{I}$ $\displaystyle=$
$\displaystyle\int_{I_{*}^{k+1}}^{I_{k+1}}\left(1-\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(i)}\right)\,di-\int_{I_{*}^{k}}^{I_{k+1}}\left(1-\frac{f_{k}(I_{*}^{k})}{f_{k}(i)}\right)\,di$
$\displaystyle=$
$\displaystyle-\int_{I_{*}^{k}}^{I_{k+1}}\left(\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(i)}-\frac{f_{k}(I_{*}^{k})}{f_{k}(i)}\right)\,di-\int_{I_{*}^{k}}^{I_{*}^{k+1}}\left(1-\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(i)}\right)\,di$
$\displaystyle=$
$\displaystyle-\int_{I_{*}^{k}}^{I_{k+1}}\left(\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(i)}-\frac{f_{k+1}(I_{*}^{k+1})}{f_{k}(i)}\right)\,di-\int_{I_{*}^{k}}^{I_{k+1}}\frac{f_{k+1}(I_{*}^{k+1})-f_{k}(I_{*}^{k})}{f_{k}(i)}\,di$
$\displaystyle-$
$\displaystyle\int_{I_{*}^{k}}^{I_{*}^{k+1}}\left(1-\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(i)}\right)\,di.$

By Lemma 5, $S_{*}^{k}\leq S_{*}^{k+1}\leq S_{k+1}$, hence

$\int_{S_{*}^{k}}^{S_{*}^{k+1}}\left(1-\frac{S_{*}^{k}}{s}\right)\,ds\geq
0,\quad\int_{S_{*}^{k+1}}^{S_{k+1}}\left(\frac{S_{*}^{k+1}-S_{*}^{k}}{s}\right)\,ds\geq
0,$ (44)

and consequently $\Delta_{S}\leq 0$. Also, the relation $R_{k+1}(I)\leq
R_{k}(I)$, $I\leq I_{k+1}$, implies that

$\int_{I_{*}^{k}}^{I_{k+1}}\left(\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(i)}-\frac{f_{k+1}(I_{*}^{k+1})}{f_{k}(i)}\right)\,di=\int_{I_{*}^{k}}^{I_{k+1}}\frac{f_{k+1}(I_{*}^{k+1})}{i}\left(\frac{1}{R_{k+1}(i)}-\frac{1}{R_{k}(i)}\right)\,di\geq
0,$ (45)

therefore

$\tilde{\Delta}V_{k}\leq\Delta_{I}\leq-\int_{I_{*}^{k}}^{I_{*}^{k+1}}\left(1-\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(i)}\right)\,di-\int_{I_{*}^{k}}^{I_{k+1}}\frac{f_{k+1}(I_{*}^{k+1})-f_{k}(I_{*}^{k})}{f_{k}(i)}\,di.$
(46)

In order to estimate the first integral in this expression, we recall that
$f_{k+1}(I)$ increases, $R_{k+1}(I)$ decreases and $I_{*}^{k+1}\leq I_{*}^{k}$
due to Lemma 5, hence

$\displaystyle 0$ $\displaystyle\leq$
$\displaystyle-\int_{I_{*}^{k}}^{I_{*}^{k+1}}\left(1-\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(i)}\right)\,di\leq(I_{*}^{k}-I_{*}^{k+1})\left(1-\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(I_{*}^{k})}\right)$
(47) $\displaystyle=$
$\displaystyle(I_{*}^{k}-I_{*}^{k+1})\frac{R_{k+1}(I_{*}^{k})I_{*}^{k}-R_{k+1}(I_{*}^{k+1})I_{*}^{k+1}}{R_{k+1}(I_{*}^{k})I_{*}^{k}}\leq(I_{*}^{k}-I_{*}^{k+1})\frac{R_{k+1}(I_{*}^{k})(I_{*}^{k}-I_{*}^{k+1})}{R_{k+1}(I_{*}^{k})I_{*}^{k}}$
$\displaystyle=$ $\displaystyle\frac{(I_{*}^{k}-I_{*}^{k+1})^{2}}{I_{*}^{k}}$

and using Lemma 5 again,

$-\int_{I_{*}^{k}}^{I_{*}^{k+1}}\left(1-\frac{f_{k+1}(I_{*}^{k+1})}{f_{k+1}(i)}\right)\,di\leq\frac{\rho^{2}q_{0}^{2}}{R_{k+1}^{2}(I_{*}^{k+1})R_{k}^{2}(I_{*}^{k})I_{*}^{k}}\,(I_{k+1}-I_{*}^{k})^{2}.$
(48)

On the other hand,

$\displaystyle-\int_{I_{*}^{k}}^{I_{k+1}}\frac{f_{k+1}(I_{*}^{k+1})-f_{k}(I_{*}^{k})}{f_{k}(i)}\,di$
$\displaystyle=-\int_{I_{*}^{k}}^{I_{k+1}}\frac{f_{k+1}(I_{*}^{k+1})-f_{k+1}(I_{*}^{k})+f_{k+1}(I_{*}^{k})-f_{k}(I_{*}^{k})}{f_{k}(i)}di$
$\displaystyle=-\int_{I_{*}^{k}}^{I_{k+1}}\frac{I_{*}^{k+1}R_{k+1}(I_{*}^{k+1})-I_{*}^{k}R_{k+1}(I_{*}^{k})+I_{*}^{k}(R_{k+1}(I_{*}^{k})-R_{k}(I_{*}^{k}))}{f_{k}(i)}di$
$\displaystyle=\int_{I_{*}^{k}}^{I_{k+1}}\frac{I_{*}^{k+1}(R_{k+1}(I_{*}^{k})-R_{k+1}(I_{*}^{k+1}))}{f_{k}(i)}di$
$\displaystyle+\int_{I_{*}^{k}}^{I_{k+1}}\frac{(I_{*}^{k}-I_{*}^{k+1})R_{k+1}(I_{*}^{k})+I_{*}^{k}(R_{k}(I_{*}^{k})-R_{k+1}(I_{*}^{k}))}{f_{k}(i)}di$
$\displaystyle\leq\frac{I_{k+1}-I_{*}^{k}}{f_{k}(I_{*}^{k})}\bigl{(}(I_{*}^{k}-I_{*}^{k+1})R_{k+1}(I_{*}^{k})+I_{*}^{k}(R_{k}(I_{*}^{k})-R_{k+1}(I_{*}^{k}))\bigr{)},$
(49)

where we used that $R_{k+1}(I_{*}^{k})-R_{k+1}(I_{*}^{k+1})\leq 0$ because
$R_{k+1}(I)$ decreases and $f_{k}(I_{*}^{k})\leq f_{k}(I)$ on
$[I_{*}^{k},I_{k+1}]$ because $f_{k}(I)$ increases. Hence, relations (40) and
(41) imply

$\displaystyle-\int_{I_{*}^{k}}^{I_{k+1}}\frac{f_{k+1}(I_{*}^{k+1})-f_{k}(I_{*}^{k})}{f_{k}(i)}\,di\leq\frac{\rho
q_{0}R_{k+1}(I_{*}^{k})(I_{k+1}-I_{*}^{k})^{2}}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})f_{k}(I_{*}^{k})}$
$\displaystyle+\frac{q_{0}I_{*}^{k}(I_{k+1}-I_{*}^{k})^{2}}{f_{k}(I_{*}^{k})}\leq\frac{q_{0}(I_{k+1}-I_{*}^{k})^{2}}{R_{k}(I_{*}^{k})}\left(\frac{\rho}{f_{k}(I_{*}^{k})}+1\right).$

Combining this estimate with (46) and (48), we obtain (43). $\Box$

Now we combine Lemmas 4 and 6 to obtain the following statement.

###### Lemma 7

If $R_{r(t_{k})}(I(t))S(t)>1$ on $(t_{k},t_{k+1})$, then

$V_{k+1}(I_{k+1},S_{k+1})-V_{k}(I_{k},S_{k})\leq-Q_{k}(I_{k+1}-I_{*}^{k})^{2}$
(50)

with

$Q_{k}=\frac{\rho}{4}\sqrt{\frac{\epsilon_{0}S_{M}}{2f_{k}(I_{k+1})}}-\frac{q_{0}}{R_{k}(I_{*}^{k})}\left(\frac{\rho^{2}q_{0}}{R_{k+1}^{2}(I_{*}^{k+1})f_{k}(I_{*}^{k})}+\frac{\rho}{f_{k}(I_{*}^{k})}+1\right),$
(51)

where $M=(I_{M},S_{M})$ is the highest point of the level line
$V_{k}(I,S)=V_{k}(I_{k+1},S_{k+1})$ of the Lyapunov function $V_{k}$.

Proof. From Lemma 1, it follows that

$\displaystyle\int_{I_{*}^{k}}^{I_{k+1}}\left(1-\frac{f_{k}(I_{*}^{k})}{f_{k}(i)}\right)di=\int_{I_{*}^{k}}^{I_{k+1}}\frac{f_{k}(i)-f_{k}(I_{*}^{k})}{f_{k}(i)}di$
(52)
$\displaystyle\geq\int_{I_{*}^{k}}^{I_{k+1}}\frac{\epsilon_{0}(i-I_{*}^{k})}{f_{k}(I_{k+1})}di\geq\frac{\epsilon_{0}(I_{k+1}-I_{*}^{k})^{2}}{2f_{k}(I_{k+1})}.$

Substituting this estimate in (33) and adding (43), we obtain (50). $\Box$

As the next step, we obtain a counterpart of Lemma 7 for the parts of the
trajectory satisfying $R_{r(t_{k})}(I(t))S(t)<1$ on $(t_{k},t_{k+1})$.

###### Lemma 8

If $R_{r(t_{k})}(I(t))S(t)<1$ on $(t_{k},t_{k+1})$, then

$V_{k+1}(I_{k+1},S_{k+1})-V_{k}(I_{k},S_{k})\leq-\tilde{Q_{k}}(I_{k+1}-I_{*}^{k})^{2}$
(53)

with

$\tilde{Q_{k}}=\frac{\rho}{4}\sqrt{\frac{\epsilon_{0}S_{m}}{2f_{k}({I_{*}^{k}})}}-q_{0}\left(\frac{\rho^{2}q_{0}}{R_{k+1}^{2}(I_{*}^{k+1})R_{k}^{2}(I_{*}^{k})I_{*}^{k}}+\frac{\rho
R_{k+1}(I_{*}^{k})}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})f_{k}(I_{k+1})}+\frac{I_{*}^{k}}{f_{k}(I_{k+1})}\right),$
(54)

where $m=(I_{m},S_{m})$ is the lowest point of the level line
$V_{k}(I,S)=V_{k}(I_{k+1},S_{k+1})$ of the Lyapunov function $V_{k}$.

Figure 4: An arc of a trajectory of system (2), (13) satisfying
$R_{k}(I(t))S(t)\leq 1$.

Proof.  We adapt the proofs of Lemmas 4 – 6 to obtain their counterparts for
the part of the trajectory satisfying $R_{r(t_{k})}(I(t))S(t)<1$ on
$(t_{k},t_{k+1})$, which implies $\dot{I}<0$. The main steps are outlined
below for completeness.

The first coordinate of the point $m=(I_{m},S_{m})$ and the value of the
Lyapunov function $V_{k}$ at this point satisfy $I_{m}=I_{k}^{*}$ and

$V_{k}(I_{*}^{k},S_{m})=S_{m}-S_{*}^{k}-S_{*}^{k}\ln{\frac{S_{m}}{S_{*}^{k}}}=V_{k}(I_{k+1},S_{k+1}),$

see Figure 4. Using Lemma 3, we obtain the estimate

$I_{*}^{k}-I_{c}\geq
I_{*}^{k}-I_{k+1}-\frac{S_{k+1}-S_{c}}{S_{k+1}-S_{m}}(I_{*}^{k}-I_{k+1})=\frac{S_{c}-S_{m}}{S_{k+1}-S_{m}}(I_{*}^{k}-I_{k+1})$
(55)

for any point $(I_{c},S_{c})$ that belongs to the arc of the level line of
$V_{k}$ connecting the points $m$ and
$(I_{k+1},S_{k+1})=(I_{k+1},1/{R_{k}(I_{k+1})})$ (cf. (34)).

Using (30) in the domain $I\leq I_{*}^{k}$, $S\leq{1}/{R_{k}(I)}$, we obtain
that along the trajectory

$\frac{dV_{k}}{dI}\geq-\frac{{\rho}(S-S_{*}^{k})^{2}}{{SS_{*}^{k}}(R_{k}(I)S-1)I}\geq
0,$

where ${1}/{R_{k}(I)}\leq{1}/{R_{k}(I_{*}^{k})}=S_{*}^{k}$ implies
$0\leq{1}/{R_{k}(I)}-S\leq S_{*}^{k}-S$, hence

$\frac{dV_{k}}{dI}\geq\frac{{\rho}(S-S_{*}^{k})^{2}}{{SS_{*}^{k}}R_{k}(I)I(S_{*}^{k}-S)}=\frac{\rho}{SS_{*}^{k}R_{k}(I)I}(S_{*}^{k}-S)$

and due to $R_{k}(I)S<1$, $S_{*}^{k}\leq 1$,

$\frac{dV_{k}}{dI}\geq\frac{\rho(S_{*}^{k}-S)}{I}.$

Since the trajectory lies below the level line of $V_{k}$, on the interval
$[I_{c},I_{*}]$ we have $S\leq S_{c}$ and
$\frac{dV_{k}}{dI}\geq\frac{\rho}{I}(S_{*}^{k}-S_{c})$, which implies

$V_{k}(I_{*}^{k},\hat{S}(I_{*}^{k}))-V_{k}(I_{c},\hat{S}(I_{c}))\geq\rho(S_{*}^{k}-S_{c})\ln{\frac{I_{*}^{k}}{I_{c}}},$

where $(I,\hat{S}(I))$, $I\in[I_{k+1},I_{k}]$ is a parameterization of the
segment of the trajectory, and using the identity $\ln{(1+x)}\geq x/(1+x)$,
$x\geq-1$, we obtain

$\displaystyle
V_{k}(I_{c},\hat{S}(I_{c}))-V_{k}(I_{*}^{k},\hat{S}(I_{*}^{k}))\leq-\rho(S_{*}^{k}-S_{c})\frac{I_{*}^{k}-I_{c}}{I_{*}^{k}}$
$\displaystyle\leq-\rho\frac{(S_{*}^{k}-S_{c})(S_{c}-S_{m})(I_{*}^{k}-I_{k+1})}{S_{k+1}-S_{m}},$

where in the last step we used $I_{*}^{k}\leq 1$ and (55). Maximizing the
right hand side of this formula over $S_{c}\in[S_{m},S_{k+1}]$ as in the proof
of Lemma 4 and taking into account that the function $V_{k}$ decreases along
the trajectory, we arrive at (cf. (37))

$\Delta
V_{k}=V_{k}(I_{k+1},S_{k+1})-V_{k}(I_{k},S_{k})\leq-\frac{\rho}{4}(S_{*}^{k}-S_{m})(I_{*}^{k}-I_{k+1}).$
(56)

In order to estimate the difference $S_{*}^{k}-S_{m}$ in (56), we equate the
values of the function $V_{k}$ at the points $(I_{*}^{k},S_{m})$ and
$(I_{k+1},S_{k+1})$ to obtain

$S_{*}^{k}\ln{\frac{S_{k+1}}{S_{m}}}-(S_{k+1}-S_{m})=\int_{I_{k+1}}^{I_{*}^{k}}\left(\frac{f_{k}(I_{*}^{k})}{f_{k}(i)}-1\right)di.$

Using the inequality $2(1+x)\ln{(1+x)}\leq{x(2+x)}$, $x\geq 0$ we get

$\ln{\frac{S_{k+1}}{S_{m}}}\leq\frac{S_{k+1}^{2}-S_{m}^{2}}{2S_{k+1}S_{m}},$

where $S_{m}<S_{k+1}<S_{*}^{k}$, therefore

$\displaystyle
S_{*}^{k}\ln{\frac{S_{k+1}}{S_{m}}}-(S_{k+1}-S_{m})\leq(S_{k+1}-S_{m})\left(\frac{S_{*}^{k}(S_{k+1}+S_{m})}{2S_{k+1}S_{m}}-1\right)$
$\displaystyle=\frac{S_{k+1}-S_{m}}{2}\left(\frac{S_{*}^{k}-S_{m}}{S_{m}}+\frac{S_{*}^{k}-S_{k+1}}{S_{k+1}}\right)\leq\frac{(S_{*}^{k}-S_{m})^{2}}{S_{m}},$

hence

$\frac{(S_{*}^{k}-S_{m})^{2}}{S_{m}}\geq\int_{I_{k+1}}^{I_{*}^{k}}\left(\frac{f_{k}(I_{*}^{k})}{f_{k}(i)}-1\right)di$

and from (56) it follows that (cf. (33))

$\Delta
V_{k}\leq-\frac{\rho}{4}(I_{*}^{k}-I_{k+1})\sqrt{S_{m}\int_{I_{k+1}}^{I_{*}^{k}}\left(\frac{f_{k}(I_{*}^{k})}{f_{k}(i)}-1\right)di}.$

Further, due to Lemma 1,

$\int_{I_{*}^{k}}^{I_{k+1}}\left(1-\frac{f_{k}(I_{*}^{k})}{f_{k}(i)}\right)di\geq\frac{\epsilon_{0}(I_{k+1}-I_{*}^{k})^{2}}{2f_{k}(I_{*}^{k})},$

consequently,

$\Delta
V_{k}\leq-\frac{\rho}{4}(I_{*}^{k}-I_{k+1})^{2}\sqrt{\frac{\epsilon_{0}S_{m}}{2f_{k}(I_{*}^{k})}}.$
(57)

As in the proof of Lemma 5,

$S_{*}^{k}-S_{*}^{k+1}=\frac{I_{*}^{k+1}-I_{*}^{k}}{\rho}=\frac{1}{R_{k}(I_{*}^{k})R_{k+1}(I_{*}^{k+1})}\bigl{(}R_{k+1}^{\prime}(\hat{I})(I_{*}^{k+1}-I_{*}^{k})+R_{k+1}(I_{*}^{k})-R_{k}(I_{*}^{k})\bigr{)}$

but this time (cf. (41))

$0\leq R_{k}(I_{*}^{k+1})-R_{k}(I_{*}^{k})\leq q_{0}(I_{*}^{k}-I_{k+1}),$ (58)

which allows us to conclude (due to $R_{k+1}^{\prime}(\hat{I})\leq 0$) that
$I_{*}^{k+1}\geq I_{*}^{k}$ and $S_{*}^{k+1}\leq S_{*}^{k}$, and the argument
used in the proof of Lemma 5 leads to the estimates (cf. (40))

$0\leq I_{*}^{k+1}-I_{*}^{k}=\rho(S_{*}^{k}-S_{*}^{k+1})\leq\frac{\rho
q_{0}}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})}\,(I_{*}^{k}-I_{k+1}).$ (59)

In particular, as opposed to the case considered in Lemmas 4 – 6,

$I_{*}^{k+1}\geq I_{*}^{k}\geq I_{k+1},\quad S_{k+1}\leq S_{*}^{k+1}\leq
S_{*}^{k},$ (60)

while

$R_{k}(I)\leq R_{k+1}(I)\quad\text{for}\quad I\geq I_{k+1}.$ (61)

Using the notation $\Delta_{S}$ and $\Delta_{I}$ introduced in the proof of
Lemma 6, from (60) we obtain (44) and $\Delta_{S}\leq 0$, while (61) leads to
(45), hence the difference (39) satisfies (46). Further, relation (47) also
holds true, hence (59) implies (48). Finally, a counterpart of (49) reads

$\displaystyle-\int_{I_{*}^{k}}^{I_{k+1}}\frac{f_{k+1}(I_{*}^{k+1})-f_{k}(I_{*}^{k})}{f_{k}(i)}\,di$
$\displaystyle\leq\frac{I_{*}^{k}-I_{k+1}}{f_{k}(I_{k+1})}\bigl{(}(I_{*}^{k+1}-I_{*}^{k})R_{k+1}(I_{*}^{k})+I_{*}^{k}(R_{k+1}(I_{*}^{k})-R_{k}(I_{*}^{k}))\bigr{)}$

and due to (58) and (59),

$-\int_{I_{*}^{k}}^{I_{k+1}}\frac{f_{k+1}(I_{*}^{k+1})-f_{k}(I_{*}^{k})}{f_{k}(i)}\,di\leq\frac{q_{0}(I_{*}^{k}-I_{k+1})^{2}}{f_{k}(I_{k+1})}\left(\frac{\rho
R_{k+1}(I_{*}^{k})}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})}+I_{*}^{k}\right).$

Combining this estimate with (46) and (48), we obtain (cf. (43))

$\tilde{\Delta}V_{k}\leq
q_{0}(I_{*}^{k}-I_{k+1})^{2}\left(\frac{\rho^{2}q_{0}}{R_{k+1}^{2}(I_{*}^{k+1})R_{k}^{2}(I_{*}^{k})I_{*}^{k}}+\frac{\rho
R_{k+1}(I_{*}^{k})}{R_{k+1}(I_{*}^{k+1})R_{k}(I_{*}^{k})f_{k}(I_{k+1})}+\frac{I_{*}^{k}}{f_{k}(I_{k+1})}\right)$

and adding (57) leads to (53). $\Box$

In order to use Lemmas 7 and 8 we need to estimate the components $I(t),S(t)$
of the trajectory from below.

###### Lemma 9

For $t\geq t_{2}$,

$I(t)\geq
i_{0}:=\rho\,\frac{R_{0}^{int}-1}{R_{0}^{int}}\,e^{-\left(\frac{2}{\rho(R_{0}^{int}-1)}+\frac{1}{R_{0}^{int}}\right)R_{0}^{nat}},\qquad
S(t)\geq s_{0}:=\frac{e^{-(1+2R_{0}^{nat})}}{R_{0}^{nat}}.$ (62)

Proof.  Consider a point $(I_{k},S_{k})=(I_{k},{1}/{R_{k}(I_{k})})$ where the
trajectory satisfies $\dot{I}=0,$ $\dot{S}<0$, see Figure 4. At this point,

$V_{k}(I_{k},S_{k})=\int_{I_{*}^{k}}^{I_{k}}\left(1-\frac{f_{k}(I_{*}^{k})}{f_{k}(i)}\right)di+\int_{S_{*}^{k}}^{S_{k}}\left(1-\frac{S_{*}^{k}}{s}\right)ds\leq
S_{k}-S_{*}^{k}+I_{k}-I_{*}^{k}\leq 2.$

Clearly, the minimal value $s_{m}^{k}=\min_{t\in[t_{k},t_{k+2}]}S(t)$ of the
$S$-component of the trajectory over the time interval $[t_{k},t_{k+2}]$ is
achieved at a time $t_{m}^{k}\in[t_{k},t_{k+1}]$, i.e.
$s_{m}^{k}=S(t_{m}^{k})$. Since $V_{k}$ decreases on the time interval
$[t_{k},t_{k+1}]$, from $V_{k}(I_{k},S_{k})\leq 2$, it follows that

$2\geq
V_{k}(I(t_{m}^{k}),S(t_{m}^{k}))\geq\int_{s_{m}^{k}}^{S_{*}^{k}}\left(\frac{S_{*}^{k}}{s}-1\right)ds=S_{*}^{k}\ln{\frac{S_{*}^{k}}{s_{m}^{k}}}-(S_{*}^{k}-s_{m}^{k})\geq
S_{*}^{k}\ln{\frac{S_{*}^{k}}{s_{m}^{k}}}-S_{*}^{k}.$

Therefore, $s_{m}^{k}\geq S_{*}^{k}e^{-\bigl{(}1+\frac{2}{S_{*}^{k}}\bigr{)}}$
and due to $S_{*}^{k}\geq 1/R_{0}^{nat}$,

$s_{m}^{k}=\min_{t\in[t_{k},t_{k+2}]}S(t)\geq\frac{e^{-(1+2R_{0}^{nat})}}{R_{0}^{nat}}.$

On the other hand, the minimal value
$i_{m}^{k}=\min_{t\in[t_{k},t_{k+2}]}I(t)$ of the $I$-component of the
trajectory over the time interval $[t_{k},t_{k+2}]$ is achieved at the point
$(I_{k+1},S_{k+1})$, i.e. $i_{m}^{k}=I_{k+1}$. Again, since $V_{k}$ decreases
on the time interval $[t_{k},t_{k+1}]$, from $V_{k}(I_{k},S_{k})\leq 2$, it
follows that

$\displaystyle 2\geq V_{k}(I_{k+1},S_{k+1})$ $\displaystyle\geq$
$\displaystyle\int_{I_{k+1}}^{I_{*}^{k}}\left(\frac{f_{k}(I_{*}^{k})}{f_{k}(i)}-1\right)di=f_{k}(I_{*}^{k})\int_{I_{k+1}}^{I_{*}^{k}}\frac{di}{iR_{k}(i)}-(I_{*}^{k}-I_{k+1})$
$\displaystyle\geq$
$\displaystyle\frac{f_{k}(I_{*}^{k})}{R_{0}^{nat}}\int_{I_{k+1}}^{I_{*}^{k}}\frac{di}{i}-I_{*}^{k}=\frac{f_{k}(I_{*}^{k})}{R_{0}^{nat}}\ln{\frac{I_{*}^{k}}{I_{k+1}}}-I_{*}^{k},$

hence

$I_{k+1}\geq
I_{*}^{k}e^{-(2+I_{*}^{k})\frac{R_{0}^{nat}}{f_{k}(I_{*}^{k})}}=I_{*}^{k}e^{-\left(\frac{2}{I_{*}^{k}R_{k}(I_{*}^{k})}+\frac{1}{R_{k}(I_{*}^{k})}\right)R_{0}^{nat}}.$

Therefore, using $I_{*}^{k}=\rho(1-S_{*}^{k})$ and
$S_{*}^{k}=1/R_{k}(I_{*}^{k})\leq 1/R_{0}^{int}$, we obtain

$I_{k+1}=\min_{t\in[t_{k},t_{k+2}]}I(t)\geq\rho(1-S_{*}^{k})e^{-\left(\frac{2}{\rho(1-S_{*}^{k})}+1\right)S_{*}^{k}R_{0}^{nat}}\geq\rho\,\frac{R_{0}^{int}-1}{R_{0}^{int}}\,e^{-\left(\frac{2}{\rho(R_{0}^{int}-1)}+\frac{1}{R_{0}^{int}}\right)R_{0}^{nat}},$

which completes the proof. $\Box$

We can now prove the convergence of the trajectories to the equilibrium set.

###### Lemma 10

Every trajectory of system (2), (13) with initial values $I(0)>0,S(0)>0$
converges to the set of endemic equilibrium states.

Proof. From Lemmas 7 – 9, it follows that the sequence
$v_{k}=V_{k}(I_{k},S_{k})$ satisfies

$v_{k+1}-v_{k}\leq-\kappa(I_{*}^{k}-I_{k+1})^{2},\qquad k=0,1,2,\ldots$ (63)

with

$\kappa:=\frac{\rho}{4}\sqrt{\frac{\epsilon_{0}s_{0}}{2R_{0}^{nat}}}-\frac{q_{0}}{i_{0}R_{0}^{int}}\left(\frac{\rho^{2}q_{0}}{(R_{0}^{int})^{3}}+\frac{\rho
R_{0}^{nat}}{(R_{0}^{int})^{2}}+1\right)>0,$ (64)

independent of $k$.

Since $f(I)=R(I)I\leq R_{0}^{nat}$ for $0\leq I\leq 1$, Lemma 1 implies that

$\int_{I_{*}^{k}}^{I}\frac{f_{k}(i)-f_{k}(I_{*}^{k})}{f_{k}(i)}di\geq\frac{\epsilon_{0}}{R_{0}^{nat}}\int_{I_{*}^{k}}^{I}(i-I_{*}^{k})di=\frac{\epsilon_{0}}{2R_{0}^{nat}}(I-I_{*}^{k})^{2}.$

On the other hand,

$f^{\prime}_{k}(I)=(IR_{k}(I))^{\prime}=R^{\prime}_{k}(I)I+R_{k}(I)\leq
R_{k}(I)\leq R^{nat}_{0},\quad f_{k}(I)=IR_{k}(I)\geq i_{0}R_{0}^{int},$

where $i_{0}$ is defined in (62), hence

$\int_{I_{*}^{k}}^{I}\frac{f_{k}(i)-f_{k}(I_{*}^{k})}{f_{k}(i)}di\leq\frac{R_{0}^{nat}}{i_{0}R_{0}^{int}}\int_{I_{*}^{k}}^{I}(i-I_{*}^{k})di=\frac{R_{0}^{nat}}{2i_{0}R_{0}^{int}}(I-I_{*}^{k})^{2},$

so we have

$\frac{\epsilon_{0}}{2R_{0}^{nat}}(I-I_{*}^{k})^{2}\leq\int_{I_{*}^{k}}^{I}\frac{f_{k}(i)-f_{k}(I_{*}^{k})}{f_{k}(i)}di\leq\frac{R_{0}^{nat}}{2i_{0}R_{0}^{int}}(I-I_{*}^{k})^{2}.$
(65)

Similarly, due to $0\leq S\leq 1$,

$\frac{(S-S_{*}^{k})^{2}}{2}\leq\int_{S_{*}^{k}}^{S}\left(1-\frac{S_{*}^{k}}{s}\right)ds\leq\frac{(S-S_{*}^{k})^{2}}{2s_{0}},$
(66)

with $s_{0}$ defined in (62). Adding (65) and (66), we obtain

$\frac{\epsilon_{0}(I-I_{*}^{k})^{2}}{2R_{0}^{nat}}+\frac{(S-S_{*}^{k})^{2}}{2}\leq
V_{k}(I,S)\leq\frac{R_{0}^{nat}(I-I_{*}^{k})^{2}}{2i_{0}R_{0}^{int}}+\frac{(S-S_{*}^{k})^{2}}{2s_{0}}.$
(67)

Relation (63) implies

$v_{n}=v_{1}+\sum_{k=1}^{n-1}(v_{k+1}-v_{k})\leq
v_{1}-\kappa\sum_{k=1}^{n-1}(I_{*}^{k}-I_{k+1})^{2},$

and due to the non-negativity of each Lyapunov function $V_{n}$,

$\sum_{k=1}^{n-1}(I_{*}^{k}-I_{k+1})^{2}\leq\frac{1}{\kappa}v_{1}$

for any $n$, hence

$\sum_{k=1}^{\infty}(I_{*}^{k}-I_{k+1})^{2}<\infty.$

Therefore $I_{k+1}-I_{*}^{k}\xrightarrow[]{}0$ (as $k\xrightarrow[]{}\infty$).
Due to (4) and (59), this implies $I_{k}-I_{*}^{k}\xrightarrow[]{}0$, and from
$S_{*}^{k}={1}/{R_{k}(I_{*}^{k})}$, $S_{k}={1}/{R_{k}(I_{k})}$ it follows that
$S_{k}-S_{*}^{k}\xrightarrow[]{}0$. Hence, (67) implies
$v_{k}=V_{k}(I_{k},S_{k})\to 0$. Since $V_{k}$ is nonnegative and decreases
along the trajectory, we conclude that

$V_{k}(I_{k},S_{k})=\max_{t\in[t_{k},t_{k+1}]}V_{k}(I(t),S(t))\xrightarrow[]{}0\quad\text{as}\quad
k\xrightarrow[]{}\infty,$ (68)

hence (67) implies that the trajectory converges to the equilibrium set.
$\Box$

It remains to show that every trajectory converges to an equilibrium. The
following lemma completes the proof of the theorem.

###### Lemma 11

Every trajectory of system (2), (13) with positive initial values $I(0),S(0)$
converges to an endemic equilibrium state.

Proof. From (67), it follows that

$V_{k}(I_{k+1},S_{k+1})\leq\frac{R_{0}^{nat}}{2i_{0}R_{0}^{int}}(I_{k+1}-I_{*}^{k})^{2}+\frac{1}{2s_{0}}\left(\frac{1}{R_{k}(I_{k+1})}-\frac{1}{R_{k}(I_{*}^{k})}\right)^{2},$

where $R_{k}(I)\geq R_{0}^{int}$ and $|R_{k}(I_{k+1})-R_{k}(I_{*}^{k})|\leq
q_{0}|I_{k+1}-I_{*}^{k}|$ due to (22), hence

$V_{k}(I_{k+1},S_{k+1})\leq a(I_{k+1}-I_{*}^{k})^{2},\qquad
a:=\frac{R_{0}^{nat}}{2i_{0}R_{0}^{int}}+\frac{q_{0}}{2s_{0}(R_{0}^{int})^{2}}.$
(69)

Lemma 6 implies that

$\tilde{\Delta}V_{k}\leq b(I_{k+1}-I_{*}^{k})^{2},\qquad
b:=\frac{q_{0}}{R_{0}^{int}}\left(\frac{\rho^{2}q_{0}}{i_{0}R_{0}^{int}(R_{0}^{int})^{2}}+\frac{\rho}{i_{0}R_{0}^{int}}+1\right),$

hence, using (69),

$V_{k+1}(I_{k+1},S_{k+1})=V_{k}(I_{k+1},S_{k+1})+\tilde{\Delta}V_{k}\leq(a+b)(I_{k+1}-I_{*}^{k})^{2}.$
(70)

If $V_{k}(I_{k},S_{k})\geq(a+b+\kappa)(I_{k+1}-I_{*}^{k})^{2}$ with $\kappa$
defined by (64), then (70) implies

$v_{k+1}=V_{k+1}(I_{k+1},S_{k+1})\leq p\,V_{k}(I_{k},S_{k})=pv_{k},\qquad
p:=\frac{a+b}{a+b+\kappa}<1.$ (71)

On the other hand, if
$V_{k}(I_{k},S_{k})\leq(a+b+\kappa)(I_{k+1}-I_{*}^{k})^{2},$ then (63) implies
(71). Therefore, in either case,

$v_{k}\leq v_{0}p^{k}.$ (72)

Since $V_{k}$ decreases along the trajectory, $V_{k}(I_{k+1},S_{k+1})\leq
V_{k}(I_{k},S_{k})=v_{k}$. But according to (67),

$V_{k}(I_{k+1},S_{k+1})\geq\frac{\epsilon_{0}}{2R_{0}^{nat}}(I_{k+1}-I_{*}^{k})^{2},$
(73)

hence

$\frac{\epsilon_{0}}{2R_{0}^{nat}}(I_{k+1}-I_{*}^{k})^{2}\leq v_{k}\leq
v_{0}p^{k}.$

Further, according to (40), (59),

$|I_{*}^{k+1}-I_{*}^{k}|\leq\frac{\rho
q_{0}}{(R_{0}^{int})^{2}}|I_{k+1}-I_{*}^{k}|,$

consequently

$|I_{*}^{k+1}-I_{*}^{k}|\leq\frac{\rho
q_{0}}{(R_{0}^{int})^{2}}\sqrt{\frac{2v_{0}R_{0}^{nat}}{\epsilon_{0}}}\,p^{\frac{k}{2}}.$

We see that $I_{*}^{k}$ is a Cauchy sequence, hence it converges to a limit
value $I_{*}$. Since $I_{k+1}-I_{*}^{k}\xrightarrow[]{}0$ as
$k\xrightarrow[]{}\infty$, we conclude that $I(t)\to I_{*}$ as $t\to\infty$,
which implies that each of the functions $R_{0}(t)$ and $S(t)$ also has a
limit as $t\to\infty$. Denoting the limit of $S(t)$ by $S_{*}$, the trajectory
converges to the equilibrium $(I_{*},S_{*})$. $\Box$

## References

  * [1] Balanov Z, Krawcewicz W, Rachinskii D, Zhezherun A (2012) Hopf bifurcation in symmetric networks of coupled oscillators with hysteresis. Journal of Dynamics and differential Equations 24(4):713-759
  * [2] Brokate M, Pokrovskii A, Rachinskii D (2006) Asymptotic stability of continuum sets of periodic solutions to systems with hysteresis. Journal of Mathematical Analysis and Applications 319(1):94-109
  * [3] Brokate M, Sprekels J (1996) Hysteresis and phase transitions. Springer, New York
  * [4] Capasso V, Serio G (1978) A generalization of the Kermack-Mckendrick deterministic epidemic model. Math Biosci 42:43-61. https://doi.org/10.1016/0025-5564(78)90006-8
  * [5] Center of Disease Control and Prevention (CDC), USA (2021) Implementation of mitigation strategies for communities with local COVID-19 transmission. https://www.cdc.gov/coronavirus/2019-ncov/community/community-mitigation.html
  * [6] Chen FH (2004) Rational behavioral response and the transmission of STDs. Theor Popul Biol 66:307-316. doi: 10.1016/j.tpb.2004.07.004
  * [7] Chen X, Fu F(2019) Imperfect vaccine and hysteresis. Proc R Soc B 286(1894):20182406. http://doi.org/10.1098/rspb.2018.2406
  * [8] Chernorutskii V, Rachinskii D (1997) On uniqueness of an initial-value problem for ODE with hysteresis. Nonlinear Differential Equations and Applications NoDEA 4 (3), 391-399
  * [9] Chladná Z, Kopfová J, Rachinskii D, Rouf SC (2020) Global dynamics of SIR model with switched transmission rate. J Math Biol 80:1209-1233.
  * [10] Debenedetto P, Ruiz B (2021) Houston region passes threshold for tougher COVID-19 restrictions. Houston Public Media.
https://www.houstonpublicmedia.org/articles/news/health-
science/coronavirus/2021/01/05/388740/houston-area-projected-to-pass-
governors-threshold-for-tougher-covid-19-restrictions-hidalgo-says/

  * [11] Department of Health and Human Services, Nebraska (2020) Phased public health restrictions tied to coronavirus hospitalization rate. https://dhhs.ne.gov/Documents/DHM-Measure-Table-ENGLISH.pdf
  * [12] Di Bernardo M, Crowther-Swanepoel D, Broderick P et al (2008) A genome-wide association study identifies six susceptibility loci for chronic lymphocytic leukemia. Nat Genet 40:1204–1210. doi: 10.1038/ng.219
  * [13] d’Onofrio A, Manfredi P (2009) Information-related changes in contact patterns may trigger oscillations in the endemic prevalence of infectious diseases. J Theor Biol 256:473-478. doi: 10.1016/j.jtbi.2008.10.005
  * [14] d’Onofrio A, Manfredi P (2020) The interplay between voluntary vaccination and reduction of risky behavior: A general behavior-implicit SIR model for vaccine preventable infections. In: Maira Aguiar et al (eds) Current Trends in Dynamical Systems in Biology and NaturalSciences. Springer, Switzerland, pp 185-203. doi:10.1007/978-3-030-41120-6k10
  * [15] Ewing JA (1885). Experimental research in magnetism. Philos Trans R Soc Lond 176: 523-640
  * [16] Fairlie R (2020) The impact of COVID-19 on small business owners: Evidence from the first three months after widespread social distancing restrictions. J Econ Manag Strategy 29:727-740. doi: 10.1111/jems.12400
  * [17] Fathian K, Rachinskii DI, Spong MW, Summers TH, Gans NR (2019) Distributed formation control via mixed barycentric coordinate and distance-based approach. American Control Conference (ACC) 51-58
  * [18] Friedman G, McCarthy S, Rachinskii D (2014) Hysteresis can grant fitness in stochastically varying environment. PLOS ONE 9(7):e103241
  * [19] Geoffard PY, Philipson T (1996) Rational epidemics and their public control. Int Econ Rev 37:603-624. doi.org/10.2307/2527443
  * [20] Gostin LO, Wiley LF (2020) Governmental public health powers during the COVID-19 pandemic: Stay-at-home orders, business closures, and travel restrictions. JAMA 323:2137-2138. doi: 10.1001/jama.2020.5460
  * [21] Guidry JPD, Laestadius LI, Vraga EK , Miller CA, Perrin PB, Burton CW, Ryan M, Fuemmeler BF, Carlyle KE (2021) Willingness to get the COVID-19 vaccine with and without emergency use authorization. Am J Infect Control 49:137-142. doi: 10.1016/j.ajic.2020.11.018
  * [22] Hadeler KP, Castillo-Chavez C (1995) A core group model for disease transmission. Math Biosci 128:41-55. doi: 10.1016/0025-5564(94)00066-9
  * [23] Hooton E, Balanov Z, Krawcewicz W, Rachinskii D (2017) Noninvasive stabilization of periodic orbits in symmetrically coupled systems near a Hopf bifurcation point. International Journal of Bifurcation and Chaos 27(06):1750087
  * [24] Hazewinkel M (1997) Kluwer encyclopedia of mathematics: Supplement 1. Kluwer Acad Publ, pp 310-384
  * [25] Hoppensteadt FC, Jäger W (1980) Pattern formation by bacteria. In: Jäger W, Rost H, Tautu P (eds) Biological growth and spread: Lecture notes in biomathematics, vol 38. Springer, Berlin, pp 68-81
  * [26] Kermack WO, McKendrick AG (1927) A contribution to the mathematical theory of epidemics. Proc R Soc Lond A 115:700–721. doi: 10.1098/rspa.1927.0118
  * [27] Kopfová J, Nábelková P, Rachinskii D, Rouf SC (2021) Dynamics of SIR model with vaccination and heterogeneous behavioral response of individuals modeled by the Preisach operator. J Math Biology, 83, 11. doi: 10.1007/s00285-021-01629-8
  * [28] Krasnosel’skii MA, Pokrovskii AV (1989) Systems with Hysteresis. Springer, Berlin, Heidelberg
  * [29] Krasnosel’skii A, Rachinskii D (2002) On a bifurcation governed by hysteresis nonlinearity. Nonlinear Differential Equations and Applications NoDEA, 9(1):93-115
  * [30] Krasnosel’skii AM, Rachinskii DI (2002) On the existence of cycles in autonomous systems. Doklady Mathematics 65(3):344-349
  * [31] Krejčí P (1996) Hysteresis, convexity and dissipation in hyperbolic equations. Gakkotosho, Tokyo
  * [32] Krejčí P, Kaltenbacher B (2019) Analysis of optimization problem for a piezoelectric energy harvester. Arch Appl Mech 89:1103-1122
  * [33] Lazarus JV, Ratzan SC, Palayew A et al (2021) A global survey of potential acceptance of a COVID-19 vaccine. Nat Med 27:225. doi: 10.1038/s41591-020-1124-9
  * [34] Leonov G, Shumafov M, Teshev V, Aleksandrov K (2017) Differential equations with hysteresis operators: Existence of solutions, stability, and oscillations. Differ Equ 53:1764-1816.
  * [35] Marquioni VM, De Aguiar M (2020) Quantifying the effects of quarantine using an IBM SEIR model on scalefree networks. Chaos Solitons Fractals 138:109999. doi:10.1016/j.chaos.2020.109999
  * [36] Matrajt L, Leung T (2020) Evaluating the effectiveness of social distancing interventions to delay or flatten the epidemic curve of coronavirus disease. Emerg Infect Dis 26:1740-1748. doi: 10.3201/eid2608.201093
  * [37] Mayergoyz ID (2003) Mathematical models of hysteresis and their applications. Academic Press, New York
  * [38] Pimenov A, Kelly TC, Korobeinikov A, O’Callaghan MJA, Pokrovskii A, Rachinskii D (2012) Memory effects in population dynamics: Spread of infectious disease as a case study. Math Model Nat Phenom 7:1-30
  * [39] Rachinskii D (2015) Realization of arbitrary hysteresis by a low-dimensional gradient flow. Discrete and Continuous Dynamical Systems B 21(1):227-243. Arxiv preprint arXiv:1506.03842
  * [40] Rachinskii D, Ruderman M (2016) Convergence of direct recursive algorithm for identification of Preisach hysteresis model with stochastic input. SIAM Journal on Applied Mathematics 76(4):1270-1295
  * [41] Reluga TC, Bauch CT, Galvani AP (2006) Evolving public perceptions and stability in vaccine uptake. Math Biosci 204:185-98. doi: 10.1016/j.mbs.2006.08.015
  * [42] Rouf S, Rachinskii D (2021) Dynamics of SIR model with heterogeneous response to intervention policy. arXiv:2103.07989
  * [43] Su Z, Wang W, Li L, Xiao J, Stanley HE (2017) Emergence of hysteresis loop in social contagions on complex networks. Sci Rep 7, 6103. https://doi.org/10.1038/s41598-017-06286-w
  * [44] Visintin A (1994) Hysteresis and semigroups: In differential models of hysteresis. Springer, Berlin, Heidelberg, pp 211-256
  * [45] Wang W, Tang M, Stanley HE, Braunstein LA (2017) Unification of theoretical approaches for epidemic spreading on complex networks. Rep Prog Phys 80:036603. https://doi.org/10.1088/1361-6633/aa5398
  * [46] Wilder-Smith A, Chiew CJ, Lee VJ (2020) Can we contain the COVID-19 outbreak with the same measures as for SARS? Lancet Infect Dis 20:102-107. doi: 10.1016/S1473-3099(20)30129-8
  * [47] Yakubovich VA, Leonov GA, Gelig AK (2004) Stability of stationary sets in control systems with discontinuous nonlinearities. World Scientific, Singapore