On the Jacobi Stability of Two SIR Epidemic Patterns With Demography

1. Introduction

Mathematical patterns that do not explicitly include births and deaths occurring in the population are called epidemic patterns without demography. They are useful for epidemic modeling on a short time scale, particularly for modeling some epidemic like influenza. Omitting population change requires that the disease develop on a much shorter time scale than the period in which significant change in the population size can occur (such as births and deaths). This is valid for fast diseases like the childhood diseases and influenza. On the other hand, there are slow diseases, such as HIV, tuberculosis, and hepatitis C, that develop for a long period of time even on an individual level. In this case, the total population does not remain constant for long periods of time, and the demography of the population cannot be ignored [1,2,3].

The main goal of the present paper is to study the Jacobi stability for two SIR models with demography: the classical SIR model and a modified SIR model with a variable transmission coefficient of the infection (linear with respect to the number I of the infected individuals). The classical stability (linear or Lyapunov stability) of different kinds of SIR (susceptible, infected and removed individuals) epidemic models without or with demography was widely completely studied in last decades [1,2,3,4,5,6]. In this work, we will study another type of stability for these systems, namely the Jacobi stability. Classical approaches of this models can be found in [1,2,3,4]. Generally, the study of mathematical deterministic models for the spread of diseases, or for the interaction between prey-predator type populations is performed only for positive values of the variables, where they has an ecological, biological or epidemiological meaning [4,5,6,7,8,9].

The Jacobi stability is a natural extension of the geometric stability of the geodesic flow, from a manifold with a Riemann or Finsler metric to a manifold with no metric [10,11,12,13,14,15]. Practically, the Jacobi stability indicates the robustness of a dynamical system defined by a system of second-order differential equations (or SODE), where this robustness measures the lack of sensitivity and adaptation to both the modifications of the internal parameters of the system and the change in the external environment. By using the Kosambi–Cartan–Chern (KCC) theory, the local behavior of dynamical systems from the point of view of the Jacobi stability has recently been studied by several authors in [11,12,16,17,18,19,20,21,22,23,24,25,26]. Hence, the local dynamics of the system is investigated by using the geometric objects that correspond to the system of the second-order differential equations (SODE), which is the system obtained from the given first-order differential system [27,28,29].

The KCC theory deals with the study of the deviation of neighboring trajectories, which allows us to to estimate the perturbation permitted around the equilibrium points of the second-order differential system. Initially, this approach was linked with the study of the variation equations (or Jacobi field equations) associated with the geometry on the differentiable manifold. More exactly, P. L. Antonelli, R. Ingarden, and M. Matsumoto started the study of the Jacobi stability for the geodesics corresponding to a Riemannian or Finslerian metric by deviating the geodesics and using the KCC-covariant derivative for the differential system in variations [10,11,12]. As a result, the second KCC-invariant appeared, also called the deviation curvature tensor, which is essential for the establishment of the Jacobi stability for geodesics and generally for the trajectories associated with a system of second-order differential equations. In differential geometry, a system of second-order differential equations (SODE) is called semi-spray. Using a semi-spray, we can define a nonlinear connection on the manifold, and conversely, using a nonlinear connection, we can define a semi-spray. Consequently, any semi-spray (or SODE) can define a geometry on the manifold by the associated geometric objects, and conversely [13,30,31,32]. Of course, these geometric objects are tensors that can check the properties of symmetry or not, depending on the particularity of the system of second-order differential equations (SODE).

The KCC theory originated from the works of D. D. Kosambi [27], E. Cartan [28], and S. S. Chern [29]; hence, the abbreviation KCC (Kosambi–Cartan–Chern). This geometric theory can be successfully applied in many research fields, from in engineering, physics, chemistry, and biology [16,19,21,23,24,33]. In addition, new approaches and results from the KCC theory in gravitation and cosmology were made in [34,35]. Moreover, in [18], a comprehensive analysis of the Jacobi stability and its relations with the Lyapunov stability for dynamical systems was made by C.G. Boehmer, T. Harko, and S.V. Sabau, who modeled the phenomena based on gravitation and astrophysics. Recently, in [36], a very interesting and complete study of the Jacobi stability for prey–predator models of Holling’s type II and III was performed.

In the second section, there will be a short, but comprehensive presentation of the classical SIR epidemic model with demography, and we will point out the main results on the local and global stability of this system. Next, in the third section, a reformulation of the classical SIR epidemic model with demography (3) as a system of second-order differential equations will be obtained, and the five geometrical invariants for this system will be computed. The obtained results for the Jacobi stability of this system near the equilibrium points will be presented in the fourth section. More precisely, we will find the necessary and sufficient conditions that would allow us to obtain the Jacobi stability of the system near the equilibrium points. Consequently, for these parameter values, it is not possible to have a chaotic behavior for the classical SIR epidemic model with demography. Furthermore, at the end of the fourth section, we will obtain the deviation equations near every equilibrium point as well as the curvature of the deviation vector, and then we perform a comparative analysis of the Jacobi stability and the Lyapunov stability in order to compare these two approaches. Using the same approach, in the fifth section we will obtain similar results for a classical SIR epidemic model with demography and vaccination.

Further, a modified SIR epidemic model with demography with a linear transmission coefficient of infection was presented in the sixth section [1]. After a second order reformulation of this modified SIR model in the seventh section, the Jacobi stability analysis of this system was done in the section eight. Conclusions and possible future research was presented in the last section.

At the end, in the appendix A was presented a comprehensive review of the basic notions and main tools of the KCC theory that are necessary for analyzing the Jacobi stability of dynamical systems. More exactly, we will present the five invariants of the KCC theory and the definition of the Jacobi stability. As usual in differential geometry, the sum over the crossed repeated indices is understood.

2. A Classical SIR Epidemic Pattern with Demography

A classical SIR epidemic model for the spread of the diseases suppose that the whole population $N\left(t\right)$ at a time t is split in three categories: $S\left(t\right)$ is the number of individuals which are susceptible at the moment t, $I\left(t\right)$ is the number of the infected individuals at t and $R\left(t\right)$ represent the number of removed individuals at the time t, where by a removed individual we understand individuals either recovered from the illness or dead after the infection [1,2,3,4,5,6]. A SIR model without demography means that the total number of individuals $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)$ is constant, i.e. no born, no immigration, no dead and no emigration. But, practically this reality it is not possible and then it was considered different SIR models with demography, for which $N\left(t\right)$ is not constant. To incorporate the demographics into a classical SIR epidemic model, we assume that all individuals are born susceptible [1]. Individuals from each category die at a per capita death rate $\mu$, so the total death rate in the susceptible category is $\mu S$, while in the infectious category, it is $\mu I$, and in the removed category, it is $\mu R$. If we denote by $\alpha$ the recovery rate, then the rate of the recoveries in the infected people is $\alpha I$. Also, as for the SIR model without demography, if we will denote by $\beta$ the transmission coefficient of the infection, then we can consider the following classical SIR model with demography:

$$\left\{\begin{array}{ccc}\hfill S^{\prime }& =& \mathrm{\Lambda }-\beta IS-\mu S\hfill \\ \hfill I^{\prime }& =& \beta IS-\alpha I-\mu I\hfill \\ \hfill R^{\prime }& =& \alpha I-\mu R\hfill \end{array}\right.$$

(1)

where $\mathrm{\Lambda }$ is the total birth rate (measured in number of people born per unit of time) and $S^{\prime }=S^{\prime }\left(t\right)$, $I^{\prime }=I^{\prime }\left(t\right)$, $R^{\prime }=R^{\prime }\left(t\right)$ denote the derivatives with respect to time t, or the rates of changes of these quantities in a short unity of time.

Let us remark that the third equation of the system (1) was added in order to obtain the next differential equation for the total population $N\left(t\right)$:

$$N^{\prime }\left(t\right)=\mathrm{\Lambda }-\mu N\left(t\right),$$

with the unique solution $N\left(t\right)=N\left(0\right)e^{-\mu t}+{\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\left(1-e^{-\mu t}\right)$. Therefore, the population size is not constant, but it is asymptotically constant, since $N\left(t\right)\to {\displaystyle \frac{\mathrm{\Lambda }}{\mu }}$ when $t\to \infty$.

Obviously, we can observe that the first two equations in (1) are independent of the third, and then we can consider the two-dimensional autonomous system of first order differential equations:

$$\left\{\begin{array}{ccc}\hfill S^{\prime }& =& \mathrm{\Lambda }-\beta IS-\mu S\hfill \\ \hfill I^{\prime }& =& \beta IS-(\alpha +\mu )I\hfill \end{array}\right.$$

(2)

where $R\left(t\right)=N\left(t\right)-S\left(t\right)-I\left(t\right)$.

Compared to $\beta$ - the transmission coefficient of the infection, which has the unit [number of people × time]${}^{-1}$, the parameters $\alpha$ - the recovery rate and $\mu$ - the dead rate, have the units [unit of time]${}^{-1}$. So, it is better to consider $\tau =(\alpha +\mu )t$ and then $\tau$ is a dimensionless quantity. If we denote $I\left(t\right)=I\left(\frac{\tau }{\alpha +\mu }\right)=\widehat{I}\left(\tau \right)$ and $S\left(t\right)=S\left(\frac{\tau }{\alpha +\mu }\right)=\widehat{S}\left(\tau \right)$, then we have $\frac{d\widehat{S}}{d\tau }=\frac{1}{\alpha +\mu }\frac{dS}{dt}$ and $\frac{d\widehat{I}}{d\tau }=\frac{1}{\alpha +\mu }\frac{dI}{dt}$, by using the chain rule. After resizing the variables $\widehat{S}$ and $\widehat{I}$ with the total limiting population size $\frac{\mathrm{\Lambda }}{\mu }$, we will obtain the new variables $x\left(\tau \right)=\frac{\mu \widehat{S}}{\mathrm{\Lambda }}$ and $y\left(\tau \right)=\frac{\mu \widehat{I}}{\mathrm{\Lambda }}$, which are also dimensionless quantities.

Consequently, the system (2) has the next form [1]:

$$\left\{\begin{array}{ccc}\hfill x^{\prime }& =& \rho (1-x)-{\mathcal{R}}_{0}xy\hfill \\ \hfill y^{\prime }& =& ({\mathcal{R}}_{0}x-1)y\hfill \end{array}\right.$$

(3)

where $\rho ={\displaystyle \frac{\mu }{\alpha +\mu }}$ and ${\mathcal{R}}_0={\displaystyle \frac{\mathrm{\Lambda }\beta }{\mu (\alpha +\mu )}}$ are both dimensionless parameters.

Therefore, we can say that the system (2) was transformed into a system with dimensionless form (3), which is equivalent to the original one, since the solutions of both systems have the same long-term behavior. Moreover, the number of parameters was reduced from four to two. The notation ${\mathcal{R}}_0$ is not choose randomly, because this dimensionless quantity is indeed the reproduction number or the basic reproduction number for this mathematical epidemiological pattern [1,5,6].

Because it is impossible to solve the first order differential system associated to this SIR pattern with demography by analytical methods, it remains only to find information about the behavior of the solutions. Especially, because we want to know what will happen to the disease in the long run: will it go out, or will it establish itself in the population and become endemic? So, the long-term behavior of the solutions is very important from an epidemiological perspective [1,2,3].

It is clear that the model has relevance when $x\ge 0$, $y\ge 0$, and then the solutions of the system (3) lie in the first quadrant ${\Sigma }_{+}^0=\left\{(x,y)\in {\mathbf{R}}^2\mid x\ge 0,y\ge 0\right\}$. Moreover, the lines $\{x=0\}$ and $\{y=0\}$ are invariant manifolds with respect to the flow of the system, that means any orbit starting from a point which belongs to ${\Sigma }_{+}=\left\{(x,y)\in {\mathbf{R}}^2\mid x>0,y>0\right\}$ remains in ${\Sigma }_{+}$. So, the orbits cannot cross any of these two invariant lines and then the study of the system where it has a epidemiological relevance is well-defined, in the sense that, an orbit starting from a zone with epidemiological relevance does not enter in a zone with epidemiological irrelevance and conversely.

However, it is very important to study the local dynamics of the system with the classical tools of dynamical systems and, also, by geometrical tools of Kosambi-Cartan-Chern (KCC) theory. First of all, we must determine the equilibrium points of the epidemic given by dynamical system (3), by solving the system

$$\left\{\begin{array}{ccc}\rho (1-x)-{\mathcal{R}}_{0}xy\hfill & =& \hfill 0\\ ({\mathcal{R}}_{0}x-1)y\hfill & =& \hfill 0\end{array}\right.$$

and so we get two equilibrium points $E_0(1,0)$, so-called the disease-free equilibrium, and $E_1\left(\frac{1}{{\mathcal{R}}_0},\rho \left(1-\frac{1}{{\mathcal{R}}_0}\right)\right)$, so-called the endemic equilibrium.

Let us remark that the endemic equilibrium exists if and only if the basic reproduction number ${\mathcal{R}}_0>1$. Else, if ${\mathcal{R}}_0=1$, then $E_1$ coincides with $E_0$, or, if ${\mathcal{R}}_0<1$, then $E_1$ is a virtual equilibrium point, i.e. $E_1$ don’t belongs to ${\Sigma }_{+}^0$ and then $E_1$ is irrelevant.

According to Hartman-Grobman Theorem, it is known that the local stability of hyperbolic equilibrium points is determined by the eigenvalues of the Jacobi matrix computed at each equilibrium point. Because the Jacobi matrix of the system (3) at the point $(x,y)$ is

$$A=\left(\begin{array}{cc}-\rho -{\mathcal{R}}_{0}y& -{\mathcal{R}}_{0}x\\ {\mathcal{R}}_{0}y& {\mathcal{R}}_{0}x-1\end{array}\right),$$

it results that:

• For the disease-free equilibrium $E_0(1,0)$, the Jacobi matrix is $A=\left(\begin{array}{cc}-\rho & -{\mathcal{R}}_0\\ 0& {\mathcal{R}}_0-1\end{array}\right)$ with eigenvalues ${\lambda }_1=-\rho$, ${\lambda }_2={\mathcal{R}}_0-1$. Then $E_0$ is unstable (saddle point) if and only if ${\mathcal{R}}_0>1$, and $E_0$ is locally asymptotically stable (stable node) if and only if ${\mathcal{R}}_0<1$. If ${\mathcal{R}}_0=1$, then $E_0$ is a non hyperbolic equilibrium point and we cannot apply Hartman-Grobman Theorem for the local behavior study.
• For the endemic equilibrium $E_1\left(\frac{1}{{\mathcal{R}}_0},\rho \left(1-\frac{1}{{\mathcal{R}}_0}\right)\right)$, the Jacobian $A=\left(\begin{array}{cc}-\rho {\mathcal{R}}_0& -1\\ \rho ({\mathcal{R}}_0-1)& 0\end{array}\right)$
  with characteristic polynomial ${\lambda }^2+\rho {\mathcal{R}}_0\lambda +\rho ({\mathcal{R}}_0-1)=0$ and eigenvalues
  ${\lambda }_{1,2}={\displaystyle \frac{1}{2}}\left(-\rho {\mathcal{R}}_0\pm \sqrt{\rho (\rho {\mathcal{R}}_0^2-4{\mathcal{R}}_0+4)}\right)$. Since ${\mathcal{R}}_0>1$, we have that ${\lambda }_1+{\lambda }_2=-\rho {\mathcal{R}}_0<0$ and ${\lambda }_1{\lambda }_2=\rho ({\mathcal{R}}_0-1)>0$, that means that, if exists, the endemic equilibrium $E_1$ is always locally asymptotically stable (stable node or stable focus). More precisely, $E_1$ is a stable node if and only if $\rho {\mathcal{R}}_0^2-4{\mathcal{R}}_0+4\ge 0$, and $E_1$ is a stable focus if and only if $\rho {\mathcal{R}}_0^2-4{\mathcal{R}}_0+4<0$.

In conclusion, the basic reproduction number ${\mathcal{R}}_0$ of the disease modelled by (3) has the next threshold role [1,5,6]:

• If ${\mathcal{R}}_0<1$, then there exists only the disease-free equilibrium point, which is an attractive equilibrium (stable node), i.e. any trajectory of the dynamical system (3) starting near to $E_0$ converges at this equilibrium when time tends to infinity, and the disease disappears from the population.
• If ${\mathcal{R}}_0>1$, then there exists two equilibrium points: the disease-free equilibrium and the endemic equilibrium. The disease-free equilibrium is not attractive (unstable, a saddle point), in the sense that there exists trajectories of the system (3) that start very close to $E_0$, but it tend to go away. Instead, the endemic equilibrium is attractive (stable node or stable focus), that means any orbit of the system (3) starting near to $E_1$ converge to $E_1$ as time goes to infinity. So, in this case, the disease remains endemic in the population.

Related to the global stability, an equilibrium point is called globally stable if it is stable for almost all initial conditions, not just those that are close to it. For this classical SIR system are known the following results [1,5] (see also the Dulac–Bendixson Theorem):

Theorem 2.1 

a) For ${\mathcal{R}}_0<1$, the disease-free equilibrium $E_0$ is globally stable.

b) For ${\mathcal{R}}_0>1$, the system (3) has no periodic orbits.

c) For ${\mathcal{R}}_0>1$, the endemic equilibrium $E_1$ is globally stable whenever $y\left(0\right)>0$.

In conclusion, we collect the obtained results for the local stability in the next Table 1:

In the next sections, we will focus only on the study of the Jacobi stability in order to clarify the behavior of the SIR system and to point out the geometric objects associated to this system of ordinary differential equations.

3. SODE Formulation of the Classical SIR Pattern with Demography

We consider the classical SIR model with demography (3). For the sake of simplicity we will denote the derivative with respect to time with "dot" over the variable and we prefer to denote t instead of $\tau$. Then the system (3) becomes

$$\left\{\begin{array}{ccc}\hfill \dot{x}& =& \rho (1-x)-{\mathcal{R}}_{0}xy\hfill \\ \hfill \dot{y}& =& ({\mathcal{R}}_{0}x-1)y\hfill \end{array}\right.$$

(4)

where $\rho ={\displaystyle \frac{\mu }{\alpha +\mu }}\in (0,1)$ and ${\mathcal{R}}_0={\displaystyle \frac{\mathrm{\Lambda }\beta }{\mu (\alpha +\mu )}}>0$.

By taking the derivative with respect to time t in both equations of this system, we obtain the following second order differential equations:

$$\left\{\begin{array}{ccc}\ddot{x}+(\rho +{\mathcal{R}}_{0}y)\dot{x}+{\mathcal{R}}_{0}x\dot{y}\hfill & =& \hfill 0\\ \ddot{y}-{\mathcal{R}}_{0}y\dot{x}+(1-{\mathcal{R}}_{0}x)\dot{y}\hfill & =& \hfill 0\end{array}\right.$$

In order to use the rule of the crossed repeated indices from differential geometry formalism, we change the notations of variables as follows:

$$x=x^1,\dot{x}=y^1,y=x^2,\dot{y}=y^2$$

Then this system of second-order differential equations (SODE) becomes:

$$\left\{\begin{array}{ccc}{\ddot{x}}^1+(\rho +{\mathcal{R}}_0{x}^2)y^1+{\mathcal{R}}_0{x}^1{y}^2\hfill & =& \hfill 0\\ {\ddot{x}}^2-{\mathcal{R}}_0{x}^2{y}^1+(1-{\mathcal{R}}_0{x}^1)y^2\hfill & =& \hfill 0\end{array}\right.$$

(5)

or, equivalently,

$$\left\{\begin{array}{ccc}\frac{d^2{x}^1}{d{t}^2}+(\rho +{\mathcal{R}}_0{x}^2)y^1+{\mathcal{R}}_0{x}^1{y}^2\hfill & =& \hfill 0\\ \frac{d^2{x}^2}{d{t}^2}-{\mathcal{R}}_0{x}^2{y}^1+(1-{\mathcal{R}}_0{x}^1)y^2\hfill & =& \hfill 0\end{array}\right.$$

(6)

where $\frac{d{x}^i}{dt}=y^i$, $i=1,2$.

This system can be written similarly to a SODE (semi-spray) from the KCC theory:

$$\left\{\begin{array}{c}\frac{d^2{x}^1}{d{t}^2}+2G^1(x^1,x^2,y^1,y^2)=0\\ \frac{d^2{x}^2}{d{t}^2}+2G^2(x^1,x^2,y^1,y^2)=0\end{array}\right.$$

(7)

where

$$\begin{array}{ccc}{G}^1(x^i,y^i)\hfill & =& \frac{1}{2}\left[(\rho +{\mathcal{R}}_0{x}^2)y^1+{\mathcal{R}}_0{x}^1{y}^2\right]\hfill \\ G^2(x^i,y^i)\hfill & =& \frac{1}{2}\left[-{\mathcal{R}}_0{x}^2{y}^1+(1-{\mathcal{R}}_0{x}^1)y^2\right]\hfill \end{array}$$

(8)

The zero-connection curvature $Z_j^i=2\frac{\partial G^i}{\partial x^j}$ is given by the next coefficients:

$$\left\{\begin{array}{ccc}\hfill Z_1^1& =& {\mathcal{R}}_0{y}^2\hfill \\ \hfill Z_2^1& =& {\mathcal{R}}_0{y}^1\hfill \\ \hfill Z_1^2& =& -{\mathcal{R}}_0{y}^2\hfill \\ \hfill Z_2^2& =& -{\mathcal{R}}_0{y}^1\hfill \end{array}\right.$$

(9)

The nonlinear connection N has the following coefficients:

$$\left\{\begin{array}{ccc}\hfill N_1^1& =& \frac{\partial G^1}{\partial y^1}=\frac{1}{2}\left(\rho +{\mathcal{R}}_0{x}^2\right)\hfill \\ \hfill N_2^1& =& \frac{\partial G^1}{\partial y^2}=\frac{1}{2}{\mathcal{R}}_0{x}^1\hfill \\ \hfill N_1^2& =& \frac{\partial G^2}{\partial y^1}=-\frac{1}{2}{\mathcal{R}}_0{x}^2\hfill \\ \hfill N_2^2& =& \frac{\partial G^2}{\partial y^2}=\frac{1}{2}\left(1-{\mathcal{R}}_0{x}^1\right)\hfill \end{array}\right.$$

(10)

Consequently, all coefficients of the Berwald connection $G_{jk}^i=\frac{\partial N_j^i}{\partial y^k}$ are null.

The first invariant of the KCC theory ${\epsilon }^i=-\left(N_j^i{y}^j-2G^i\right)$ has the following components:

$$\left\{\begin{array}{ccc}\hfill {\epsilon }^1& =& \frac{1}{2}(\rho +{\mathcal{R}}_0{x}^2)y^1+\frac{1}{2}{\mathcal{R}}_0{x}^1{y}^2\hfill \\ \hfill {\epsilon }^2& =& -\frac{1}{2}{\mathcal{R}}_0{x}^2{y}^1+\frac{1}{2}(1-{\mathcal{R}}_0{x}^1)y^2\hfill \end{array}\right.$$

(11)

Let us remark that ${\epsilon }^i=G^i$ for $i=1,2$, which means that $\frac{\partial G^i}{\partial y^j}{y}^j=1·G^i$ for $i=1,2$, or equivalently, the functions $G^i$ are 1-homogeneous with respect to $y^i$.

Next, by using (A10), we obtain the components of the second invariant of the KCC theory, which means that the deviation curvature tensor of the classical SIR model with demography (4) is as follows:

$$\begin{array}{ccc}\hfill P_1^1& =& -\frac{1}{2}{\mathcal{R}}_0{y}^2+\frac{1}{4}{\left(\rho +{\mathcal{R}}_0{x}^2\right)}^2-\frac{1}{4}{\mathcal{R}}_0^2{x}^1{x}^2\hfill \\ \hfill P_2^1& =& -\frac{1}{2}{\mathcal{R}}_0{y}^1+\frac{1}{4}{\mathcal{R}}_0{x}^1\left(\rho +1+{\mathcal{R}}_0{x}^2-{\mathcal{R}}_0{x}^1\right)\hfill \\ \hfill P_1^2& =& \frac{1}{2}{\mathcal{R}}_0{y}^2-\frac{1}{4}{\mathcal{R}}_0{x}^2\left(\rho +1+{\mathcal{R}}_0{x}^2-{\mathcal{R}}_0{x}^1\right)\hfill \\ \hfill P_2^2& =& \frac{1}{2}{\mathcal{R}}_0{y}^1+\frac{1}{4}{\left(1-{\mathcal{R}}_0{x}^1\right)}^2-\frac{1}{4}{\mathcal{R}}_0^2{x}^1{x}^2\hfill \end{array}$$

(12)

Then, the trace and the determinant of the following deviation curvature matrix:

$$P=\left(\begin{array}{cc}{P}_1^1& P_2^1\\ P_1^2& P_2^2\end{array}\right)$$

are tr $P=P_1^1+P_2^2$ and $detP=P_1^1{P}_2^2-P_1^2{P}_2^1$.

Following Theorem A2, we have:

Theorem 3.1 

All the roots of the characteristic polynomial of P are negative or have negative real parts (i.e., Jacobi stability) if and only if

$$trP=P_1^1+P_2^2<0\mathrm{and}detP=P_1^1{P}_2^2-P_1^2{P}_2^1>0.$$

Taking into account that $P_{jk}^i=\frac{1}{3}\left(\frac{\partial P_j^i}{\partial y^k}-\frac{\partial P_k^i}{\partial y^j}\right)$, $P_{jkl}^i=\frac{\partial P_{jk}^i}{\partial y^l}$, $D_{jkl}^i=\frac{\partial G_{jk}^i}{\partial y^l}$, we obtained the third, fourth, and fifth invariants of the classical SIR model with demography (4) as follows:

Theorem 3.2 

The third KCC invariant $P_{jk}^i$, called the torsion tensor, has all eight components null, i.e.,

$$P_{jk}^i=0\mathrm{for}\mathrm{all}i,j,k.$$

(13)

The fourth KCC invariant $P_{jkl}^i$, called the Riemann–Christoffel curvature tensor, has all sixteen components null, i.e.,

$$P_{jkl}^i=0\mathrm{for}\mathrm{all}i,j,k,l.$$

(14)

The fifth KCC invariant $D_{jkl}^i$, called the Douglas tensor, has all sixteen components null, which means that

$$D_{jkl}^i=0\mathrm{for}\mathrm{all}i,j,k,l.$$

(15)

4. Jacobi Stability Analysis of the Classical SIR System with Demography

In this section, we will determine the first two invariants at the equilibrium points of the SIR model with demography (4), and we will analyze the Jacobi stability of the system near each equilibrium point.

Furthermore, for equilibrium points $E_0(1,0)$ and $E_1\left(\frac{1}{{\mathcal{R}}_0},\rho \left(1-\frac{1}{{\mathcal{R}}_0}\right)\right)$ of the initial SIR model with demography (3), we have the corresponding equilibrium points $E_0(1,0,0,0)$ and $E_1\left(\frac{1}{{\mathcal{R}}_0},\rho \left(1-\frac{1}{{\mathcal{R}}_0}\right),0,0\right)$ for SODE (6).

For $E_0(1,0,0,0)$, the first invariant of the KCC theory has the components ${\epsilon }^1={\epsilon }^2=0$, and the matrix with the components of the second KCC invariant is as follows:

$$P=\left(\begin{array}{cc}\frac{1}{4}{\rho }^2& \frac{1}{4}{\mathcal{R}}_0\left(\rho +1-{\mathcal{R}}_0\right)\\ 0& \frac{1}{4}{\left(1-{\mathcal{R}}_0\right)}^2\end{array}\right).$$

Since $trP=\frac{1}{4}{\rho }^2+\frac{1}{4}{\left(1-{\mathcal{R}}_0\right)}^2>0$ and $detP=\frac{1}{16}{\rho }^2{\left(1-{\mathcal{R}}_0\right)}^2>0$, by using Theorem 3.1, we obtain the following result:

Theorem 4.1 

The disease-free equilibrium point $E_0$ is always Jacobi unstable.

For $E_1\left(\frac{1}{{\mathcal{R}}_0},\rho \left(1-\frac{1}{{\mathcal{R}}_0}\right),0,0\right)$ the first invariant of the KCC theory has the components ${\epsilon }^1=0$, ${\epsilon }^2=0$, and the matrix with the components of the second KCC invariant is as follows:

$$P=\left(\begin{array}{cc}\frac{1}{4}\rho \left(\rho {\mathcal{R}}_0^2-{\mathcal{R}}_0+1\right)& \frac{1}{4}\rho {\mathcal{R}}_0\\ -\frac{1}{4}{\rho }^2{\mathcal{R}}_0\left({\mathcal{R}}_0-1\right)& -\frac{1}{4}\rho \left({\mathcal{R}}_0-1\right)\end{array}\right).$$

Since $trP=\frac{1}{4}\rho \left(\rho {\mathcal{R}}_0^2-2{\mathcal{R}}_0+2\right)$ and $detP=\frac{1}{16}{\rho }^2{\left({\mathcal{R}}_0-1\right)}^2>0$, by using Theorem 3.1, we obtain the following result:

Theorem 4.2 

The endemic equilibrium point $E_1$ is Jacobi stable if and only if $\rho {\mathcal{R}}_0^2-2{\mathcal{R}}_0+2<0$.

Remark 4.3. 

Whenever $E_1$ exists and is Jacobi-stable, a chaotic behavior of the SIR system in a small enough neighborhood of this point is not possible.

In order to more clearly represent the relationship between linear (or Lyapunov) stability and Jacobi stability for this SIR system, we will consider the following diagram with respect to the parameters $\rho =x$ and ${\mathcal{R}}_0=y$ of the system (see Figure 1):

Taking into account that $\rho ={\displaystyle \frac{\mu }{\alpha +\mu }}$, ${\mathcal{R}}_0={\displaystyle \frac{\mathrm{\Lambda }\beta }{\mu (\alpha +\mu )}}$ and according with Theorem 4.2, we obtain the following result:

Theorem 4.4 

If the endemic equilibrium point $E_1$ exists, then it is Jacobi-stable if and only if the reproduction number ${\mathcal{R}}_0$ satisfies the conditions:

$$1+\frac{\alpha }{\mu }-\sqrt{\frac{{\alpha }^2}{{\mu }^2}-1}<{\mathcal{R}}_0<1+\frac{\alpha }{\mu }+\sqrt{\frac{{\alpha }^2}{{\mu }^2}-1}.$$

Proof. 

By using the Jacobi stability condition $\rho {\mathcal{R}}_0^2-2{\mathcal{R}}_0+2<0$ from Theorem 4.2 and the study of the sign of the second order function in ${\mathcal{R}}_0$, $\mu {\mathcal{R}}_0^2-2(\alpha +\mu ){\mathcal{R}}_0+2(\alpha +\mu )$, with the discriminant equal with $4({\alpha }^2-{\mu }^2)$, the theorem is proved. □

According to the expression of the characteristic polynomial at the endemic equilibrium point $E_1$, we get the next result about the local linear stability of the endemic equilibrium.

Theorem 4.5 

If exists, the endemic equilibrium $E_1$, is a stable focus if and only if the reproduction number ${\mathcal{R}}_0$ satisfies the conditions:

$$2\left(1+\frac{\alpha }{\mu }-\sqrt{\frac{{\alpha }^2}{{\mu }^2}+\frac{\alpha }{\mu }}\right)<{\mathcal{R}}_0<2\left(1+\frac{\alpha }{\mu }+\sqrt{\frac{{\alpha }^2}{{\mu }^2}+\frac{\alpha }{\mu }}\right).$$

Otherwise, the endemic equilibrium $E_1$ is a stable node.

Proof. 

Because the discriminant associated to the characteristic polynomial at $E_1$ is equal with $\rho (\rho {\mathcal{R}}_0^2-4{\mathcal{R}}_0+4)$, is necessary to study the sign of the second order function in ${\mathcal{R}}_0$, $\mu {\mathcal{R}}_0^2-4(\alpha +\mu ){\mathcal{R}}_0+4(\alpha +\mu )$. □

Let us remark that apart from the reproduction number ${\mathcal{R}}_0$, the parameter $\alpha$ is the most important parameter of this epidemic model. Only if $\alpha >\mu$ (i.e. $\rho <\frac{1}{2}$), the endemic equilibrium (if exists) can be stable from the Jacobi point of view. Therefore, recovery rate from illness $\alpha$ plays an unexpectedly crucial role both for Lyapunov stability and for Jacobi stability of this classical SIR system.

In order to obtain characterizations of the local stability of the endemic equilibrium point $E_1$ related to the parameter $\beta$ - the transmission coefficient of the infection, the next results are obtained:

Theorem 4.6 

If the endemic equilibrium point $E_1$ exists, then it is Jacobi-stable if and only if the transmission coefficient β satisfies the conditions:

$${\displaystyle \frac{\alpha +\mu }{{\mathrm{\Lambda }}_{\infty }}}\left(1+\frac{\alpha }{\mu }-\sqrt{\frac{{\alpha }^2}{{\mu }^2}-1}\right)<\beta <{\displaystyle \frac{\alpha +\mu }{{\mathrm{\Lambda }}_{\infty }}}\left(1+\frac{\alpha }{\mu }+\sqrt{\frac{{\alpha }^2}{{\mu }^2}-1}\right),$$

where ${\mathrm{\Lambda }}_{\infty }={\displaystyle \frac{\mathsf{\Lambda }}{\mu }}$ is the limit size of the entire population.

Proof. 

By replacing both dimensionless parameters $\rho$ and ${\mathcal{R}}_0$ with the four parameters of the initial system (1), and using the Jacobi stability condition from Theorem 4.2, the theorem is proved by following the sign of the second order function in $\beta$, ${\mathrm{\Lambda }}^2{\beta }^2-2\mathrm{\Lambda }{(\alpha +\mu )}^2\beta +2\mu {(\alpha +\mu )}^3$, with the discriminant equal with $4{\mathrm{\Lambda }}^2{(\alpha +\mu )}^2({\alpha }^2-{\mu }^2)$. □

Theorem 4.7 

If exists, the endemic equilibrium $E_1$, is a stable focus if and only if the transmission coefficient β satisfies the conditions:

$${\displaystyle \frac{2(\alpha +\mu )}{{\mathrm{\Lambda }}_{\infty }}}\left(1+\frac{\alpha }{\mu }-\sqrt{\frac{{\alpha }^2}{{\mu }^2}+\frac{\alpha }{\mu }}\right)<\beta <{\displaystyle \frac{2(\alpha +\mu )}{{\mathrm{\Lambda }}_{\infty }}}\left(1+\frac{\alpha }{\mu }+2\sqrt{\frac{{\alpha }^2}{{\mu }^2}+\frac{\alpha }{\mu }}\right),$$

where ${\mathrm{\Lambda }}_{\infty }={\displaystyle \frac{\mathrm{\Lambda }}{\mu }}$ is the limit size of the entire population.

Otherwise, the endemic equilibrium $E_1$ is a stable node.

Proof. 

Because the discriminant associated to the characteristic polynomial at $E_1$ is equal with $\rho (\rho {\mathcal{R}}_0^2-4{\mathcal{R}}_0+4)$, is necessary to study the sign of the second order function in $\beta$, ${\mathrm{\Lambda }}^2{\beta }^2-4\mathrm{\Lambda }{(\alpha +\mu )}^2\beta +4\mu {(\alpha +\mu )}^3$. □

If the mortality rate is known at a value ${\mu }_0$, in order to obtain threshold values of the parameter $\alpha$, the recovery rate, it is enough to replace $\rho$ with ${\displaystyle \frac{{\mu }_0}{\alpha +{\mu }_0}}$ in Theorem 4.2 and it results:

Theorem 4.8 

a) If the endemic equilibrium point $E_1$ exists, then it is Jacobi-stable if and only if the recovery rate α fulfills the condition:

$$\alpha >{\alpha }_J={\mu }_0{\displaystyle \frac{{({\mathcal{R}}_0-1)}^2+1}{2({\mathcal{R}}_0-1)}}.$$

b) If exists, the endemic equilibrium $E_1$, is a stable focus if and only if the recovery rate α fulfills the condition:

$$\alpha >{\alpha }_F={\mu }_0{\displaystyle \frac{{({\mathcal{R}}_0-2)}^2}{4({\mathcal{R}}_0-1)}}.$$

Otherwise, the endemic equilibrium $E_1$ is a stable node.

Let us remark that ${\alpha }_J=2{\alpha }_F+{\mu }_0$, i.e. ${\alpha }_J>{\alpha }_F$ because the Jacobi stability of an equilibrium point implies that this equilibrium is a focus.

Conversely, if we fixed the value of the recovery rate $\alpha$ at the value ${\alpha }_0$, then we obtain the following threshold values for the mortality rate $\mu$:

Theorem 4.9 

a) If the endemic equilibrium point $E_1$ exists, then it is Jacobi-stable if and only if the mortality rate μ fulfills the condition:

$$\mu <{\mu }_J={\alpha }_0{\displaystyle \frac{2({\mathcal{R}}_0-1)}{{({\mathcal{R}}_0-1)}^2+1}}.$$

b) If exists, the endemic equilibrium $E_1$, is a stable focus if and only if the mortality rate μ fulfills the condition:

$$\alpha <{\mu }_F={\alpha }_0{\displaystyle \frac{4({\mathcal{R}}_0-1)}{{({\mathcal{R}}_0-2)}^2}}.$$

Otherwise, the endemic equilibrium $E_1$ is a stable node.

Of course, ${\mu }_J<{\mu }_F$ because ${\mu }_F-{\mu }_J=2{\alpha }_0{\displaystyle \frac{{\mathcal{R}}_0^2({\mathcal{R}}_0-1)}{{({\mathcal{R}}_0-2)}^2\left({({\mathcal{R}}_0-1)}^2+1\right)}}>0$.

4.1. Dynamics of the Deviation Vector for the Classical SIR System with Demography

The behavior of the deviation vector ${\xi }^i$, $i=1,2$, giving the trajectory behavior of the dynamical system near an equilibrium point is described by the deviation equations (or Jacobi equations) (A8), or in a covariant form such as equations (A9).

For this classical SIR system with demography, the deviation equations become:

$$\left\{\begin{array}{ccc}\frac{d^2{\xi }^1}{d{t}^2}+\left(\rho +{\mathcal{R}}_0{x}^2\right)\frac{d{\xi }^1}{dt}+{\mathcal{R}}_0{x}^1\frac{d{\xi }^2}{dt}+{\mathcal{R}}_0{y}^2{\xi }^1+{\mathcal{R}}_0{y}^1{\xi }^2\hfill & =& \hfill 0\\ \frac{d^2{\xi }^2}{d{t}^2}-{\mathcal{R}}_0{x}^2\frac{d{\xi }^1}{dt}+\left(1-{\mathcal{R}}_0{x}^1\right)\frac{d{\xi }^2}{dt}-{\mathcal{R}}_0{y}^2{\xi }^1-{\mathcal{R}}_0{y}^1{\xi }^2\hfill & =& \hfill 0\end{array}\right.$$

(16)

The length of the deviation vector $\xi \left(t\right)=\left({\xi }^1\left(t\right),{\xi }^2\left(t\right)\right)$ is given by

$$∥\xi \left(t\right)∥=\sqrt{{\left({\xi }^1\left(t\right)\right)}^2+{\left({\xi }^2\left(t\right)\right)}^2}.$$

Next, we will write the deviation equations near to each equilibrium point for the classical SIR system with demography. Then, the dynamics of the deviation vector near to the disease-free equilibrium point $E_0(1,0,0,0)$ is described by the following SODE:

$$\left\{\begin{array}{ccc}\frac{d^2{\xi }^1}{d{t}^2}+\rho \frac{d{\xi }^1}{dt}+{\mathcal{R}}_0\frac{d{\xi }^2}{dt}\hfill & =& \hfill 0\\ \frac{d^2{\xi }^2}{d{t}^2}+\left(1-{\mathcal{R}}_0\right)\frac{d{\xi }^2}{dt}\hfill & =& \hfill 0\end{array}\right.$$

(17)

The dynamics of the deviation vector near to the endemic equilibrium point

$E_1\left(\frac{1}{{\mathcal{R}}_0},\rho \left(1-\frac{1}{{\mathcal{R}}_0}\right),0,0\right)$ is described by the deviation equations:

$$\left\{\begin{array}{ccc}\frac{d^2{\xi }^1}{d{t}^2}+\rho {\mathcal{R}}_0\frac{d{\xi }^1}{dt}+\frac{d{\xi }^2}{dt}\hfill & =& \hfill 0\\ \frac{d^2{\xi }^2}{d{t}^2}+\rho \left(1-{\mathcal{R}}_0\right)\frac{d{\xi }^1}{dt}\hfill & =& \hfill 0\end{array}\right.$$

(18)

According to the standard approach used in differential geometry of plane curves [19], the curvature $\kappa \left(t\right)$ of the integral curve $\xi \left(t\right)=\left({\xi }^1\left(t\right),{\xi }^2\left(t\right)\right)$ corresponding to deviation Equation (16) represents a quantitative description of the behavior of the deviation vector and is given by the following:

$$\kappa \left(t\right)=\frac{{\dot{\xi }}^1\left(t\right){\ddot{\xi }}^2\left(t\right)-{\ddot{\xi }}^1\left(t\right){\dot{\xi }}^2\left(t\right)}{{\left[{\left({\dot{\xi }}^1\left(t\right)\right)}^2+{\left({\dot{\xi }}^2\left(t\right)\right)}^2\right]}^{3/2}}$$

(19)

where ${\dot{\xi }}^i\left(t\right)=\frac{d{\xi }^i}{dt}$, ${\ddot{\xi }}^i\left(t\right)=\frac{d^2{\xi }^i}{d{t}^2}$, $i=1,2$.

5. A Simple SIR Epidemic Pattern with Demography and Vaccination

In this section we propose a classical and simple SIR model with vaccination as in [37]. Also, a very interesting four dimensional model with vaccination was studied recently in [38]. Next, we will reconsider the classical SIR model with demography (1) and we will assume that a part of the susceptible individuals are vaccinated with the rate of vaccination p. According to the classical SIR model this vaccination rate is measured in [unit of time]${}^{-1}$. If we will suppose that all vaccinated individuals will not be infected and more that, they can be considered to belongs to the catogory of removed individuals, then the model (1) becomes

$$\left\{\begin{array}{ccc}\hfill S^{\prime }& =& \mathrm{\Lambda }-\beta IS-\mu S-pS\hfill \\ \hfill I^{\prime }& =& \beta IS-\alpha I-\mu I\hfill \\ \hfill R^{\prime }& =& \alpha I-\mu R+pS\hfill \end{array}\right.$$

(20)

Let us observe that the first two equations of the system (20) are independent of the third equation, because we can obtain the removed population by using $R\left(t\right)=N\left(t\right)-S\left(t\right)-I\left(t\right)$ and due to the fact that the total population $N\left(t\right)$ can be obtained from the differential equation

$$N^{\prime }\left(t\right)=\mathrm{\Lambda }-\mu N\left(t\right),$$

with the unique solution $N\left(t\right)=N\left(0\right)e^{-\mu t}+{\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\left(1-e^{-\mu t}\right)$.

Obviously, as for the classical model (1), the population size is not constant, but it is asymptotically constant, since $N\left(t\right)\to {\displaystyle \frac{\mathrm{\Lambda }}{\mu }}$ when $t\to \infty$.

Taking into account all these proposals, we can consider the two-dimensional autonomous system of first order differential equations:

$$\left\{\begin{array}{ccc}\hfill S^{\prime }& =& \mathrm{\Lambda }-\beta IS-(\mu +p)S\hfill \\ \hfill I^{\prime }& =& \beta IS-(\alpha +\mu )I\hfill \end{array}\right.$$

(21)

By applying similar techniques as in Section 2 and by changing of time and variables, $\tau =(\alpha +\mu )t$, $x\left(\tau \right)=\frac{p+\mu }{\mathrm{\Lambda }}\widehat{S}\left(\tau \right)$, $y\left(\tau \right)=\frac{p+\mu }{\mathrm{\Lambda }}\widehat{I}\left(\tau \right)$, where $\widehat{S}\left(\tau \right)=S\left(\frac{\tau }{\alpha +\mu }\right)$, $\widehat{I}\left(\tau \right)=I\left(\frac{\tau }{\alpha +\mu }\right)$, we obtain the equivalent dimensionless form of (21):

$$\left\{\begin{array}{ccc}\hfill x^{\prime }& =& \rho (1-x)-{\mathcal{R}}_{0}xy\hfill \\ \hfill y^{\prime }& =& {\mathcal{R}}_{0}xy-y\hfill \end{array}\right.$$

(22)

where $\rho ={\displaystyle \frac{p+\mu }{\alpha +\mu }}$ and ${\mathcal{R}}_0={\displaystyle \frac{\mathrm{\Lambda }\beta }{(p+\mu )(\alpha +\mu )}}$ are both dimensionless parameters.

Let us remark that system (22) is exactly the dimensionless system (3) obtained for the classical SIR model in Section 2. Therefore, all obtained results remain available also for this SIR epidemic model with demography and vaccination. Only interpretations can be different because the dimensionless parameters $\rho$ and ${\mathcal{R}}_0$ are different. Of course, ${\mathcal{R}}_0$ is the basic reproduction number for this model. Moreover, the number of parameters was reduced from five to two.

According to the expressions of $\rho >0$ and ${\mathcal{R}}_0>0$, it is clear that $\rho <1$ if and only if $\alpha >p$ and $\rho <\frac{1}{2}$ if and only if $\alpha >2p+\mu$. So, if we take account of the results presented in Figure 1, for this model the endemic equilibrium $E_1$ can be a stable focus only if the recovery rate $\alpha$ is greater than the vaccination rate p and $E_1$ can be Jacobi stable only if the recovery rate $\alpha$ is greater than $2p+\mu$.

6. A Modified SIR Epidemic Pattern with Demography

Taking into account that for the classical SIR model (1) we have no periodic orbits, and then Hopf bifurcations don’t occur, next we will consider a modified SIR system. More exactly, we will assume that the transmission coefficient of infection $\beta$ is not constant, but is a linear function of the number of infected individuals, that means $\beta$ is replaced by $\beta (1+\nu I)$, where $\nu >0$ is a parameter [1,39]. This imply that the transmission coefficient of infection and, also, the contact rate increase with the number of infectious individuals and new infections occur much faster in comparison with the classical model [1].

Therefore, the model becomes:

$$\left\{\begin{array}{ccc}\hfill S^{\prime }& =& \mathrm{\Lambda }-\beta (1+\nu I)IS-\mu S\hfill \\ \hfill I^{\prime }& =& \beta (1+\nu I)IS-(\alpha +\mu )I\hfill \\ \hfill R^{\prime }& =& \alpha I-\mu R\hfill \end{array}\right.$$

(23)

Let us remark that the total population size $N=S+I+R$ satisfies $N^{\prime }\left(t\right)=\mathrm{\Lambda }-\mu N\left(t\right)$ and then we can omit the third equation for recovered individuals R. It follow that it is necessary and sufficient to study the two-dimensional autonomous system of first order differential equations:

$$\left\{\begin{array}{ccc}\hfill S^{\prime }& =& \mathrm{\Lambda }-\beta (1+\nu I)IS-\mu S\hfill \\ \hfill I^{\prime }& =& \beta (1+\nu I)IS-(\alpha +\mu )I\hfill \end{array}\right.$$

(24)

where $R\left(t\right)=N\left(t\right)-S\left(t\right)-I\left(t\right)$.

In order to obtain the dimensionless version of the system, using similar techniques as for the classical system (2), after changing the time variable t by $\tau =(\alpha +\mu )t$ and rescaling the variables $\widehat{S}=S\left(\frac{\tau }{\alpha +\mu }\right)$, $\widehat{I}=I\left(\frac{\tau }{\alpha +\mu }\right)$ with the total limiting population size ${\mathrm{\Lambda }}_{\infty }=\frac{\mathrm{\Lambda }}{\mu }$, we obtain the new variables $x\left(\tau \right)=\frac{\widehat{S}}{{\mathrm{\Lambda }}_{\infty }}$, $y\left(\tau \right)=\frac{\widehat{I}}{{\mathrm{\Lambda }}_{\infty }}$, which are dimensionless quantities. Consequently, if we denote $\omega =\nu {\mathrm{\Lambda }}_{\infty }$, then the system (24) becomes:

$$\left\{\begin{array}{ccc}\hfill x^{\prime }& =& \rho (1-x)-{\mathcal{R}}_0(1+\omega y)xy\hfill \\ \hfill y^{\prime }& =& {\mathcal{R}}_0(1+\omega y)xy-y\hfill \end{array}\right.$$

(25)

where $\rho ={\displaystyle \frac{\mu }{\alpha +\mu }}$ and ${\mathcal{R}}_0={\displaystyle \frac{\mathrm{\Lambda }\beta }{\mu (\alpha +\mu )}}$ are both dimensionless parameters.

Therefore, we can say that the system (24) was transformed the system into a dimensionless form (25), which is equivalent to the original one, since the solutions of both systems have the same long-term behavior. Moreover, the number of parameters was reduced from five ($\mathrm{\Lambda }$, $\beta$, $\nu$, $\mu$, $\alpha$) to three ($\rho$, ${\mathcal{R}}_0$, $\omega$), all being strictly positive real numbers. Of course, ${\mathcal{R}}_0$ represent the reproduction number for this modified SIR model and $0<\rho <1$.

Further, our study for this modified SIR model has relevance when $x\ge 0$, $y\ge 0$, i.e. on the first quadrant ${\Sigma }_{+}^0=\left\{(x,y)\in {\mathbf{R}}^2\mid x\ge 0,y\ge 0\right\}$ or even on the open first quadrant ${\Sigma }_{+}=\left\{(x,y)\in {\mathbf{R}}^2\mid x>0,y>0\right\}$.

The Jacobian matrix of the system (25) at a point $(x,y)$ is

$$A=\left(\begin{array}{cc}-\rho -{\mathcal{R}}_0(1+\omega y)y& -{\mathcal{R}}_{0}x-2{\mathcal{R}}_0\omega xy\\ {\mathcal{R}}_0(1+\omega y)y& {\mathcal{R}}_{0}x+2{\mathcal{R}}_0\omega xy-1\end{array}\right)$$

In order to find the equilibrium points of (25), by analyzing the system,

$$\left\{\begin{array}{ccc}\rho (1-x)-{\mathcal{R}}_0(1+\omega y)xy\hfill & =& \hfill 0\\ {\mathcal{R}}_0(1+\omega y)xy-y\hfill & =& \hfill 0\end{array}\right.$$

we obtain at most three equilibrium, as follows:

For $y=0$, it results $E_0(1,0)$, the disease-free equilibrium point, with Jacobian $A=\left(\begin{array}{cc}-\rho & -{\mathcal{R}}_0\\ 0& {\mathcal{R}}_0-1\end{array}\right)$ and eigenvalues ${\lambda }_1=-\rho$, ${\lambda }_2={\mathcal{R}}_0-1$. Then $E_0$ is local asymptotically stable (stable node) if and only if ${\mathcal{R}}_0<1$, or $E_0$ is unstable (saddle point) if and only if ${\mathcal{R}}_0>1$. For ${\mathcal{R}}_0=1$, $E_0$ is a non-hyperbolic equilibrium point because ${\lambda }_2=0$.

For $y\ne 0$, we have ${\mathcal{R}}_0(1+\omega y)x=1$, and then $\rho \left(1-\frac{1}{{\mathcal{R}}_0\left(1+\omega y\right)}\right)-y=0$. Because $\rho \left(1-\frac{1}{{\mathcal{R}}_0\left(1+\omega y\right)}\right)-y=\frac{\rho {\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega y-\rho -{\mathcal{R}}_{0}y-{\mathcal{R}}_0{y}^2\omega }{R\left(1+\omega y\right)}$, it results the roots of second order algebraic equation in y

$$-\omega {\mathcal{R}}_0{y}^2+{\mathcal{R}}_0(\omega \rho -1)y+\rho ({\mathcal{R}}_0-1)=0,$$

(26)

$$y_1=\frac{1}{2{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0(\omega \rho -1)-\sqrt{\mathrm{\Delta }}\right),y_2=\frac{1}{2{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right),$$

where $\mathrm{\Delta }={\mathcal{R}}_0^2{(\omega \rho +1)}^2-4{\mathcal{R}}_0\omega \rho$ is the discriminant of the second order equation (26).

If we denote by $f\left(y\right)=-\omega {\mathcal{R}}_0{y}^2+{\mathcal{R}}_0(\omega \rho -1)y+\rho ({\mathcal{R}}_0-1)$, the second order function associated to the equation (26), then we have the following three cases:

Case 1. If $f\left(0\right)=\rho ({\mathcal{R}}_0-1)>0$, i.e. ${\mathcal{R}}_0>1$, then only the second root $y_2$ is strictly positive and we have only one endemic equilibrium point $E\left(x,y\right)$, with coordinates

$$x=\frac{2}{{\mathcal{R}}_0(1+\omega \rho )+\sqrt{\mathrm{\Delta }}},y=\frac{1}{2{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right).$$

(27)

Let us remark that $\mathrm{\Delta }>0$, because $\mathrm{\Delta }={\mathcal{R}}_0^2{(\omega \rho +1)}^2-4{\mathcal{R}}_0\omega \rho >{\mathcal{R}}_0^2\left({(\omega \rho +1)}^2-4\omega \rho \right)>{(\omega \rho +1)}^2-4\omega \rho ={(\omega \rho -1)}^2$, for any positive parameters $\omega$ and $\rho$. Then, taking into account that ${\mathcal{R}}_0(1+\omega y)x=1$, the Jacobian at the unique endemic equilibrium $E\left(x,y\right)$ is

$$A=\left(\begin{array}{cc}-\rho -{\mathcal{R}}_0\left(1+\omega y\right)y& -{\mathcal{R}}_0\omega xy-1\\ {\mathcal{R}}_0\left(1+\omega y\right)y& {\mathcal{R}}_0\omega xy\end{array}\right),$$

(28)

where the trace is $\mathrm{tr}A={\lambda }_1+{\lambda }_2=-\rho -{\mathcal{R}}_0\left(1+\omega y\right)y+{\mathcal{R}}_0\omega xy$ and the determinant is $detA={\lambda }_1{\lambda }_2=-\rho {\mathcal{R}}_0\omega xy+{\mathcal{R}}_0\left(1+\omega y\right)y$, where ${\lambda }_1$, ${\lambda }_2$ are eigenvalues of A.

By replacing x and y from (27), it results

$$\mathrm{tr}A=-\frac{1}{4}\frac{\left({\left(\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}^2+{\mathcal{R}}_0(4\rho \omega -{\mathcal{R}}_0)\right)\left({\mathcal{R}}_0+\rho R\omega +\sqrt{\mathrm{\Delta }}\right)-4R\omega \left(\rho {\mathcal{R}}_0\omega -{\mathcal{R}}_0+\sqrt{\mathrm{\Delta }}\right)}{{\mathcal{R}}_0\omega \left({\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}$$

and

$$detA=\frac{1}{4}\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right)\frac{{\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)+\sqrt{\mathrm{\Delta }}\right)}{{\mathcal{R}}_0\omega \left({\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}.$$

Due to $y>0$ and ${\mathcal{R}}_0>1$, it results $detA>0$ for any positive values of parameters. Therefore, the endemic equilibrium $E\left(x,y\right)$ is local asymptotically stable (stable node or stable focus) if and only if $\mathrm{tr}A<0$, i.e.

$$\left({\left(\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}^2+{\mathcal{R}}_0(4\rho \omega -{\mathcal{R}}_0)\right)\left({\mathcal{R}}_0+\rho R\omega +\sqrt{\mathrm{\Delta }}\right)-4R\omega \left(\rho {\mathcal{R}}_0\omega -{\mathcal{R}}_0+\sqrt{\mathrm{\Delta }}\right)>0.$$

Otherwise, the endemic equilibrium E is unstable (unstable node or unstable focus) if and only if $\mathrm{tr}A>0$, i.e.

$$\left({\left(\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}^2+{\mathcal{R}}_0(4\rho \omega -{\mathcal{R}}_0)\right)\left({\mathcal{R}}_0+\rho R\omega +\sqrt{\mathrm{\Delta }}\right)-4R\omega \left(\rho {\mathcal{R}}_0\omega -{\mathcal{R}}_0+\sqrt{\mathrm{\Delta }}\right)<0.$$

Let us point out that, if $\mathrm{tr}A=0$ (i.e. ${\lambda }_1=-{\lambda }_2$), then ${\lambda }_{1,2}$ are pure imaginary roots of the characteristic polynomial at E, with Re ${\lambda }_{1,2}=0$. In this case it is possible to have Hopf bifurcations along the curve tr $A=0$, i.e

$$\left({\left(\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}^2+{\mathcal{R}}_0(4\rho \omega -{\mathcal{R}}_0)\right)\left({\mathcal{R}}_0+\rho R\omega +\sqrt{\mathrm{\Delta }}\right)-4R\omega \left(\rho {\mathcal{R}}_0\omega -{\mathcal{R}}_0+\sqrt{\mathrm{\Delta }}\right)=0.$$

Case 2. If $f\left(0\right)=\rho ({\mathcal{R}}_0-1)=0$, i.e. ${\mathcal{R}}_0=1$, then $\mathrm{\Delta }={(\omega \rho -1)}^2\ge 0$ and the roots of the equation (26) are $y_1=0$ and $y_2=\frac{\omega \rho -1}{\omega }$. Then we have that the first endemic equilibrium $E_1\left(x_1,y_1\right)$ coincides with $E_0(1,0)$, and so, this equilibrium point is non hyperbolic, but the second endemic equilibrium becomes $E_2\left(x_2,y_2\right)$, with coordinates $x_2=\frac{1}{\omega \rho }$, $y_2=\frac{\omega \rho -1}{\omega }$.

Obviously, $E_2\in {\Sigma }_{+}$ if and only if $\omega \rho -1>0$. The Jacobian at $E_2$ is $A=\left(\begin{array}{cc}-{\rho }^2\omega & -\frac{2\rho \omega -1}{\rho \omega }\\ \left(\rho \omega -1\right)\rho & \frac{\rho \omega -1}{\rho \omega }\end{array}\right)$, with eigenvalues ${\lambda }_{1,2}=\frac{1}{2\rho \omega }\left(\rho \omega -1-{\rho }^3{\omega }^2\pm \sqrt{{\mathrm{\Delta }}_1}\right)$, where ${\mathrm{\Delta }}_1={\left({\rho }^3{\omega }^2-\rho \omega +1\right)}^2-4{\rho }^2\omega {\left(\rho \omega -1\right)}^2$.

Tacking into account that $\mathrm{tr}A={\lambda }_1+{\lambda }_2=-\frac{{\rho }^3{\omega }^2-\rho \omega +1}{\rho \omega }$ and $detA={\lambda }_1{\lambda }_2=\frac{{\left(\rho \omega -1\right)}^2}{\omega }>0$, it results that $E_2$ is local asymptotically stable (stable node or stable focus) if and only if $\mathrm{tr}A<0$, i.e. ${\rho }^3{\omega }^2-\rho \omega +1>0$. Else, $E_2$ is unstable (unstable node or unstable focus) if and only if $\mathrm{tr}A>0$, i.e. ${\rho }^3{\omega }^2-\rho \omega +1<0$.

Let us point out that, if $\mathrm{tr}A=0$ (i.e. ${\lambda }_1=-{\lambda }_2$), then ${\mathrm{\Delta }}_1=-4{\rho }^2\omega {\left(\rho \omega -1\right)}^2<0$, i.e. ${\lambda }_{1,2}$ are pure imaginary roots with Re ${\lambda }_{1,2}=0$. In this case it is possible to have Hopf bifurcations along the curve ${\rho }^3{\omega }^2-\rho \omega +1=0$. More exactly, for $\omega =\frac{1\pm \sqrt{1-4\rho }}{2{\rho }^2}$ and $0<\rho <\frac{1}{4}$.

If $\frac{1}{4}<\rho <1$, then ${\rho }^3{\omega }^2-\rho \omega +1>0$. For $\rho =\frac{1}{4}$, ${\rho }^3{\omega }^2-\rho \omega +1>0$ for any $\omega \ne \frac{1}{8}$ and ${\rho }^3{\omega }^2-\rho \omega +1=0$ only for $\omega =\frac{1}{8}$.

Case 3. If $f\left(0\right)=\rho ({\mathcal{R}}_0-1)<0$, i.e. ${\mathcal{R}}_0<1$, then we have two endemic equilibrium points $E_1\left(x_1,y_1\right)$ and $E_2\left(x_2,y_2\right)$, with positive coordinates given by

$$x=\frac{1}{{\mathcal{R}}_0(1+\omega y)}and-\omega {\mathcal{R}}_0{y}^2+{\mathcal{R}}_0(\omega \rho -1)y+\rho ({\mathcal{R}}_0-1)=0$$

if and only if $\mathrm{\Delta }={\mathcal{R}}_0^2{(\omega \rho +1)}^2-4{\mathcal{R}}_0\omega \rho >0$, that means $\frac{4\omega \rho }{{(\omega \rho +1)}^2}<{\mathcal{R}}_0<1$.

More exactly, this two equilibrium points has the following coordinates:

$$x_1=\frac{2}{{\mathcal{R}}_0(1+\omega \rho )-\sqrt{\mathrm{\Delta }}},y_1=\frac{1}{2{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0(\omega \rho -1)-\sqrt{\mathrm{\Delta }}\right)$$

(29)

and, respectively,

$$x_2=\frac{2}{{\mathcal{R}}_0(1+\omega \rho )+\sqrt{\mathrm{\Delta }}},y_2=\frac{1}{2{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right).$$

(30)

Let us remark that $\omega \rho -1>0$, because $\frac{y_1+y_2}{2}=\frac{\omega \rho -1}{2\omega }>0$.

For the first endemic equilibrium $E_1\left(x_1,y_1\right)$, taking into account that ${\mathcal{R}}_0(1+\omega y_1)x_1=1$, it results that the Jacobian matrix is

$$A=\left(\begin{array}{cc}-\rho -{\mathcal{R}}_0\left(1+\omega y_1\right)y_1& -{\mathcal{R}}_0\omega x_1{y}_1-1\\ {\mathcal{R}}_0\left(1+\omega y_1\right)y_1& {\mathcal{R}}_0\omega x_1{y}_1\end{array}\right)$$

with $\mathrm{tr}A={\lambda }_1+{\lambda }_2=-\rho -{\mathcal{R}}_0\left(1+\omega y_1\right)y_1+{\mathcal{R}}_0\omega x_1{y}_1$ and $detA={\lambda }_1{\lambda }_2=-{\mathcal{R}}_0\omega x_1{y}_1\rho +{\mathcal{R}}_0\left(1+\omega y_1\right)y_1$, where ${\lambda }_1$, ${\lambda }_2$ are the eigenvalues of A.

By replacing $x_1$ and $y_1$ from (29), it results that

$$\mathrm{tr}A=-\frac{1}{4}\frac{\left({\left(\rho {\mathcal{R}}_0\omega -\sqrt{\mathrm{\Delta }}\right)}^2+{\mathcal{R}}_0(4\rho \omega -{\mathcal{R}}_0)\right)\left({\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega -\sqrt{\mathrm{\Delta }}\right)-4{\mathcal{R}}_0\omega \left(\rho \omega {\mathcal{R}}_0-{\mathcal{R}}_0-\sqrt{\mathrm{\Delta }}\right)}{{\mathcal{R}}_0\omega \left({\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega -\sqrt{\mathrm{\Delta }}\right)}$$

and

$$detA=\frac{1}{4}\left({\mathcal{R}}_0(\omega \rho -1)-\sqrt{\mathrm{\Delta }}\right)\frac{{\left({\mathcal{R}}_0(\omega \rho -1)-\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)-\sqrt{\mathrm{\Delta }}\right)}{{\mathcal{R}}_0\omega \left({\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega -\sqrt{\mathrm{\Delta }}\right)}.$$

Because $y_1=\frac{1}{2{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0(\omega \rho -1)-\sqrt{\mathrm{\Delta }}\right)>0$ and ${\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega -\sqrt{\mathrm{\Delta }}>0$, it results that the sign of $detA$ is given by ${\left({\mathcal{R}}_0(\omega \rho -1)-\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)-\sqrt{\mathrm{\Delta }}\right)$. Since $\frac{4\omega \rho }{{(\omega \rho +1)}^2}<{\mathcal{R}}_0<1$, we have $\rho \omega ({\mathcal{R}}_0-1)-\sqrt{\mathrm{\Delta }}<0$ and then the expression ${\left({\mathcal{R}}_0(\omega \rho -1)-\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)-\sqrt{\mathrm{\Delta }}\right)$ can be negative. Indeed, if we use the identity

$$\begin{array}{ccc}\hfill {\left({\mathcal{R}}_0(\omega \rho -1)-\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)-\sqrt{\mathrm{\Delta }}\right)& =& {\left(1+\rho \omega \right)}^2\left({\mathcal{R}}_0-\frac{4\omega \rho }{{(\omega \rho +1)}^2}\right)\left({\mathcal{R}}_0-\frac{2\sqrt{\mathrm{\Delta }}}{\omega \rho +1}\right)\hfill \\ & +& \mathrm{\Delta }-\frac{8\rho \omega \sqrt{\mathrm{\Delta }}}{1+\omega \rho }\hfill \end{array}$$

and we observe that ${\mathcal{R}}_0-\frac{2\sqrt{\mathrm{\Delta }}}{\omega \rho +1}<0$$\left(dueto{\mathcal{R}}_0>\frac{4\omega \rho }{{(\omega \rho +1)}^2}\right)$, $\mathrm{\Delta }-\frac{8\rho \omega \sqrt{\mathrm{\Delta }}}{1+\omega \rho }<0$ for ${\mathcal{R}}_0<1$ and $\omega \rho >1$ (because ${(\omega \rho +1)}^2<64{\omega }^2{\rho }^2+4\omega \rho$ for $\omega \rho >1$), it result that

$${\left({\mathcal{R}}_0(\omega \rho -1)-\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)-\sqrt{\mathrm{\Delta }}\right)<0\mathrm{for}\mathrm{any}{\mathcal{R}}_0\in \left(\frac{4\omega \rho }{{(\omega \rho +1)}^2},1\right).$$

Therefore, the first endemic equilibrium $E_1$ is unstable (a saddle point) because $detA<0$.

For the second endemic equilibrium $E_2\left(x_2,y_2\right)$, taking into account that ${\mathcal{R}}_0(1+\omega y_2)x_2=1$, the Jacobian matrix is

$$A=\left(\begin{array}{cc}-\rho -{\mathcal{R}}_0\left(1+\omega y_2\right)y_2& -{\mathcal{R}}_0\omega x_2{y}_2-1\\ {\mathcal{R}}_0\left(1+\omega y_2\right)y_2& {\mathcal{R}}_0\omega x_2{y}_2\end{array}\right)$$

with $\mathrm{tr}A={\lambda }_1+{\lambda }_2=-\rho -{\mathcal{R}}_0\left(1+\omega y_2\right)y_2+{\mathcal{R}}_0\omega x_2{y}_2$ and $detA={\lambda }_1{\lambda }_2=-{\mathcal{R}}_0\omega x_2{y}_2\rho +{\mathcal{R}}_0\left(1+\omega y_2\right)y_2$, where ${\lambda }_1$, ${\lambda }_2$ are the eigenvalues of A.

By replacing $x_2$ and $y_2$ from (30), it results that

$$\mathrm{tr}A=-\frac{1}{4}\frac{\left({\left(\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}^2+{\mathcal{R}}_0(4\rho \omega -{\mathcal{R}}_0)\right)\left({\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)-4{\mathcal{R}}_0\omega \left(\rho {\mathcal{R}}_0\omega -{\mathcal{R}}_0+\sqrt{\mathrm{\Delta }}\right)}{{\mathcal{R}}_0\omega \left({\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}$$

and

$$detA=\frac{1}{4}\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right)\frac{{\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)+\sqrt{\mathrm{\Delta }}\right)}{{\mathcal{R}}_0\omega \left({\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}.$$

Since $y_2=\frac{1}{2{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right)>0$ and ${\mathcal{R}}_0+\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}>0$, we have that the sign of $detA$ is given by ${\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)+\sqrt{\mathrm{\Delta }}\right)$. From $\frac{4\omega \rho }{{(\omega \rho +1)}^2}<{\mathcal{R}}_0<1$, we have $\rho \omega ({\mathcal{R}}_0-1)+\sqrt{\mathrm{\Delta }}$ can be negative and also the expression ${\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)+\sqrt{\mathrm{\Delta }}\right)$ can be negative. But, if we use the identity

$$\begin{array}{ccc}\hfill {\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)+\sqrt{\mathrm{\Delta }}\right)& =& {\left(1+\rho \omega \right)}^2\left({\mathcal{R}}_0-\frac{4\omega \rho }{{(\omega \rho +1)}^2}\right)\left({\mathcal{R}}_0+\frac{2\sqrt{\mathrm{\Delta }}}{\omega \rho +1}\right)\hfill \\ & +& \mathrm{\Delta }+\frac{8\rho \omega \sqrt{\mathrm{\Delta }}}{1+\omega \rho }\hfill \end{array}$$

we obtain that ${\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right)}^2+4{\mathcal{R}}_0\left(\rho \omega ({\mathcal{R}}_0-1)+\sqrt{\mathrm{\Delta }}\right)>0$ for any ${\mathcal{R}}_0\in \left(\frac{4\omega \rho }{{(\omega \rho +1)}^2},1\right)$.

Then $detA>0$ and the second endemic equilibrium $E_2$ is local asymptotically stable (stable node or stable focus) if and only of $\mathrm{tr}A<0$, i.e.

$$\left({\left(\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}^2+{\mathcal{R}}_0(4\rho \omega -{\mathcal{R}}_0)\right)\left({\mathcal{R}}_0+\rho R\omega +\sqrt{\mathrm{\Delta }}\right)-4R\omega \left(\rho {\mathcal{R}}_0\omega -{\mathcal{R}}_0+\sqrt{\mathrm{\Delta }}\right)>0.$$

Otherwise, the endemic equilibrium $E_2$ is unstable (unstable node or unstable focus) if and only if $\mathrm{tr}A<0$, i.e.

$$\left({\left(\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}^2+{\mathcal{R}}_0(4\rho \omega -{\mathcal{R}}_0)\right)\left({\mathcal{R}}_0+\rho R\omega +\sqrt{\mathrm{\Delta }}\right)-4R\omega \left(\rho {\mathcal{R}}_0\omega -{\mathcal{R}}_0+\sqrt{\mathrm{\Delta }}\right)<0.$$

Let us point out that, if $\mathrm{tr}A=0$ (i.e. ${\lambda }_1=-{\lambda }_2$), then ${\lambda }_{1,2}$ are pure imaginary roots of the characteristic polynomial at $E_2$, with Re ${\lambda }_{1,2}=0$. In this case it is possible to have Hopf bifurcations along the curve tr $A=0$, i.e

$$\left({\left(\rho {\mathcal{R}}_0\omega +\sqrt{\mathrm{\Delta }}\right)}^2+{\mathcal{R}}_0(4\rho \omega -{\mathcal{R}}_0)\right)\left({\mathcal{R}}_0+\rho R\omega +\sqrt{\mathrm{\Delta }}\right)-4R\omega \left(\rho {\mathcal{R}}_0\omega -{\mathcal{R}}_0+\sqrt{\mathrm{\Delta }}\right)=0.$$

If $\mathrm{\Delta }=0$, i.e. ${\mathcal{R}}_0=\frac{4\omega \rho }{{(\omega \rho +1)}^2}$, then there is an unique endemic equilibrium point $E\left(x,y\right)$, with $x=\frac{\rho \omega +1}{2\omega \rho }$ and $y=\frac{\omega \rho -1}{2\omega }$. The Jacobian at this $E\left(x,y\right)$ is

$$A=\left(\begin{array}{cc}-2{\rho }^2\frac{\omega }{\rho \omega +1}& -2\rho \frac{\omega }{\rho \omega +1}\\ \frac{\rho }{\rho \omega +1}\left(\rho \omega -1\right)& \frac{1}{\rho \omega +1}\left(\rho \omega -1\right)\end{array}\right),$$

with eigenvalues ${\lambda }_1=0,{\lambda }_2=\frac{-2{\rho }^2\omega +\rho \omega -1}{\rho \omega +1}$, that means that this endemic equilibrium point E is a non hyperbolic equilibrium.

If $\mathrm{\Delta }<0$, i.e. ${\mathcal{R}}_0<\frac{4\omega \rho }{{(\omega \rho +1)}^2}$, then there is no endemic equilibrium.

In conclusion, we collect the obtained results in the next Table 2:

Remark 6.1 

Following the ideas and the same techniques as in Section 5, if we will consider the modified SIR model with demography and vaccination,

$$\left\{\begin{array}{ccc}\hfill S^{\prime }& =& \mathrm{\Lambda }-\beta (1+\nu I)IS-\mu S-pS\hfill \\ \hfill I^{\prime }& =& \beta (1+\nu I)IS-(\alpha +\mu )I\hfill \\ \hfill R^{\prime }& =& \alpha I-\mu R+pS\hfill \end{array}\right.$$

(31)

then, after changing of time and variables, we will obtain the same dimensionless system (25), where $\rho ={\displaystyle \frac{p+\mu }{\alpha +\mu }}$, ${\mathcal{R}}_0={\displaystyle \frac{\mathrm{\Lambda }\beta }{(p+\mu )(\alpha +\mu )}}$ and $\omega =\nu \frac{\mathrm{\Lambda }}{p+\mu }$. All results obtained remains available also for this SIR epidemic model with demography and vaccination. Only the interpretations can be different because the dimensionless parameters ρ and ${\mathcal{R}}_0$ are different.

7. SODE Formulation of the Modified SIR Pattern with Demography

We consider the modified SIR model with demography (25). For the sake of simplicity we will denote the derivative with respect to time with "dot" over the variable and we prefer to denote t instead of $\tau$. Then the system (25) becomes

$$\left\{\begin{array}{ccc}\hfill \dot{x}& =& \rho (1-x)-{\mathcal{R}}_0(1+\omega y)xy\hfill \\ \hfill \dot{y}& =& {\mathcal{R}}_0(1+\omega y)xy-y\hfill \end{array}\right.$$

(32)

where $\rho ={\displaystyle \frac{\mu }{\alpha +\mu }}\in (0,1)$, $\omega =\nu \frac{\mathrm{\Lambda }}{\mu }>0$, ${\mathcal{R}}_0={\displaystyle \frac{\mathrm{\Lambda }\beta }{\mu (\alpha +\mu )}}>0$.

By taking the derivative with respect to time t in both equations of this system, we obtain the following second order differential equations:

$$\left\{\begin{array}{ccc}\ddot{x}+\left[\rho +{\mathcal{R}}_0(1+\omega y)y\right]\dot{x}+\left[{\mathcal{R}}_0\omega xy+{\mathcal{R}}_0(1+\omega y)x\right]\dot{y}\hfill & =& \hfill 0\\ \ddot{y}-{\mathcal{R}}_0(1+\omega y)y\dot{x}+\left[1-{\mathcal{R}}_0\omega xy-{\mathcal{R}}_0(1+\omega y)x\right]\dot{y}\hfill & =& \hfill 0\end{array}\right.$$

In order to use the rule of the crossed repeated indices from differential geometry formalism, we change the notations of variables as follows:

$$x=x^1,\dot{x}=y^1,y=x^2,\dot{y}=y^2$$

Then this system of second-order differential equations (SODEs) becomes:

$$\left\{\begin{array}{ccc}{\ddot{x}}^1+\left[\rho +{\mathcal{R}}_0(1+\omega x^2)x^2\right]y^1+\left[{\mathcal{R}}_0\omega x^1{x}^2+{\mathcal{R}}_0(1+\omega x^2)x^1\right]y^2\hfill & =& \hfill 0\\ {\ddot{x}}^2-{\mathcal{R}}_0(1+\omega x^2)x^2{y}^1+\left[1-{\mathcal{R}}_0\omega x^1{x}^2-{\mathcal{R}}_0(1+\omega x^2)x^1\right]y^2\hfill & =& \hfill 0\end{array}\right.$$

(33)

or, equivalently,

$$\left\{\begin{array}{ccc}\frac{d^2{x}^1}{d{t}^2}+\left[\rho +{\mathcal{R}}_0(1+\omega x^2)x^2\right]y^1+\left[{\mathcal{R}}_0\omega x^1{x}^2+{\mathcal{R}}_0(1+\omega x^2)x^1\right]y^2\hfill & =& \hfill 0\\ \frac{d^2{x}^2}{d{t}^2}-{\mathcal{R}}_0(1+\omega x^2)x^2{y}^1+\left[1-{\mathcal{R}}_0\omega x^1{x}^2-{\mathcal{R}}_0(1+\omega x^2)x^1\right]y^2\hfill & =& \hfill 0\end{array}\right.$$

(34)

where $\frac{d{x}^i}{dt}=y^i$, $i=1,2$.

This system can be written similarly to a SODE (semispray) from the KCC theory:

$$\left\{\begin{array}{c}\frac{d^2{x}^1}{d{t}^2}+2G^1(x^1,x^2,y^1,y^2)=0\\ \frac{d^2{x}^2}{d{t}^2}+2G^2(x^1,x^2,y^1,y^2)=0\end{array}\right.$$

(35)

where

$$\begin{array}{ccc}{G}^1(x^i,y^i)\hfill & =& \frac{1}{2}\left[\rho y^1+{\mathcal{R}}_0\left(1+\omega x^2\right)x^2{y}^1+{\mathcal{R}}_0\left(1+2\omega x^2\right)x^1{y}^2\right]\hfill \\ G^2(x^i,y^i)\hfill & =& \frac{1}{2}\left[-{\mathcal{R}}_0(1+\omega x^2)x^2{y}^1+y^2-{\mathcal{R}}_0\left(1+2\omega x^2\right)x^1{y}^2\right]\hfill \end{array}$$

(36)

The zero-connection curvature $Z_j^i=2\frac{\partial G^i}{\partial x^j}$ is given by the next coefficients:

$$\begin{array}{ccc}\hfill Z_1^1& =& {\mathcal{R}}_0{y}^2+2{\mathcal{R}}_0\omega x^2{y}^2\hfill \\ \hfill Z_2^1& =& {\mathcal{R}}_0{y}^1+2{\mathcal{R}}_0\omega (x^2{y}^1+x^1{y}^2)\hfill \\ \hfill Z_1^2& =& -{\mathcal{R}}_0{y}^2-2{\mathcal{R}}_0\omega x^2{y}^2\hfill \\ \hfill Z_2^2& =& -{\mathcal{R}}_0{y}^1-2{\mathcal{R}}_0\omega (x^2{y}^1+x^1{y}^2)\hfill \end{array}$$

The nonlinear connection N has the coefficients $N_j^i=\frac{\partial G^i}{\partial y^j}$:

$$\left\{\begin{array}{ccc}\hfill N_1^1& =& \frac{1}{2}\left[\rho +{\mathcal{R}}_0\left(1+\omega x^2\right)x^2\right]\hfill \\ \hfill N_2^1& =& \frac{1}{2}{\mathcal{R}}_0\left(1+2\omega x^2\right)x^1\hfill \\ \hfill N_1^2& =& -\frac{1}{2}{\mathcal{R}}_0\left(1+\omega x^2\right)x^2\hfill \\ \hfill N_2^2& =& \frac{1}{2}\left[1-{\mathcal{R}}_0\left(1+\omega x^2\right)x^1\right]\hfill \end{array}\right.$$

(37)

Consequently, all coefficients of the Berwald connection $G_{jk}^i=\frac{\partial N_j^i}{\partial y^k}$ are equal with zero.

The first invariant of the KCC theory ${\epsilon }^i=-\left(N_j^i{y}^j-2G^i\right)$ has the following components:

$$\left\{\begin{array}{ccc}\hfill {\epsilon }^1& =& \frac{1}{2}\rho y^1+\frac{1}{2}{\mathcal{R}}_0\left(1+\omega x^2\right)x^2{y}^1+\frac{1}{2}{\mathcal{R}}_0\left(1+2\omega x^2\right)x^1{y}^2\hfill \\ \hfill {\epsilon }^2& =& \frac{1}{2}{y}^2-\frac{1}{2}{\mathcal{R}}_0(1+\omega x^2)x^2{y}^1-\frac{1}{2}{\mathcal{R}}_0\left(1+2\omega x^2\right)x^1{y}^2\hfill \end{array}\right.$$

(38)

Let us remark that ${\epsilon }^i=G^i$ for $i=1,2$, which means that $\frac{\partial G^i}{\partial y^j}{y}^j=1·G^i$ for $i=1,2$, or equivalently, that the functions $G^i$ are 1-homogeneous with respect to $y^i$.

Next, by using (A10), we obtain the components of the second invariant of the KCC theory, which means that the deviation curvature tensor of the modified SIR model with demography (32) is as follows:

$$\begin{array}{ccc}\hfill P_1^1& =& -\frac{1}{2}{\mathcal{R}}_0\left(1+2\omega x^2\right)y^2+\frac{1}{4}{\left[\rho +{\mathcal{R}}_0\left(1+\omega x^2\right)x^2\right]}^2-\frac{1}{4}{\mathcal{R}}_0^2\left(1+\omega x^2\right)\left(1+2\omega x^2\right)x^1{x}^2\hfill \\ \hfill P_2^1& =& -\frac{1}{2}{\mathcal{R}}_0{y}^1-{\mathcal{R}}_0\omega (x^2{y}^1+x^1{y}^2)+\frac{1}{4}{\mathcal{R}}_0\left(1+2\omega x^2\right)x^1\left[\rho +1+{\mathcal{R}}_0\left(x^2-x^1+\omega {\left(x^2\right)}^2-2\omega x^1{x}^2\right)\right]\hfill \\ \hfill P_1^2& =& \frac{1}{2}{\mathcal{R}}_0\left(1+2\omega x^2\right)y^2-\frac{1}{4}{\mathcal{R}}_0\left(1+\omega x^2\right)x^2\left[\rho +1+{\mathcal{R}}_0\left(x^2-x^1+\omega {\left(x^2\right)}^2-2\omega x^1{x}^2\right)\right]\hfill \\ \hfill P_2^2& =& \frac{1}{2}{\mathcal{R}}_0{y}^1+{\mathcal{R}}_0\omega (x^2{y}^1+x^1{y}^2)-\frac{1}{4}{\mathcal{R}}_0^2\left(1+\omega x^2\right)\left(1+2\omega x^2\right)x^1{x}^2+\frac{1}{4}{\left[1-{\mathcal{R}}_0\left(1+2\omega x^2\right)x^1\right]}^2\hfill \end{array}$$

(39)

Then, the trace and the determinant of the following deviation curvature matrix:

$$P=\left(\begin{array}{cc}{P}_1^1& P_2^1\\ P_1^2& P_2^2\end{array}\right)$$

are tr $P=P_1^1+P_2^2$ and $detP=P_1^1{P}_2^2-P_1^2{P}_2^1$.

Therefore, following Theorem A2, we obtain:

Theorem 7.1 

All the roots of the characteristic polynomial of P are negative or have negative real parts (i.e., Jacobi stability) if and only if

$$trP=P_1^1+P_2^2<0\mathrm{and}detP=P_1^1{P}_2^2-P_1^2{P}_2^1>0.$$

Taking into account that $P_{jk}^i=\frac{1}{3}\left(\frac{\partial P_j^i}{\partial y^k}-\frac{\partial P_k^i}{\partial y^j}\right)$, $P_{jkl}^i=\frac{\partial P_{jk}^i}{\partial y^l}$, $D_{jkl}^i=\frac{\partial G_{jk}^i}{\partial y^l}$, we obtained the third, fourth, and fifth invariants of the modified SIR model with demography (32) as follows:

Theorem 7.2 

The third KCC invariant $P_{jk}^i$, called the torsion tensor, has all eight components null, i.e.,

$$P_{jk}^i=0\mathrm{for}\mathrm{all}i,j,k.$$

(40)

The fourth KCC invariant $P_{jkl}^i$, called the Riemann–Christoffel curvature tensor, has all sixteen components null, i.e.,

$$P_{jkl}^i=0\mathrm{for}\mathrm{all}i,j,k,l.$$

(41)

The fifth KCC invariant $D_{jkl}^i$, called the Douglas tensor, has all sixteen components null, which means that

$$D_{jkl}^i=0\mathrm{for}\mathrm{all}i,j,k,l.$$

(42)

8. Jacobi Stability Analysis of the Modified SIR System with Demography

In this section, we will determine the first two invariants at the equilibrium points of the modified SIR model with demography (32), and we will analyze the Jacobi stability of the system near each equilibrium point.

For the disease-free equilibrium point $E_0(1,0)$, we have the corresponding equilibrium point $E_0(1,0,0,0)$ of SODE (34) and the results are the same as for the classical model. So, the first invariant of the KCC theory has the components ${\epsilon }^1={\epsilon }^2=0$, and the matrix with the components of the second KCC invariant is as follows:

$$P=\left(\begin{array}{cc}\frac{1}{4}{\rho }^2& \frac{1}{4}{\mathcal{R}}_0\left(\rho +1-{\mathcal{R}}_0\right)\\ 0& \frac{1}{4}{\left(1-{\mathcal{R}}_0\right)}^2\end{array}\right).$$

Since $trP=\frac{1}{4}{\rho }^2+\frac{1}{4}{\left(1-{\mathcal{R}}_0\right)}^2>0$ and $detP=\frac{1}{16}{\rho }^2{\left(1-{\mathcal{R}}_0\right)}^2>0$, we obtain the following result:

Theorem 8.1 

The disease-free equilibrium point $E_0$ of (32) is always Jacobi unstable.

Using (38) and tacking into account that an endemic equilibrium point $E(x,y)$ of (25) have the corresponding equilibrium point $E_0(x,y,0,0)$ of SODE (34), it results that the first invariant of the KCC theory has all components null, i.e. ${\epsilon }^1={\epsilon }^2=0$ for any endemic equilibrium.

Further, we will study the Jacobi stability for endemic equilibrium points following the three cases: ${\mathcal{R}}_0>1$, ${\mathcal{R}}_0=1$ and $0<{\mathcal{R}}_0<1$.

8.1. Jacobi Stability Near to Endemic Equilibrium for ${\mathcal{R}}_0>1$

If ${\mathcal{R}}_0>1$ then there is only one endemic equilibrium $E\left(x,y\right)$ in ${\Sigma }_{+}$, where

$$x=\frac{2}{{\mathcal{R}}_0(1+\omega \rho )+\sqrt{\mathrm{\Delta }}}\mathrm{and}y=\frac{1}{2{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0(\omega \rho -1)+\sqrt{\mathrm{\Delta }}\right),$$

where $\mathrm{\Delta }={\mathcal{R}}_0^2{(\omega \rho +1)}^2-4{\mathcal{R}}_0\omega \rho$.

Using (39), with $x^1=x$, $x^2=y$, $y^1=0$, $y^2=0$, and tacking into account that ${\mathcal{R}}_0(1+\omega y)x=1$, it results that the coefficients of second invariant of KCC-theory (the curvature deviation tensor) at E are the following:

$$\begin{array}{ccc}\hfill P_1^1& =& \frac{1}{4}{\left(\rho +\frac{y}{x}\right)}^2-\frac{1}{4}{\mathcal{R}}_0(1+2\omega y)y\hfill \\ \hfill P_2^1& =& \frac{1}{4}{\mathcal{R}}_0\left(\rho +1\right)(1+2\omega y)x+\frac{1}{4}{\mathcal{R}}_0(1+2\omega y)y-\frac{1}{4}{\mathcal{R}}_0^2{(1+2\omega y)}^2{x}^2\hfill \\ \hfill P_1^2& =& -\frac{1}{4}\left(\rho +1\right)\frac{y}{x}-\frac{1}{4}\frac{y^2}{x^2}+\frac{1}{4}{\mathcal{R}}_0(1+2\omega y)y\hfill \\ \hfill P_2^2& =& -\frac{1}{4}{\mathcal{R}}_0(1+2\omega y)y+\frac{1}{4}{\left(1-{\mathcal{R}}_0(1+2\omega y)x\right)}^2\hfill \end{array}$$

Then

$$\mathrm{tr}P=P_1^1+P_2^2=\frac{1}{4}{\left(\rho +\frac{y}{x}\right)}^2-\frac{1}{2}{\mathcal{R}}_0(1+2\omega y)y+\frac{1}{4}{\left(1-{\mathcal{R}}_0(1+2\omega y)x\right)}^2$$

or

$$4x^2\mathrm{tr}P=x^4{\left(1+2\omega y\right)}^2{\mathcal{R}}_0^2-2x^2\left(y+x\right)\left(1+2\omega y\right){\mathcal{R}}_0+{\rho }^2{x}^2+2\rho xy+y^2+x^2,$$

and

$$detP=P_1^1{P}_2^2-P_2^1{P}_1^2=\frac{{\left(-\rho x+\rho x^2{\mathcal{R}}_0+2\rho x^2{\mathcal{R}}_0\omega y-y\right)}^2}{16x^2}>0.$$

Therefore, following Theorem 7.1, we obtain:

Theorem 8.2 

The endemic equilibrium point E of (32) is Jacobi stable if and only if $\mathrm{tr}P<0$.

Using the identity

$$\mathrm{tr}P={\left(\frac{1+2\omega y}{2}\right)}^2\left({\mathcal{R}}_0-\frac{y+x-\sqrt{2xy-{\rho }^2{x}^2-2\rho xy}}{x^2\left(1+2\omega y\right)}\right)\left({\mathcal{R}}_0-\frac{y+x+\sqrt{2xy-{\rho }^2{x}^2-2\rho xy}}{x^2\left(1+2\omega y\right)}\right),$$

it results

Corollary 8.3 

If $2xy-{\rho }^2{x}^2-2\rho xy>0$, then the endemic equilibrium point E of (32) is Jacobi stable if and only if

$$\frac{y+x-\sqrt{2xy-{\rho }^2{x}^2-2\rho xy}}{x^2\left(1+2\omega y\right)}<{\mathcal{R}}_0<\frac{y+x+\sqrt{2xy-{\rho }^2{x}^2-2\rho xy}}{x^2\left(1+2\omega y\right)}$$

or, equivalently,

$$y+x-\sqrt{2xy-{\rho }^2{x}^2-2\rho xy}<{\mathcal{R}}_0{x}^2\left(1+2\omega y\right)<y+x+\sqrt{2xy-{\rho }^2{x}^2-2\rho xy}.$$

Otherwise, if $2xy-{\rho }^2{x}^2-2\rho xy\le 0$, then the endemic equilibrium E is Jacobi unstable.

In order to deep the study of the sign of $\mathrm{tr}P$, we replace x, y from (27), and it results

$$\begin{array}{ccc}\hfill P_1^1& =& \frac{1}{4}{\left(\rho +\frac{1}{4{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\left({\mathcal{R}}_0\left(1+\omega \rho \right)+\sqrt{\mathrm{\Delta }}\right)\right)}^2\hfill \\ \hfill & & -\frac{1}{8\omega }\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\left(1+\frac{1}{{\mathcal{R}}_0}\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\right)\hfill \\ \hfill P_2^1& =& \frac{1}{2}{\mathcal{R}}_0\left(\rho +1\right)\frac{1+\frac{1}{{\mathcal{R}}_0}\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)}{{\mathcal{R}}_0\left(1+\omega \rho \right)+\sqrt{\mathrm{\Delta }}}\hfill \\ \hfill & & +\frac{1}{8\omega }\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\left(1+\frac{1}{{\mathcal{R}}_0}\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\right)-{\mathcal{R}}_0^2\frac{{\left(1+\frac{1}{{\mathcal{R}}_0}\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\right)}^2}{{\left({\mathcal{R}}_0\left(1+\omega \rho \right)+\sqrt{\mathrm{\Delta }}\right)}^2}\hfill \\ \hfill P_1^2& =& -\frac{1}{16}\frac{\rho +1}{{\mathcal{R}}_0\omega }\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\left({\mathcal{R}}_0\left(1+\omega \rho \right)+\sqrt{\mathrm{\Delta }}\right)\hfill \\ \hfill & & -\frac{1}{64{\mathcal{R}}_0^2{\omega }^2}{\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)}^2{\left({\mathcal{R}}_0\left(1+\omega \rho \right)+\sqrt{\mathrm{\Delta }}\right)}^2\hfill \\ & & +\frac{1}{8\omega }\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\left(1+\frac{1}{{\mathcal{R}}_0}\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\right)\hfill \\ \hfill P_2^2& =& -\frac{1}{8\omega }\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\left(1+\frac{1}{{\mathcal{R}}_0}\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)\right)\hfill \\ & & +\frac{1}{4}{\left(1-2{\mathcal{R}}_0\frac{1+\frac{1}{{\mathcal{R}}_0}\left({\mathcal{R}}_0\left(\omega \rho -1\right)+\sqrt{\mathrm{\Delta }}\right)}{{\mathcal{R}}_0\left(1+\omega \rho \right)+\sqrt{\mathrm{\Delta }}}\right)}^2\hfill \end{array}$$

(43)

Then

$$\begin{array}{c}16{\mathcal{R}}_0\omega {\left({\mathcal{R}}_0(1+\rho \omega )+\sqrt{\mathrm{\Delta }}\right)}^2\mathrm{tr}P={\mathcal{R}}_0\omega {\rho }^2{\left({\mathcal{R}}_0(1+\rho \omega )+\sqrt{\mathrm{\Delta }}\right)}^4\\ -4{\left({\mathcal{R}}_0(1+\rho \omega )+\sqrt{\mathrm{\Delta }}\right)}^2\left({\mathcal{R}}_0(\rho \omega -1)+\sqrt{\mathrm{\Delta }}\right)\left({\mathcal{R}}_0\rho \omega +\sqrt{\mathrm{\Delta }}\right)+4{\mathcal{R}}_0\omega {\left({\mathcal{R}}_0(\rho \omega -1)+\sqrt{\mathrm{\Delta }}\right)}^2\end{array}$$

and we obtain the following result:

Theorem 8.4 

The endemic equilibrium point E of (32) is Jacobi stable if and only if

$$\begin{array}{c}{\mathcal{R}}_0\omega {\rho }^2{\left({\mathcal{R}}_0(1+\rho \omega )+\sqrt{\mathrm{\Delta }}\right)}^4+4{\mathcal{R}}_0\omega {\left({\mathcal{R}}_0(\rho \omega -1)+\sqrt{\mathrm{\Delta }}\right)}^2\\ <4{\left({\mathcal{R}}_0(1+\rho \omega )+\sqrt{\mathrm{\Delta }}\right)}^2\left({\mathcal{R}}_0(\rho \omega -1)+\sqrt{\mathrm{\Delta }}\right)\left({\mathcal{R}}_0\rho \omega +\sqrt{\mathrm{\Delta }}\right).\end{array}$$

8.2. Jacobi Stability Near to Endemic Equilibrium for ${\mathcal{R}}_0=1$

For ${\mathcal{R}}_0=1$ and $\omega \rho -1>0$, we have two so called endemic equilibrium points: $E_1(1,0)$ and $E_2(x,y)$ in ${\Sigma }_{+}$, with coordinates $x=\frac{1}{\omega \rho }$, $y=\frac{\omega \rho -1}{\omega }$. But, because $E_1$ coincides with $E_0$, remains to study the Jacobi stability only for $E_2$. Then, we have computed the coefficients of second invariant of KCC-theory (the curvature deviation tensor) at $E_2$:

$$\begin{array}{ccc}\hfill P_1^1& =& \frac{1}{4}\frac{{\rho }^4{\omega }^3+3\omega \rho -2{\rho }^2{\omega }^2-1}{\omega }\hfill \\ \hfill P_2^1& =& \frac{1}{4}\left(-1+2\omega \rho \right)\frac{-\omega \rho +{\rho }^3{\omega }^2+1}{{\rho }^2{\omega }^2}\hfill \\ \hfill P_1^2& =& -\frac{1}{4}\left(\omega \rho -1\right)\frac{-\omega \rho +{\rho }^3{\omega }^2+1}{\omega }\hfill \\ \hfill P_2^2& =& -\frac{1}{4}\left(\omega \rho -1\right)\frac{-{\rho }^2\omega +2{\rho }^3{\omega }^2-\omega \rho +1}{{\rho }^2{\omega }^2}\hfill \end{array}$$

Then it results that

$$\mathrm{tr}P=P_1^1+P_2^2=\frac{1}{4}\frac{{\left({\rho }^3{\omega }^2-\rho \omega +1\right)}^2-2{\rho }^2\omega {\left(\rho \omega -1\right)}^2}{{\rho }^2{\omega }^2}$$

and

$$detP=P_1^1{P}_2^2-P_2^1{P}_1^2=\frac{1}{16}\frac{{\left(\omega \rho -1\right)}^4}{{\omega }^2}>0,\mathrm{for}\mathrm{any}\mathrm{values}\mathrm{of}\mathrm{parameters}.$$

Theorem 8.5 

The endemic equilibrium $E_2$ of (32) is Jacobi stable if and only if the following condition is satisfied

$${\left({\rho }^3{\omega }^2-\rho \omega +1\right)}^2<2{\rho }^2\omega {\left(\rho \omega -1\right)}^2,$$

or, equivalently,

$$-\rho \sqrt{2\omega }\left(\rho \omega -1\right)\le {\rho }^3{\omega }^2-\rho \omega +1\le \rho \sqrt{2\omega }\left(\rho \omega -1\right).$$

Let us remark that $E_2$ is focus (stable or unstable) if and only if the following condition is satisfied

$${\left({\rho }^3{\omega }^2-\rho \omega +1\right)}^2<4{\rho }^2\omega {\left(\rho \omega -1\right)}^2,$$

or, equivalently,

$$-2\rho \sqrt{\omega }\left(\rho \omega -1\right)\le {\rho }^3{\omega }^2-\rho \omega +1\le 2\rho \sqrt{\omega }\left(\rho \omega -1\right).$$

If we denote by $S=\mathrm{tr}A$ and $P=detA$, then

$$S^2-4P=\frac{{\left({\rho }^3{\omega }^2-\rho \omega +1\right)}^2-4{\rho }^2\omega {\left(\rho \omega -1\right)}^2}{{\rho }^2{\omega }^2}$$

and

$$S^2-2P=\frac{{\left({\rho }^3{\omega }^2-\rho \omega +1\right)}^2-2{\rho }^2\omega {\left(\rho \omega -1\right)}^2}{{\rho }^2{\omega }^2}.$$

Consequently, we have $S^2-4P<0$ if and only if ${\left({\rho }^3{\omega }^2-\rho \omega +1\right)}^2-4{\rho }^2\omega {\left(\rho \omega -1\right)}^2<0$ and $S^2-2P<0$ if and only if ${\left({\rho }^3{\omega }^2-\rho \omega +1\right)}^2-2{\rho }^2\omega {\left(\rho \omega -1\right)}^2<0$.

Therefore, it is confirm once again that the Jacobi stability of an equilibrium point implies that this equilibrium is a focus (stable or unstable), see Figure A1. More that, we can claim that the Jacobi stability near to an equilibrium point not allows a chaotic behavior of the dynamical system around this equilibrium point.

8.3. Jacobi Stability Near to Endemic Equilibrium for ${\mathcal{R}}_0<1$

In this case, and for $\mathrm{\Delta }={\mathcal{R}}_0^2{(\omega \rho +1)}^2-4{\mathcal{R}}_0\omega \rho >0$, that means $\frac{4\omega \rho }{{(\omega \rho +1)}^2}<{\mathcal{R}}_0<1$, we have two endemic equilibrium points $E_1\left(x_1,y_1\right)$ and $E_2\left(x_2,y_2\right)$, with positive coordinates given by (29), and respectively (30). Because the first endemic equilibrium $E_1\left(x_1,y_1\right)$ is a saddle point, it results that we have no Jacobi stability at $E_1$. Instead, for the second endemic equilibrium $E_2\left(x_2,y_2\right)$, it obtained the same results about Jacobi stability as for the endemic equilibrium from the case ${\mathcal{R}}_0>1$, see subsection 8.1.

Obviously, when $\mathrm{\Delta }=0$, i.e. ${\mathcal{R}}_0=\frac{4\omega \rho }{{(\omega \rho +1)}^2}$, we have no Jacobi stability at this non hyperbolic equilibrium point.

Remark 8.6 

Similarly with subsection 4.1, near each equilibrium point, the dynamics of the deviation vector for this modified SIR system with demography can be obtained from the deviation equations (A8) as follow:

$$\left\{\begin{array}{ccc}\frac{d^2{\xi }^1}{d{t}^2}+\left[\rho +{\mathcal{R}}_0\left(1+\omega x^2\right)x^2\right]\frac{d{\xi }^1}{dt}+{\mathcal{R}}_0\left(1+2\omega x^2\right)x^1\frac{d{\xi }^2}{dt}\hfill & +& \\ \left[{\mathcal{R}}_0{y}^2+2{\mathcal{R}}_0\omega x^2{y}^2\right]{\xi }^1+\left[{\mathcal{R}}_0{y}^1+2{\mathcal{R}}_0\omega (x^2{y}^1+x^1{y}^2)\right]{\xi }^2\hfill & =& \hfill 0\\ \frac{d^2{\xi }^2}{d{t}^2}-{\mathcal{R}}_0\left(1+\omega x^2\right)x^2\frac{d{\xi }^1}{dt}+\left[1-{\mathcal{R}}_0\left(1+\omega x^2\right)x^1\right]\frac{d{\xi }^2}{dt}\hfill & -& \\ \left[{\mathcal{R}}_0{y}^2+2{\mathcal{R}}_0\omega x^2{y}^2\right]{\xi }^1-\left[{\mathcal{R}}_0{y}^1+2{\mathcal{R}}_0\omega (x^2{y}^1+x^1{y}^2)\right]{\xi }^2\hfill & =& \hfill 0\end{array}\right.$$

(44)

9. Conclusions

In this work, we conducted a study on the Jacobi stability of two SIR models with demography for the spread of diseases by using the geometric tools of the KCC theory. It’s about of the classical SIR model with demography (and vaccination) and a modified SIR model with demography (and vaccination), but with a linear transmission coefficient of infection. First of all, for each of this patterns, we presented the classical local dynamics around equilibrium points, and then we reformulated each first-order nonlinear differential system into a system of second-order differential equations (SODE) in order to determine the five geometrical invariants of the KCC theory. We calculated the first and the second invariant of the Kosambi–Cartan–Chern theory, and we were able to confirm that the third, fourth, and fifth invariants are all with null components, and that the Berwald connection vanishes. Furthermore, we determined the zero-connection curvature tensor, the nonlinear connection associated with the semispray (SODE), and we computed the deviation curvature tensor at each equilibrium point in order to determine the Jacobi stability conditions.

Moreover, in order to compare these two approaches, a comparative analysis of the Jacobi stability and the Lyapunov (linear) stability near equilibrium points was made. Additionally, the deviation equations near every equilibrium point were determined. A next approach could be the performance of a computational study on the time variation of the deviation vector and its curvature in order to highlight the behavior of the SIR system near the equilibrium points.

Appendix A.1. Comparison between Lyapunov Stability and Jacobi Stability for Two-Dimensional Systems

Following (A10) and according to [11,18,40], the matrix of the curvature deviation tensor $P_j^i$ at the equilibrium point $E(x_1,x_2,0,0)$ has the following expression:

$${\left.\left(P_j^i\right)\right|}_{(x_1,x_2,0,0)}=-2{\left.\left(\frac{\partial G^i}{\partial x^j}\right)\right|}_{(x_1,x_2,0,0)}+{\left.\left(N_l^i{N}_j^l\right)\right|}_{(x_1,x_2,0,0)}=\frac{1}{4}{A}^2$$

(A15)

where A is the Jacobian matrix at $(x_1,x_2)$, and $E(x_1,x_2)$ is the equilibrium point of the initial first-order system from which the system of the second-order differential equations (A1) will be obtained.

If ${\lambda }_1$, ${\lambda }_2$ are the eigenvalues of A, then $\frac{1}{4}{\lambda }_1^2$, $\frac{1}{4}{\lambda }_2^2$ are the eigenvalues of ${\left.\left(P_j^i\right)\right|}_{(x_1,x_2,0,0)}$.

Since ${\lambda }_1$, ${\lambda }_2$ are the roots of the following characteristic equation:

$${\lambda }^2-trA\lambda +detA=0,$$

then we have ${\lambda }_{1,2}=\frac{trA\pm \sqrt{\mathrm{\Delta }}}{2}$, where $\mathrm{\Delta }={\left(trA\right)}^2-4detA$.

Therefore, the equilibrium point $E(x_1,x_2,0,0)$ is Jacobi-stable if and only if the real parts of the eigenvalues of P are negative, i.e.,

$$\mathrm{\Delta }<0andRe{\lambda }_{1,2}^2=\frac{{\left(trA\right)}^2+\mathrm{\Delta }}{4}=\frac{{\left(trA\right)}^2-2detA}{2}<0,$$

because ${\lambda }_{1,2}^2=\frac{1}{4}\left({\left(trA\right)}^2+\mathrm{\Delta }\pm 2itrA\sqrt{-\mathrm{\Delta }}\right)$.

Therefore, the equilibrium point E is Jacobi-stable if and only if ${\left(trA\right)}^2-4detA<0$ and ${\left(trA\right)}^2-2detA<0$.

As in [25], in order to more clearly represent the relationship between linear (or Lyapunov) stability and Jacobi stability for two-dimensional systems, we will consider the following diagram with respect to $S={\lambda }_1+{\lambda }_2=trA$ and $P={\lambda }_1{\lambda }_2=detA$ (see Figure A1):

Figure A1. Relationschip between Lyapunov stability and Jacobi stability for 2D systems.

Figure A1. Relationschip between Lyapunov stability and Jacobi stability for 2D systems.