SIR Model with Vaccination: Bifurcation Analysis

Abstract

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.

1 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 Double-zero (DZ) singularity and the dynamics it may unfold [25]. Sensitivity analysis may be particularly useful in this context.

1.1 State of Art

There are several studies involving singularities in epidemic models.¹

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.

1.2 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 $(\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. Writing the parameters (of the model) as function of $\mathcal {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 disease-free equilibria and is asymptotically unstable (repelling).

For the sake of completeness, in Sect. 2, we describe some terminology that are going to be useful throughout this article.

2 Preliminaries

In this section, we introduce some terminology for vector fields acting on a $\mathbb {R}^n$, $n\in {\mathbb {N}}$, that will be used in the remaining sections. Let f be a smooth vector field on $\mathbb {R}^n$ with flow given by the unique solution $x(t)=\varphi (t, x)\in \mathbb {R}^n$ of the two-parameter family

$$\begin{aligned} \dot{x}=f_{(\eta _1, \eta _2)} (x), \qquad x(0)=x_0, \end{aligned}$$

(1)

where $(\eta _1, \eta _2)\in \mathbb {R}^2$.

Definition 1

We say that $\mathcal {K} \subset (\mathbb {R}^+_0)^2$ is a positively flow-invariant set for (1) if for all $x \in \mathcal {K}$ the trajectory of $\varphi (t, x)$ is contained in $\mathcal {K}$ for $t \ge 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)\cap 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\in \mathbb {R}^n$, the set of its accumulation points, as t goes to $+\infty$, is the $\omega$-limit set of x. More formally, if $\overline{A}$ is the topological closure of $A\subset \mathbb {R}^n$, then:

Definition 3

If $x\in \mathbb {R}^n$, the $\omega$-limit of x is:

$$\begin{aligned} \omega (x)=\bigcap _{T=0}^{+\infty } \overline{\left( \bigcup _{t>T}\varphi (t, x)\right) }. \end{aligned}$$

It is well known that $\omega (x)$ is closed and flow-invariant, and if the $\varphi$–trajectory of x is contained in a compact set, then $\omega (x)$ is non-empty. If E is an invariant set of (1), we say that E is a global attractor if $\omega (x) \subset E$, for Lebesgue almost all points x in $\mathbb {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_{(\eta _1, \eta _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$, $\theta <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.

3 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:

$$\begin{aligned} \begin{array}{lcl} \dot{X} = \mathcal {F} (X) \quad \Leftrightarrow \quad {\left\{ \begin{array}{ll} &{}\dot{S} = S(A-S) - \beta I S - pm \\ \\ &{}\dot{I} = \beta IS - (\mu +d) I - gI \\ \\ &{}\dot{R} = p m + gI - \mu R, \end{array}\right. } \end{array} \end{aligned}$$

(2)

where

$$\begin{aligned} \begin{array}{lcl} X(t) &{}=&{} \left( S(t), I(t), R(t) \right) \in (\mathbb {R}^+)^3, \\ \\ \dot{X} &{}=&{} (\dot{S}, \dot{I}, \dot{R}) \,\,\, = \,\,\, \displaystyle \left( \frac{\textrm{d}S}{\textrm{d}t},\frac{\textrm{d}I}{\textrm{d}t},\frac{\textrm{d}R}{\textrm{d}t}\right) . \\ \end{array} \end{aligned}$$

Remark 1

When $S=0$, we may extend non-smoothly $\mathcal {F}$ to the vector field

$$\begin{aligned}(0, - (\mu +d) I - gI, pm + gI - \mu R).\end{aligned}$$

The parameters of (2) may be interpreted as follows:

A:

carrying capacity of susceptible individuals when $\beta =0$ (in the absence of disease) and $p = 0$ (in the absence of vaccination);

$\beta$:

transmission rate of the disease;

m:

birth rate;

p:

proportion of susceptible individuals successfully vaccinated at birth, for $p \in [0,1]$;

$\mu$:

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 $\mathcal {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$, $\mathcal {R}_0$ may be explicitly computed as [9, 39, 40]:

$$\begin{aligned} \begin{array}{lcl} \mathcal {R}_0 =\displaystyle \lim _{T\rightarrow +\infty } \frac{1}{T} \int _0^T \dfrac{A \beta }{\mu +d+ g} \, \textrm{d}t = \dfrac{A \beta }{\mu +d+ g} > 0. \end{array} \end{aligned}$$

(3)

3.1 Hypotheses and Motivation

With respect to system (2), we also assume the following conditions:

(C1):

All parameters are positive;

(C2):

For all $t \in \mathbb {R}^+_0, S(t) \le A$;

(C3):

Vaccination is considered only when $\mathcal {R}_0 > 1$, i.e. $p>0$ if and only if $\mathcal {R}_0 > 1$ (in other words, for $\mathcal {R}_0 < 1$ no preventive measures involving vaccination will be considered).

The phase space associated to (2) is a subset of $(\mathbb {R}_0^+)^3$ induced with the usual topology, and the set of parameters is:

$$\begin{aligned} \Omega = \left\{ \omega = (A,\beta ,m,\mu ,d,g) \in (\mathbb {R}^+)^6 \right\} \quad \text {and} \quad p \in \mathbb {R}_0^+. \end{aligned}$$

Fig. 1

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

System (2) has been inspired in the classical SIR model [6] with slight modifications as we proceed to explain:

• 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), $\dot{S}$ and $\dot{I}$, are independent of R. Therefore, we may consider the system:

$$\begin{aligned} \left\{ \begin{array}{lcl} &{}\dot{S} = S(A-S) - \beta I S - pm \\ \\ &{}\dot{I} = \beta IS - \left( \mu +d \right) I - g I, \end{array}\right. \end{aligned}$$

(4)

with $x=(S,I)\in (\mathbb {R}^+)^2$. Because of Remark 1, we may extend (4) to the line $S=0$:

$$\begin{aligned} \left\{ \begin{array}{lcl} &{}\dot{S} = 0 \\ \\ &{}\dot{I} = - \left( \mu +d \right) I - g I. \end{array}\right. \end{aligned}$$

(5)

The vector field associated to (4) and (5) will be denoted by f and its flow is

$$\begin{aligned}\varphi (t, (S_0, I_0)), \quad t \in \mathbb {R}^+_0,\quad (S_0, I_0) \in (\mathbb {R}^+_0)^2.\end{aligned}$$

This model will be the object of study of the present article; in order to shorten the notation, we denote by $\sigma$ the sum $\mu + d$.

Lemma 1

The region defined by:

$$\begin{aligned} \mathcal {M} = \left\{ (S,I) \in (\mathbb {R}_0^+)^2: \quad 0 \le S \le A, \quad 0 \le S+I \le \dfrac{A (\sigma + g + A)}{\sigma + g}, \quad S, I\ge 0 \right\} , \end{aligned}$$

is positively flow-invariant for (4) and (5).

Proof

It is easy to check that $(\mathbb {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)\in \mathcal {M}$, then $\varphi _0(t, (S_0, I_0))$, $t \in \mathbb {R}_0^+$, is contained in $\mathcal {M}$. Let us define $\phi (t) = S(t) + I(t)$ associated to the trajectory $\varphi (t, (S_0, I_0))$. Omitting the dependence of the variables on t (when there is no risk of ambiguity), one knows that:

$$\begin{aligned} \begin{array}{lcl} \dot{\phi } &{} = &{} \dot{S} + \dot{I} \\ \\ &{} = &{} S(A-S) - \beta I S - pm + \beta I S - (\sigma + g) I \\ \\ &{} = &{} S(A-S) - pm - (\sigma + g)I,\\ \end{array} \end{aligned}$$

from where we deduce that:

$$\begin{aligned} \begin{array}{lcl} \dot{\phi } + (\sigma + g) \phi &{} = &{} S(A-S) - pm - (\sigma + g)I + (\sigma + g)S + (\sigma + g)I \\ \\ &{} = &{} S(A-S) + (\sigma + g)S \\ \\ &{} \le &{} (\sigma + g + A) S.\\ \end{array} \end{aligned}$$

If $\beta =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 $\beta ,p > 0$. In particular, we may conclude that

$$\begin{aligned} \begin{array}{lcl} \dot{\phi } + (\sigma + g) \phi \le (\sigma + g + A) A. \end{array} \end{aligned}$$

The classical differential version of the Gronwall’s Lemma² says that for all $t \in \mathbb {R}^+_0$, we have:

$$\begin{aligned} \phi (t)\le \phi (0) e^{-(\sigma +g) t} - \frac{(\sigma + g + A) A}{(\sigma +g)}\left( e^{-(\sigma +g) t}-1\right) . \end{aligned}$$

Taking the limit when $t\rightarrow +\infty$, we get:

$$\begin{aligned} 0\le \lim _{t\rightarrow +\infty } \phi (t)\le & {} \lim _{t\rightarrow +\infty } \left[ \phi (0) e^{-(\sigma +g) t} - \frac{(\sigma + g + A) A}{(\sigma +g)}\left( e^{-(\sigma +g) t}-1\right) \right] \\= & {} \frac{(\sigma + g + A) A}{(\sigma +g)}\\ \end{aligned}$$

Since $\displaystyle \lim \nolimits _{t\rightarrow +\infty } \phi (t)= \displaystyle \lim \nolimits _{t\rightarrow +\infty }\left( S(t) + I(t)\right)$, the result follows. $\square$

4 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³: two disease-free equilibria and one endemic equilibrium, when they exist. The disease-free equilibria of system (4) are:

$$\begin{aligned} E_0^p =(S_0^p, I_0^p)= \left( \dfrac{A-\sqrt{A^2 - 4pm}}{2},0 \right) \end{aligned}$$

and

$$\begin{aligned} E_1^p = (S_1^p,I_1^p)= \left( \dfrac{A+\sqrt{A^2 - 4pm}}{2},0 \right) , \end{aligned}$$

where $p \le \frac{A^2}{4m} \le 1$, which implies $A \le 2\sqrt{m}$. The endemic formal equilibrium of (4) is:

$$\begin{aligned} E_2^p =\left( S_2^p,I_2^p\right) = \left( \dfrac{\sigma + g}{\beta },\dfrac{-pm\beta ^2 + A\left( \sigma +g\right) \beta -\left( \sigma +g\right) ^2}{\beta ^2 \left( \sigma +g\right) }\right) , \end{aligned}$$

where $\displaystyle -pm\beta ^2 + A\left( \sigma +g\right) \beta -\left( \sigma +g\right) ^2> 0 \, \overset{(3)}{\Leftrightarrow }\,\mathcal {R}_0 > 1 + \frac{pm\beta ^2}{\left( \sigma + g \right) }$.

We are able to study the map f as a two-parameter family depending on the basic reproduction number $\mathcal {R}_0$ and the proportion of vaccination p. Our main result says that, in the bifurcation parameter $(\mathcal {R}_0, p)$, $E_2^p$ is a DZ singularity for the vector field $f \mapsto f_{(\mathcal {R}_0, p)}$.

Theorem A

The endemic equilibrium $E_2^p$ of (4) undergoes a DZ bifurcation at $(\mathcal {R}_0^{\star },p^\star ) = \left( 2, \frac{A^2}{4\,m}\right)$. The local representations of the bifurcation curves in the space of parameters $(\mathcal {R}_0,p) \in (\mathbb {R}_0^+)^2$ are as follows:

• (i) Saddle-node bifurcation curve:

  $$\begin{aligned} \varvec{SN} = \left\{ (\mathcal {R}_0,p) \in (\mathbb {R}_0^+)^2 \,: \, p = \frac{A^2}{4m} \right\} \end{aligned}$$

• (ii) Transcritical bifurcation curve:

  $$\begin{aligned} \varvec{T} = \left\{ (\mathcal {R}_0,p) \in (\mathbb {R}_0^+)^2 \,: \, p = \dfrac{A^2}{m} \dfrac{\left( \mathcal {R}_0-1\right) }{\mathcal {R}_0^2} \right\} \end{aligned}$$

• (iii) Belyakov transition curve:

  $$\begin{aligned} \varvec{B_t} = \left\{ (\mathcal {R}_0,p) \in (\mathbb {R}_0^+)^2 \,: \, p = \dfrac{\left( -2\beta +1+2\displaystyle \sqrt{\beta \left( \mathcal {R}_0 + \beta - 2\right) }\right) A^2}{m\mathcal {R}_0^2}\right\} \end{aligned}$$

• (iv) Heteroclinic cycle bifurcation curve:

  $$\begin{aligned} \mathbf{Het.} = \left\{ (\mathcal {R}_0,p) \in (\mathbb {R}_0^+)^2 \,: \, p \approx \dfrac{4.495}{\mathcal {R}_0^{2.313}} - 0.039 \right\} \end{aligned}$$

• (v) Hopf bifurcation curve:

  $$\begin{aligned} \varvec{H} = \left\{ (\mathcal {R}_0,p) \in (\mathbb {R}_0^+)^2 \,: \, p = \dfrac{A^2}{m\mathcal {R}_0^2} \right\} . \end{aligned}$$

The Belyakov transition is a curve in the bifurcation space where the eigenvalues of $Df_{(\mathcal {R}_0, p)}$ at $E_2^p$ change from non-real to real or vice-versa.

Fig. 2

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

4.2 Proof of Theorem A and the Structure of the Article

The jacobian matrix of the vector field associated to (4) at $E_2^p$ when $(\mathcal {R}_0^{\star },p^\star ) = \left( 2, \frac{A^2}{4\,m}\right)$ is given by:

$$\begin{aligned}{} & {} \left( \begin{array}{cc} \dfrac{p^\star m\beta ^2-\left( \sigma +g\right) ^2}{\beta \left( \sigma +g\right) } &{} - \left( \sigma +g\right) \\ \\ \dfrac{-p^\star m\beta ^2 + A\left( \sigma +g\right) \beta - \left( \sigma +g\right) ^2}{\beta \left( \sigma +g\right) } &{} 0 \end{array}\right) \nonumber \\{} & {} \qquad =\left( \begin{array}{cc} \dfrac{p^\star m}{A}\mathcal {R}_0^{\star } - \dfrac{A}{\mathcal {R}_0^{\star }} &{} - \left( \sigma +g\right) \\ \\ -\dfrac{p^\star m}{A}\mathcal {R}_0^{\star } + A - \dfrac{A}{\mathcal {R}_0^{\star }} &{} 0 \end{array}\right) \nonumber \\{} & {} \qquad =\left( \begin{array}{cc} 0 &{} - \left( \sigma +g\right) \\ \\ 0 &{} 0 \end{array}\right) . \end{aligned}$$

(6)

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 $(\mathcal {R}_0^{\star },p^\star ) = \left( 2, \frac{A^2}{4\,m}\right)$, the point $E_2^p$ is a singularity of codimension 2. The associated DZ bifurcations exist if the vector field $f_{(\mathcal {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 $(\mathcal {R}_0^{\star },p^\star ) =\left( 2, \frac{A^2}{4\,m}\right)$, as pointed out in Table 1.

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

For the sake of completeness, we have added in Sect. 5 a study of model (4) for $p=0$. The curve representing the heteroclinic cycle, denoted by $\mathcal {H}$, in the space of parameters $(\mathcal {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 Sect. 6.6 and “Appendix 1”. We simulate the dynamics in all hyperbolic regions associated to the DZ singularity in Fig. 5. Section 8 finishes this article.

4.3 Biological Consequences

As a direct consequence of Theorem A, we are able to locate three regions in the parameter space $(\mathcal {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 $\omega$-limit is the endemic equilibrium. In regions C and D, trajectories starting “below” $W^s(E_0^p)$ tend to $E_2^p$. In region E, there is a repelling periodic solution $\mathcal {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 $\beta$, the associated flow exhibits a strange repeller and hyperbolic horseshoes, as a consequence of the works [9, 41]. The heteroclinic cycle $\mathcal {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.

5 The Case Without Vaccination $(p=0)$

For the sake of completeness, we analyze model (4) considering $p = 0$:

$$\begin{aligned} \dot{x} = f_{(\mathcal {R}_0, 0)}(x) \quad \Leftrightarrow \quad \left\{ \begin{array}{ll} &{}\dot{S} = S(A-S) - \beta I S \\ \\ &{}\dot{I} = \beta IS - \left( \sigma +g\right) I. \end{array}\right. \end{aligned}$$

(7)

The disease-free equilibria of system (7) are obtained by imposing $I=0$. Then we get:

$$\begin{aligned} E_0 \equiv E_0^0= \left( 0,0 \right) \quad \text {and} \quad E_1 \equiv E_1^0= \left( A,0 \right) . \end{aligned}$$

With respect to the endemic equilibrium, we know that:

$$\begin{aligned} E_2 = \left( \dfrac{\sigma + g}{\beta },\dfrac{A\beta -\left( \sigma +g\right) }{\beta ^2}\right) \overset{(3)}{=} \left( \dfrac{A}{\mathcal {R}_0}, \dfrac{A}{\beta } \left( 1-\dfrac{1}{\mathcal {R}_0} \right) \right) . \end{aligned}$$

The formal equilibrium $E_2\equiv E_2^0$ lies in the interior of the first quadrant when $\mathcal {R}_0 > 1$.

The jacobian matrix of the vector field associated to (7) at a general point $E = (S, I) \in (\mathbb {R}_0^+)^2$ is given by:

$$\begin{aligned} \begin{array}{lcl} J(E)=\left( \begin{array}{cc} -\beta I + A - 2S &{} - \beta S \\ \\ \beta I &{} \beta S - \left( \sigma +g\right) \end{array}\right) \end{array}. \end{aligned}$$

(8)

Evaluating J(E) at $E_0$, $E_1$ and $E_2$ we have:

$$\begin{aligned} J(E_0)=\left( \begin{array}{cc} A &{} 0 \\ \\ 0 &{} - \left( \sigma +g\right) \end{array}\right) \qquad J(E_1)=\left( \begin{array}{cc} -A &{} -A\beta \\ \\ 0 &{} A\beta - \left( \sigma +g\right) \end{array}\right) . \end{aligned}$$

and

$$\begin{aligned} J(E_2)=\left( \begin{array}{cc} -\dfrac{\sigma +g}{\beta } &{} - \left( \sigma +g\right) \\ \\ \dfrac{A\beta - \left( \sigma +g\right) }{\beta } &{} 0 \end{array}\right) , \end{aligned}$$

respectively.

Lemma 2

System (7) exhibits:

• (1) two disease-free equilibria, $E_0$ and  $E_1$, such that:

  • (a) for all $\mathcal {R}_0 \in \mathbb {R}^+$, $E_0$ is a saddle;

  • (b) if  $\mathcal {R}_0 < 1$, then  $E_1$ is a sink. If  $\mathcal {R}_0 > 1$, then  $E_1$ is a saddle;

  • (c) at  $\mathcal {R}_0 = 1$,  $E_1$ undergoes a transcritical bifurcation with $E_2$.

• (2) an endemic equilibrium $E_2$ such that:

  • (a) if  $1<\mathcal {R}_0 < 1+\frac{1}{4\beta }$, then  $E_2$ is a stable node;

  • (b) if  $\mathcal {R}_0 > 1 + \frac{1}{4\beta }$, then  $E_2$ is a stable focus;

  • (c) at  $\mathcal {R}_0 = 1+\frac{1}{4\beta }$,  $E_2$ undergoes a Belyakov transition.

Proof

• (1) The eigenvalues of $J(E_0)$ are  $-\left( \sigma +g\right) < 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 $\mathcal {R}_0 \in \mathbb {R}^+$. The eigenvalues of $J(E_1)$ are $A \beta - \left( \sigma +g\right)$ and $-A<0$. If

  $$\begin{aligned} A\beta - \left( \sigma +g\right)< 0 \quad \Leftrightarrow \quad \dfrac{A\beta }{\sigma + g}< 1 \quad \overset{(3)}{\Leftrightarrow } \quad \mathcal {R}_0 <1, \end{aligned}$$

  then $E_1$ is a sink. If  $\mathcal {R}_0 > 1$, then $E_1$ is a saddle. Therefore, in the vertical direction, $E_1$ interchanges its stability with $E_2$ from stable to unstable at $\mathcal {R}_0 = 1$ (transcritical bifurcation).

• (2) The eigenvalues of $J(E_2)$ are given by:

  $$\begin{aligned}\lambda _1 :=\dfrac{- \left( \sigma +g\right) - \sqrt{\Delta }}{2\beta } \qquad \text {and} \qquad \lambda _2 :=\dfrac{- \left( \sigma +g\right) + \sqrt{\Delta }}{2\beta },\end{aligned}$$

  where $\Delta = -4\left( \sigma +g\right) \left[ \left( -\beta - \dfrac{1}{4} \right) g + A\beta ^2 - \sigma \beta - \dfrac{\sigma }{4} \right]$. Since $\beta >0$ and $\sigma +g>0$, it is easy to verify that $\lambda _1$ has negative real part. The eigenvalue $\lambda _2$ has negative real part if and only if

  $$\begin{aligned}{} & {} -\left( \sigma + g \right) + \sqrt{\Delta }< 0 \\{} & {} \quad \Leftrightarrow \sqrt{\Delta }< \sigma + g \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{}\Delta \ge 0 \\ \\ &{}\Delta< \left( \sigma + g \right) ^2 \end{array}\right. } \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{} -4\left( \sigma +g\right) \left[ \left( -\beta - \dfrac{1}{4} \right) g + A\beta ^2 - \sigma \beta - \dfrac{\sigma }{4} \right] \ge 0 \\ \\ &{}-4\left( \sigma +g\right) \left[ \left( -\beta - \dfrac{1}{4} \right) g + A\beta ^2 - \sigma \beta - \dfrac{\sigma }{4} \right]< \left( \sigma + g \right) ^2 \end{array}\right. } \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{} \left( -\beta - \dfrac{1}{4} \right) g + A\beta ^2 - \sigma \beta - \dfrac{\sigma }{4} \le 0 \\ \\ &{} -4\left[ \left( -\beta - \dfrac{1}{4} \right) g + A\beta ^2 - \sigma \beta - \dfrac{\sigma }{4} \right]< \sigma + g \end{array}\right. } \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{} 4A\beta ^2 - \left( 4 \sigma \beta + \sigma \right) - \left( 4\beta + 1 \right) g \le 0 \\ \\ &{} g - A\beta +\sigma < 0 \end{array}\right. } \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{} 4\beta \mathcal {R}_0 - \dfrac{\sigma \left( 4 \beta + 1 \right) }{\sigma + g} - \dfrac{g\left( 4 \beta + 1 \right) }{\sigma + g} \le 0 \\ \\ &{} \mathcal {R}_0> 1 \end{array}\right. } \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{} \mathcal {R}_0 \le 1 + \dfrac{1}{4\beta } \\ \\ &{} \mathcal {R}_0 > 1. \end{array}\right. } \end{aligned}$$

  If $1< \mathcal {R}_0 < 1+\frac{1}{4\beta }$, then $\lambda _2 < 0$ and $E_2$ is a stable node. In an analogous way, if $\mathcal {R}_0 > 1 + \frac{1}{4\beta }$, then $\Delta < 0$ and $E_2$ is a stable focus. Hence, $E_2$ evolves from stable node to stable focus at $\mathcal {R}_0 = 1+\frac{1}{4\beta }$ (Belyakov transition).

$\square$

In Fig. 3, one observes the scheme of the equilibria stability of system (7) for different values of $\mathcal {R}_0$. The transcritical bifurcation at $\mathcal {R}_0 = 1$ represents the threshold for the existence of the endemic disease in the population, agreeing well with the empirical belief.

Fig. 3

Phase diagram of (7) for different values of $\mathcal {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 $\mathcal {R}_0 = 1$. From II. to III. $E_2$ evolves a Belyakov transition at $\mathcal {R}_0 = 1+\frac{1}{4\beta }$

6 The Case with Constant Vaccination $(p>0)$

We analyze (4) by considering $p > 0$ and $\mathcal {R}_0 > 1$ (see (C3)):

$$\begin{aligned} \dot{x} = f_{(\mathcal {R}_0, 0)}(x) \quad \Leftrightarrow \quad \left\{ \begin{array}{lll} &{}\dot{S} = S(A-S) - \beta I S - pm \\ \\ &{}\dot{I} = \beta IS - \left( \sigma +g\right) I. \end{array}\right. \end{aligned}$$

(9)

For A and m fixed, we settle the following maps that will be used throughout this text:

$$\begin{aligned} p_{SN}(\mathcal {R}_0) :=\dfrac{A^2}{4m} \quad \text {,} \quad p_T(\mathcal {R}_0) :=\dfrac{A^2}{m} \dfrac{\left( \mathcal {R}_0-1\right) }{\mathcal {R}_0^2} \quad \text {and} \quad p_H(\mathcal {R}_0) :=\dfrac{A^2}{m\mathcal {R}_0^2}, \end{aligned}$$

(10)

for $\mathcal {R}_0 > 1$. It is easy to verify that:

• $p_{SN}$ is a constant map and

• if $1< \mathcal {R}_0 < 2$, then $p_H (\mathcal {R}_0)>p_{SN}(\mathcal {R}_0)$, which does not make sense (because $E_2^p$ does not exist in the first quadrant as a consequence of Lemmas 3, 5). This is the reason why we consider the set $[2, +\infty [$ 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):

• (1) If  $\mathcal {R}_0 > 1$ and  $\mathcal {R}_0 \ne 2$, then  $p_T (\mathcal {R}_0)< p_{SN} (\mathcal {R}_0)$;

• (2)  $p_H (2)= p_T(2) = p_{SN}(2)$;

• (3) If  $\mathcal {R}_0 > 2$, then  $p_H (\mathcal {R}_0)< p_T(\mathcal {R}_0)<p_{SN}(\mathcal {R}_0)$.

Proof

• (1) From (10), it follows that:

  $$\begin{aligned}{} & {} p_T(\mathcal {R}_0)< p_{SN} (\mathcal {R}_0) \nonumber \\{} & {} \quad \Leftrightarrow \dfrac{A^2}{m} \dfrac{\left( \mathcal {R}_0-1\right) }{\mathcal {R}_0^2}< \dfrac{A^2}{4m}\nonumber \\{} & {} \quad \Leftrightarrow \dfrac{\mathcal {R}_0-1}{\mathcal {R}_0^2} < \dfrac{1}{4} \nonumber \\{} & {} \quad \Leftrightarrow \mathcal {R}_0^2 - 4\mathcal {R}_0 + 4> 0\nonumber \\{} & {} \quad \Leftrightarrow (\mathcal {R}_0-2)^2 > 0. \end{aligned}$$

  (11)

• (2) This item follows if we evaluate $p_H$, $p_T$ and $p_{SN}$ at $\mathcal {R}_0 = 2$.

• (3) From (10), we have:

  $$\begin{aligned}{} & {} p_H (\mathcal {R}_0)< p_T (\mathcal {R}_0) \nonumber \\{} & {} \quad \Leftrightarrow \dfrac{A^2}{m\mathcal {R}_0^2} < \dfrac{A^2}{m} \dfrac{\left( \mathcal {R}_0-1\right) }{\mathcal {R}_0^2} \nonumber \\{} & {} \quad \Leftrightarrow \mathcal {R}_0-1> 1 \nonumber \\{} & {} \quad \Leftrightarrow \mathcal {R}_0 > 2. \end{aligned}$$

  (12)

  The remaining inequality follows from the previous item and transitivity.

$\square$

6.2 Saddle-Node Bifurcation

The general jacobian matrix of (9) at $E^p = (S^p, I^p) \in (\mathbb {R}_0^+)^2$ coincides with that of (8). At $E_0^p$, $E_1^p$ and $E_2^p$, the matrix (8) takes the form:

$$\begin{aligned} J(E_0^p)=\left( \begin{array}{cc} A - 2S_0^p &{} -\beta S_0^p \\ \\ 0 &{} \beta S_0^p - \left( \sigma +g\right) \end{array}\right) , \qquad J(E_1^p)=\left( \begin{array}{cc} A - 2S_1^p &{} -\beta S_1^p \\ \\ 0 &{} \beta S_1^p - \left( \sigma +g\right) \end{array}\right) \end{aligned}$$

and

$$\begin{aligned} J(E_2^p)=\left( \begin{array}{cc} \dfrac{pm\beta ^2-\left( \sigma +g\right) ^2}{\beta \left( \sigma +g\right) } &{} - \left( \sigma +g\right) \\ \\ \dfrac{-pm\beta ^2 + A\left( \sigma +g\right) \beta - \left( \sigma +g\right) ^2}{\beta \left( \sigma +g\right) } &{} 0 \end{array}\right) , \end{aligned}$$

(13)

respectively.

Lemma 4

The following statements hold for model (9):

• (1) If  $p > p_{SN}(\mathcal {R}_0 )$, then $E_0^p$ and $E_1^p$ do not exist in the first quadrant of (S, I).

• (2) If  $p_T (\mathcal {R}_0 )< p < p_{SN}(\mathcal {R}_0 )$, then:

  • (a) $E_0^p$ is a saddle and $E_1^p$ is a sink for  $1< \mathcal {R}_0 < 2$;

  • (b) $E_0^p$ is a source and $E_1^p$ is a saddle for  $\mathcal {R}_0 > 2$.

• (3) If  $p<p_T(\mathcal {R}_0 )$ and  $\mathcal {R}_0 > 1$, then  $E_0^p$ and  $E_1^p$ are saddles.

Proof

The eigenvalues of $J(E_0^p)$ are $A-2S_0^p$ and $\beta S_0^p - (\sigma +g)$. We know that

$$\begin{aligned}{} & {} A-2S_0^p> 0 \\{} & {} \quad \Leftrightarrow \displaystyle \dfrac{2A - 2A + 2\sqrt{A^2 - 4pm}}{2}> 0 \\{} & {} \quad \Leftrightarrow \displaystyle \sqrt{A^2 - 4pm} > 0 \quad \hbox { if}\ p \le p_{SN}(\mathcal {R}_0 ) \end{aligned}$$

• (1) If $p> p_{SN}(\mathcal {R}_0 )$, then $E_0^p$ and $E_1^p$ do not exist.

• (2) We may write:

  $$\begin{aligned}{} & {} \beta S_0^p - (\sigma +g)< 0 \nonumber \\{} & {} \quad \Leftrightarrow S_0^p< \dfrac{\sigma + g}{\beta } \nonumber \\{} & {} \quad \overset{(6)}{\Leftrightarrow } \displaystyle \dfrac{A-\sqrt{A^2 - 4pm}}{2}< \dfrac{\sigma + g}{\beta } \nonumber \\{} & {} \quad \Leftrightarrow \displaystyle \sqrt{A^2 - 4pm}> A - \dfrac{2\left( \sigma +g\right) }{\beta } \nonumber \\{} & {} \quad \Leftrightarrow \displaystyle \sqrt{A^2 - 4pm}> \dfrac{A\mathcal {R}_0}{\mathcal {R}_0} - \dfrac{2A}{\mathcal {R}_0} \nonumber \\{} & {} \quad \Leftrightarrow \displaystyle \sqrt{A^2 - 4pm}> \dfrac{A\left( \mathcal {R}_0 - 2\right) }{\mathcal {R}_0} \nonumber \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{}A^2 - 4pm \ge 0, \qquad 1< \mathcal {R}_0< 2 \\ \\ &{}A^2 - 4pm> \dfrac{A^2\left( \mathcal {R}_0 - 2\right) ^2}{\mathcal {R}_0^2}, \qquad \mathcal {R}_0 \ge 2 \end{array}\right. } \nonumber \\{} & {} \quad \overset{(10)}{\Leftrightarrow } {\left\{ \begin{array}{ll} &{}p \le p_{SN}, \qquad 1< \mathcal {R}_0< 2 \\ \\ &{} p< \dfrac{A^2\left( \mathcal {R}_0 - 1\right) }{m\mathcal {R}_0^2} , \qquad \mathcal {R}_0 \ge 2 \end{array}\right. } \nonumber \\{} & {} \quad \overset{(10)}{\Rightarrow } {\left\{ \begin{array}{ll} &{}p< p_{SN}(\mathcal {R}_0), \qquad 1< \mathcal {R}_0< 2 \\ \\ &{}p < p_T(\mathcal {R}_0), \qquad \mathcal {R}_0 > 2 \end{array}\right. }, \end{aligned}$$

  (14)

  then the eigenvalues of $J(E_0^p)$ have real part with opposite signs and $E_0^p$ is a saddle. On the other hand, if

  $$\begin{aligned}{} & {} \beta S_0 - (\sigma +g)> 0 \nonumber \\{} & {} \quad \Leftrightarrow S_0> \dfrac{\sigma + g}{\beta } \nonumber \\{} & {} \quad \Leftrightarrow \sqrt{A^2 - 4pm}< \dfrac{A\left( \mathcal {R}_0 -2 \right) }{\mathcal {R}_0} \nonumber \\{} & {} \quad \Leftrightarrow 0 \le A^2 - 4pm< \dfrac{A^2 \left( \mathcal {R}_0 - 2 \right) ^2}{\mathcal {R}_0^2}, \qquad \mathcal {R}_0> 2 \nonumber \\{} & {} \quad \overset{(10)}{\Rightarrow } p_T(\mathcal {R}_0)< p < p_{SN}(\mathcal {R}_0), \qquad \mathcal {R}_0 > 2, \end{aligned}$$

  (15)

  then the eigenvalues of $J(E_0^p)$ have real part with positive signs, and $E_0^p$ is a source. Moreover, the eigenvalues of $J(E_1^p)$ are $A - 2S_1^p < 0$ and $\beta S_1^p - \left( \sigma +g\right)$. Hence, if

  $$\begin{aligned}{} & {} \beta S_1^p - (\sigma +g)< 0 \nonumber \\{} & {} \quad \Leftrightarrow S_1^p< \dfrac{\sigma + g}{\beta } \nonumber \\{} & {} \quad \Leftrightarrow \sqrt{A^2 - 4pm}< \dfrac{A\left( 2 - \mathcal {R}_0 \right) }{\mathcal {R}_0} \nonumber \\{} & {} \quad \Leftrightarrow 0 \le A^2 - 4pm< \dfrac{A^2 \left( 2 - \mathcal {R}_0 \right) ^2}{\mathcal {R}_0^2}, \qquad 1< \mathcal {R}_0< 2 \nonumber \\{} & {} \quad \overset{(10)}{\Rightarrow } p_T(\mathcal {R}_0)< p< p_{SN}(\mathcal {R}_0), \qquad 1< \mathcal {R}_0 < 2, \end{aligned}$$

  (16)

  then both eigenvalues of $J(E_1^p)$ have negative real part and $E_1^p$ is a sink. If

  $$\begin{aligned}{} & {} \beta S_1^p - (\sigma +g)> 0 \nonumber \\{} & {} \quad \Leftrightarrow S_1^p> \dfrac{\sigma + g}{\beta } \nonumber \\{} & {} \quad \Leftrightarrow \sqrt{A^2 - 4pm}> \dfrac{A\left( 2 - \mathcal {R}_0 \right) }{\mathcal {R}_0} \nonumber \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{}A^2 - 4pm> \dfrac{A^2\left( 2 - \mathcal {R}_0 \right) ^2}{\mathcal {R}_0^2}, \qquad 1< \mathcal {R}_0 \le 2 \\ \\ &{}A^2 - 4pm \ge 0, \qquad \mathcal {R}_0> 2 \end{array}\right. } \nonumber \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{}p< \dfrac{A^2}{m}\left( \dfrac{\mathcal {R}_0 - 1}{\mathcal {R}_0^2}\right) , \qquad 1< \mathcal {R}_0 \le 2 \\ \\ &{}p \le \dfrac{A^2}{4m}, \qquad \mathcal {R}_0> 2 \end{array}\right. } \nonumber \\{} & {} \quad \overset{(10)}{\Rightarrow } {\left\{ \begin{array}{ll} &{}p< p_T(\mathcal {R}_0), \qquad 1< \mathcal {R}_0< 2 \\ \\ &{}p < p_{SN}(\mathcal {R}_0), \qquad \mathcal {R}_0 > 2, \end{array}\right. } \end{aligned}$$

  (17)

  then the eigenvalues of $J(E_1^p)$ have real part with opposite signs and $E_1^p$ is a saddle. Therefore, from (14) and (16), we conclude that:

  • (a) if  $p_T (\mathcal {R}_0)< p < p_{SN}(\mathcal {R}_0)$ and  $1< \mathcal {R}_0 < 2$, then $E_0^p$ is a saddle and $E_1^p$ is a sink;

  • (b) if  $p_T(\mathcal {R}_0)< p < p_{SN}(\mathcal {R}_0)$, then $E_0^p$ is a source and $E_1^p$ is a saddle, for  $\mathcal {R}_0 > 2$, under conditions (15) and (17).

• (3) From (14) and (17), we conclude that if  $p < p_T(\mathcal {R}_0)$, then $E_0^p$ and $E_1^p$ are saddles.

$\square$

6.3 Transcritical Bifurcation

Lemma 5

The endemic equilibrium $E_2^p$ undergoes a transcritical bifurcation at $p = p_T(\mathcal {R}_0)$ and lies in the interior of the first quadrant if  $p <p_T(\mathcal {R}_0)$.

Proof

The equilibrium $E_2^p$ lies in the interior of the first quadrant if both $S_2^p$ and $I_2^p$ are positive. It is clear that $S_2^p =\dfrac{\sigma + g}{\beta }> 0$. Since $\beta ^2(\sigma + g) > 0$, we have

$$\begin{aligned}I_2^p> 0 \quad \Leftrightarrow \quad -pm\beta ^2 + A\left( \sigma +g\right) \beta -\left( \sigma +g\right) ^2 > 0.\end{aligned}$$

Solving the previous inequality in order to p, one gets:

$$\begin{aligned}{} & {} -pm\beta ^2 + A\left( \sigma +g\right) \beta -\left( \sigma +g\right) ^2 > 0 \\{} & {} \quad \Leftrightarrow p< \dfrac{A\left( \sigma +g\right) \beta -\left( \sigma +g\right) ^2}{m\beta ^2} \\{} & {} \quad \Leftrightarrow p< \dfrac{A\left( \sigma +g\right) ^2\beta }{m\left( \sigma +g\right) \beta ^2}-\dfrac{\left( \sigma +g\right) ^2A^2}{A^2m\beta ^2} \\{} & {} \quad \Leftrightarrow p< \dfrac{\left( \sigma +g\right) ^2 A^2}{m\beta ^2A^2}\mathcal {R}_0 - \dfrac{1}{\mathcal {R}_0^2}\dfrac{A^2}{m} \\{} & {} \quad \Leftrightarrow p< \dfrac{1}{\mathcal {R}_0} \dfrac{A^2}{m} - \dfrac{1}{\mathcal {R}_0^2}\dfrac{A^2}{m} \\{} & {} \quad \Leftrightarrow p< \dfrac{A^2}{m} \dfrac{\left( \mathcal {R}_0 - 1\right) }{\mathcal {R}_0^2} \\{} & {} \quad \overset{(10)}{\Leftrightarrow } p < p_T(\mathcal {R}_0) \end{aligned}$$

since $\mathcal {R}_0 > 1$. Hence, $E_2^p$ undergoes a trancritical bifurcation along the line $p = p_T(\mathcal {R}_0)$ and interchanges its stability with $E_2^p$ (if $\mathcal {R}_0<2$) and $E_1^p$ (if $\mathcal {R}_0>2$). $\square$

6.4 Belyakov Transition

We settle the following functions that will be used throughout this text:

$$\begin{aligned} \Delta _2(p):= & {} m^2 p^2 \beta ^4 + 4pm\left( \sigma +g\right) ^2\beta ^3 - \left[ 4A\left( \sigma +g\right) +2pm\right] \left( \sigma +g\right) ^2\beta ^2 \nonumber \\{} & {} +\, 4\left( \sigma +g\right) ^4\beta +\left( \sigma +g\right) ^4 \nonumber \\= & {} m^2 \beta ^4 p^2 + 2m\beta ^2\left( 2\beta -1 \right) \left( \sigma + g \right) ^2 p + \left( 4\beta +1\right) \left( \sigma +g\right) ^4 \nonumber \\{} & {} -\, 4A\left( \sigma +g\right) ^3\beta ^2, \end{aligned}$$

(18)

$$\begin{aligned} p_{B_t}^{(1)}\left( \mathcal {R}_0 \right):= & {} \dfrac{\left( -2\beta +1-2\displaystyle \sqrt{\beta \left( \mathcal {R}_0 + \beta - 2\right) }\right) A^2}{m\mathcal {R}_0^2} , \end{aligned}$$

(19)

$$\begin{aligned} p_{B_t}^{(2)}\left( \mathcal {R}_0 \right):= & {} \dfrac{\left( -2\beta +1+2\displaystyle \sqrt{\beta \left( \mathcal {R}_0 + \beta - 2\right) }\right) A^2}{m\mathcal {R}_0^2} , \end{aligned}$$

(20)

for  $\beta \left( \mathcal {R}_0 + \beta - 2 \right)> 0 \Leftrightarrow \beta > 2 - \mathcal {R}_0$, where  $p_{B_t}^{(1)}(\mathcal {R}_0)$ and  $p_{B_t}^{(2)}(\mathcal {R}_0)$ are the square roots of  $\Delta _2(p)$ written as function of $\mathcal {R}_0$. It is easy to verify that  $p_{B_t}^{(1)} (\mathcal {R}_0)< p_{B_t}^{(2)}(\mathcal {R}_0)$.

Lemma 6

If  $\beta \ge \frac{1}{2}$ and  $\mathcal {R}_0 > 1 + \frac{1}{4\beta }$, then  $p_{B_t}^{(2)} (\mathcal {R}_0)> 0$.

Proof

From (20) we know that if

$$\begin{aligned}{} & {} \dfrac{\left( -2\beta +1+2\displaystyle \sqrt{\beta \left( \mathcal {R}_0 + \beta - 2\right) }\right) A^2}{m\mathcal {R}_0^2}> 0 \\{} & {} \quad \Leftrightarrow \left( -2\beta +1+2\displaystyle \sqrt{\beta \left( \mathcal {R}_0 + \beta - 2\right) }\right) A^2> 0 \\{} & {} \quad \Leftrightarrow 2\displaystyle \sqrt{\beta \left( \mathcal {R}_0 + \beta - 2\right) }> 2\beta - 1\\{} & {} \quad \Leftrightarrow \displaystyle \sqrt{\beta \left( \mathcal {R}_0 + \beta - 2\right) }> \beta - \dfrac{1}{2} \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{} \beta - \dfrac{1}{2} \ge 0, \qquad \beta \left( \mathcal {R}_0 + \beta - 2 \right)> \left[ \beta -\dfrac{1}{2}\right] ^2 \\ \\ &{} \beta - \dfrac{1}{2}< 0, \qquad \beta \left( \mathcal {R}_0 + \beta - 2 \right) \ge 0 \end{array}\right. } \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{} \beta \ge \dfrac{1}{2}, \qquad \beta \mathcal {R}_0 - 2\beta> -\beta + \dfrac{1}{4} \\ \\ &{}\beta< \dfrac{1}{2}, \qquad \beta ^2 + \mathcal {R}_0 \beta - 2\beta \ge 0 \end{array}\right. } \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{} \beta \ge \dfrac{1}{2}, \qquad \beta \left( \mathcal {R}_0 - 1 \right)> \dfrac{1}{4} \\ \\ &{}\beta< \dfrac{1}{2}, \qquad \beta \ge 2 - \mathcal {R}_0 \end{array}\right. } \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{} \beta \ge \dfrac{1}{2}, \qquad \mathcal {R}_0 > 1+\dfrac{1}{4\beta } \\ \\ &{}\beta < \dfrac{1}{2}, \qquad \mathcal {R}_0 \ge 2 - \beta \end{array}\right. }, \end{aligned}$$

then $p_{B_t}^{(2)} (\mathcal {R}_0)> 0$. If  $\mathcal {R}_0 = 1+\frac{1}{4\beta }$, then  $p_{B_t}^{(2)}(\mathcal {R}_0) = 0$ and if  $\mathcal {R}_0 = 2$, then  $p_{B_t}^{(2)}(\mathcal {R}_0) = p_{SN}(\mathcal {R}_0)$. Hence we conclude that if  $\beta \ge \frac{1}{2}$ and $\mathcal {R}_0 > 1 + \frac{1}{4\beta }$, then  $p_{B_t}^{(2)}(\mathcal {R}_0) > 0$. $\square$

Before we analyze the stability of $E_2^p$, we remind the readers that $E_2^p$ exists in the first quadrant when $p < p_T(\mathcal {R}_0)$ (by Lemma 5).

Lemma 7

The following statements hold for model (9):

• (1) If  $p_{B_t}^{(2)}(\mathcal {R}_0)<p<p_T(\mathcal {R}_0)$, then  $\Delta _2(p) > 0$ and $E_2^p$ is:

  • (a) stable node for  $\mathcal {R}_0<2$;

  • (b) unstable node for  $p > p_H(\mathcal {R}_0)$.

• (2) If  $p<p_{B_t}^{(2)}(\mathcal {R}_0)$, then  $\Delta _2(p) < 0$ and  $E_2^p$ is:

  • (a) stable focus for  $p < p_H(\mathcal {R}_0)$;

  • (b) unstable focus for  $p > p_H(\mathcal {R}_0)$.

Hence, $E_2^p$ undergoes a Belyakov transition at  $p=p_{B_t}^{(2)}(\mathcal {R}_0)$.

Proof

With respect to (9), the eigenvalues of the jacobian matrix of $J(E_2^p)$ are:

$$\begin{aligned}\lambda _1 = \dfrac{p m \beta ^2 - \left( \sigma + g \right) ^2 - \displaystyle \sqrt{\Delta _2(p)}}{2\beta \left( \sigma + g\right) } \qquad \text {and} \qquad \lambda _2 = \dfrac{p m \beta ^2 - \left( \sigma + g \right) ^2 + \displaystyle \sqrt{\Delta _2(p)}}{2\beta \left( \sigma + g\right) }.\end{aligned}$$

where $\lambda _1 < \lambda _2$. Hence, if $\lambda _1>0$, then $\lambda _2>0$. If $\lambda _2<0$, then $\lambda _1<0$. We know from (18), (19), (20) and Lemma 6 that:

$$\begin{aligned} {\left\{ \begin{array}{ll}\,\, \Delta _2 > 0 \,, \qquad \text {if} \quad \,0<p<p_{B_t}^{(1)}{} or \,p_{B_t}^{(2)}<p<+\infty \\ \\ \,\, \Delta _2< 0 \,, \qquad \text {if} \quad \,p_{B_t}^{(1)}<p<p_{B_t}^{(2)} \end{array}\right. }, \end{aligned}$$

where $p_{B_t}^{(2)} (\mathcal {R}_0)< p_T(\mathcal {R}_0)$.

(1) Therefore, if

$$\begin{aligned}{} & {} \lambda _1> 0 \\{} & {} \quad \Leftrightarrow \dfrac{p m \beta ^2 - \left( \sigma + g \right) ^2 - \displaystyle \sqrt{\Delta _2(p)}}{2\beta \left( \sigma + g\right) }> 0 \\{} & {} \quad \Leftrightarrow \displaystyle \sqrt{\Delta _2(p)}< p m \beta ^2 - \left( \sigma + g \right) ^2 \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{}\Delta _2(p) \ge 0 \\ \\ &{}\Delta _2(p)< \left[ p m \beta ^2 - \left( \sigma + g \right) ^2 \right] ^2 \end{array}\right. } \qquad , \quad p m \beta ^2 - \left( \sigma + g \right) ^2> 0 \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{}p \ge p_{B_t}^{(2)} (\mathcal {R}_0)\\ \\ &{} p< p_T(\mathcal {R}_0) \end{array}\right. } \qquad , \quad p> p_H (\mathcal {R}_0) \\{} & {} \quad \Rightarrow p_{B_t}^{(2)}(\mathcal {R}_0)< p < p_T(\mathcal {R}_0) , \qquad \text {for} \quad p > p_H (\mathcal {R}_0) \end{aligned}$$

then $\lambda _1> 0 \Rightarrow \lambda _2 > 0$ and $E_2^p$ is an unstable node. Analogously, we proceed in the same way for $\lambda _2 < 0$. If

$$\begin{aligned}{} & {} \lambda _2< 0 \\{} & {} \quad \Leftrightarrow \dfrac{p m \beta ^2 - \left( \sigma + g \right) ^2 + \displaystyle \sqrt{\Delta _2(p)}}{2\beta \left( \sigma + g\right) }< 0 \\{} & {} \quad \Leftrightarrow \displaystyle \sqrt{\Delta _2(p)}< \left( \sigma + g \right) ^2 - p m \beta ^2 \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{}\Delta _2(p) \ge 0 \\ \\ &{}\Delta _2(p)< \left[ \left( \sigma + g \right) ^2 - p m \beta ^2 \right] ^2 \end{array}\right. } \qquad , \quad \left( \sigma + g \right) ^2 - p m \beta ^2 > 0 \\{} & {} \quad \Leftrightarrow {\left\{ \begin{array}{ll} &{}p \ge p_{B_t}^{(2)} (\mathcal {R}_0)\\ \\ &{} p< p_T(\mathcal {R}_0) \end{array}\right. } \qquad , \quad p< p_H(\mathcal {R}_0) \\{} & {} \quad \Rightarrow p_{B_t}^{(2)} (\mathcal {R}_0)< p<p_T (\mathcal {R}_0) , \qquad \text {for} \quad p < p_H (\mathcal {R}_0) \end{aligned}$$

then $\lambda _2< 0 \Rightarrow \lambda _1 < 0$ and $E_2^p$ is a stable node.

(2) Now, let us assume that  $\Delta _2(p) < 0$, which means that  $p<p_{B_t}^{(2)}(\mathcal {R}_0)$, the eigenvalues are complex and  $E_2^p$ is a focus. If the real part of the eigenvalues of  $J(E_2^p)$ is negative (and  $\beta$ is positive), i.e. if

$$\begin{aligned}{} & {} pm\beta ^2 - \left( \sigma +g\right) ^2< 0 \\{} & {} \quad \Leftrightarrow pm\beta ^2< \left( \sigma +g\right) ^2 \\{} & {} \quad \Leftrightarrow p< \dfrac{\left( \sigma +g\right) ^2}{m\beta ^2} \\{} & {} \quad \Leftrightarrow p< \dfrac{A^2}{m\mathcal {R}_0^2} \\ \overset{(10)}{\Leftrightarrow }{} & {} p < p_H(\mathcal {R}_0), \end{aligned}$$

then  $E_2^p$ is a stable focus. Otherwise  $E_2^p$ is an unstable focus. Also, we conclude that  $E_2^p$ undergoes a Belyakov transition at  $p=p_{B_t}^{(2)}(\mathcal {R}_0)$, since $E_2^p$ changes its stability from a node to a focus. $\square$

6.5 Hopf Bifurcation

In this section we will find a line in the parameter space $(p, \mathcal {R}_0)$ where a Hopf bifurcation occurs.

Lemma 8

If  $p = p_H(\mathcal {R}_0)$, then $E_2^p$ undergoes a subcritical Hopf bifurcation associated to an unstable periodic solution $\mathcal {C}$ (as $\mathcal {R}_0$ increases).

Proof

Let the eigenvalues of $J(E_2^p)$ given by

$$\begin{aligned}\lambda _1 = \dfrac{p m \beta ^2 - \left( \sigma + g \right) ^2 - \displaystyle \sqrt{\Delta _2(p)}}{2\beta \left( \sigma + g\right) } \qquad \text {and} \qquad \lambda _2 = \dfrac{p m \beta ^2 - \left( \sigma + g \right) ^2 + \displaystyle \sqrt{\Delta _2(p)}}{2\beta \left( \sigma + g\right) }.\end{aligned}$$

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_2^p)$ (see (13)) have the form $\pm i\omega$ ($\omega >0$), which is equivalent to $\textrm{Tr} \,J(E_2^p) = 0$, where $\textrm{Tr}$ represents the usual trace operator. Indeed,

  $$\begin{aligned}{} & {} \textrm{Tr} \,J(E_2^p) = 0 \\{} & {} \quad \Leftrightarrow \dfrac{pm\beta ^2-\left( \sigma +g\right) ^2}{\beta \left( \sigma +g\right) } = 0\\{} & {} \quad \Leftrightarrow pm\beta ^2-\left( \sigma +g\right) ^2 = 0 \\{} & {} \quad \Leftrightarrow p = \dfrac{\left( \sigma +g\right) ^2}{m\beta ^2} \\{} & {} \quad \Leftrightarrow p = \dfrac{A^2}{m\mathcal {R}_0^2} \\{} & {} \quad \overset{(10)}{\Leftrightarrow } p = p_H (\mathcal {R}_0). \end{aligned}$$

• (2) $\dfrac{\textrm{d}}{\textrm{d}\mathcal {R}_0} \left( {\text {Re}} \lambda _j(\mathcal {R}_0)\right) \big |_{\mathcal {R}_0 \in p_H} \ne 0$,    for $j \in \{1,2\}$:

  Hence, for $p=p_H(\mathcal {R}_0)$ we get

  $$\begin{aligned} \dfrac{\textrm{d}}{\textrm{d}\mathcal {R}_0} \left( {\text {Re}} \lambda _j(\mathcal {R}_0)\right) \big |_{\mathcal {R}_0 = \mathcal {R}_0^{\star }}= & {} \dfrac{\textrm{d}}{\textrm{d}\mathcal {R}_0} \left( \dfrac{pm\beta ^2}{2\beta \left( \sigma +g\right) } - \dfrac{\left( \sigma +g\right) ^2}{2\beta \left( \sigma +g\right) }\right) \Big |_{\mathcal {R}_0 = \mathcal {R}_0^{\star }} \nonumber \\= & {} \dfrac{\textrm{d}}{\textrm{d}\mathcal {R}_0} \left( \dfrac{\sigma + g}{2\beta } - \dfrac{A}{2\mathcal {R}_0}\right) \Big |_{\mathcal {R}_0 = \mathcal {R}_0^{\star }} \nonumber \\= & {} \dfrac{A}{2{\mathcal {R}_0^{\star }}^2} \ne 0. \end{aligned}$$

  (21)

Hence, $E_2^p$ undergoes a subcritical Hopf bifurcation (as $\mathcal {R}_0$ increases; the periodic solution appears for $\mathcal {R}_0<\mathcal {R}_0^{\star }$ as illustrated in Fig. 4. Figure 2 provides a better perspective of all bifurcations under consideration. $\square$

Fig. 4

Phase portraits of system (9) when $E_2^p$ undergoes a subcritical Hopf bifurcation as $\mathcal {R}_0$ increases

Table 2 13 interpolating points of the heteroclinic cycle function in the space of parameters $(\mathcal {R}_0,p)$

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 1”.

Using MATLAB_R2018a, we locate thirteen pairs $({\mathcal {R}_0}_i,p_i)$ in the parameter space $(\mathcal {R}_0,p)$ where one observes a heteroclinic cycle associated to the disease-free equilibria $E_0^{p_i}$ and $E_1^{p_i}$ (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 (Fig. 5):

$$\begin{aligned}y = ax^b + c \quad \text {for}\quad a,b,c \in \mathbb {R}.\end{aligned}$$

Using again MATLAB_R2018a, we obtain $a=4.495$, $b= -2.313$ and $c= - 0.039$, and the graph of $p_{\text {Het.}}\left( \mathcal {R}_0 \right) = \dfrac{4.495}{\mathcal {R}_0^{2.313}} - 0.039$ is plotted in Fig. 6.

The heteroclinic bifurcation organizes the dynamics, stressing a geometric configuration in its unfolding.

7 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 Figs. 2 and 5.

Fig. 5

Numerical simulations of DZ bifurcation diagram with $A = 1.1$, $\beta = 1.3$, $m = \sigma = g = 0.35$ for regions A, D, Het., E, F, G and H. For region B we use $A = 1$, $\beta = 1.3$, $m = 0.35$ and $\sigma = g = 0.50$, and for region C we use $A = 1.1$, $\beta = 0.91$ and $m = \sigma = g = 0.35$. For all regions we consider $p \in [0, 1]$. Compare with the theoretical description in Fig. 2

A: :

The vaccination coverage is so large that Infectious individuals tend to disappear.

B: :

The flow exhibits two disease-free equilibria, a saddle $E_0^p$ and a sink $E_1^p$. 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_2^p$ 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_2^p$, so the disease remains persistently in the population. Initial conditions lying “below” $W^s(E_0^p)$ have $E_2^p$ as $\omega$-limit.

E: :

The unstable heteroclinic cycle $\mathcal {H}$ is broken giving rise to an unstable periodic solution $\mathcal {C}$. The region bounded by $\mathcal {C}$ is conducive to the persistence of the disease since all initial conditions tend toward $E_2^p$. All initial conditions lying in the unbounded region defined by $\mathcal {C}$ are propitious to the disappearance of the disease. The period of $\mathcal {C}$ is very large when $(\mathcal {R}_0, p)$ is close to the curve Het.

F and G: :

The equilibrium $E_2^p$ 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_0^p$ and $E_1^p$ repel Lebesgue almost all trajectories towards the origin.

8 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.

Fig. 6

Function obtained by interpolation of 13 points (see Table 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 $(\mathcal {R}_0,p)$ is $p\left( \mathcal {R}_0 \right) \approx 4.495 \mathcal {R}_0^{-\,2.313} - 0.039$

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_2^p$) and we have described all associated bifurcation curves passing through the bifurcation point $\left( \mathcal {R}_0^{\star },p^\star \right) = \left( 2,\frac{A^2}{4\,m}\right)$ in the parameter space. In Fig. 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 $\mathcal {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_{ \left( 2, \frac{A^2}{4m}\right) }$ at $E_2^p$ 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(\mathcal {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 $\mathcal {R}_0 \ge 1$:

• $1<\mathcal {R}_0 <2$: the increase of the vaccination rate plays an essencial role to decrease the number of Infectious;

• $2\le \mathcal {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 Fig. 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) in the following way:

$$\begin{aligned} \dot{S} = S(A-S) - \beta IS - p \displaystyle \sum _{n=0}^{+\infty } S(nT^-) \delta (t-nT)\end{aligned}$$

where $n \in {\mathbb {N}}$ and:

• T is the period between two consecutive “shots” (of vaccination);

• $\displaystyle S(nT^-)=\lim \nolimits _{t\rightarrow n T^-} S(t)$ is the number of Susceptible individuals at the instant immediately before being vaccinated and $\delta$ is the usual Dirac function.

The analysis of an impulsive differential equation involving pulse vaccination is an ongoing work.

MATLAB R2018a: Code and Procedure

We set the code of MATLAB_R2018a to obtain the interpolation function that best fits the $({\mathcal {R}_0}_i, p_i)$ points in the parameter space $(\mathcal {R}_0, p)$.

f=fit(x,y,‘power2’)

plot(f, x, y)

xlim([2 3.8])

ylim([0.1 0.75])

ax = gca;

ax.FontSize = 30;

xlabel(‘$$\mathcal {R}\_0$$’,

‘Interpreter’,‘latex’)

ylabel(‘$p$’,

‘Interpreter’,‘latex’)

leg1 = legend(‘($$\mathcal {R}\_0$_i,p_i$)

points’,‘Het. cycle function’)

lgd.FontSize = 20;

legend boxoff

set(leg1,‘Interpreter’,‘latex’);

where x and y are the vectors of ${\mathcal {R}_0}_i$ and $p_i$ values, respectively. The interpolation function given by the MATLAB_R2018a has a correlation coefficient equal to 1 (precision: $10^{-5}$; confidence bound: 95%), as we can confirm in Fig. 7.

Fig. 7

MATLAB_R2018a output for the two-term power series model used in $({\mathcal {R}_0}_i,p_i)$ points interpolation