Oscillatory behaviour on a non-autonomous hybrid SIR-model

. We study the impact of some abstract agent intervention on the disease spread modelled by a SIR-model with linear growth infectivity. The intervention is meant to decrease the infectivity, which are activated by a threshold on the number of infected individuals. The coupled model is represented as a nonlinear non-autonomous hybrid sys-tem. Stability and reduction results are obtained using the notions of non-autonomous attractors, Bohl exponents, and dichotomy spectrum. Numerical examples are given where the number of infected individuals can oscillate around a equilibrium point or be a succession of bump functions, which are validated with a tool based on the notion of δ -complete decision procedures for solving satisﬁability modulo theories problems over the real numbers and bounded δ -reachability. These ﬁndings seem to show that hybrid SIR-models are more ﬂexible than standard models and generate a vast set of solution proﬁles. It also raises questions regarding the possibility of the agent intervention been somehow responsible for the shape and intensity of future outbreaks.


Introduction
Mathematically, the model of choice to represent the dynamics of the epidemic is the SIR-model and its variants; introduced by Kermack and McKendrick [11].Since then, the literature on the subject is quite vast.However, one of the key issues in the subject is that the simplicity of (autonomous) SIR-models do not produce solutions with complicate oscillatory behaviour.
Although less common, non-autonomous SIR-models have been introduced and studied in the literature and may overcome partially such limitations, e.g.see [2,18,23,4,3,16,21]. Usually these models introduce some type of seasonality behaviour, for example, through a periodic infectivity function.Indeed, Bacaër et al. [2] introduced a generalization of the basic reproduction number, and Boatto et.all [4] considered a SIR-model with birth and death terms and time-varying infectivity as a sinusoidal, showing that the (average) basic reproduction number, the initial phase, the amplitude and the period are all relevant issues.Moreover, they show the existence of a periodic orbit.Bai et.al. [3] studied a model with a seasonal contact rate and a staged treatment strategy, showing two different bistable behaviours under certain conditions: the stable disease-free state coexists with a stable endemic periodic solution, and three endemic periodic solutions coexist with two of them being stable.Kuniya [16] deals with an age-structured SIR epidemic model with time periodic coefficients, obtaining the basic reproduction number as the spectral radius of the next generation operator and showing that it plays the role of a threshold value for the existence of a nontrivial periodic solution; based on a Krasnoselskii fixed point theorem argument.Another approach to produce non-autonomous systems is to couple different SIR-models by a non-autonomous function, e.g.Rocha et al. [21] introduced a tuberculosis (TB) mathematical model, with 25 state-space variables where 15 are evolution disease states (EDSs), which takes into account the (seasonal) flux of populations between a high incidence TB country (A) and a host country (B) with low TB incidence, where (B) is divided into a community (G) with high percentage of people from (A) plus the rest of the population (C).
In this work, we consider an infectivity function which grows linearly (i.e. the most simple non-autonomous function), but we also study the effect of an agent intervention on the model in the form of some action policies.The policies are meant to reduce to decrease the infectivity, which are activated by a threshold on the number of infected individuals.Such approach turns the full model into a nonlinear non-autonomous hybrid system, see Section 2. Stability and reduction results are obtained using the notions of non-autonomous attractors, Bohl exponents, and dichotomy spectrum, which are presented in Section 3. In Section 4, we give some numerical examples, where the number of infected individuals can oscillate around a equilibrium point or be a succession of bump functions.The last example is quite interesting and raises questions regarding the possibility of the agent intervention time been somehow responsible for the shape and intensity of some future outbreaks.We end this work with some brief concluding remarks.

The class of non-autonomous ODEs
Consider the basic SIR epidemic model together with a piecewise linear continuous infection coefficient β ξ , described by (a)

and (b)
where γ > 0, α ≥ 0, ζ ≥ −α, and ξ ∈ R is a bifurcation parameter; e.g. for ξ = 0 the model is autonomous.The values S(t), I(t), R(t) are, respectively, the number of healthy individuals (susceptible), infected individuals and recovered individuals; and α is a parameter of birth and death, γ is a recovery rate without possibility of re-infection, and ζ accounts for the rate of individuals that become healthy but may be re-infected in the future.We assume a (normalized) constant population S + I + R = 1, so the system (1)(a) evolves on the simplex defined by Σ 1 = {(S, I, R) ∈ R 3 : S, I, R ≥ 0, S + I + R = 1}, meaning that system (1) may be written as (a) and (b) For mathematical reasons, which will be clear in what follow, e.g.use of pullback limits and Bohl exponents, we work with an unbounded from below time interval T = (−∞, T ], for T > 0, with an initial time t 0 ∈ (0, T ).For convenience, from now on, we use the notations T t0 = [t 0 , T ], SIR(ξ) to describe the set of equations ( 2) for a given parameter ξ ∈ R, which account for the (linear) increase/decrease ratio of the disease.
Clearly, equation (2)(b) may be extended by using other growth functions, e.g.accounting for saturation phenomena, instead of the simple linear change in the infection coefficient.However, for the purpose of this work, such is enough in order to discover the main differences from the standard (autonomous) SIRmodel, vastly used in the literature.

Non-autonomous hybrid SIR-models generated by simple action policies
In the model under study, we have two main entities, i.e. the natural evolution of the disease (nature) versus the evolution of the disease together with some abstract agent action with the purpose of reducing the transmission rate.Each will have a on/off-state, but makes sense to suppose nature is always in the on-state when the agent action is in the off-state, and vice-versa.Since they are complementary, we consider states in the viewpoint of the agent action.Further, to model the action from the agent, we assume that it depends on the current number of infected individuals I(t) and has a maximum fixed time of intervention T * > 0. For that, we establish two threshold values as triggers to the on/offstates, namely, I b ∈ (0, 1] and I s ∈ [0, 1).Then, the agent strategies considered are: (S 0 ) the action starts at time t ∈ T t0 , if it was in the off-state and I( t) = I b , then stops at time t = t + T * (i.e.I b = 1); (S 1 ) the action starts at time t ∈ T t0 , if it was in the off-state and I( t) = I b , then stops at the first time t > t with I(t) = I s (i.e.T * = +∞); (S 2 ) the action starts at time t ∈ T t0 , if it was in the off-state and I( t) = I b , then stops when (S 0 ) or (S 1 ) are satisfied.
Although in general, in each on/off-state, we may have different behaviours, e.g.applying different techniques to reduce the (time dependent) transmission rate, for this work we assume that there is only one behaviour in the on-state.Such corresponds to restrict the values of the parameter ξ to the set {β − , β + }, for given constants β − < 0 < β + .The value β + accounts for the natural increasing effect of the disease (i.e.no agent intervention) and β − accounts for the result of an agent action for controlling/reducing the transmission rate.Thus, in this model, ξ ≡ ξ(t) turns now to be a piecewise function defined on T t0 with values on {β − , β + }, where the discontinuity instances are precisely the switching times generated by the application of one of the agent strategies (S 0 )−(S 2 ).Moreover, ξ(t 0 ) = β + , I(t 0 ) < I b , and the system may alternate (none or some bounded number of times) between the values β − and β + .Hence, it makes sense to define ξ − ⊆ T t0 as the support where ξ(t) = β − , ξ + ⊆ T t0 as the support where ξ(t) = β + , N ξ ∈ N 0 the number of switchings, and t 0 , t 1 , . . ., t N ξ , with t 0 < t 1 < • • • < t N ξ < T , the corresponding times of switchings.
Realistic constraints impose further that, in equation ( 2)(b), we assume β ξ (t) > 0 on T t0 ≡ [0, T ] and I(t) > 0 on T t0 (i.e. in the time window there are always infected individuals), otherwise the problem is not interesting or meaningful.Additionally, for mathematical reasons, we require that, besides β ξ being a continuous integrable bounded function on T t0 , to be defined also on R\T t0 .In particular, we will have the following structure where χ S (t) is the characteristic function of the set S. Therefore, there are positive constants β * and β * , such that β ξ (t) ∈ [β * , β * ] for t ∈ R.
Regarding the triggers values there are two situations: I b > I s and I b < I s .The most natural situations is I b > I s , but I b < I s makes sense in specific and limit situations.In either cases, because of the agent action, β ξ (t) ≡ β ξ (t, I(t)) and there is a memory effect, not present in equations (2), which controls in which state the system is running.In general, the model under study is neither an ordinary differential equation or a differential inclusion, but can be treated in the setting of (generic) hybrid systems, e.g.see [10,20] for definitions and properties.
[off] SIR(β + ) [on] SIR(β − ) The hybrid model is generally described in Fig. 1 and in more detail in Fig. 2, when expanding the invariant sets and dealing with the situations I b > I s and I b < I s .

Stability and bifurcation in each node
Restricting to a hybrid system node and for mathematical reasons, we may assume that the infectivity function β ξ , satisfying (2)(b) on [t 0 , T ], is defined in all R with the structure (4) where β 0 = β λ (t 0 ) > max{0, −ξT } and χ S (t) is the characteristic function of the set S. Therefore, there are positive constants β * and β * , such that β λ (t) ∈ [β * , β * ] for t ∈ R. In fact, β * = min{β 0 , β 0 + ξT } and β * = max{β 0 , β 0 + ξT }.Further, each node is a nonlinear non-autonomous ordinary differential equation, so standard results of (autonomous) SIR-models do not apply.Such happen even for non-autonomous linear system as x = A(t)x.For instance, the wellknown fact that the origin is globally asymptotically stable if the real part of all eigenvalues of the matrix A are negative turn out to be wrong, as the following Nemytskii-Vinograd counterexample shows which has the constant eigenvalues λ 1 = −1 and λ 2 = −1, but have the fundamental matrix X(t) ≡ X(t, 0) = e t sin(2t) e t cos(2t) , so solutions are unstable.To overcome these difficulties, we will use tools as the Chueshov's notion of nonautonomous equilibrium solution in a random dynamical system through a pullback limit, and the notions of Bohl exponents and exponential dichotomy.We start by briefly collecting stability results for the autonomous SIR-model in the next section.

Stability of the autonomous SIR-model
The corresponding autonomous SIR-model of equations (2)(a), i.e. when ξ = 0, may be written as where R 0 is the so-called basic reproduction number.When R 0 ≤ 1, the disease-free equilibrium (S * , I * ) = (1, 0) is a globally asymptotic equilibrium point, proved by using Lyapunov-LaSalle function V (S, I) = I and the LaSalle's Invariance Principle in the compact positively invariant set Σ 1 .The endemic equilibrium only belongs to the simplex In fact, for R 0 > 1, the disease-free equilibrium is unstable and the endemic equilibrium is asymptotically stable on Σ 1 \M 0 , where M 0 = [0, 1] × {0} is the stable manifold of the disease-free equilibrium, see [1].The stability is obtained by using the Lyapunov function For the linear growth infectivity setting, i.e. ξ = 0, the disease-free equilibrium (S * (t), I * (t)) = (1, 0) is still valid (in some regimes) but the endemic equilibrium ( S, Ī) do not make sense a priori as an equilibrium point, since the ODE is non-autonomous.Hence, we introduce the notion of nontrivial nonautonomous attractor, obtained by a pullback limit mechanism when using the skew product flow formalism, in the next two sections.

The skew product flow formalism
From Lemma 1, for each initial condition x 0 ∈ Σ 1 , there exists a unique solution x(t; t 0 , x 0 ), so we may define the flow φ t,t0 (x 0 ) = x(t), which satisfies the following set of conditions: is differentiable and (t, t 0 , x 0 ) → φ t,t0 (x 0 ) is continuous.
Let (X, d X ) and (P, d P ) be metric spaces.A skew product flow (, ϕ) is defined in terms of a cocycle mapping ϕ : R + 0 × P × X → X which is driven by an autonomous dynamical system ψ : R × P → P acting on a base or parameter space P and the time set R. For convenience, we write ϕ p t (x) to denote ϕ(t, p, x).The driving system ψ on P is a group of homeomorphisms (ψ) t∈R under the composition on P with the properties that x for all (p, x) ∈ P × X; (b) ϕ p t+s = ϕ q t • ϕ p s with q = ψ s (p) for all s, t ∈ R + 0 , and p ∈ P ; (c) the mapping (t, p, x) → ϕ(t, p, x) is continuous.
In particular, system (2) can be seen as and (b) So, we have One of the advantages of this formalism is that, in our case, X and P are both compact metric spaces.Natural extensions to random dynamical systems are obtained when replacing P by a probability space and the continuity property in (P 1 )(c) by measurability.The solutions are then generated by solutions of the corresponding Itô stochastic differential equation, see [8] as an introduction to the subject.

Non-autonomous attractors
An entire solution is a continuous function u : R → Σ 1 such that u(t + s) = ϕ p t (u(s)) for all s ∈ R and t ∈ R + 0 .We say that A is a non-autonomous attractor p-family if it is a set of nonempty compact subsets A p ⊆ Σ 1 such that ϕ p t (A p ) = A ψt(p) for all t ∈ R + 0 and p ∈ P .Such sets are made up of entire solutions.Let dim S (A, B) denote the Hausdorff semidistance between the nonempty compact subsets A and B of S. We call a non-autonomous attractor p-family A a pullback attractor p-family if it holds the pullback convergence and a forward attractor p-family if it holds the forward convergence If the convergence is uniform in p, the pullback and forward convergences coincide.A pullback absorbing set B is a nonempty subset of Σ 1 such that, for all p ∈ P and (bounded) D ⊆ Σ 1 , there exists a time T ≡ T (p, D) > 0, and If there exists a pullback absorbing set B then A is a pullback attractor p-family if we define There is a corresponding formulation for t-families.Recall the system representation (2.3) with flow φ t,t0 .We say that A is a non-autonomous attractor t-family if it is a set of nonempty compact subsets A t ⊆ Σ 1 such that φ t,t0 (A t0 ) = A t for all t ≥ t 0 .Then it is a pullback attractor t-family if This notion will play an important role in finding (nontrivial) non-autonomous equilibrium points.
Considering the complexity of the model ( 2), in a first step, we study a sub-case proposed in [13], where the equation for R(t) will not appear and the system evolves on the simplex Thus, it have the form (a) and will be denoted by SI(ξ)-model, although other variations exist on the literature with similar notation.For studying the stability properties of (10), it will be relevant the following results which may be checked computationally.
, where erf is the error function and .
Lemma 3. The general Bernoulli differential equation for some arbitrary functions a and b, has the unique solution with explicit solution (see Lemma 3) For the t → β ξ (t) function, ξ 0 ∈ R is the so-called shovel bifurcation parameter, due to a change in the range, and ξ is a transcritical bifurcation parameter, due to a change in the amplitude (see [12]).
Note that no information is given in the case β * ≤ α + ζ < β * .A direct use of the above Lemma gives the following result.Lemma 6.The globally asymptotically stability of E 0 w.r.t.Σ 0 , satisfies: it is stable when ξ 0 ≤ −T β + and unstable when ξ 0 > 0.
We now define the so-called Bohl exponents, introduced by Piers Bohl [5], which give information on the uniform exponential growth, whereas Lyapunov exponents, introduced by Aleksandr M. Lyapunov [17], only measure the exponential growth.It is well-known from ODEs that the Bohl exponent compared with the Lyapunov exponent is the appropriate concept in the setting of nonautonomous systems.For a summary of the history of Lyapunov and Bohl exponents see [7].In detail, the upper Bohl exponent of a locally integrable function f : T → R is defined by and the lower Bohl exponent by The exponents are finite when f is integrally bounded, besides other properties (e.g.see [7]).From now on, we use the convection a + [b, c] = [b + a, c + a] for any a, b, c ∈ R.
We say that a linear system x = A(t)x has an exponential dichotomy (for short, E.D.) on R, if there exists a projection P : R n → R n and positive constants C, α, β such that The dichotomy spectrum [22] is the set which is considered in the literature as the appropriate counterpart to eigenvalues in the non-autonomous setting.Then dichotomy spectrum is related with the Bohl exponents in the following way.
Lemma 9 is false for SI(ξ), because sign(S ) = −sign(I ), and also shows that the flow, associated with the hybrid system of Fig. 2, can be quite complex since each node can (generally) initiate in any monotonicity situation.
Poincaré, in 1892, started the theory of normal forms as a technique to simplifying a nonlinear system in the neighborhood of a reference solution by a smooth change of coordinates.Let us summarize the ideas.Consider the autonomous system x = Ax + f (x), where A is a constant matrix and f (x) = O( x 2 ) as x → 0. Then by a formal coordinate transformation x = y + +∞ i=2 h i y i the above system can be changed into the system y = By + +∞ i=1 g i y i where B us the complex Jordan form of A, g i = (g 1 i , . . ., g n i ) and for ξ i eigenvalues of A, p ∈ Z n + and n i=1 p i ≥ 2. In addition, if f is analytic in the origin, we have the so-called Poincaré-Dulac's analytic normal forms.Consider the systems (a) A change of variables x = P (t)y is said to be a Lyapunov-Perron transformation (for short, L.P.) if P (t) is nonsingular for all t ∈ R and P, P −1 , P are uniform bounded in t ∈ R. The system (15)(a) is locally analytically equivalent to the system y = G(t, y) if there exists a coordinate substitution x = P (t)y + h(t, y) which transforms one to the other, where f, G, P, h are analytic in Bρ (0) × R, for some ρ > 0, f (t, 0) = G(t, 0) = h(t, 0) = 0, P is a LP transformation and h(t, y) = O( y 2 ) as y → 0. Assume the dichotomy spectrum of (15 Lemma 10. (see [24]) We have for type-II, system (15)(a) is locally analytically equivalent to its linear part (15)(b).
In this section, we always assume x = (S, I, R), and γ > 0. Under the linear transformation the system (2)(a) is transformed into (the Jordan canonical linear form) .

Three illustrative examples
We present numerical examples that illustrate occurrences that may not appear in SIR-models without agent actions: (E 1 ) the agent action is not able to decrease the number of infected individuals and they tend to a non-autonomous attractor (see section 3.3 for a precise definition); (E 2 ) the agent action introduces an oscillatory behaviour in the number of infected individuals around some nonautonomous attractor; and (E 3 ) the agent action are able, in each period, to significantly decrease the number of infected individuals but such mechanism introduces a succession of bumps along time.
Since hybrid systems return solutions selected by discrete jump events its validation and error control is a key issue, requiring tailored tools based on first order logics.The numerical calculations of the given examples below were produced using two packages dReal and dReach [9,14,15].dReal is an automated reasoning tool focused on solving problems that can be encoded as first-order logic formulas over the real numbers by implementing the framework of δ-complete decision procedures.dReach deals with the bounded δ-reachability problem.For a hybrid system H =< X, Q, f low, jump, inv, init >, where flow, jump, inv, init are SMT formulas that dReal can handle and specifying a numerical error bound δ, any formula φ can have its δ-perturbation counterpart φ δ .Then, a δ-perturbation of H is defined as H δ =< X, Q, f low δ , jump δ , inv δ , init δ >, by relaxing the logic formulas in H. Now, choosing n ∈ N to be a bound on the number of discrete mode changes, T ∈ R + an upper bound on the time duration, and unsafe to encode a subset of X × Q, the bounded δ-reachability problem asks for one of the following answers: (a) "safe" if H cannot reach unsafe in n steps within time T ; (b) "δ-unsafe" if H δ can reach unsafe δ in n steps within time T .In this way, we ensure that our examples are numerically correct, since the SMT tool produce a logical proof of reachability.Although, there is no stantard basic reproduction number in the non-autonomous setting (but there exist generalizations as the notion [2]), we would still look to the value of R 0 (t) =

Conclusions
A non-autonomous hybrid SIR-model was introduced as the result of agent action policies on diseases modelled by a SIR-model with linear infectivity growth.The coupled system shows a great variety of profile solutions and extends the standard SIR-model (i.e. when max I(t) < I b ).Two ingredients make the problem difficult: (a) its non-autonomous nature; and (b) the jumps between ODEs (i.e. the hybrid system nodes) are controlled by the values of the state variable I(t).This work is a first step to study the properties of such hybrid SIR-models, since several issues are still to be clear, e.g.complete scheme of stability of the hybrid system, existence of nontrivial periodic solutions crossing several nodes, behaviour of the system under the assumption I s > I b , etc.. Nevertheless, the stability results obtained and examples provided already show the richness and potential application of the model to better fit oscillatory real data.Considering that the number of infected individuals is the most observable variable in reality, Fig. 5 turns out to be quite interesting, since the choices of I b , I s are critical to determine the solution profile.Hence, it raises several questions: (a) Are human disease control strategies somehow responsible for the oscillatory behaviour of some diseases?(b) For each choice of parameters and agent action, are there optimal values for I b , I s such the maximum of I(t) is reduced?(c) How different is the solution if the agent action activates by a stochastic process?

Appendix
Proof.(Lemma 1) Let h = y − x for x, y ∈ Σ 1 .Recall that the mean value inequality, for a vector value function F : R × Σ 1 → R 3 , says that when the Jacobian matrix of F at w = x + τ h, i.e.J F (w), is uniformly bounded by some constant L > 0 for any τ ∈ [0, 1] and t ∈ R, then Hence, the function F in equation (2.3) is locally Lipschitz continuous in the second variable since and w(t) is bounded for any τ ∈ [0, 1], by definition.So, the Picard-Lindelöf theorem ensures the existence and uniqueness of solution in each node of the hybrid system.Its not difficult to see that the solution is globally defined.Further, the hybrid system Fig. 2 is deterministic and has only one jump condition in each node, so we conclude the proof.

Lemma 4 .
Let t ∈ [0, T ].Define the complex value function and, since I(t) > 0 and (ψ I +ψ R )I = β ξ SI −ζI −α(I +R) = β ξ SI −ζI −α(1−S), we have that SIR(ξ) can be written as (a)    S = −(ψ I + ψ R )I, I = ψ I I, R = ψ R I, and (b) b p .We say that system (15)(b) is of type: (type-I) when a 1 b p > 0 (i.e. it is in the Poincaré domain); (type-II) when a 1 b p > 0 and it is non-resonant, i.e. 0 ∈ p i=1 a i m i − a j , p i=1 b i m i − b j with m ∈ N p ; (type-III) when a 1 b p > 0 and A(t) is block diagonal w.r.t. the spectral interval [a i , b i ] (i.e. of Poincaré-Dulac type).