key: cord-0623262-8h0gnflv
authors: Guan, Ruofei; Kopfov'a, Jana; Rachinskii, Dmitrii
title: Global stability of SIR model with heterogeneous transmission rate modeled by the Preisach operator
date: 2022-01-15
journal: nan
DOI: nan
sha: 0a3e229c91c31c5f74b7ef593cd4c0f5498e6ef3
doc_id: 623262
cord_uid: 8h0gnflv

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.

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 .

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.1. The systems considered below adapt the standard dimensionless SIR model

with simple immigration demography, where I and S are the densities of the infected and susceptible populations, respectively; β is the transmission coefficient; γ is the recovery rate; and µ, 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 ≥ 0, I + S ≤ 1 is positively invariant for system (1).

Re-scaling the time, one obtains the normalized systeṁ

with the dimensionless parameters

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

is globally stable in the open positive quadrant. If

then the endemic equilibrium is of focus type. and R 0 = R int 0 , was considered in Chladná et al., 2020; Rouf and Rachinskii, 2021. Following this work, we assume that the parameter ρ 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 = α 1 and I = α 2 , satisfying 0 < α 1 < α 2 < 1. More precisely, it is postulated that the basic reproduction number instantaneously switches from the value R nat 0 to the value R int 0 when the variable I = I(t) reaches the upper threshold value α 2 . On the other hand, R 0 switches back from the value R int 0 to the value R nat 0 as I(t) reaches the lower threshold value α 1 , see Figure 1 .

This two-threshold switching rule can be formalized using the standard nonideal 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 α 1 , α 2 , and we use the notation α = (α 1 , α 2 ). The function I(t) : R + → R is called the input of the relay. The state of the relay, denoted r α (t), equals either 0 or 1 at any moment t ∈ R + . Specifically, given any continuous input I(t) : R + → R and an initial value of the state, r α (0) = r 0 α , which satisfies the constraints

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

Function (7) will be denoted by

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) : R + → R rather than an instantaneous value of I as an input. Moreover, the function r α (t) : R + → {0, 1} depends both on the input I(t) : R + → R and the initial state r α (0) = r 0 α of the relay. By definition (7) of the input-to-state map (8) , the state satisfies the constraints r α (t) = 1 whenever I(t) ≥ α 2 ; r α (t) = 0 whenever I(t) ≤ α 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 ≤ t ≤ t 1 .

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

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 α 2 and is revoked when the density of infection drops below the lower threshold value α 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 ≥ 0, I + S ≤ 1 converges to either an equilibrium state or a periodic orbit. In particular, depending on the positioning of the thresholds α 1 , α 2 , the system can have either one attractor (an equilibrium state or a periodic orbit), or two attractors (two coexisting 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 α 1 and α 2 , the point (I(t), R 0 (t)) moves clockwise along the rectangular hysteresis loop shown in Figure 1 . (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 α, 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 α parameterized by points α of a subset Π ⊂ {α = (α 1 , α 2 ) : α 1 < α 2 } of the α-plane, assuming a homogeneous response within each group S α . 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 α is given by (10) and the contact structure is separable. Then, the average basic reproduction number for the entire population at time t equals

where the cumulative distribution function F (α) describes the distribution of the susceptible population S over the index set Π (the set of threshold pairs α = (α 1 , α 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 (α) is independent of time, i.e. F (α) does not change with variations of I. In this case, the mapping of the space of continuous inputs I(t) : R + → R to the space of outputs R 0 (t) : (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, . . . , N ) in place, each decreasing the transmission coefficient by a certain amount ∆β 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 α n = (α n 1 , α n 2 ). The response of the population is assumed to be homogeneous. Under these assumptions, the basic reproduction number of system (2) is given by

whereq n = ∆β n /(γ + µ); hence, R 0 is set to change at multiple thresholds α n 1 , α n 2 . 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 (α) (see Appendix 5.1 for details). The corresponding operator (11) , which is called the continuous Preisach model, can be written in the equivalent form

where 0 < R int 0 < R nat 0 and q = q(α) : Π → R + is a strictly positive measurable and bounded probability density function, i.e.

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

and

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(τ ) achieved on the time interval τ ∈ [0, t]. Fundamental properties of the history-dependence, which are typical of hysteresis operators, are briefly discussed in Appendix 5.1. 

3.1. Trajectories of system (2) coupled with the operator relationship (13) lie in the infinite-dimensional phase space U of triplets (I, S, r α ), where the measurable function r α = r(α) : Π → {0, 1} of the variable α = (α 1 , α 2 ) describes the states of the relays R α 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 α 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 α = r(α) : Π → {0, 1} by the equation

Due to assumption (15) and the compatibility constraint (9), the inclusion α ∈ Π implies that all the relays are in state r α = 0 when I = 0. Therefore, system (2), (13) has a unique infection free equilibrium state

in which the function r α : Π → {0, 1} is the identical zero and the basic reproduction number equals R 0 = R nat 0 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

which are characterized by different proportions of the infected, susceptible and recovered populations and different values R θ 0 of the basic reproduction number (Rachinskii and Rouf 2021).

The phase space U of system (2), (13) is naturally embedded into the space R 2 × L 1 (Π; 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

Theorem 1 For any given ρ > 0, R nat 0 > 1 and a strictly positive probability density function q = q(α) : Π → R + of the Preisach operator, there is an ε ∈ (0, R nat 0 − 1) such that if 1 + ε ≤ R int 0 < R nat 0 , then the set of endemic equilibrium states of system (2), (13) is globally asymptotically stable in U 0 . Moreover, each trajectory from U 0 converges to an endemic equilibrium state. In particular, the components I(t), S(t) of any trajectory from U 0 satisfy (

The proof presented in the Appendix provides an explicit estimate for the interval of values of R int 0 ∈ [1 + ε, R nat 0 ), for which the set of endemic equilibrium states is globally stable. In particular, the global stability is guaranteed if the quantities 0 and κ defined by (27) , (64) using (23), (62) are both positive.

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 α 1 and ends the intervention when the number of active cases drops below another threshold α 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, σ 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 σ 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 σ 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. It involves a collection of non-ideal relays R α , which respond to the same continuous input I(t) : R + → R independently according to formula (7) . The relays contributing to the system have different pairs of thresholds α = (α 1 , α 2 ) ∈ Π, where the subset Π of the half-plane {α = (α 1 , α 2 ) : α 1 < α 2 } is assumed to be measurable and bounded; the α-plane is called the Preisach plane. The output of the continuous Preisach model is the scalar-valued function R 0 (t) : R + → R defined by (13) , where q(α) : Π → R + is a positive bounded measurable function (probability density function) representing the weights of the relays; and, r 0 α is the initial state of the relay R α for any given α ∈ Π. The function r 0 α = r 0 (α) : Π → {0, 1} of the variable α = (α 1 , α 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 ≥ 0 and, furthermore, the output R 0 (t) : R + → R of the Preisach model is a continuous function of time.

The function (7) with a fixed t ≥ 0 and varying α ∈ Π 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 α (t) = r(t; α) : Π → {0, 1} is an element of the space L 1 (Π; R) for each t ≥ 0. In this sense, formula (7) defines the evolution of the state function in the space L 1 (Π; R) for a given input provided that the initial state-input pair is compatible. Moreover, the state function r α (t) = r(t; α) : Π → {0, 1} and the concurrent input value I(t) are compatible for each t ≥ 0 and the semigroup property holds:

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

where both the input I(t) : R + → R and the initial state function r 0 α = r 0 (α) : Π → {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) : R + → R. An important property of the Preisach operator (13) is that it is Lipschitz continuous if q(α) : Π → R + is bounded (Krasnosel'skii et al. 1989 ). More precisely, the relations

and 0 ≤ I 1 (t), I 2 (t) ≤ 1 (t ≥ 0) imply

for any τ ≥ 0 with the Lipschitz constant The set 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 α ) ∈ 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 α of the Preisach opeartor (see e.g. the survey Leonov et al. 2017 ). This solution satisfies (I(t), S(t), r α (t)) ∈ U for each t in its domain. Further, the positive invariance of D implies that each solution is extendable to the whole semi-axis t ∈ R + . Due to the semigroup property (20) , these solutions induce a continuous semi-flow in the phase space U, which is endowed with a metric by the natural embedding into the space R 2 × L 1 (Π; 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 α ) ∈ U and a periodic solution is a periodic function (I(t), S(t), r α (t)) : R + → U where the last component, viewed as a function r α (t) = r(t; α) : R + × Π → {0, 1} of two variables t ∈ R + and α ∈ Π, 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).

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 θ(t) : R + → R + (transformation of time) inputs I(t) and I(θ(t)) produce identical input-output diagrams such as in Figure ? ?. Due to the rate-independence property of the Preisach operator, given an initial state function r 0 = r 0 α : Π → {0, 1}, there is a continuous function R r 0 : R → R such that every monotone input I(t) : R + → 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)), t ≥ 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 α parameterizes the infinite set of all possible branches. Due to the assumption that the density q = q(α) of the Preisach operator is positive, every branch R r 0 (I) is a decreasing function on the interval 0 ≤ I ≤ 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)) : R + → 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 systeṁ

with

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 α 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İ = 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)) : R + → U will be referred to as a solution of system (2), (13) whenever this doesn't cause confusion.

Throughout the proofs we assume that R nat 0 −R int 0 > 0 is sufficiently small to ensure that the quantities 0 , κ defined below by (27) , (64) using (23), (62) satisfy 0 , κ > 0.

We first recall well-known stability properties of the planar ordinary differential system (25) , which are ensured by the following lemmas. Proof. By assumption, 1/R(0) = 1/R nat < 1 and 1/R(ρ) > 0, hence equation (28) has a solution I * in the interval [0, ρ] ⊂ [0, 1]. The measure density q(α) 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/ρ decreases, this solution is unique.

The Jacobi matrix of system (25) at the infection-free equilibrium (0, 1) has eigenvalues λ 1 = R nat − 1, λ 2 = −R nat − ρ, 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

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

is a Lyapunov function for system (25) 

where we used that I * = ρ(1 − S * ) = ρ(1 − I * /f * ) due to Lemma 2. Equivalently,

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

Hence, the equilibrium (I * , S * ) is globally stable in the positive quadrant.

It turns out that the function V is convex. 

and that f (I) > 0 due to Lemma 1, we obtain κ < 0, which completes the proof.

Let us consider the trajectory (I(·), S(·), r(·)) : R + → U of system (2), (13) and the corresponding planar curve (I(·), S(·)) : R + → R 2 . Consider the sequence T = {t k } of successive switching moments 0 = t 0 < t 1 < t 2 < · · · such that • The restriction of (I(·), S(·)) : R + → R 2 to each interval (t k , t k+1 ), t k , t k+1 ∈ T is a solution of system (25) with R(I) = R r(t k ) (I) satisfying R r(t k ) (I(t))S(t) = 1, t ∈ (t k , t k+1 ); R(I(t k+1 ))S(t k+1 ) = 1.

• In addition, if the set T = {t 0 , . . . , t N } is finite, then the restriction of (I(·), S(·)) :

If the set T is finite, the trajectory of system (2), (13) follows the flow (25), (26) on the interval [t N , ∞) 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 T = {t k } is infinite, i.e. there is an infinite sequence of switches along the trajectory. Let

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

of this function along the trajectory on the time interval [t k , t k+1 ]. (31)), then the increment (32) of the Lyapunov function V k along the trajectory satisfies

where 

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

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, Figure 3 ). Equivalently,

This implies the estimate

In the domain I ≥ I k * , S ≥ 1/R k (I), we have 1/R k (I) ≥ 1/R k (I k * ) = S k * , which implies 0 ≤ S − 1/R k (I) ≤ S − S k * . Hence, relation (30) implies that along the trajectory,

Since V k decreases along the trajectory, the segment of the trajectory with t ∈ (t k , t k+1 ) lies above the level line of V k , hence for this segment I ∈ [I k * ,

The relationİ(t) > 0, t ∈ (t k , t k+1 ), implies that the arc of the trajectory for this time interval admits a parameterization (I,Ŝ(I)), I ∈ [I k , I k+1 ]. Therefore, integrating (35) over the interval [I * k , I b ], we obtain

Further, using the identity ln(1 + x) ≥ x/(1 + x), x ≥ −1, we arrive at

and due to (34) ,

The maximum of the function

Therefore, from (36) it follows that

Since V k decreases along the trajectory, this estimate implies

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

We estimate the left hand side of this equation using ln (1 + x) ≥ x/1 + x, x ≥ −1:

Since S k * < S k+1 , this implies

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

Next, we estimate the differencẽ

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,

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:

.

These equations imply

for an intermediate value ofÎ between I k * and I k+1 * . If we assume that I k+1 * > I k * , then we arrive at a contradiction:

because R k+1 (I) decreases and (22) implies that

Hence, I k+1 * ≤ I k * , and consequently S k+1 * ≥ S k * as well as

Combining (41) and (42), we obtain (40).

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

Proof. By the definition of the Lyapunov function,∆V k = ∆ S + ∆ I , where

and consequently ∆ S ≤ 0. Also, the relation R k+1 (I) ≤ R k (I), I ≤ I k+1 , implies that

thereforẽ

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 * ≤ I k * due to Lemma 5, hence

and using Lemma 5 again,

On the other hand,

where we used that R k+1 (I k

Hence, relations (40) and (41) 

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

Now we combine Lemmas 4 and 6 to obtain the following statement. Proof. From Lemma 1, it follows that

.

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

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 ). 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 implieṡ 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

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

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 ≤ I k * , S ≤ 1/R k (I), we obtain that along the trajectory

and due to R k (I)S < 1, S k * ≤ 1,

Since the trajectory lies below the level line of V k , on the interval [I c , I * ] we have S ≤ S c and dV k dI ≥ ρ I (S k * − S c ), which implies

where (I,Ŝ(I)), I ∈ [I k+1 , I k ] is a parameterization of the segment of the trajectory, and using the identity ln (1

where in the last step we used I k * ≤ 1 and (55). Maximizing the right hand side of this formula over S c ∈ [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))

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

Using the inequality 2(1 + x) ln (1 + x) ≤ x(2 + x), x ≥ 0 we get

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

Further, due to Lemma 1,

, consequently,

.

As in the proof of Lemma 5,

which allows us to conclude (due to R k+1 (Î) ≤ 0) that I k+1 * ≥ I k * and S k+1 * ≤ S k * , and the argument used in the proof of Lemma 5 leads to the estimates (cf. (40) )

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

while

Using the notation ∆ S and ∆ I introduced in the proof of Lemma 6, from (60) we obtain (44) and ∆ S ≤ 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

and due to (58) and (59),

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

and adding (57) leads to (53).

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

For t ≥ t 2 ,

Proof. Consider a point (I k , S k ) = (I k , 1/R k (I k )) where the trajectory satisfieṡ I = 0,Ṡ < 0, see Figure 4 . At this point,

Clearly, the minimal value s k m = min t∈[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 k m ∈ [t k , t k+1 ], i.e. s k m = S(t k m ). Since V k decreases on the time interval [t k , t k+1 ], from V k (I k , S k ) ≤ 2, it follows that

S k * and due to S k * ≥ 1/R nat 0 ,

On the other hand, the minimal value i k m = min t∈[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 k m = I k+1 . Again, since V k decreases on the time interval [t k , t k+1 ], from V k (I k , S k ) ≤ 2, it follows that

Therefore, using I k * = ρ(1 − S k * ) and S k * = 1/R k (I k * ) ≤ 1/R int 0 , we obtain

which completes the proof.

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 ≤ −κ(I k * − I k+1 ) 2 , k = 0, 1, 2, . . .

with κ := ρ 4

independent of k.

Since f (I) = R(I)I ≤ R nat 0 for 0 ≤ I ≤ 1, Lemma 1 implies that

On the other hand, f k (I) = (IR k (I)) = R k (I)I + R k (I) ≤ R k (I) ≤ R nat 0 , f k (I) = IR k (I) ≥ i 0 R int 0 , where i 0 is defined in (62), hence

Similarly, due to 0 ≤ S ≤ 1,

with s 0 defined in (62). Adding (65) and (66), we obtain

Relation (63) implies

and due to the non-negativity of each Lyapunov function V n , Therefore I k+1 −I k * − → 0 (as k − → ∞). Due to (4) and (59), this implies I k −I k * − → 0, and from S k * = 1/R k (I k * ), S k = 1/R k (I k ) it follows that S k − S k * − → 0. Hence, (67) implies v k = V k (I k , S k ) → 0. Since V k is nonnegative and decreases along the trajectory, we conclude that

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

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

where R k (I) ≥ R int 0 and |R k (I k+1 ) − R k (I k * )| ≤ q 0 |I k+1 − I k * | due to (22) , hence

Lemma 6 implies that

hence, using (69),

If V k (I k , S k ) ≥ (a + b + κ)(I k+1 − I k * ) 2 with κ defined by (64), then (70) implies

On the other hand, if V k (I k , S k ) ≤ (a + b + κ)(I k+1 − I k * ) 2 , then (63) implies (71). Therefore, in either case, v k ≤ v 0 p k .

Since V k decreases along the trajectory, V k (I k+1 , S k+1 ) ≤ V k (I k , S k ) = v k . But according to (67),

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

We see that I k * is a Cauchy sequence, hence it converges to a limit value I * . Since I k+1 − I k * − → 0 as k − → ∞, we conclude that I(t) → I * as t → ∞, which implies that each of the functions R 0 (t) and S(t) also has a limit as t → ∞. Denoting the limit of S(t) by S * , the trajectory converges to the equilibrium (I * , S * ).

Hopf bifurcation in symmetric networks of coupled oscillators with hysteresis

Asymptotic stability of continuum sets of periodic solutions to systems with hysteresis

Hysteresis and phase transitions

A generalization of the Kermack-Mckendrick deterministic epidemic model

USA (2021) Implementation of mitigation strategies for communities with local COVID-19 transmission

Rational behavioral response and the transmission of STDs

Imperfect vaccine and hysteresis

On uniqueness of an initial-value problem for ODE with hysteresis

Global dynamics of SIR model with switched transmission rate

Houston region passes threshold for tougher COVID-19 restrictions

Nebraska (2020) Phased public health restrictions tied to coronavirus hospitalization rate

A genomewide association study identifies six susceptibility loci for chronic lymphocytic leukemia

Information-related changes in contact patterns may trigger oscillations in the endemic prevalence of infectious diseases

The interplay between voluntary vaccination and reduction of risky behavior: A general behavior-implicit SIR model for vaccine preventable infections

Experimental research in magnetism

The impact of COVID-19 on small business owners: Evidence from the first three months after widespread social distancing restrictions

Distributed formation control via mixed barycentric coordinate and distancebased approach

Hysteresis can grant fitness in stochastically varying environment

Rational epidemics and their public control

Governmental public health powers during the COVID-19 pandemic: Stay-at-home orders, business closures, and travel restrictions

Willingness to get the COVID-19 vaccine with and without emergency use authorization

A core group model for disease transmission

Noninvasive stabilization of periodic orbits in symmetrically coupled systems near a Hopf bifurcation point

Kluwer encyclopedia of mathematics: Supplement 1

Pattern formation by bacteria

A contribution to the mathematical theory of epidemics

Dynamics of SIR model with vaccination and heterogeneous behavioral response of individuals modeled by the Preisach operator

Systems with Hysteresis

On a bifurcation governed by hysteresis nonlinearity

On the existence of cycles in autonomous systems

Hysteresis, convexity and dissipation in hyperbolic equations

Analysis of optimization problem for a piezoelectric energy harvester

A global survey of potential acceptance of a COVID-19 vaccine

Differential equations with hysteresis operators: Existence of solutions, stability, and oscillations

Quantifying the effects of quarantine using an IBM SEIR model on scalefree networks

Evaluating the effectiveness of social distancing interventions to delay or flatten the epidemic curve of coronavirus disease

Mathematical models of hysteresis and their applications

Memory effects in population dynamics: Spread of infectious disease as a case study

Realization of arbitrary hysteresis by a low-dimensional gradient flow

Convergence of direct recursive algorithm for identification of Preisach hysteresis model with stochastic input

Evolving public perceptions and stability in vaccine uptake

Dynamics of SIR model with heterogeneous response to intervention policy

Emergence of hysteresis loop in social contagions on complex networks. Sci Rep

Hysteresis and semigroups: In differential models of hysteresis

Unification of theoretical approaches for epidemic spreading on complex networks

Can we contain the COVID-19 outbreak with the same measures as for SARS?

Stability of stationary sets in control systems with discontinuous nonlinearities