The study of the singular bifurcation of a 3D dynamical sys-tem of interest in ecology

. The classical bifurcations that occur in the SIR system were already studied in many papers. Our attention is focused on the singular bifurcations, when the system has in inﬁnite number of periodic points at the bifurcation limit. The dynamics of the system is analyzed and the results are interpreted from a practical point of view.


Introduction
Compartmental models of epidemic spread are widely used in practice in order to understand the mechanism of the propagation of some contagious diseases. Despite the fact that they rely on various simplifying assumptions (like the homogeneity of the population inside each compartment, neglecting the stochastic effects or the use of well defined parameters), these models have a main advantage: the impact of the modification of each parameter on some qualitative dynamical properties of the system can be easily studied. The analytical study of such models, in particular of their bifurcations, is important not only from theoretical point of view, but also from a practical one: it is useful to have an overview on all possible behaviors of the system because it is sometimes difficult to obtain the precise values of some parameters from statistical data.
In this paper we will study the bifurcations that occur in a SIR-type model, our attention being focused on the singular bifurcation, when the system has in infinite number of periodic points at the bifurcation limit.
The SIR model (susceptible, infected, recovered) was introduced by Kermack and McKendrick in the beginning of the last century ( [1]): where S = S (t) , I = I (t), R = R (t) represent the susceptible, infected, respectively recovered population and N (t) = S (t) + I (t) + R (t).
The parameters of the system are the transmission rate β (the product of the average number of contacts per person and per unit time by the probability of disease transmission in a contact between a susceptible and an infectious individual), the transition rate γ ( 1 γ is the average duration of the infection of an individual).
From (1) it results that N (t) = 0, so N (t) = cst = N (0). One can consider the fractions s = S /N, i = I/N, r = R/N of the susceptible, infected, and recovered individuals among the whole population.
It is a closed system, with no incoming persons and no persons who leave the population. In this paper we will study a slightly different system, which takes into account the incoming constant population "a" (newborns, immigrants, travelers etc.) and the persons who leave the population with the rate µ (deaths, emigrants, travelers).
A simpler model, corresponding to a = µ was proposed and analyzed in [2]. We are interested to study the singular bifurcation that occurs for a = µ = 0. It will be shown that the open system (3) has a completely different asymptotic behavior comparing with the closed system (2).
We study (3) for small parameters a and µ using a specific tool, namely the entry-exit function, and we explain the transition between the two systems in a long but finite time interval.
The paper is organized as follows: section 1 is a short presentation of the problem we will study, section 2 is devoted to the analysis of the closed system, in section 3 the main results concerning the open system are presented, in section 4 we study the dynamics of the open system for small parameters a and µ and some conclusions are presented in section 5.

The closed system
Some properties of the closed system (2) were already pointed out in many papers, but we re-interpret part of them.
An important and classical observation is that the domain is invariant under the action of the flow of (2). Because n(t) = s (t) + i (t) + r (t) = 1, one can conclude that r (t) = 1 − s (t) − i (t), so the reduced system perfectly describes the dynamics of (2) Remark 2.1 R 0 = β/γ is called the basic reproductive number and can be interpreted as the average number of persons infected by one case in a totally susceptible population, in absence of interventions aimed to control the infection. This is one of the most important parameters in modelling of any epidemics.

Remark 2.2
The corresponding equilibria of (2) are E = (s, 0, 1 − s), s ∈ [0, 1]. They are called disease free equilibria because finally there are not infected people in the population, but only susceptible and recovered people.

Remark 2.3
The existence of the heteroclinic orbits is in agreement with the fact that the equilibrium point E 1/R 0 = γ β , 0 has λ 1 = λ 2 = 0 so its central manifold is two dimensional. In figure 1 is presented the phase portrait of SIR system with β = 1/10, γ = 1/6, i.e. R 0 = 0.6 < 1. It can be observed that all equilibrium points E s , s ∈ [0, 1] are normally attractive.

The open system
The open system takes into consideration incoming and people who leave the population, with the rate "a", respectively "µ": It is a slight generalization of the system proposed and studied in [2] , where a = µ, i.e. the incoming rate is equal with the outgoing rate.
Because (s + i + r) = a − µ (s + i + r), one may observe that the total population n (t) = . Using the natural assumption n 0 < a/µ, we observe that n is an increasing function and lim t→∞ n (t) = a µ . It results that the set gives the main information on the behavior of (5). It will be studied in what follows. We define the basic reproduction number of (5) by R 0 = aβ µ(γ+µ) . It has the same empirical interpretation as for the closed system: Proposition 3.1 a) The system (6) has the equilibrium point E 0 = a µ , 0 (the free disease equilibrium point) b) For R 0 > 1 the system (6) has the equilibrium point E * = γ+µ β , µ β R 0 − 1 (the endemic equilibrium point).

It results that E 0 is asymptotically stable if and only if
More of this, it can be shown that the equilibrium point E 0 , respectively E * are globally stable. In [3], [4] the Lyapunov functions are obtained for the SIR system and other epidemiological models. The next proposition is based on these results. Proposition 3.3 a) for R 0 < 1 the disease free equilibrium E 0 is globally asymptotically stable b) for R 0 > 1 the endemic equilibrium E * is globally asymptotically stable Even if the global stability is stronger that the local stability, we considered important to point out the last one because some results are useful in the following analysis.

Fast-slow oscillations in the open system
We consider the system (6) for small a > 0 and µ > 0. For a → 0, µ → 0 the system (6) approach the system (4), however their dynamics are completely different, which gives rise to a singular bifurcation. The aim of this section is to explain the behavior of (6) for very small parameters and to give an impression on the bifurcation mechanism in the frame of the singular perturbation theory. A singular perturbation problem is one for which the flow of the perturbed system is qualitatively different from the flow of the unperturbed one and the aim of the singular perturbation theory is to determine the behavior of the solution of the perturbed system when the perturbation is going to 0. The theory that is widely used for explaining this behavior is the geometric singular perturbation theory initiated by N. Fenichel in 1979 [5].
It is addressed to the so called "fast-slow systems", i.e. systems described by ẋ = f (x, y, ε) y = ε · g (x, y, ε) where x ∈ R n is the fast variable and y ∈ R m is the slow variable (with a slow variation due to the fact that y << 1 when 0 < ε << 1). The fast time is t, i.e. x = x (t) and y = y (t). The slow time "τ" is obtained by the change of variable τ = t/ε. The corresponding system is The systems (7) and (8) are equivalent for ε > 0, but they give rise to systems with very different behavior for ε = 0.
From (7) is obtained the fast subsystem (the layer problem) ẋ = f (x, y, 0) y = 0 that describes the evolution of the fast variable "x" far from the critical manifold that describes the evolution of the slow variable "y" on C 0 . Fenichel's theory requires certain assumption on C 0 : there is a m-dimensional manifold M 0 ⊂ C 0 which is normally hyperbolic (i.e. all eigenvalues of D x ( f (x, y, 0)) have non-zero real part on M 0 ) and locally invariant (i.e. the flow can leave M 0 only through its boundary).
In practical term, Fenichel's theorems show that, for ε > 0 sufficiently small, the systems (7), (8) are regular perturbations of (9), (11) and the orbits of (8) can be obtained by the concatenation of orbits of the fast and slow subsystems In our study, for fixed parameters a, µ, β, γ we will consider the system where ε > 0 is a small parameter and R 0 > 1. Obviously lim The system (12) is a fast-slow system, for the moment not being written in the standard form. To obtain the standard form we perform a change of variables: In (13) the fast variable is "x" and the slow variable is "y", so "i", respectively "s" are the fast, respectively slow variables in (12) The critical manifold is It is not normally hyperbolic in (0, 0) , which correspond to 1 R 0 , 0 in the system (4). It means that the Fenichel's theory may not be applied. There are more sophisticated tools that could be applied in such situations. One of them, applicable in our case, is the "entry-exit function" [6], [7], [8] that gives an estimate, in terms of a Poincare map, of the behavior of the orbits near the point in which the critical manifold changes stability. The following analysis is based on the results presented in [6], which addresses to systems described by with (x, y) ∈ R 2 and sgn ( f (0, y, 0)) = sgn (y) for all y 0, g (0, y, 0) > 0 for all y ∈ R.
In this case The assumption sgn ( f (0, y, 0)) = sgn (y) for all y 0 shows that the axis x = 0 is formed by equilibria that are normally attracting for y < 0 and normally repelling for y > 0 .
For ε > 0 the axis x = 0 remains invariant and the assumption g (0, y, 0) > 0 for all y ∈ R shows that the flow is going from the left to the right. For small ε > 0 , the orbit starting from (y 0 , x 0 ) with x 0 > 0 enough small and y 0 < 0 is attracted quickly toward the axis x = 0, then drifts to the right along the axis, passes beyond the line y = 0 and finally is repelled from the axis. It re-intersects the line x = x 0 at a point (p ε (y 0 ) , x 0 ). As ε → 0, the return map p ε (y 0 ) approaches a function p 0 (y 0 ) .
In other words, the orbit does not leave the y-axis as soon as it becomes unstable at y = 0, but stay near the y-axis until the repulsion balances the attraction that occurred before y = 0. This phenomenon has been called "bifurcation delay" [9] and comprehensive analysis is made in [10].
Some results presented in [6] can be summarized in the following theorem.
The system (13) fulfills the hypothesis of the previous theorem.
In this case f (x, y, ε) = γ · y − εµ , g (x, y, ε) = aR 0 − (µ + βx) (y + 1) and the physically relevant interval for y is −1 < y < a µ R 0 − 1.  x = x 0 in the point (p ε (y 0 ) , x 0 ) that is close to (p 0 (y 0 ) , x 0 ), where p 0 (y 0 ) > 1 R 0 is the unique solution of the equation In order to construct (in a semi-analytical manner) this orbit of (12) starting from (s(0), i(0)) with s(0) > 1/R 0 , we come back to the old variables i and s and use the following algorithm: 1) determine (q 0 , i * ) the intersection of the heteroclinic orbit of (4) passing through (s(0), i(0)) with the line i = i * . q 0 < 1/R 0 is the solution of the equation repeat the steps 1-4 for the initial value (s 1 , i * ) Because the maximum of i on the heteroclinic orbits is obtained for s = 1 R 0 it results that the two sequences of points (s n ) n∈N and (q n ) n∈N are related through the inequalities q 0 < q 1 < ... < q n < ... < s n < ... < s 1 .
It is reasonable to suppose that the two sequences have the same limit s ∞ = q ∞ ≈ 1 R 0 . This result is in agreement with the existence of the endemic equilibrium point E * = γ+εµ β , εµ β R 0 − 1 which is an asymptotically stable focus. The point E * approaches the non-hyperbolic point on the critical manifold 1 R 0 , 0 when ε → 0 For ε > 0 enough small, the orbit of (s (0) , i (0)) is approximated by the concatenation of fast orbits, that are heteroclinic orbits of (12) for ε = 0 (i.e. SIR system) and slow orbits obtained from the system (13) through the previous presented algorithm.
The intersections of the orbit with the line i = i * = 0.0012 are denoted with s 1 , s 2 ,... respectively q 0 , q 1 , q 2 , q 3 .... . They are 0.0649 < 0.2595 < 0.2801 < 0.2939 < ... < 0.3824 < 0.4005 < 0.4299 For the same values of the main parameters of the system and ε = 0, 001, the intersection of the orbit with the line i = i * is 0, 0602 < 0, 2776 < 0, 2896 < 0, 2985 < ... < 0, 3750 < 0, 3859 < 0, 4020 In figure 5 the evolution of the susceptible and infected populations is presented. The fast oscillations of the infected people correspond to a short period when the disease spreads rapidly and they are followed by relaxation periods when the disease is almost eradicated (because i(t) ≈ 0). It can be observed that the two populations exhibit slow, respectively fast oscillations in the proximity of the endemic equilibrium point (which confirm the analytical results). The time t spent (ε) spent near the the axis i = 0 after the first intersection (q 0 , i * ) with the line i = i * is larger when ε is smaller: t spent (0.01) = 1640, t spent (0.005) = 2990, t spent (0.001) = 13211 and t spent (0.0005) > 35000.
In the same time the variable s increases more slowly when ε decreases. It can be concluded that the open system has a completely different asymptotic behavior compared to that of the closed system. However, for small enough ε, the orbit approaches the heteroclinic orbit passing through (s (0) , i (0)) and remains a long time near its intersection with the axis i = 0. In such a way, it mimics the behavior of the closed system in a finite time interval (longer when ε is smaller)

Conclusions
In this work we studied the bifurcation that occurs in an open SIR system, when the incoming and outgoing rates are 0, so the system becomes closed.
We found that the asymptotic behavior of the two system is completely different, so the bifurcation is singular.
It was also found that the open system approaches, for very small values of the incoming and outgoing rates, the dynamics of the closed system in a finite time interval.