SIR model with vaccination: bifurcation analysis

There are few adapted SIR models in the literature that combine vaccination and logistic growth. In this article, we study bifurcations of a SIR model where the class of Susceptible individuals grows logistically and has been subject to constant vaccination. We explicitly prove that the endemic equilibrium is a codimension two singularity in the parameter space $(\mathcal{R}_0, p)$, where $\mathcal{R}_0$ is the basic reproduction number and $p$ is the proportion of Susceptible individuals successfully vaccinated at birth. We exhibit explicitly the Hopf, transcritical, Belyakov, heteroclinic and saddle-node bifurcation curves unfolding the singularity. The two parameters $(\mathcal{R}_0, p)$ are written in a useful way to evaluate the proportion of vaccinated individuals necessary to eliminate the disease and to conclude how the vaccination may affect the outcome of the epidemic. We also exhibit the region in the parameter space where the disease persists and we illustrate our main result with numerical simulations, emphasizing the role of the parameters.


Introduction
Mathematical models applied to epidemiology play an important role to understand the dynamics of infectious diseases.Understanding how a disease evolves in a population can be beneficial, not only to predict epidemic outbreaks, but also to improve approaches against its spread [1,2,3,4,5].One of the most widely used models to study population interactions is the SIR model based on the division of the population in three classes of individuals: Susceptible (S), Infectious (I) and Recovered (R) [6,7].
When a disease enters a population, the medical community aims (as quick as possible) to find ways to stop its evolution.There are studies in the literature that investigate the behaviour of populations in the presence of infectious diseases by using a wide range of techniques [8,9,10,11,12].One of the most effective ways to prevent the disease's progression is via a vaccination policy [13, 14,15,16,17].As far as we know, there are few articles in the literature that combine logistic growth in the Susceptible population and the effect of vaccination.
In mathematical modelling, one may identify three vaccination strategies: (i) Constant vaccination [15,18,19] -it consists of vaccinating a prescribed ratio of newborns; (ii) Pulse vaccination [18,20] -it implies vaccinating a percentage of the susceptible population periodically; (iii) Mixed vaccination [18] -combination of constant and pulse vaccination.
Finding the most appropriate strategy is a challenge because we need to combine the effect of several factors, namely the efficiency of the vaccination and its cost to the public health policies.The application of bifurcation analysis to epidemiology may give clues about the evolution of a given disease in the presence of several factors namely vaccination, seasonality, sub-optimal immunity and nonlinear incidence [21,22,23,24].In the present article, we are interested in a modified SIR model exhibiting a Doublezero (DZ) singularity and the dynamics it may unfold [25].Sensitivity analysis may be particularly useful in this context.

State of art.
There are several studies involving singularities in epidemic models 1 .
In Jin et al. [26], the authors studied global dynamics of a SIRS model with nonlinear incidence rate.They established a threshold for a disease to be extinct or endemic, analyzed the existence and stability of equilibria and verified the existence of bistable states.Using normal forms and the Dulac criterion, they investigated backward bifurcation and obtained all curves associated to the Bogdanov-Takens bifurcation.
Zhang and Qiao [27] analyzed bifurcations in a SIR model including bilinear incidence rate, vaccination and hospital resources (number of available beds).They exhibited conditions to ensure the existence of a codimension three singularity for the system.They concluded that a high value of the vaccination parameter allows the disappearance of the disease in the population.See also [28] where the authors considered a SIR model with a standard incidence rate and a nonlinear recovery rate, formulated to consider the impact of available resources of the public health system (essentially the number of hospital beds).
Alexander and Moghadas [29] studied a SIRS epidemic model with a generalized nonlinear incidence as a function of the number of infected individuals.It is assumed that the natural immunity acquired by the infection is not permanent but wanes with time.Normal forms have been derived for the different types of bifurcation that the model undergoes.The Bogdanov-Takens normal form has been used to formulate the local bifurcation curves for a family of homoclinic orbits arising when a Hopf and a saddle-node bifurcation merge.They provided conditions for the occurrence of Hopf bifurcations in terms of two parameters: the basic reproductive number and the rate of loss of immunity acquired by the infection.
In 2022, Pan et al. [30] proposed a SIRS model undergoing a degenerate codimension three bifurcation inducing intermittency into the dynamics.The authors provided sufficient conditions to guarantee the global asymptotic stability of the unique endemic equilibrium.Vaccination and its boosted version have not been taken into account.
In 2018, Li and Teng [31] presented a SIRS model with generalized non-monotone incidence rate and qualitatively proved the existence of Bogdanov-Takens bifurcations.They located regions where the disease either persists or disappears.This phenomenon indicates that the initial conditions of an epidemic may determine the final states of an epidemic to go extinct or not.We also address the reader to [32,33] for other variations of the SIR/SIRS models.The reference [33] includes numerical simulations and data-fitting of the influenza data in China.
Novelty.Our work contributes to the mathematical understanding of the modified SIR model, where the class of Susceptible individuals (subject to a constant vaccination) grows logistically.We find a DZ singularity in the bifurcation space (R 0 , p), where R 0 is the basic reproduction number and p is the proportion of Susceptible individuals successfully vaccinated at birth.Writing the parameters (of the model) as function of R 0 and p is our first breakthrough.
We exhibit explicit expressions for the saddle-node, transcritical, Hopf and heteroclinic bifurcation curves associated to the DZ bifurcation.These unfolding curves have the same qualitative properties to the truncated amplitude system associated to the Hopf-zero normal form (8.81) of [25].The heteroclinic cycle is associated to two diseasefree equilibria and is asymptotically unstable (repelling).
For the sake of completeness, in Section 2, we describe some terminology that are going to be useful throughout this article.

Preliminaries
In this section, we introduce some terminology for vector fields acting on a R n , n ∈ N, that will be used in the remaining sections.Let f be a smooth vector field on R n with flow given by the unique solution x(t) = ϕ(t, x) ∈ R n of the two-parameter family where (η 1 , η 2 ) ∈ R 2 .
Definition 1.We say that K ⊂ (R + 0 ) 2 is a positively flow-invariant set for (1) if for all x ∈ K the trajectory of ϕ(t, x) is contained in K for t ≥ 0. Definition 2. Given two hyperbolic equilibria A and B of (1), a heteroclinic connection from A to B, is a solution of (1) contained in W u (A) ∩ W s (B), the intersection of the unstable manifold of A and the stable manifold of B.
For a solution of (1) passing through x ∈ R n , the set of its accumulation points, as t goes to +∞, is the ω-limit set of x.More formally, if A is the topological closure of A ⊂ R n , then: It is well known that ω(x) is closed and flow-invariant, and if the ϕ-trajectory of x is contained in a compact set, then ω(x) is non-empty.If E is an invariant set of (1), we say that E is a global attractor if ω(x) ⊂ E, for Lebesgue almost all points x in R n .
The center manifold of a non-hyperbolic equilibrium is the set of solutions whose behaviour around the equilibrium point is not controlled neither by the attraction of the stable manifold nor by the repulsion of the unstable manifold.If the linearized part of Df (η 1 ,η 2 ) (at a given equilibrium) has an eigenvalue with zero real part, the center manifold plays an important goal and it is the right set where bifurcations occur.
Throughout this article, we study the DZ singularity of codimension two for a family of differential equations corresponding to the case s = 1, θ < 0 in Equation (8.81) of [25].The unfolding of this singularity involves lines of saddle-node, Belyakov transition, Hopf and homo/heteroclinic bifurcations.We suggest the reading of [25] for a complete understanding of these bifurcations as well as the sufficient conditions that prompt their existence.

The model
Inspired by the classical SIR model, we are going to divide the individuals of a given (human) population into three classes of individuals [7,34]: • Susceptible (S): proportion of healthy individuals who are susceptible to the disease; • Infectious (I): proportion of infected individuals who can transmit the disease to susceptible individuals; • Recovered (R): proportion of individuals who recovered naturally from the disease or through immunity conferred by the vaccine.This comprises individuals who have definitive immunity and can not transmit the disease.
We assume that Susceptible individuals have never been in contact with the disease, but may become infected when they are in contact with the population of the Infectious, and then become part of this class.Whereas in the class of Infectious, these individuals can recover naturally and become part of the class of Recovered ones.The Susceptible individuals may also have been successfully vaccinated at birth, thus becoming immune to the disease [18,35,36].Inspired by [7,18,37], the nonlinear system of ODE in variables S, I and R (depending on time t), is given by: where Remark 1.When S = 0, we may extend non-smoothly F to the vector field The parameters of (2) may be interpreted as follows: A: carrying capacity of susceptible individuals when β = 0 (in the absence of disease) and p = 0 (in the absence of vaccination); β: transmission rate of the disease; m: birth rate; p: proportion of susceptible individuals successfully vaccinated at birth, for p ∈ [0, 1]; µ: natural death rate of infected and recovered individuals; d: death rate of infected individuals due to the disease; g: natural recovery rate.
Figure 1 illustrates the interaction between the classes of Susceptible, Infectious and Recovered individuals in model (2).The basic reproduction number, denoted by R 0 , may be seen as the number of secondary infections caused by a single infected person in a susceptible population [38].For model (2) with p = 0, R 0 may be explicitly computed as [9,39,40]: (3) 3.1.Hypotheses and motivation.With respect to system (2), we also assume the following conditions: (C1) All parameters are positive; (in other words, for R 0 < 1 no preventive measures involving vaccination will be considered).
The phase space associated to (2) is a subset of (R + 0 ) 3 induced with the usual topology, and the set of parameters is: 6  and p ∈ R + 0 .System (2) has been inspired in the classical SIR model [6] with slight modifications as we proceed to explain:

S
• We consider logistic growth in the susceptible population due to the crowding and natural competition for resources [9,34,37], instead of linear or exponential growth; • Instead of a mixed or pulse vaccination strategy [18, Subsection 2.1], we have assumed constant vaccination [15,18,20].
3.2.Two-dimensional system.The first two equations of (2), Ṡ and İ, are independent of R. Therefore, we may consider the system: with x = (S, I) ∈ (R + ) 2 .Because of Remark 1, we may extend (4) to the line S = 0: The vector field associated to ( 4) and ( 5) will be denoted by f and its flow is . This model will be the object of study of the present article; in order to shorten the notation, we denote by σ the sum µ + d.Lemma 1.The region defined by: is positively flow-invariant for (4) and (5).
Proof.It is easy to check that (R + 0 ) 2 is flow invariant (note that (5) leaves invariant the vertical lines S = 0 and I = 0).Now, we show that if (S 0 , I 0 ) ∈ M, then ϕ 0 (t, (S 0 , I 0 )), t ∈ R + 0 , is contained in M. Let us define φ(t) = S(t) + I(t) associated to the trajectory ϕ(t, (S 0 , I 0 )).Omitting the dependence of the variables on t (when there is no risk of ambiguity), one knows that: from where we deduce that: If β = p = 0, then the first component of ( 4) would be the logistic growth and thus its solution is limited by A (see (C2)), a property which remains for β, p > 0. In particular, we may conclude that φ + (σ + g)φ ≤ (σ + g + A)A.
The classical differential version of the Gronwall's Lemma 2 says that for all t ∈ R + 0 , we have: Taking the limit when t → +∞, we get: ), the result follows.

Main result and consequences
We state the main results of the article, as well as its structure.We also discuss some consequences.
4.1.Main result.System (4) may have three formal equilibria 3 : two disease-free equilibria and one endemic equilibrium, when they exist.The disease-free equilibria of system (4) are: The endemic formal equilibrium of ( 4) is: where We are able to study the map f as a two-parameter family depending on the basic reproduction number R 0 and the proportion of vaccination p.Our main result says that, in the bifurcation parameter (R 0 , p), E p 2 is a DZ singularity for the vector field a (e at − 1). 3 The term "formal equilibria" means that they are equilibria of the system regardless it makes sense or not in the context of the epidemic problem.
Theorem A. The endemic equilibrium E p 2 of (4) undergoes a DZ bifurcation at (R 0 , p ) = 2, A 2 4m .The local representations of the bifurcation curves in the space of parameters (R 0 , p) ∈ (R + 0 ) 2 are as follows: (i) Saddle-node bifurcation curve: (ii) Transcritical bifurcation curve: (iii) Belyakov transition curve: (iv) Heteroclinic cycle bifurcation curve: The Belyakov transition is a curve in the bifurcation space where the eigenvalues of Df (R 0 ,p) at E p 2 change from non-real to real or vice-versa.4.2.Proof of Theorem A and the structure of the article.The jacobian matrix of the vector field associated to (4) at E p 2 when (R 0 , p ) = 2, A 2 4m is given by: It is easy to verify that the matrix ( 6) is non-hyperbolic and has a double zero eigenvalue.This is why we say that, for (R 0 , p ) = 2, A 2 4m , the point E p 2 is a singularity of codimension 2. The associated DZ bifurcations exist if the vector field f (R 0 ,p) , at the bifurcation point, satisfies the nondegeneracy conditions described in [25, pp. 316-322].Instead of verifying these additional conditions, we check analytically the existence of all bifurcation curves that pass through the point (R 0 , p ) = 2, A 2 4m , as pointed out in Table 1.
For the sake of completeness, we have added in Section 5 a study of model ( 4) for p = 0.The curve representing the heteroclinic cycle, denoted by H, in the space of parameters (R 0 , p) was obtained by interpolation.The estimated correlation between the two parameters where one observes a repelling heteroclinic cycle is 1 (precision: 10 −5 ), as the reader may check in Section 6.6 and Appendix A. We simulate the dynamics in all hyperbolic regions associated to the DZ singularity in Figure 5. Section 8 finishes this article.

Biological consequences.
As a direct consequence of Theorem A, we are able to locate three regions in the parameter space (R 0 , p) where the disease persists (C, D and E), i.e. there is a set with positive Lebesgue measure (in the phase space) whose ω-limit is the endemic equilibrium.In regions C and D, trajectories starting "below" W s (E p 0 ) tend to E p 2 .In region E, there is a repelling periodic solution C (arising from the Hopf bifurcation) which is responsible for the maintenance of an endemic region where the disease persists.
In region E, if we assume seasonality in the transmission rate β, the associated flow exhibits a strange repeller and hyperbolic horseshoes, as a consequence of the works [9,41].The heteroclinic cycle H is a repeller.In the regions F, G and H, the high vaccination does not play a major role in the elimination of the disease because the carrying capacity A is small compared with the number of Infected individuals that are being generated.

The case without vaccination (p = 0)
For the sake of completeness, we analyze model (4) considering p = 0: The disease-free equilibria of system (7) are obtained by imposing I = 0. Then we get: With respect to the endemic equilibrium, we know that: (3) The formal equilibrium E 2 ≡ E 0 2 lies in the interior of the first quadrant when R 0 > 1.
The jacobian matrix of the vector field associated to (7) at a general point E = (S, I) ∈ (R + 0 ) 2 is given by: Evaluating J(E) at E 0 , E 1 and E 2 we have: and Lemma 2. System (7) exhibits: (1) two disease-free equilibria, E 0 and E 1 , such that: (2) an endemic equilibrium E 2 such that: Proof.
(1) The eigenvalues of J(E 0 ) are − (σ + g) < 0 and A > 0. Since the eigenvalues of J(E 0 ) have real part with opposite signs, then E 0 is a saddle for all R 0 ∈ R + .The eigenvalues of J(E 1 ) are Aβ − (σ + g) and −A < 0. If Therefore, in the vertical direction, E 1 interchanges its stability with E 2 from stable to unstable at R 0 = 1 (transcritical bifurcation).
(2) The eigenvalues of J(E 2 ) are given by: and , . Since β > 0 and σ +g > 0, it is easy to verify that λ 1 has negative real part.The eigenvalue λ 2 has negative real part if and only if , then λ 2 < 0 and E 2 is a stable node.In an analogous way, if R 0 > 1 + 1 4β , then ∆ < 0 and E 2 is a stable focus.Hence, E 2 evolves from stable node to stable focus at R 0 = 1 + 1 4β (Belyakov transition).
In Figure 3, one observes the scheme of the equilibria stability of system (7) for different values of R 0 .The transcritical bifurcation at R 0 = 1 represents the threshold for the existence of the endemic disease in the population, agreeing well with the empirical belief.Besides the E 0 and E 1 , there exists a stable node E 2 .III.Besides the E 0 and E 1 , there exists a stable focus E 2 .From I. to II.E 1 undergoes a transcritical bifurcation at R 0 = 1.From II. to III.E 2 evolves a Belyakov transition at R 0 = 1 + 1 4β .

The case with constant vaccination (p > 0)
We analyze (4) by considering p > 0 and R 0 > 1 (see (C3)): For A and m fixed, we settle the following maps that will be used throughout this text: for R 0 > 1.It is easy to verify that: • p SN is a constant map and , which does not make sense (because E p 2 does not exist in the first quadrant as a consequence of Lemmas 3 and 5).This is the reason why we consider the set [2, +∞[ as the domain of p H . 6.1.Preparatory section.We state a preliminary result that will be used in the sequel.
Lemma 3. The following statements hold for the maps given in (10):

Proof.
(1) From ( 10), it follows that: (2) This item follows if we evaluate p H , p T and p SN at R 0 = 2.
(3) From ( 10), we have: The remaining inequality follows from the previous item and transitivity.
Lemma 4. The following statements hold for model ( 9): (1) If p > p SN (R 0 ), then E p 0 and E p 1 do not exist in the first quadrant of (S, I).
(2) If p T (R 0 ) < p < p SN (R 0 ), then: (a) E p 0 is a saddle and E p 1 is a sink for 1 < R 0 < 2; (b) E p 0 is a source and E p 1 is a saddle for R 0 > 2.
(3) If p < p T (R 0 ) and R 0 > 1, then E p 0 and E p 1 are saddles.
Proof.The eigenvalues of J(E p 0 ) are A − 2S p 0 and βS p 0 − (σ + g).We know that , then E p 0 and E p 1 do not exist.
(2) We may write: then the eigenvalues of J(E p 0 ) have real part with opposite signs and E p 0 is a saddle.On the other hand, if then the eigenvalues of J(E p 0 ) have real part with positive signs, and E p 0 is a source.Moreover, the eigenvalues of J(E p 1 ) are A − 2S p 1 < 0 and βS p 1 − (σ + g).Hence, if then both eigenvalues of J(E p 1 ) have negative real part and E p 1 is a sink.If then the eigenvalues of J(E p 1 ) have real part with opposite signs and E p 1 is a saddle.Therefore, from ( 14) and ( 16), we conclude that: 0 is a source and E p 1 is a saddle, for R 0 > 2, under conditions ( 15) and ( 17).

Transcritical bifurcation.
Lemma 5.The endemic equilibrium E p 2 undergoes a transcritical bifurcation at p = p T (R 0 ) and lies in the interior of the first quadrant if p < p T (R 0 ).
Proof.The equilibrium E p 2 lies in the interior of the first quadrant if both S p 2 and I p 2 are positive.It is clear that S p 2 = σ + g β > 0. Since β 2 (σ + g) > 0, we have Solving the previous inequality in order to p, one gets: Hence, E p 2 undergoes a trancritical bifurcation along the line p = p T (R 0 ) and interchanges its stability with 6.4.Belyakov transition.We settle the following functions that will be used throughout this text: for Bt (R 0 ) and p Bt (R 0 ) are the square roots of ∆ 2 (p) written as function of R 0 .It is easy to verify that p (1) Bt (R 0 ) > 0. Proof.From (20) we know that if Bt (R 0 ) > 0. Before we analyze the stability of E p 2 , we remind the readers that E p 2 exists in the first quadrant when p < p T (R 0 ) (by Lemma 5).
Hence, E p 2 undergoes a Belyakov transition at p = p Bt (R 0 ).
(1) Therefore, if then λ 1 > 0 ⇒ λ 2 > 0 and E p 2 is an unstable node.Analogously, we proceed in the same way for λ 2 < 0. If Bt (R 0 ) (2) Now, let us assume that ∆ 2 (p) < 0, which means that p < p Bt (R 0 ), the eigenvalues are complex and E p 2 is a focus.If the real part of the eigenvalues of J(E p 2 ) is negative (and β is positive), i.e. if then E p 2 is a stable focus.Otherwise E p 2 is an unstable focus.Also, we conclude that E p 2 undergoes a Belyakov transition at p = p Bt (R 0 ), since E p 2 changes its stability from a node to a focus.6.5.Hopf bifurcation.In this section we will find a line in the parameter space (p, R 0 ) where a Hopf bifurcation occurs.  2 ) given by From [42, pp.151-153, Theorem 3.4.2(adapted)] we know that the Hopf bifurcation occurs when the following conditions hold: (1) the eigenvalues of the map Df (E p 2 ) (see ( 13)) have the form ±iω (ω > 0), which is equivalent to Tr J(E p 2 ) = 0, where Tr represents the usual trace operator.Indeed, (2) d dR 0 (Re λ j (R 0 )) R 0 ∈p H = 0, for j ∈ {1, 2}: Hence, for p = p H (R 0 ) we get Hence, E p 2 undergoes a subcritical Hopf bifurcation (as R 0 increases; the periodic solution appears for R 0 < R 0 as illustrated in Figure 4. Figure 2 provides a better perspective of all bifurcations under consideration.6.6.Heteroclinic cycle.Finding explicitly a non-robust heteroclinic cycle (in the phase space) to two hyperbolic saddles and its bifurcating curve is a difficult task.In this section, we use polynomial interpolation (of rational degree) to find the latter curve.We proceed to explain the method that we have used to compute the map.All steps are described in Appendix A.
Using MATLAB_R2018a, we locate thirteen pairs (R 0i , p i ) in the parameter space (R 0 , p) where one observes a heteroclinic cycle associated to the disease-free equilibria E p i 0 and E p i 1 (see Table 2).The type of function for the interpolation should be appropriately chosen to fit the curve to our finite set of points [43, Part I, Section 1].
As suggested by [23], points in the parameter space that correspond to the heteroclinic cycle associated to the disease-free equilibria should lie on the graph of: The heteroclinic bifurcation organizes the dynamics, stressing a geometric configuration in its unfolding.

Numerics
In this section we describe the dynamics of the bifurcations of model (9) and we perform numerical simulations associated to parameters in all hyperbolic (open) regions of Figures 2 and 5.
A: The vaccination coverage is so large that Infectious individuals tend to disappear.

B:
The flow exhibits two disease-free equilibria, a saddle E p 0 and a sink E p 1 .All trajectories starting in the first quadrant evolve in such a way that their I component tend to disappear.No endemic equilibria exist.

C and D:
The endemic equilibrium E p 2 lies in the interior of the first quadrant and it is Lyapunov stable.There is a positive Lebesgue measure set of initial conditions converging to E p 2 , so the disease remains persistently in the population.Initial conditions lying "below" W s (E p 0 ) have E p 2 as ω-limit.

E:
The unstable heteroclinic cycle H is broken giving rise to an unstable periodic solution C. The region bounded by C is conducive to the persistence of the disease since all initial conditions tend toward E p 2 .All initial conditions lying in the unbounded region defined by C are propitious to the disappearance of the disease.The period of C is very large when (R 0 , p) is close to the curve Het.
F and G: The equilibrium E p 2 is unstable and all initial conditions are in a region conducive to the disappearance of the disease.

H:
No endemic equilibria exist.The equilibria E p 0 and E p 1 repel Lebesgue almost all trajectories towards the origin.

Concluding remarks
In this article, we have performed a bifurcation analysis of a modified SIR model accomodating constant vaccination and logistic growth in the Susceptible population, which may be seen as a contribution towards the study of strategies to control infectious diseases.
8.1.Conclusions.Model (2) exhibits an endemic Double-zero (DZ) singularity.To prove the main result, we have checked that the vector field associated to (2) has a double zero eigenvalue (at the endemic equilibrium E p 2 ) and we have described all associated bifurcation curves passing through the bifurcation point (R 0 , p ) = 2, A 2 4m in the parameter space.In Figure 2, we have exhibited explicitly the regions in the space of parameters where the disease persists in the population.
We have established a threshold for the disease to be extinct or endemic and we have analyzed the existence and asymptotic stability of equilibria.Writing the bifurcation curves as a function R 0 and p is useful since these variables are key parameters in many epidemiological studies.
The inclusion of a constant vaccination and a logistic function in the Susceptible population increases the complexity of the dynamics of the classical SIR model.The system may exhibit two disease-free equilibria and none of them corresponds to the equilibrium (0, 0).Similar dynamics occurs in the context of prey-predator models [44,45].
Although Df 2, A 2 4m at E p 2 has a double zero eigenvalue, the associated bifurcation curves do not correspond to those of the classic Bogdanov-Takens bifurcation.It is a variant of the truncated amplitude system associated to the Hopf-zero normal form (8.81) of [25].This implies chaos and suspended horseshoes in the presence of seasonal parameters near p H (R 0 ), as a consequence of the theory developed in [9,46].
In the presence of the endemic equilibrium (in the first quadrant of (S, I)), we have detected two different regimes for R 0 ≥ 1: • 1 < R 0 < 2: the increase of the vaccination rate plays an essencial role to decrease the number of Infectious; • 2 ≤ R 0 : although vaccination policies are beneficial, the control of the endemic disease is mostly due to the saturation of the class of Susceptible individuals.
With constant vaccination, the transcritical and the Hopf bifurcation curves induce a change in the stability of the endemic equilibrium point.At the epidemiological level, this means that these curves determine if the disease either disappears or remains.The saddle-node bifurcation curves indicate the moment where disease-free equilibria appear.
Although the model under analysis has limitations in terms of biologic validity, the existence of a DZ bifurcation organizes the dynamics and stresses the role played by the vaccination; the model agrees well with the empirical beliefs.8.2.Future work.As already referred, our work generalizes that of Shulgin et al. [18] (with constant vaccination) who modelled individual growth just with birth rate.
The strategy of constant vaccination may not be the most effective if the proportion of successfully vaccinated newborns is low.It may be necessary to adopt more effective strategies such as the pulse vaccination.This technique aims to find an optimal period between "shots" of the vaccine, allowing a proportion of infected individuals lower than a given threshold.It should be possible to apply this technique to the regions where the disease persists (C, D and E of Figure 2).We would like to derive the period of the vaccination "shot" that (efficiently) allows to control the number of infected individuals.
Pulse vaccination can have different effects on the dynamics of infectious diseases, depending on our starting model.Our next goal is to modify the first component of ( 2 where n ∈ N: • T is the period between two consecutive "shots" (of vaccination);    2) for which it is possible to observe the heteroclinic cycle in phase space (S, I).The function obtained by interpolation for which we can find the heteroclinic cycle in the space of parameters (R 0 , p) is p (R 0 ) ≈ 4.495R −2.313 0 − 0.039.

Figure 1 .
Figure 1.Schematic diagram of model (2).Boxes represent compartments, and arrows indicate the flow between S, I and R.

Figure 2 .
Figure 2. DZ bifurcation diagram associated to (4).In the regions C, D and E, the disease persists in the population.Compare with numerics of Figure 5.

Figure 3 .
Figure 3. Phase diagram of(7) for different values of R 0 .I. E 1 is a global attractor when restricted to the closure of the first quadrant.II.Besides the E 0 and E 1 , there exists a stable node E 2 .III.Besides the E 0 and E 1 , there exists a stable focus E 2 .From I. to II.E 1 undergoes a transcritical bifurcation at R 0 = 1.From II. to III.E 2 evolves a Belyakov transition at R 0 = 1 + 1 4β .

Lemma 8 .
If p = p H (R 0 ), then E p 2 undergoes a subcritical Hopf bifurcation associated to an unstable periodic solution C (as R 0 increases).

Figure 4 .
Figure 4. Phase portraits of system (9) when E p 2 undergoes a subcritical Hopf bifurcation as R 0 increases.

Figure 6 .
Figure 6.Function obtained by interpolation of 13 points (see Table2) for which it is possible to observe the heteroclinic cycle in phase space (S, I).The function obtained by interpolation for which we can find the heteroclinic cycle in the space of parameters (R 0 , p) is p (R 0 ) ≈ 4.495R−2.313

Figure 7 .
Figure 7. MATLAB_R2018a output for the two-term power series model used in (R 0i , p i ) points interpolation.

Table 1 .
Codimension 1 bifurcations that characterizes the DZ bifurcation and structure of the proof of Theorem A.

Table 2 .
13 interpolating points of the heteroclinic cycle function in the space of parameters (R 0 , p).These points were obtained in MATLAB_R2018a.