# Stability, Bifurcation, and a Pair of Conserved Quantities in a Simple Epidemic System with Reinfection for the Spread of Diseases Caused by Coronaviruses

Authors: Jorge Fernando Camacho; Cruz Vargas-De-León
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

In this paper, we study a modified SIRI model without vital dynamics, based on a system of nonlinear ordinary differential equations, for epidemics that exhibit partial immunity after infection, reinfection, and disease-induced death. This model can be applied to study epidemics caused by SARS-CoV, MERS-CoV, and SARS-CoV-2 coronaviruses, since there is the possibility that, in diseases caused by these pathogens, individuals recovered from the infection have a decrease in their immunity and can be reinfected. On the other hand, it is known that, in populations infected by these coronaviruses, individuals with comorbidities or older people have significant mortality rates or deaths induced by the disease. By means of qualitative methods, we prove that such system has an endemic equilibrium and an infinite line of nonhyperbolic disease-free equilibria, we determine the local and global stability of these equilibria, and we also show that it has no periodic orbits. Furthermore, we calculate the basic reproductive numberR0 and find that the system exhibits a forward bifurcation: disease-free equilibria are stable when R0<1/σ and unstable when R0>1/σ, while the endemic equilibrium consist of an asymptotically stable upper branch that appears from R0>1/σ, σ being the rate that quantifies reinfection. We also show that this system has two conserved quantities. Additionally, we show some of the most representative numerical solutions of this system.


---

## Body

## 1. Introduction

Coronaviruses (CoVs), first discovered in a domestic poultry in the 1930s, are a large family of viruses that can cause respiratory, gastrointestinal, liver, and neurological diseases in animals; of all of them, only seven coronaviruses are known to cause disease in humans. Four of these seven coronaviruses most frequently cause symptoms of the common cold; in rare cases, they can produce severe lower respiratory tract infections, including pneumonia, mainly in infants, older people, and the immunocompromised. The other three of the seven led to much more serious respiratory illnesses in humans, sometimes with fatal consequences, than other coronaviruses and have caused major outbreaks of deadly pneumonia in this century: SARS-CoV (also called as SARS-CoV-1), identified as the cause of the severe acute respiratory syndrome (SARS); MERS-CoV, responsible for the Middle East respiratory syndrome (MERS); and the novel coronavirus SARS-CoV-2 identified as the cause of coronavirus disease 2019 (COVID-19) [1, 2].SARS-CoV was first detected in Guangdong Province of China in November 2002 and subsequently spread to more than 30 countries [3]; MERS-CoV was first identified in September 2012 in Saudi Arabia, although there is a documented outbreak in April 2012 in Jordan [4]; and the new strain of coronavirus, SARS-CoV-2, was first reported in the city of Wuhan China in late 2019, and since then, it has spread widely around the world becoming an ongoing pandemic [1]. These coronaviruses that cause severe respiratory infections in humans are zoonotic pathogens (they can be transmitted from animals to people). Based on extensive studies, today, it is known that SARS-CoV was acquired from civet and that contagions of MERS-CoV were transferred by dromedary; however, the origin of the animal transmission of SARS-CoV-2 has not yet been clarified. Regarding SARS-CoV-2, it should be noted that it presents a high level of person-to-person contagion with respect to the other two. Besides, person-to-person spread not only occurs through contact with infected secretions but also occurs via contact with a surface contaminated by respiratory droplets. It is known that symptomatic, as well as asymptomatic and presymptomatic patients can transmit the virus, which, by the way, appears more transmissible than SARS-CoV [1].SARS-CoV, MERS-CoV, and by far, currently SARS-CoV-2, have been responsible for serious epidemic respiratory infections with great international repercussion due to their morbidity and mortality. At the outbreak of SARS-CoV, more than 8000 cases were reported worldwide, with 774 deaths, whereas through 2019, worldwide, nearly 2500 cases of MERS-CoV infection (with at least 850 related deaths) have been reported from 27 countries [1]. However, as for SARS-CoV-2, the numbers increased dramatically for just over a year, and it has spread to practically every country in the world; until July 21 of this year, this pandemic has caused 191,148,056 confirmed cases and 4,109,303 deaths [5]. Since the great impact of the abovementioned data, it is not unnecessary to say that, to overcome this global health crisis, it is essential to study and understand the nature of coronaviruses: their transmission, characteristics, spatial and temporal evolution, preventive and therapeutic measures, and, of course, their complications.It should be mentioned that infection by SARS-CoV-2 not only affects older people but also those with some type of comorbidity [6, 7], as occurring in some countries of the American continent, such as Mexico [8], whose populations show high levels of overweight, obesity, diabetes, and hypertension, to name a few cases [9]. The latent possibility that hospital services will collapse makes urgent the need for effective medical treatments. Instead of creating compounds from zero that can take years to develop, it has been chosen to reuse drugs that are already approved for other diseases and known to be safe. Unapproved drugs that have worked well in animal studies are also being tested with the other two fatal infections caused by SARS-CoV and MERS-CoV [10]. Since there are currently no specific medications available to inhibit the viral replication of SARS-CoV-2, government measures each country has taken to limit its spread are based on nonpharmaceutical interventions: cancellation of mass events, social distancing, confinement measures, and even curfews [11]. Although different types of SARS-CoV-2 vaccines have been rapidly developed today and started to be applied (to date, only 26% of the world’s population has received at least one dose and 13.5% is fully vaccinated [12]), the impact of vaccination on the COVID-19 pandemic is far from significant since there are new outbreaks and it continues to expand worldwide.Diseases caused by respiratory coronaviruses in humans have the peculiarity of reinfecting them throughout their lives [13]. Although, for example, the persistence of antibodies against MERS-CoV has been found in infected individuals almost three years after the outbreak of this disease in Jordan in 2012 [14], for the case of SARS, the duration of immunity has been estimated at more than 3 years [15]. In a cohort study, 470 infected subjects with human coronavirus (HKU1, OC43, NL63, and 229E) were followed up for some years. In this cohort, 8.5% of infected subjects had reinfection with new or the same coronavirus occurring 1.55 and 1.85 years, respectively [16]. A key unresolved question is the duration of acquired immunity in people recovered from COVID-19. The recent appearance of SARS-CoV-2 prevents a direct study on this virus. However, there have been some studies that suggest that recovered COVID-19 patients lose their immunity within a few months [17, 18]. The first cases of two instances of SARS-CoV-2 infection in the same individual have been reported in the United States [19] and in Hong Kong [20]. The first confirmed cases of reinfection in Mexico were reported by Murillo-Zamora et al. [11]. Several systematic reviews have documented the different cases of reinfection reported in the literature [21–23]. To date, no cohort studies have been reported on the follow-up of infected subjects who present reinfection with SARS-CoV-2 [24].Previous evidence indicates that a coronavirus reinfection is likely. This will be precisely one of our working hypotheses in this research. We will describe the temporal evolution of epidemics with these characteristics, using a compartmental mathematical model based on a set of nonlinear ordinary differential equations, whose variables are the different groups, individuals, or compartments of the community affected by the infection according to their state with respect to the disease [25, 26]. Mathematical models, especially compartmental ones, play a fundamental role in understanding and predicting epidemics caused by these deadly coronaviruses, such as the recent COVID-19 pandemic. Additionally, they have also become one of the elements used by governments in the design of effective public policies that allow satisfying the growing demand for medical attention that the population requires and the implementation of social confinement measures to mitigate its growth.Since there is the possibility that individuals recovered from these diseases present partial immunity and reinfection, we will use to describe them a SIRI (susceptible-infectious-recovered-infectious) model that takes both characteristics into account; in this model, a population of infectious individuals can only recover temporarily and once it recovers cannot remain that condition and become infected again [25]. The SIRI model has already been applied to study these types of epidemics (with partial immunity and reinfection), from practical situations such as the effects of vaccination based on the information available on a coronavirus outbreak [27], to relevant analytical considerations of the model itself [28]. This model has also been used to describe epidemics of another type, such as that caused by overweight and obesity [29, 30].It has already been mentioned that SARS-CoV, MERS-CoV, and SARS-Cov-2 are viruses that cause epidemics with high death rates. The description of this relevant aspect is also part of the purposes of this research work. One way to include this situation in our model is to expand it with an additional variable, which is the population of individuals who died due to the disease; however, this procedure has the great drawback, since our system now has one more variable, of hindering its qualitative analysis. An alternative way that avoids the previous problem is to introduce a loss term that we identify as a disease-induced death term directly into the equation for the evolution of infected individuals; the details of how this is proposed are given in the next section.Several transmission models have been developed for the dynamics of novel SARS-CoV infection [31–33]. These models have considered the natural history of the disease, as well as the control measures such as quarantine and isolation, but these have not considered the reinfection phenomenon.In sum up, in this work, we propose a variant of the SIRI model to describe epidemics in human populations caused by the previously described coronaviruses, which have the possibility of granting temporary immunity and reinfection, as well as causing high lethality. The objective of this work is to analyze this system using the methods of the qualitative theory of ordinary differential equations and show some of the most representative numerical solutions.The rest of this paper is organized as follows. In Section2, we present our modified SIRI model and reduce it to a two-dimensional system. Section 3 focuses on qualitative analysis: here, we carry out a local analysis to establish equilibrium points and their corresponding local stabilities, as well as bifurcations; we also carry out a global analysis using a suitable Lyapunov function to establish the global stability of the endemic equilibrium and, using the Dulac criterion, the absence of periodic orbits; additionally, we show that our system presents two first integrals. In Section 4, some of the most representative numerical solutions of our model are shown. Finally, in Section 5, we discuss the results obtained and present our conclusions.


## 2. Model

From a population or “macroscopic” point of view, in these epidemic outbreaks caused by the aforementioned coronaviruses, three types of populations of individuals can be identified: those infectedI, who, by interacting with the susceptible S, get infected, causing the eventual spread of the disease throughout the population, and recovered R, who apparently acquire immunity (cannot be infected again) and do not affect the subsequent development of the epidemic. One of the most successful and common mathematical tools that describe the abovementioned situation is the epidemiological model SIR [34, 35]. In this, St, It, and Rt denote the numbers of individuals in the groups or compartments S, I, and R, at time t, respectively; the total population is considered to be of fixed size N, so that St+It+Rt=N. In addition, these epidemics have the characteristic of developing rapidly; that is, they spread in a relatively short time scale. For this reason, it is irrelevant to take into account in its description the natural births and deaths of the population, whose significant changes in size occur over a much longer period. In other words, vital dynamics factors or explicit demography can be neglected in our model [25, 26].Furthermore, in our description of these types of epidemic outbreaks, we must take into account that individuals recovered from the infection have temporary immunity and the possibility of reinfection, so considering aSIRI model (a variant of the SIR that incorporates these features) is more convenient [25, 27, 28]. Since we will also consider that these outbreaks can present a high lethality, with the purpose of modeling them, we propose a modified SIRI model with a disease-induced death rate and without vital dynamics, written in terms of the following set of ordinary differential equations:(1)dSdt=−βSI,dIdt=βSI−δI+σβRI−αI,dRdt=δI−σβRI.In equation (1) β, δ,α, and σ are positive parameters that represent the disease transmission, recovery, disease-induced death, and recurrence rates, respectively. In this description, note that the value of parameter β depends on whether the disease is more or less contagious, 1/δ is the mean recovery period, and σ is the reduction in susceptibility due to previous infection that takes values between zero and one 0<σ<1. It should also be noted that system (1) is in effect a modified SIRI model: the typical term σβRI represents the reinfection; that is, the recovered people, after some time, become infected again; the change comes from the disease-induced death term αI, which quantifies, in a first approximation, the losses attributed to unnatural deaths of individuals infected by the disease. It can also be see that, in the total absence of reinfection and disease-induced death terms, system (1) is reduced to the corresponding well-known SIR model without vital dynamics. Figure 1 shows schematically with a flowchart the assumptions made in the deduction of this system. From system (1), we have(2)dNdt=dS+I+Rdt=−αI,so the rate of change of the total population N is not constant, it decreases proportionally with the number of infected individuals I.Figure 1
Flowchart of the variant of modelSIRI considered. This is the well-known SIR model without vital dynamics, with its term of reinfection σβRI and modified by the term of disease-induced death αI.


### 2.1. Adimensionalization and Reduction of a Two-Dimensional System

It is much more convenient to handle system (1), simplifying its writing by dimensionless. With this purpose, we propose to change the variables of number of individuals S,I,R and the temporal t, by the dimensionless quantities [25].(3)x=SN,y=IN,z=RN,(4)τ=δ+αt.The latter form is motivated by the fact that the parametersα,δ have dimensions of (units of time)−1. From (3), we have that Nx˙=S˙−xN˙,Ny˙=I˙−yN˙ and Nz˙=R˙−zN˙; note that the point above the quantity in question indicates its derivative with respect to time t. On the other hand, from (4), using the chain rule, we obtain that dx/dt=δ+αdx/dτ,dy/dt=δ+αdy/dτ,dz/dt=δ+αdz/dτ. Therefore, system (1), written in terms of these new set of dimensionless variables, takes the form(5)dxdτ=−R0−ρ1xy,dydτ=−y+R0xy+σR0yz+ρ1y2,dzdτ=ρ2y−σR0−ρ1yz,where we have also defined the dimensionless quantities(6)ρ1=αδ+α,ρ2=δδ+α,(7)R0=β′δ+α,which is called the basic reproductive number, and that represents the number of new infections produced by a single infected person if the entire population is susceptible. R0 is commonly considered as a key parameter in the evolution of the epidemic: the smaller it is, the slower the epidemic will be (in practice, observing the epidemic makes it possible to measure R0 and consequently estimate β′). In (7), we have defined the quantity β′=βN. It should be noted that since the parameters involved δ,β′ are positive and α is nonnegative, then 0≤ρ1<1, 0<ρ2≤1 and 0<R0; besides, we have that ρ1+ρ2=1. Note that the new variables x,y,z meet the condition(8)x+y+z=1.Taking into account precisely this condition, it can be shown that the set of equations (5) can be simplified, for the variables y,z, in the two-dimensional system(9)dydτ=R0−1y−1−σR0yz−R0−ρ1y2,dzdτ=ρ2y−σR0−ρ1yz.It may be shown easily that the region of biological sense is(10)Λ=y,z∈ℝ2|y≥0,z≥0,y+z<1,that is, all solutions starting in Λ remain there for all τ≥0. Clearly, the set Λ is positively invariant with respect to system (9).


### 2.2. Basic Reproductive Number

The procedure in whichR0 was previously obtained is a consequence of rewriting system (1), by means of the dimensionless variables (3) and (4), in form (5). An alternative way of finding R0 is to consider that initially, S0=S0,I0>0 with I0≪S0 and R0=0. In this case, second equation (1) takes the form(11)dI0dt=βS0−δ+σβR0−αI0.Hence, a necessary and sufficient condition for an initial increase in the number of infected (sinceR0=0) is βS0−δ−α>0, or equivalently R0>1, where(12)R0=βS0δ+α=β′x0δ+α,where x0=S0/N. On the other hand, R0<1 is the sufficient condition to have an initial decrease of infected and that the epidemic does not spread. Thus, equation (12) is the expression for R0 that we are looking for and is in complete agreement, when x0=1, with equation (7) found previously.


## 2.1. Adimensionalization and Reduction of a Two-Dimensional System

It is much more convenient to handle system (1), simplifying its writing by dimensionless. With this purpose, we propose to change the variables of number of individuals S,I,R and the temporal t, by the dimensionless quantities [25].(3)x=SN,y=IN,z=RN,(4)τ=δ+αt.The latter form is motivated by the fact that the parametersα,δ have dimensions of (units of time)−1. From (3), we have that Nx˙=S˙−xN˙,Ny˙=I˙−yN˙ and Nz˙=R˙−zN˙; note that the point above the quantity in question indicates its derivative with respect to time t. On the other hand, from (4), using the chain rule, we obtain that dx/dt=δ+αdx/dτ,dy/dt=δ+αdy/dτ,dz/dt=δ+αdz/dτ. Therefore, system (1), written in terms of these new set of dimensionless variables, takes the form(5)dxdτ=−R0−ρ1xy,dydτ=−y+R0xy+σR0yz+ρ1y2,dzdτ=ρ2y−σR0−ρ1yz,where we have also defined the dimensionless quantities(6)ρ1=αδ+α,ρ2=δδ+α,(7)R0=β′δ+α,which is called the basic reproductive number, and that represents the number of new infections produced by a single infected person if the entire population is susceptible. R0 is commonly considered as a key parameter in the evolution of the epidemic: the smaller it is, the slower the epidemic will be (in practice, observing the epidemic makes it possible to measure R0 and consequently estimate β′). In (7), we have defined the quantity β′=βN. It should be noted that since the parameters involved δ,β′ are positive and α is nonnegative, then 0≤ρ1<1, 0<ρ2≤1 and 0<R0; besides, we have that ρ1+ρ2=1. Note that the new variables x,y,z meet the condition(8)x+y+z=1.Taking into account precisely this condition, it can be shown that the set of equations (5) can be simplified, for the variables y,z, in the two-dimensional system(9)dydτ=R0−1y−1−σR0yz−R0−ρ1y2,dzdτ=ρ2y−σR0−ρ1yz.It may be shown easily that the region of biological sense is(10)Λ=y,z∈ℝ2|y≥0,z≥0,y+z<1,that is, all solutions starting in Λ remain there for all τ≥0. Clearly, the set Λ is positively invariant with respect to system (9).


## 2.2. Basic Reproductive Number

The procedure in whichR0 was previously obtained is a consequence of rewriting system (1), by means of the dimensionless variables (3) and (4), in form (5). An alternative way of finding R0 is to consider that initially, S0=S0,I0>0 with I0≪S0 and R0=0. In this case, second equation (1) takes the form(11)dI0dt=βS0−δ+σβR0−αI0.Hence, a necessary and sufficient condition for an initial increase in the number of infected (sinceR0=0) is βS0−δ−α>0, or equivalently R0>1, where(12)R0=βS0δ+α=β′x0δ+α,where x0=S0/N. On the other hand, R0<1 is the sufficient condition to have an initial decrease of infected and that the epidemic does not spread. Thus, equation (12) is the expression for R0 that we are looking for and is in complete agreement, when x0=1, with equation (7) found previously.


## 3. Analysis

### 3.1. Equilibria and Nullclines

In order to find the equilibrium points of system (9), we follow the usual way of equating the left members of its equations with zero to obtain the values of the pair of variables y,z that satisfy this equality. As a result of this procedure, we obtain two equilibrium points. One is disease-free equilibria of the form(13)E00,z,with z∈ℝ; that is, the straight line y=0 with an infinite number of points that constitute the vertical coordinate axis of the yz plane. The second is the endemic equilibrium(14)E∗y∗,z∗=σR0−1σR0−ρ1,ρ2σR0−ρ1.So, the dynamics develop in regionΛ, and equilibrium points E∗ must be in the first quadrant of plane yz, that is, y∗≥0 and z∗>0; in other words, σR0−ρ1>0 and σR0−1≥0. Just when σR0−1=0, E∗ is located on the vetical z-axis; that is, it is part of the points E0.An alternative way to find the equilibrium points of a system is through the intersections of its nullclines [36, 37]. Since we will use these ideas later, we will begin this discussion below. According to system (9), we have the nullclines dy/dτ=0 (vertical directions), given by the functions(15)z=my+b,y=0,which consist of two straight lines: one of slope m=−R0−ρ1/1−σR0 and ordered to origin b=R0−1/1−σR0 and another that coincides with the vertical axis z. In addition, we have the nullclines dz/dτ=0 (horizontal directions), determined by(16)z=z∗,y=0,which are again two lines: one horizontal that cuts the vertical axis at point 0,ρ2/σR0−ρ1 and another vertical that also coincides with the z-axis. It should be noted that, as 0<σ<1 and R0>0, the quantity 1−σR0 is always positive; consequently, the slope of the vertical nullcline is only determined by the sign of R0−ρ1 and its ordinate to the origin by that of R0−1. As we are interested in the points E∗ located in region Λ, it must happen that σR0−ρ1>0,σR0−1≥0, and σR0−ρ1>0, as already mentioned previously; therefore, note that R0−ρ1=1−σR0+σR0−ρ1>0 and R0−1=1−σR0+σR0−1>0. These last two results imply that, in region Λ, the first vertical nullcline (15) must have a negative slope and its ordinate to the origin positive. In addition, it was also mentioned that, for E∗ to be in Λ, σR0−ρ1>0; therefore, the first horizontal nullcline (16) must be located in the upper half plane of plane yz.The intersections of different nullclines coincide with the equilibrium points. Note that, aty=0, nullclines overlap in both vertical and horizontal directions; therefore y=0 consists of an infinite line of equilibrium points, that is, the equilibria E0. On the other hand, another intersection between the other different nullclines gives rise to point E∗. Figure 2 illustrates the nullclines (15) and (16) in a configuration compatible with region Λ; it indicates horizontal nullclines in blue and verticals in red, vertical axis z (of infinite equilibrium points), as well as equilibrium points with black circles. This figure also shows the different directions that the vector field has along all the nullclines.Figure 2
Representation of nullclines of system (9): horizontal direction z=z∗ (in blue) and vertical direction z=my+b (in red); the vertical axis zy=0 is a combination of both. Equilibrium points are indicated by black circles. (a) σR0−ρ1>0,σR0−1>0,R0−ρ1>0, and R0−1>0; (b) σR0−ρ1>0,σR0−1<0,R0−ρ1>0, and R0−1>0.
(a)(b)The directions (left or right) of the horizontal nullclines are determined by substitutingz=z∗, first equation (16), in the first equation (9), which gives rise to(17)dydτ=R0−ρ1yy∗−y.To determine the positive or negative sign of (17), it is required to know the signs of y∗ (which depends in turn on the signs of σR0−1 and σR0−ρ1), the factor R0−ρ1, the values of y, and their order relation with y∗. Since Λ is satisfied that σR0−ρ1>0,σR0−1≥0 (that is, y∗≥0) and R0−ρ1>0(z∗ is on the nonnegative z-axis), the horizontal nullcline z=z∗ cuts the nonnegative vertical axis at point 0,z∗ and the nullcline of vertical directions, first equation (15), at y∗,z∗, which is an equilibrium point. Therefore, according to (16), if y∗<y<∞, then dy/dτ<0; that is, in this interval, the vector field is directed to the left if we are to the right of y∗,z∗; when 0<y<y∗, then dy/dτ>0 so that the vector field points to the right if we are to the left of y∗,z∗ and to the right of 0,z∗, while if −∞<y<0, then dy/dτ<0; that is, the direction of the field is to the left if we are to the left of 0,z∗. This result indicates that, on horizontal nullcline z=z∗, the vector field leaves point 0,z∗ and enters y∗,z∗.The abovementioned analysis can be repeated for all combinations of the quantitiesσR0−1,σR0−ρ1 (which determine the sign of y∗) and R0−ρ1. It can be shown that there are four possible cases. For y∗≥0, if R0−ρ1>0, then vector field on the horizontal nullcline z=z∗ leaves 0,z∗ and enters y∗,z∗, and when R0−ρ1<0, it enters 0,z∗ and leaves y∗,z∗; Figure 2(a) illustrates the first case. Also, for y∗<0, if R0−ρ1>0, then field enters 0,z∗ and leaves y∗,z∗, and when R0−ρ1<0, it leaves 0,z∗ and enters y∗,z∗; Figure 2(b) shows the third case. Note that only in the first case the equilibrium points are in region Λ.On the other hand, the directions (up or down) of the vertical nullclines are determined by replacing the vertical nullcline, first equation (15), in the second equation (9), which allows obtaining(18)dzdτ=−R0−ρ1σR0−ρ11−σR0yy∗−y.In order to determine the positive or negative sign of (18), we need to know the signs of y∗ (that depends on the signs of σR0−1 and σR0−ρ1), the factor R0−ρ1, the values of y, and their relation with y∗. Since in Λ, we have that σR0−ρ1>0,σR0−1≥0,R0−ρ1>0 and R0−1>0, the first vertical nullcline (16) cuts horizontal nullcline z=z∗ at the equilibrium point y∗,z∗ and the positive vertical axis at the point 0,b and has a negative slope. Therefore, according to (18), if y∗<y<∞, then dz/dτ>0; that is, in this interval, the direction of the vector field is up, while if 0<y<y∗, then dz/dτ<0 so that, in this interval, directions are down, and when −∞<y<0,dz/dτ>0, in this interval, directions are up.A complete analysis of the directions (up or down) of this vertical nullcline requires taking into account all combinations of signs of the quantitiesσR0−1,σR0−ρ1 and R0−ρ1. It can also be shown that there are four possible cases. For y∗>0, if R0−ρ1 and σR0−ρ1 are of opposite signs, then the direction of vector field on this vertical nullcline is down in the interval −∞<y<0, is up in 0<y<y∗, and is down in y∗<y<∞, whereas if R0−ρ1 and σR0−ρ1 are of the same signs, then the vector field is up in −∞<y<0, down in 0<y<y∗, and up in y∗<y<∞. Figure 2(a) illustrates the second case. Besides, for y∗<0, if R0−ρ1 and σR0−ρ1 are of opposite signs, then the vector field is down in −∞<y<y∗, up in y∗<y<0, and down in 0<y<∞, whereas if R0−ρ1 and σR0−ρ1 are of the same signs, then the vector field is up in −∞<y<y∗, down in y∗<y<0, and up in 0<y<∞. Figure 2(b) shows the fourth case. Note that only in the second case the equilibrium points are in region Λ.


### 3.2. Local Stability

In order to determine the stability of these families of equilibrium points, we need to calculate the Jacobian matrix of system (9) and evaluate it on them. This matrix overall is given as(19)Jy,z=−1−σR0z−2R0−ρ1y+R0−1−1−σR0yρ2−σR0−ρ1z−σR0−ρ1y.To perform the local stability analysis, we will consider separately the disease-free equilibria and the endemic equilibrium cases.


#### 3.2.1. The Disease-Free Equilibria,E0=0,z

The Jacobian matrix (19) evaluated in the disease-free equilibria E0 takes the simple form(20)J0,z=−1−σR0z+R0−10ρ2−σR0−ρ1z0,whose eigenvalues are λ1=0 and λ2=−1−σR0z+R0−1. Note that these equilibrium points are nonhyperbolic. Due to this nature of the points E0 and the infinite number of them, it is quite difficult to determine their stability. One way to visualize the dynamics around them is through a geometric analysis of the nullclines of system (9). This allows us to find the directions left or right of the vector field along the horizontal nullclines and its directions up or down on the vertical nullclines; in particular, this gives a glimpse of how it is around the equilibrium points. Of particular interest, here are the results of the analysis performed in the previous section of the directions of the vector field on the horizontal nullclines z=z∗ around their point of intersection 0,z∗ with the vertical axis z. Those results, rewritten in terms of the ideas of stability, can be written in the following way: for y∗>0, if R0−ρ1>0, then point 0,z∗ is unstable, and when R0−ρ1<0, it is stable, whereas for y∗<0, if R0−ρ1>0, then point 0,z∗ is stable, and when R0−ρ1<0, it is unstable, see Figure 2. Since we have an infinite number of points along the vertical axis z, the stability along a horizontal line that cuts one of them, say 0,z∗, must be identical to the stability found at its intersection with the horizontal nullcline z=z∗. It should be noted that the analysis based on the method of the nullclines reveals that points 0,z∗ are stable if z>b and unstable if z<b, see again Figure 2.Since0<σ<1 and 0≤ρ1<1, we can establish the relations 0≤ρ1<ρ1/σ and ρ1/σ<1/σ; besides that, we know that 1<1/σ; therefore, the order relation between these quantities is 0≤ρ1<ρ1/σ<1/σ<∞, which will be called segment L. On the other hand, it has just been mentioned that the stability of the nonhyperbolic points 0,z∗ depends on the signs of σR0−ρ1,σR0−1 and R0−ρ1, or in other words, it depends of the location that R0 has in L.Thus, whenσR0−ρ1>0,σR0−1≥0 (which determine that y∗>0) and R0−ρ1>0, according to segment L, the intersection of these three subintervals is 1/σ<R0<∞. Since the previous stability result establishes that, in this situation, points 0,z∗ are unstable, we conclude that the equilibrium points E0 are unstable when 1/σ<R0<∞. On other hand, when σR0−ρ1>0,σR0−1<0 and R0−ρ1>0, according to L, intersection is ρ1/σ<R0<1/σ; since points 0,z∗ are stable here, we have that points E0 are unstable when ρ1/σ<R0<1/σ. By means of an analysis similar to the previous one, it can be shown that points E0 are unstable in ρ1<R0≤ρ1/σ and stable in 0≤R0<ρ1. These obtained results can be summarized as follows.Theorem 1.
Disease-free equilibrium pointsE0, which constitute the vertical axis z, are nonhyperbolic. These points are(i)
unstable, if1/σ<R0<∞ and ρ1<R0<ρ1/σ(ii)
stable, ifρ1/σ<R0<1/σ and 0<R0<ρ1Since whenα=0, our model is reduced to the conventional SIRI model, in this situation ρ1=0, so from Theorem 1, it follows that the nonhyperbolic equilibrium points E0 are stable if 0<R0<1/σ and unstable if 1/σ<R0<∞.


#### 3.2.2. The Endemic Equilibrium,E∗=y∗,z∗

By substituting the endemic equilibriumE∗=y∗,z∗ in system (9), the identities are obtained.(21)−1−σR0z∗=R0−ρ1y∗−R0−1,ρ2−σR0−ρ1z∗=0,which allows writing the Jacobian matrix (19), evaluated in such a point, as(22)Jy∗,z∗=−R0−ρ1y∗−1−σR0y∗0−σR0−ρ1y∗.The trace of (22) is given by(23)TrJ=2ρ1−σ+1R0y∗=2ρ1−σ+1R0σR0−1σR0−ρ1,and its determinant is(24)DetJ=R0−ρ1σR0−ρ1y∗2=R0−ρ1σR0−12σR0−ρ1,whereas it can be shown that its discriminant, defined as Δ=TrJ2−4DetJ, is(25)Δ=−1−σR0y∗2=1−σ2R02σR0−12σR0−ρ12,which is a nonnegative quantity. According to the trace-determinant plane method [36–38], the analysis of matrix (22) indicates that points E∗, since Δ≥0, can be saddle points (if DetJ<0) and degenerate equilibrium (if DetJ=0) or nodes of two “straight-line” solutions (if DetJ>0), the latter pair being stable or unstable if TrJ<0 or TrJ>0, respectively. Furthermore, the eigenvalues of (22) are given by(26)λ1,2=12TrJ±Δ,which can be explicitly written as(27)λ1=−R0−ρ1y∗=−R0−ρ1σR0−1σR0−ρ1,(28)λ2=−σR0−ρ1y∗=−σR0−1.Note that eigenvalues (27) and (28) are functions of R0 and their signs will depend on the location of this quantity in segment L already indicated above. In general, they are different from zero, so the equilibrium points E∗ are hyperbolic.Next, we will analyze the signs thatλ1 and λ2 could have in L. According to (27) and (28), λ1<0 and λ2<0, if σR0−1>0,R0−ρ1>0, and σR0−ρ1>0; the intersection of these three subintervals occurs at 1/σ<R0<∞, so in this region, the equilibrium points E∗ correspond to nodes of two “straight-line” solutions that are locally asymptotically stable. In this case, in which both eigenvalues are negative, points E∗ lie within region Λ; however, for other combinations of signs of the eigenvalues λ1 and λ2, points E∗ lie outside of Λ. Thus, it can be shown that λ1>0 and λ2>0 when ρ1/σ<R0<1/σ and 0≤R0<ρ1 so that, in these regions, points E∗ correspond to unstable nodes of two “straight lines.” Additionally, λ1>0 and λ2<0, if ρ1<R0<ρ1/σ, so there points E∗ are saddle points; the case λ1<0 and λ2>0 does not exist because the intersection of the corresponding subintervals gives rise to the empty set. Moreover, λ1=0 and λ2>0, when R0=ρ1, so here, E∗ is a singly degenerate equilibrium. Lastly, λ1=λ2=0, when R0=1/σ, so E∗ is a doubly degenerate equilibrium. These last two cases correspond to nonhyperbolic equilibrium points. The previous results are summarized in Theorem 2.Theorem 2.
The endemic equilibrium pointsE∗ are in general hyperbolic. Such points are(i)
nodes (of two “straight-line” solutions) locally asymptotically stable, if1/σ<R0<∞(ii)
a doubly degenerate equilibrium, ifR0=1/σ (in which case it is nonhyperbolic)(iii)
unstable nodes (of two “straight-line” solutions), if0≤R0<ρ1 and ρ1/σ<R0<1/σ(iv)
saddle points, ifρ1<R0<ρ1/σ(v)
a singly degenerate equilibrium, ifR0=ρ1 (in which case it is nonhyperbolic)(vi)
unstable nodes (of two “straight-line” solutions), if0≤R0<ρ1It is important to note that, in Theorem2, subinterval 1/σ<R0<∞ (that of condition i)) is the only one that guarantees that points E∗ are in region Λ. Again, since when α=0, our model is reduced to the conventional SIRI model, in this situation, ρ1=0, so for this case, from Theorem 2, it follows that the equilibrium points E∗ are nodes (of two “straight-line” solutions) locally asymptotically stable if 1/σ<R0<∞ and unstable if 0<R0<1/σ.


### 3.3. Presence of a Forward Bifurcation

The analysis that allowed to formulate the two previous theorems shows that the coordinates of both pointsE0 and E∗ are a function of the reproductive number R0. In particular, it is of great interest to know how the family of equilibrium points y∗ (which is the abscissa of points E∗) varies with R0 and to visualize such dependence in a diagram R0y∗; this is important because R0 could play the role of a bifurcation parameter. For our study problem, the y∗ families consist of two curves: one horizontal line y∗=0 of disease-free equilibrium points located on the R0 axis and another determined by(29)y∗R0=σR0−1σR0−ρ1,of endemic equilibrium points. By inspection, it can be seen that equation (29) has a discontinuity at R0=ρ1/σ, which gives rise to one left branch that cuts the vertical axis at y∗=1/ρ1 when R0=0 and moves upward indefinitely as R0 approaches to R0=ρ1/σ, as well as another right branch for values greater than R0=ρ1/σ that comes from below cutting to the horizontal axis at R0=1/σ and asymptotically approaches to the horizontal line y∗=1 as R0 tends to infinity, see Figure 3.Figure 3
Bifurcation diagram of system (9). Locally asymptotically stable and unstable points have been indicated, respectively, with continuous and discontinuous lines. The region of biological sense Λ is indicated with gray color.Figure3, in fact, is the schematic representation of the results of the two previous theorems: both illustrate the behavior of the different families of equilibrium points as a function of R0 according to their positions in segment L defined above. Thus, for example, this figure shows how these families exhibit an abrupt change at the point 0,1/σ; in other words, in this nonhyperbolic point, a bifurcation occurs. Besides, it also indicates, in gray, the region of biological sense Λ. This bifurcation diagram corresponds to what is called as forward bifurcation [26]. It should be noted that this type of bifurcation is really a kind of transcritical bifurcation [37, 38], with bifurcation point 0,1/σ; that is, it consists of a horizontal branch that coincides with the R0 axis, which is stable before the bifurcation point and unstable after it, and another branch that is unstable below the bifurcation point and stable above it, see again Figure 3.


### 3.4. Nonexistence of Periodic Orbits and Global Stability

We will now discuss two results directly related to the global behavior of the trajectories of system (9). These are the nonexistence of periodic orbits in the region of biological sense Λ and the global stability of the endemic equilibrium point E∗ in the said region.


#### 3.4.1. Nonexistence of Periodic Orbits

The nonexistence of periodic orbits in the two-dimensional system (9) implies, of course, that neither system (5) has this type of trajectories. The first statement is shown in Theorem 3.Theorem 3.
System (9) does not have periodic orbits in the interior of Λ.Proof.
Let us identify the right-hand parts of the system of equations (9) as(30)Py,z≡R0−1y−1−σR0yz−R0−ρ1y2,Qy,z≡ρ2y−σR0−ρ1yz.
We propose the Dulac function(31)Φy,z=1y.
In this way, it can be easily verified that(32)∂ΦP∂y+∂ΦQ∂z=−σR0−ρ1+R0−ρ1<0,since as mentioned in Section 3, in Λ, it is true that σR0−ρ1>0 and as a consequence, in this region, R0−ρ1>0. Therefore, according to the Dulac criterion [37, 38], there are no periodic orbits in the interior of Λ.


#### 3.4.2. Global Stability of the Endemic Equilibrium Point

We will establish the conditions that the endemic equilibrium pointE∗ must meet to guarantee its global stability. This will be carried out by following Lyapunov’s second method or Lyapunov stability theorem [37, 38], as is shown in Theorem 4.Theorem 4.
The unique endemic equilibrium,E∗, is globally asymptotically stable in the interior of Λ.Proof.
We consider the common quadratic Lyapunov function(33)Vy,z=12y−y∗2+R01−σ2σR0−ρ1z−z∗2,which is continuous and a positive definite function, Vy,z>0, in a region that contains E∗=y∗,z∗ that is in turn located in Λ. It can be verified by inspection that Vy,z takes the value Vy,z=0 at E∗ so that the global minimum value of Vy,z is located at that point. Calculating the total derivative of (32) and using the two-dimensional system (9), we obtain(34)ddτVy,z=∂V∂ydydτ+∂V∂zdzdτ,=y−y∗R0−1y−1−σR0yz−R0−ρ1y2+R01−σσR0−ρ1z−z∗ρ2y−σR0−ρ1yz.
Taking into account the relations (21), the previous expression (34) takes the form(35)ddτVy,z=−y143+σR0−ρ1y−y∗2−y1−σR012y−y∗+z−z∗2<0.
Note in (35) that 1/43+σR0−ρ1=R0−ρ1−1/41−σR0>0 in region Λ; this is true since, as already mentioned, R0−ρ1=1−σR0+σR0−ρ1>0 in such a region so that R0−ρ1>1−σR0, and as 1−σR0>1/41−σR0, we finally obtain R0−ρ1−1/41−σR0>0. Thus, since also 1−σR0>0, the negative inequality in (35) holds for all y>0 and z>0, that is, for every point E∗ located in Λ (except the positive coordinate axes and the origin). Therefore, according to Lyapunov stability theorem, the endemic equilibrium E∗ is globally asymptotically stable in the interior of Λ.


### 3.5. Existence of Two First Integrals

In this section, we will refer to the three-dimensional system (5) instead of the two-dimensional one (9). It was found that system (5) has an interesting characteristic: it has two quantities that are invariant in time along their trajectories; that is, it has two first integrals.To show this, first note that starting with (2) and using (3), (4), and (8), we obtain(36)dx+y+zdτ=0,so that (8) is in itself worth the redundancy, a relation among the variables x, y, and z that remains constant. Therefore, this is a first integral, as shown below.Theorem 5.
Dynamic system (5) has a first integral, given by(37)H1x,y,z=x+y+z.Proof.
Calculating the total derivative of (37) and using the three equations (5), we have(38)ddτH1x,y,z=∂H1∂xdxdτ+∂H1∂ydydτ+∂H1∂zdzdτ,=−R0−ρ1xy−y+R0xy+σR0yz+ρ1y2+ρ2y−σR0−ρ1yz,=yρ2+ρ1x+y+z−1,=yρ1+ρ2−1=0,since (8) is true and also ρ1+ρ2=1. The proof is complete.On the other hand, note that the first and third equations (5) are only coupled by the variable y; therefore, by forming the quotient between them, we obtain(39)dxdz=−R0−ρ1xρ2−σR0−ρ1z.This expression can be solved, in principle, by the method of separating variables, giving rise to(40)1R0−ρ1lnx=1σR0−ρ1lnρ2−σR0−ρ1z+C,where C is a integration constant. It can be shown that, after some algebraic steps, it gives rise to(41)ρ2σR0−ρ1−zx−σR0−ρ1/R0−ρ1=K,K being a constant quantity. Consequently, we have obtained a quantity, dependent on the variables x and z, which remains constant for all τ. It turns out that this quantity is also a first integral, as shown below.Theorem 6.
Dynamic system (5) has another first integral, given by(42)H2x,z=ρ2σR0−ρ1−zx−σR0−ρ1/R0−ρ1.Proof.
Calculating the total derivative of (39) and using the first and third equations (5), we have(43)ddτH2x,z=∂H2∂xdxdτ+∂H2∂zdzdτ,=σR0−ρ1R0−ρ1ρ2σR0−ρ1−zx−σR0−ρ1/R0−ρ1−1R0−ρ1xy−x−σR0−ρ1/R0−ρ1ρ2y−σR0−ρ1yz,=ρ2yx−σR0−ρ1/R0−ρ1−σR0−ρ1yzx−σR0−ρ1/R0−ρ1−ρ2yx−σR0−ρ1/R0−ρ1+σR0−ρ1yzx−σR0−ρ1/R0−ρ1=0.
The proof is complete.It should be noted that expression (42) is reduced, when α=0, to the first integral of Lemma 1 reported in [28]. In this sense, the former is a generalization of the latter.


#### 3.5.1. Calculation of the Maximum ofyτ

The calculation of the first previous integrals allows us to estimate, for example, without having to know the corresponding exact solution, the maximum value that the infectious proportionyτ reaches. To do this, we first consider that we have the initial conditions x0=x0 and z0=z0. Thus, it can be shown that, from (41), we obtain(44)z=ρ2σR0−ρ1−ρ2σR0−ρ1−z0xx0σR0−ρ1/R0−ρ1.Therefore, using (8) and (44), we have the following expression for the variable y written only in terms of x:(45)y=1−x−ρ2σR0−ρ1+ρ2σR0−ρ1−z0xx0σR0−ρ1/R0−ρ1.Equating the derivative ofyτ with respect to τ to zero allows us to find its critical point; by deriving (45) in this way, it can be shown that this point is given by(46)x^=x0R0−ρ1x0ρ2−z0σR0−ρ1−R0−ρ1/1−σR0and also that it is a maximum. Therefore, substitution of (46) in (45) allows us to find the maximum of the infectious proportion yτ at any value of τ and in the presence of the parameter α, which we will denote as ymaxα; that is,(47)ymaxα=1−ρ2σR0−ρ1+R0−ρ1σR0−ρ1−1ρ2/R0−ρ1−z0σR0−ρ1/R0−ρ1x0σR0−ρ1/R0−ρ1R0−ρ1/1−σR0.It should be noted that equation (47) reduces completely, when the limit case α⟶0 is taken (that is, ρ1⟶0 and ρ2⟶1), to expression (ii) (a) of Theorem 1 in [28]. In this sense, our expression represents a generalization of that one for the case in which the parameter α of disease-induced death is taken into account.On the other hand, by means of equation (47), an estimate can be made of the difference between the maximum values of y in the presence and absence of the parameter α, that is, ymaxα−ymaxα=0. If we consider α small α≪1, we can approximate definitions (6) by ρ1≃α/δ and ρ2≃1. Therefore, for this situation, it can be shown that the said difference is given by(48)ymaxα−ymaxα=0≃1σ−1−αδσR01/R0−z0σ+z0α/δR0x0σ−α/δR01/1−σ−α/δ1−σR0−1σ−11/R0−z0σx0σ1/1−σ.A slight simplification of expressions (47) and (48) can be obtained, if we consider z0=0 in them; that is, there are no infected individuals at the beginning of the epidemic.


## 3.1. Equilibria and Nullclines

In order to find the equilibrium points of system (9), we follow the usual way of equating the left members of its equations with zero to obtain the values of the pair of variables y,z that satisfy this equality. As a result of this procedure, we obtain two equilibrium points. One is disease-free equilibria of the form(13)E00,z,with z∈ℝ; that is, the straight line y=0 with an infinite number of points that constitute the vertical coordinate axis of the yz plane. The second is the endemic equilibrium(14)E∗y∗,z∗=σR0−1σR0−ρ1,ρ2σR0−ρ1.So, the dynamics develop in regionΛ, and equilibrium points E∗ must be in the first quadrant of plane yz, that is, y∗≥0 and z∗>0; in other words, σR0−ρ1>0 and σR0−1≥0. Just when σR0−1=0, E∗ is located on the vetical z-axis; that is, it is part of the points E0.An alternative way to find the equilibrium points of a system is through the intersections of its nullclines [36, 37]. Since we will use these ideas later, we will begin this discussion below. According to system (9), we have the nullclines dy/dτ=0 (vertical directions), given by the functions(15)z=my+b,y=0,which consist of two straight lines: one of slope m=−R0−ρ1/1−σR0 and ordered to origin b=R0−1/1−σR0 and another that coincides with the vertical axis z. In addition, we have the nullclines dz/dτ=0 (horizontal directions), determined by(16)z=z∗,y=0,which are again two lines: one horizontal that cuts the vertical axis at point 0,ρ2/σR0−ρ1 and another vertical that also coincides with the z-axis. It should be noted that, as 0<σ<1 and R0>0, the quantity 1−σR0 is always positive; consequently, the slope of the vertical nullcline is only determined by the sign of R0−ρ1 and its ordinate to the origin by that of R0−1. As we are interested in the points E∗ located in region Λ, it must happen that σR0−ρ1>0,σR0−1≥0, and σR0−ρ1>0, as already mentioned previously; therefore, note that R0−ρ1=1−σR0+σR0−ρ1>0 and R0−1=1−σR0+σR0−1>0. These last two results imply that, in region Λ, the first vertical nullcline (15) must have a negative slope and its ordinate to the origin positive. In addition, it was also mentioned that, for E∗ to be in Λ, σR0−ρ1>0; therefore, the first horizontal nullcline (16) must be located in the upper half plane of plane yz.The intersections of different nullclines coincide with the equilibrium points. Note that, aty=0, nullclines overlap in both vertical and horizontal directions; therefore y=0 consists of an infinite line of equilibrium points, that is, the equilibria E0. On the other hand, another intersection between the other different nullclines gives rise to point E∗. Figure 2 illustrates the nullclines (15) and (16) in a configuration compatible with region Λ; it indicates horizontal nullclines in blue and verticals in red, vertical axis z (of infinite equilibrium points), as well as equilibrium points with black circles. This figure also shows the different directions that the vector field has along all the nullclines.Figure 2
Representation of nullclines of system (9): horizontal direction z=z∗ (in blue) and vertical direction z=my+b (in red); the vertical axis zy=0 is a combination of both. Equilibrium points are indicated by black circles. (a) σR0−ρ1>0,σR0−1>0,R0−ρ1>0, and R0−1>0; (b) σR0−ρ1>0,σR0−1<0,R0−ρ1>0, and R0−1>0.
(a)(b)The directions (left or right) of the horizontal nullclines are determined by substitutingz=z∗, first equation (16), in the first equation (9), which gives rise to(17)dydτ=R0−ρ1yy∗−y.To determine the positive or negative sign of (17), it is required to know the signs of y∗ (which depends in turn on the signs of σR0−1 and σR0−ρ1), the factor R0−ρ1, the values of y, and their order relation with y∗. Since Λ is satisfied that σR0−ρ1>0,σR0−1≥0 (that is, y∗≥0) and R0−ρ1>0(z∗ is on the nonnegative z-axis), the horizontal nullcline z=z∗ cuts the nonnegative vertical axis at point 0,z∗ and the nullcline of vertical directions, first equation (15), at y∗,z∗, which is an equilibrium point. Therefore, according to (16), if y∗<y<∞, then dy/dτ<0; that is, in this interval, the vector field is directed to the left if we are to the right of y∗,z∗; when 0<y<y∗, then dy/dτ>0 so that the vector field points to the right if we are to the left of y∗,z∗ and to the right of 0,z∗, while if −∞<y<0, then dy/dτ<0; that is, the direction of the field is to the left if we are to the left of 0,z∗. This result indicates that, on horizontal nullcline z=z∗, the vector field leaves point 0,z∗ and enters y∗,z∗.The abovementioned analysis can be repeated for all combinations of the quantitiesσR0−1,σR0−ρ1 (which determine the sign of y∗) and R0−ρ1. It can be shown that there are four possible cases. For y∗≥0, if R0−ρ1>0, then vector field on the horizontal nullcline z=z∗ leaves 0,z∗ and enters y∗,z∗, and when R0−ρ1<0, it enters 0,z∗ and leaves y∗,z∗; Figure 2(a) illustrates the first case. Also, for y∗<0, if R0−ρ1>0, then field enters 0,z∗ and leaves y∗,z∗, and when R0−ρ1<0, it leaves 0,z∗ and enters y∗,z∗; Figure 2(b) shows the third case. Note that only in the first case the equilibrium points are in region Λ.On the other hand, the directions (up or down) of the vertical nullclines are determined by replacing the vertical nullcline, first equation (15), in the second equation (9), which allows obtaining(18)dzdτ=−R0−ρ1σR0−ρ11−σR0yy∗−y.In order to determine the positive or negative sign of (18), we need to know the signs of y∗ (that depends on the signs of σR0−1 and σR0−ρ1), the factor R0−ρ1, the values of y, and their relation with y∗. Since in Λ, we have that σR0−ρ1>0,σR0−1≥0,R0−ρ1>0 and R0−1>0, the first vertical nullcline (16) cuts horizontal nullcline z=z∗ at the equilibrium point y∗,z∗ and the positive vertical axis at the point 0,b and has a negative slope. Therefore, according to (18), if y∗<y<∞, then dz/dτ>0; that is, in this interval, the direction of the vector field is up, while if 0<y<y∗, then dz/dτ<0 so that, in this interval, directions are down, and when −∞<y<0,dz/dτ>0, in this interval, directions are up.A complete analysis of the directions (up or down) of this vertical nullcline requires taking into account all combinations of signs of the quantitiesσR0−1,σR0−ρ1 and R0−ρ1. It can also be shown that there are four possible cases. For y∗>0, if R0−ρ1 and σR0−ρ1 are of opposite signs, then the direction of vector field on this vertical nullcline is down in the interval −∞<y<0, is up in 0<y<y∗, and is down in y∗<y<∞, whereas if R0−ρ1 and σR0−ρ1 are of the same signs, then the vector field is up in −∞<y<0, down in 0<y<y∗, and up in y∗<y<∞. Figure 2(a) illustrates the second case. Besides, for y∗<0, if R0−ρ1 and σR0−ρ1 are of opposite signs, then the vector field is down in −∞<y<y∗, up in y∗<y<0, and down in 0<y<∞, whereas if R0−ρ1 and σR0−ρ1 are of the same signs, then the vector field is up in −∞<y<y∗, down in y∗<y<0, and up in 0<y<∞. Figure 2(b) shows the fourth case. Note that only in the second case the equilibrium points are in region Λ.


## 3.2. Local Stability

In order to determine the stability of these families of equilibrium points, we need to calculate the Jacobian matrix of system (9) and evaluate it on them. This matrix overall is given as(19)Jy,z=−1−σR0z−2R0−ρ1y+R0−1−1−σR0yρ2−σR0−ρ1z−σR0−ρ1y.To perform the local stability analysis, we will consider separately the disease-free equilibria and the endemic equilibrium cases.


### 3.2.1. The Disease-Free Equilibria,E0=0,z

The Jacobian matrix (19) evaluated in the disease-free equilibria E0 takes the simple form(20)J0,z=−1−σR0z+R0−10ρ2−σR0−ρ1z0,whose eigenvalues are λ1=0 and λ2=−1−σR0z+R0−1. Note that these equilibrium points are nonhyperbolic. Due to this nature of the points E0 and the infinite number of them, it is quite difficult to determine their stability. One way to visualize the dynamics around them is through a geometric analysis of the nullclines of system (9). This allows us to find the directions left or right of the vector field along the horizontal nullclines and its directions up or down on the vertical nullclines; in particular, this gives a glimpse of how it is around the equilibrium points. Of particular interest, here are the results of the analysis performed in the previous section of the directions of the vector field on the horizontal nullclines z=z∗ around their point of intersection 0,z∗ with the vertical axis z. Those results, rewritten in terms of the ideas of stability, can be written in the following way: for y∗>0, if R0−ρ1>0, then point 0,z∗ is unstable, and when R0−ρ1<0, it is stable, whereas for y∗<0, if R0−ρ1>0, then point 0,z∗ is stable, and when R0−ρ1<0, it is unstable, see Figure 2. Since we have an infinite number of points along the vertical axis z, the stability along a horizontal line that cuts one of them, say 0,z∗, must be identical to the stability found at its intersection with the horizontal nullcline z=z∗. It should be noted that the analysis based on the method of the nullclines reveals that points 0,z∗ are stable if z>b and unstable if z<b, see again Figure 2.Since0<σ<1 and 0≤ρ1<1, we can establish the relations 0≤ρ1<ρ1/σ and ρ1/σ<1/σ; besides that, we know that 1<1/σ; therefore, the order relation between these quantities is 0≤ρ1<ρ1/σ<1/σ<∞, which will be called segment L. On the other hand, it has just been mentioned that the stability of the nonhyperbolic points 0,z∗ depends on the signs of σR0−ρ1,σR0−1 and R0−ρ1, or in other words, it depends of the location that R0 has in L.Thus, whenσR0−ρ1>0,σR0−1≥0 (which determine that y∗>0) and R0−ρ1>0, according to segment L, the intersection of these three subintervals is 1/σ<R0<∞. Since the previous stability result establishes that, in this situation, points 0,z∗ are unstable, we conclude that the equilibrium points E0 are unstable when 1/σ<R0<∞. On other hand, when σR0−ρ1>0,σR0−1<0 and R0−ρ1>0, according to L, intersection is ρ1/σ<R0<1/σ; since points 0,z∗ are stable here, we have that points E0 are unstable when ρ1/σ<R0<1/σ. By means of an analysis similar to the previous one, it can be shown that points E0 are unstable in ρ1<R0≤ρ1/σ and stable in 0≤R0<ρ1. These obtained results can be summarized as follows.Theorem 1.
Disease-free equilibrium pointsE0, which constitute the vertical axis z, are nonhyperbolic. These points are(i)
unstable, if1/σ<R0<∞ and ρ1<R0<ρ1/σ(ii)
stable, ifρ1/σ<R0<1/σ and 0<R0<ρ1Since whenα=0, our model is reduced to the conventional SIRI model, in this situation ρ1=0, so from Theorem 1, it follows that the nonhyperbolic equilibrium points E0 are stable if 0<R0<1/σ and unstable if 1/σ<R0<∞.


### 3.2.2. The Endemic Equilibrium,E∗=y∗,z∗

By substituting the endemic equilibriumE∗=y∗,z∗ in system (9), the identities are obtained.(21)−1−σR0z∗=R0−ρ1y∗−R0−1,ρ2−σR0−ρ1z∗=0,which allows writing the Jacobian matrix (19), evaluated in such a point, as(22)Jy∗,z∗=−R0−ρ1y∗−1−σR0y∗0−σR0−ρ1y∗.The trace of (22) is given by(23)TrJ=2ρ1−σ+1R0y∗=2ρ1−σ+1R0σR0−1σR0−ρ1,and its determinant is(24)DetJ=R0−ρ1σR0−ρ1y∗2=R0−ρ1σR0−12σR0−ρ1,whereas it can be shown that its discriminant, defined as Δ=TrJ2−4DetJ, is(25)Δ=−1−σR0y∗2=1−σ2R02σR0−12σR0−ρ12,which is a nonnegative quantity. According to the trace-determinant plane method [36–38], the analysis of matrix (22) indicates that points E∗, since Δ≥0, can be saddle points (if DetJ<0) and degenerate equilibrium (if DetJ=0) or nodes of two “straight-line” solutions (if DetJ>0), the latter pair being stable or unstable if TrJ<0 or TrJ>0, respectively. Furthermore, the eigenvalues of (22) are given by(26)λ1,2=12TrJ±Δ,which can be explicitly written as(27)λ1=−R0−ρ1y∗=−R0−ρ1σR0−1σR0−ρ1,(28)λ2=−σR0−ρ1y∗=−σR0−1.Note that eigenvalues (27) and (28) are functions of R0 and their signs will depend on the location of this quantity in segment L already indicated above. In general, they are different from zero, so the equilibrium points E∗ are hyperbolic.Next, we will analyze the signs thatλ1 and λ2 could have in L. According to (27) and (28), λ1<0 and λ2<0, if σR0−1>0,R0−ρ1>0, and σR0−ρ1>0; the intersection of these three subintervals occurs at 1/σ<R0<∞, so in this region, the equilibrium points E∗ correspond to nodes of two “straight-line” solutions that are locally asymptotically stable. In this case, in which both eigenvalues are negative, points E∗ lie within region Λ; however, for other combinations of signs of the eigenvalues λ1 and λ2, points E∗ lie outside of Λ. Thus, it can be shown that λ1>0 and λ2>0 when ρ1/σ<R0<1/σ and 0≤R0<ρ1 so that, in these regions, points E∗ correspond to unstable nodes of two “straight lines.” Additionally, λ1>0 and λ2<0, if ρ1<R0<ρ1/σ, so there points E∗ are saddle points; the case λ1<0 and λ2>0 does not exist because the intersection of the corresponding subintervals gives rise to the empty set. Moreover, λ1=0 and λ2>0, when R0=ρ1, so here, E∗ is a singly degenerate equilibrium. Lastly, λ1=λ2=0, when R0=1/σ, so E∗ is a doubly degenerate equilibrium. These last two cases correspond to nonhyperbolic equilibrium points. The previous results are summarized in Theorem 2.Theorem 2.
The endemic equilibrium pointsE∗ are in general hyperbolic. Such points are(i)
nodes (of two “straight-line” solutions) locally asymptotically stable, if1/σ<R0<∞(ii)
a doubly degenerate equilibrium, ifR0=1/σ (in which case it is nonhyperbolic)(iii)
unstable nodes (of two “straight-line” solutions), if0≤R0<ρ1 and ρ1/σ<R0<1/σ(iv)
saddle points, ifρ1<R0<ρ1/σ(v)
a singly degenerate equilibrium, ifR0=ρ1 (in which case it is nonhyperbolic)(vi)
unstable nodes (of two “straight-line” solutions), if0≤R0<ρ1It is important to note that, in Theorem2, subinterval 1/σ<R0<∞ (that of condition i)) is the only one that guarantees that points E∗ are in region Λ. Again, since when α=0, our model is reduced to the conventional SIRI model, in this situation, ρ1=0, so for this case, from Theorem 2, it follows that the equilibrium points E∗ are nodes (of two “straight-line” solutions) locally asymptotically stable if 1/σ<R0<∞ and unstable if 0<R0<1/σ.


## 3.2.1. The Disease-Free Equilibria,E0=0,z

The Jacobian matrix (19) evaluated in the disease-free equilibria E0 takes the simple form(20)J0,z=−1−σR0z+R0−10ρ2−σR0−ρ1z0,whose eigenvalues are λ1=0 and λ2=−1−σR0z+R0−1. Note that these equilibrium points are nonhyperbolic. Due to this nature of the points E0 and the infinite number of them, it is quite difficult to determine their stability. One way to visualize the dynamics around them is through a geometric analysis of the nullclines of system (9). This allows us to find the directions left or right of the vector field along the horizontal nullclines and its directions up or down on the vertical nullclines; in particular, this gives a glimpse of how it is around the equilibrium points. Of particular interest, here are the results of the analysis performed in the previous section of the directions of the vector field on the horizontal nullclines z=z∗ around their point of intersection 0,z∗ with the vertical axis z. Those results, rewritten in terms of the ideas of stability, can be written in the following way: for y∗>0, if R0−ρ1>0, then point 0,z∗ is unstable, and when R0−ρ1<0, it is stable, whereas for y∗<0, if R0−ρ1>0, then point 0,z∗ is stable, and when R0−ρ1<0, it is unstable, see Figure 2. Since we have an infinite number of points along the vertical axis z, the stability along a horizontal line that cuts one of them, say 0,z∗, must be identical to the stability found at its intersection with the horizontal nullcline z=z∗. It should be noted that the analysis based on the method of the nullclines reveals that points 0,z∗ are stable if z>b and unstable if z<b, see again Figure 2.Since0<σ<1 and 0≤ρ1<1, we can establish the relations 0≤ρ1<ρ1/σ and ρ1/σ<1/σ; besides that, we know that 1<1/σ; therefore, the order relation between these quantities is 0≤ρ1<ρ1/σ<1/σ<∞, which will be called segment L. On the other hand, it has just been mentioned that the stability of the nonhyperbolic points 0,z∗ depends on the signs of σR0−ρ1,σR0−1 and R0−ρ1, or in other words, it depends of the location that R0 has in L.Thus, whenσR0−ρ1>0,σR0−1≥0 (which determine that y∗>0) and R0−ρ1>0, according to segment L, the intersection of these three subintervals is 1/σ<R0<∞. Since the previous stability result establishes that, in this situation, points 0,z∗ are unstable, we conclude that the equilibrium points E0 are unstable when 1/σ<R0<∞. On other hand, when σR0−ρ1>0,σR0−1<0 and R0−ρ1>0, according to L, intersection is ρ1/σ<R0<1/σ; since points 0,z∗ are stable here, we have that points E0 are unstable when ρ1/σ<R0<1/σ. By means of an analysis similar to the previous one, it can be shown that points E0 are unstable in ρ1<R0≤ρ1/σ and stable in 0≤R0<ρ1. These obtained results can be summarized as follows.Theorem 1.
Disease-free equilibrium pointsE0, which constitute the vertical axis z, are nonhyperbolic. These points are(i)
unstable, if1/σ<R0<∞ and ρ1<R0<ρ1/σ(ii)
stable, ifρ1/σ<R0<1/σ and 0<R0<ρ1Since whenα=0, our model is reduced to the conventional SIRI model, in this situation ρ1=0, so from Theorem 1, it follows that the nonhyperbolic equilibrium points E0 are stable if 0<R0<1/σ and unstable if 1/σ<R0<∞.


## 3.2.2. The Endemic Equilibrium,E∗=y∗,z∗

By substituting the endemic equilibriumE∗=y∗,z∗ in system (9), the identities are obtained.(21)−1−σR0z∗=R0−ρ1y∗−R0−1,ρ2−σR0−ρ1z∗=0,which allows writing the Jacobian matrix (19), evaluated in such a point, as(22)Jy∗,z∗=−R0−ρ1y∗−1−σR0y∗0−σR0−ρ1y∗.The trace of (22) is given by(23)TrJ=2ρ1−σ+1R0y∗=2ρ1−σ+1R0σR0−1σR0−ρ1,and its determinant is(24)DetJ=R0−ρ1σR0−ρ1y∗2=R0−ρ1σR0−12σR0−ρ1,whereas it can be shown that its discriminant, defined as Δ=TrJ2−4DetJ, is(25)Δ=−1−σR0y∗2=1−σ2R02σR0−12σR0−ρ12,which is a nonnegative quantity. According to the trace-determinant plane method [36–38], the analysis of matrix (22) indicates that points E∗, since Δ≥0, can be saddle points (if DetJ<0) and degenerate equilibrium (if DetJ=0) or nodes of two “straight-line” solutions (if DetJ>0), the latter pair being stable or unstable if TrJ<0 or TrJ>0, respectively. Furthermore, the eigenvalues of (22) are given by(26)λ1,2=12TrJ±Δ,which can be explicitly written as(27)λ1=−R0−ρ1y∗=−R0−ρ1σR0−1σR0−ρ1,(28)λ2=−σR0−ρ1y∗=−σR0−1.Note that eigenvalues (27) and (28) are functions of R0 and their signs will depend on the location of this quantity in segment L already indicated above. In general, they are different from zero, so the equilibrium points E∗ are hyperbolic.Next, we will analyze the signs thatλ1 and λ2 could have in L. According to (27) and (28), λ1<0 and λ2<0, if σR0−1>0,R0−ρ1>0, and σR0−ρ1>0; the intersection of these three subintervals occurs at 1/σ<R0<∞, so in this region, the equilibrium points E∗ correspond to nodes of two “straight-line” solutions that are locally asymptotically stable. In this case, in which both eigenvalues are negative, points E∗ lie within region Λ; however, for other combinations of signs of the eigenvalues λ1 and λ2, points E∗ lie outside of Λ. Thus, it can be shown that λ1>0 and λ2>0 when ρ1/σ<R0<1/σ and 0≤R0<ρ1 so that, in these regions, points E∗ correspond to unstable nodes of two “straight lines.” Additionally, λ1>0 and λ2<0, if ρ1<R0<ρ1/σ, so there points E∗ are saddle points; the case λ1<0 and λ2>0 does not exist because the intersection of the corresponding subintervals gives rise to the empty set. Moreover, λ1=0 and λ2>0, when R0=ρ1, so here, E∗ is a singly degenerate equilibrium. Lastly, λ1=λ2=0, when R0=1/σ, so E∗ is a doubly degenerate equilibrium. These last two cases correspond to nonhyperbolic equilibrium points. The previous results are summarized in Theorem 2.Theorem 2.
The endemic equilibrium pointsE∗ are in general hyperbolic. Such points are(i)
nodes (of two “straight-line” solutions) locally asymptotically stable, if1/σ<R0<∞(ii)
a doubly degenerate equilibrium, ifR0=1/σ (in which case it is nonhyperbolic)(iii)
unstable nodes (of two “straight-line” solutions), if0≤R0<ρ1 and ρ1/σ<R0<1/σ(iv)
saddle points, ifρ1<R0<ρ1/σ(v)
a singly degenerate equilibrium, ifR0=ρ1 (in which case it is nonhyperbolic)(vi)
unstable nodes (of two “straight-line” solutions), if0≤R0<ρ1It is important to note that, in Theorem2, subinterval 1/σ<R0<∞ (that of condition i)) is the only one that guarantees that points E∗ are in region Λ. Again, since when α=0, our model is reduced to the conventional SIRI model, in this situation, ρ1=0, so for this case, from Theorem 2, it follows that the equilibrium points E∗ are nodes (of two “straight-line” solutions) locally asymptotically stable if 1/σ<R0<∞ and unstable if 0<R0<1/σ.


## 3.3. Presence of a Forward Bifurcation

The analysis that allowed to formulate the two previous theorems shows that the coordinates of both pointsE0 and E∗ are a function of the reproductive number R0. In particular, it is of great interest to know how the family of equilibrium points y∗ (which is the abscissa of points E∗) varies with R0 and to visualize such dependence in a diagram R0y∗; this is important because R0 could play the role of a bifurcation parameter. For our study problem, the y∗ families consist of two curves: one horizontal line y∗=0 of disease-free equilibrium points located on the R0 axis and another determined by(29)y∗R0=σR0−1σR0−ρ1,of endemic equilibrium points. By inspection, it can be seen that equation (29) has a discontinuity at R0=ρ1/σ, which gives rise to one left branch that cuts the vertical axis at y∗=1/ρ1 when R0=0 and moves upward indefinitely as R0 approaches to R0=ρ1/σ, as well as another right branch for values greater than R0=ρ1/σ that comes from below cutting to the horizontal axis at R0=1/σ and asymptotically approaches to the horizontal line y∗=1 as R0 tends to infinity, see Figure 3.Figure 3
Bifurcation diagram of system (9). Locally asymptotically stable and unstable points have been indicated, respectively, with continuous and discontinuous lines. The region of biological sense Λ is indicated with gray color.Figure3, in fact, is the schematic representation of the results of the two previous theorems: both illustrate the behavior of the different families of equilibrium points as a function of R0 according to their positions in segment L defined above. Thus, for example, this figure shows how these families exhibit an abrupt change at the point 0,1/σ; in other words, in this nonhyperbolic point, a bifurcation occurs. Besides, it also indicates, in gray, the region of biological sense Λ. This bifurcation diagram corresponds to what is called as forward bifurcation [26]. It should be noted that this type of bifurcation is really a kind of transcritical bifurcation [37, 38], with bifurcation point 0,1/σ; that is, it consists of a horizontal branch that coincides with the R0 axis, which is stable before the bifurcation point and unstable after it, and another branch that is unstable below the bifurcation point and stable above it, see again Figure 3.


## 3.4. Nonexistence of Periodic Orbits and Global Stability

We will now discuss two results directly related to the global behavior of the trajectories of system (9). These are the nonexistence of periodic orbits in the region of biological sense Λ and the global stability of the endemic equilibrium point E∗ in the said region.


### 3.4.1. Nonexistence of Periodic Orbits

The nonexistence of periodic orbits in the two-dimensional system (9) implies, of course, that neither system (5) has this type of trajectories. The first statement is shown in Theorem 3.Theorem 3.
System (9) does not have periodic orbits in the interior of Λ.Proof.
Let us identify the right-hand parts of the system of equations (9) as(30)Py,z≡R0−1y−1−σR0yz−R0−ρ1y2,Qy,z≡ρ2y−σR0−ρ1yz.
We propose the Dulac function(31)Φy,z=1y.
In this way, it can be easily verified that(32)∂ΦP∂y+∂ΦQ∂z=−σR0−ρ1+R0−ρ1<0,since as mentioned in Section 3, in Λ, it is true that σR0−ρ1>0 and as a consequence, in this region, R0−ρ1>0. Therefore, according to the Dulac criterion [37, 38], there are no periodic orbits in the interior of Λ.


### 3.4.2. Global Stability of the Endemic Equilibrium Point

We will establish the conditions that the endemic equilibrium pointE∗ must meet to guarantee its global stability. This will be carried out by following Lyapunov’s second method or Lyapunov stability theorem [37, 38], as is shown in Theorem 4.Theorem 4.
The unique endemic equilibrium,E∗, is globally asymptotically stable in the interior of Λ.Proof.
We consider the common quadratic Lyapunov function(33)Vy,z=12y−y∗2+R01−σ2σR0−ρ1z−z∗2,which is continuous and a positive definite function, Vy,z>0, in a region that contains E∗=y∗,z∗ that is in turn located in Λ. It can be verified by inspection that Vy,z takes the value Vy,z=0 at E∗ so that the global minimum value of Vy,z is located at that point. Calculating the total derivative of (32) and using the two-dimensional system (9), we obtain(34)ddτVy,z=∂V∂ydydτ+∂V∂zdzdτ,=y−y∗R0−1y−1−σR0yz−R0−ρ1y2+R01−σσR0−ρ1z−z∗ρ2y−σR0−ρ1yz.
Taking into account the relations (21), the previous expression (34) takes the form(35)ddτVy,z=−y143+σR0−ρ1y−y∗2−y1−σR012y−y∗+z−z∗2<0.
Note in (35) that 1/43+σR0−ρ1=R0−ρ1−1/41−σR0>0 in region Λ; this is true since, as already mentioned, R0−ρ1=1−σR0+σR0−ρ1>0 in such a region so that R0−ρ1>1−σR0, and as 1−σR0>1/41−σR0, we finally obtain R0−ρ1−1/41−σR0>0. Thus, since also 1−σR0>0, the negative inequality in (35) holds for all y>0 and z>0, that is, for every point E∗ located in Λ (except the positive coordinate axes and the origin). Therefore, according to Lyapunov stability theorem, the endemic equilibrium E∗ is globally asymptotically stable in the interior of Λ.


## 3.4.1. Nonexistence of Periodic Orbits

The nonexistence of periodic orbits in the two-dimensional system (9) implies, of course, that neither system (5) has this type of trajectories. The first statement is shown in Theorem 3.Theorem 3.
System (9) does not have periodic orbits in the interior of Λ.Proof.
Let us identify the right-hand parts of the system of equations (9) as(30)Py,z≡R0−1y−1−σR0yz−R0−ρ1y2,Qy,z≡ρ2y−σR0−ρ1yz.
We propose the Dulac function(31)Φy,z=1y.
In this way, it can be easily verified that(32)∂ΦP∂y+∂ΦQ∂z=−σR0−ρ1+R0−ρ1<0,since as mentioned in Section 3, in Λ, it is true that σR0−ρ1>0 and as a consequence, in this region, R0−ρ1>0. Therefore, according to the Dulac criterion [37, 38], there are no periodic orbits in the interior of Λ.


## 3.4.2. Global Stability of the Endemic Equilibrium Point

We will establish the conditions that the endemic equilibrium pointE∗ must meet to guarantee its global stability. This will be carried out by following Lyapunov’s second method or Lyapunov stability theorem [37, 38], as is shown in Theorem 4.Theorem 4.
The unique endemic equilibrium,E∗, is globally asymptotically stable in the interior of Λ.Proof.
We consider the common quadratic Lyapunov function(33)Vy,z=12y−y∗2+R01−σ2σR0−ρ1z−z∗2,which is continuous and a positive definite function, Vy,z>0, in a region that contains E∗=y∗,z∗ that is in turn located in Λ. It can be verified by inspection that Vy,z takes the value Vy,z=0 at E∗ so that the global minimum value of Vy,z is located at that point. Calculating the total derivative of (32) and using the two-dimensional system (9), we obtain(34)ddτVy,z=∂V∂ydydτ+∂V∂zdzdτ,=y−y∗R0−1y−1−σR0yz−R0−ρ1y2+R01−σσR0−ρ1z−z∗ρ2y−σR0−ρ1yz.
Taking into account the relations (21), the previous expression (34) takes the form(35)ddτVy,z=−y143+σR0−ρ1y−y∗2−y1−σR012y−y∗+z−z∗2<0.
Note in (35) that 1/43+σR0−ρ1=R0−ρ1−1/41−σR0>0 in region Λ; this is true since, as already mentioned, R0−ρ1=1−σR0+σR0−ρ1>0 in such a region so that R0−ρ1>1−σR0, and as 1−σR0>1/41−σR0, we finally obtain R0−ρ1−1/41−σR0>0. Thus, since also 1−σR0>0, the negative inequality in (35) holds for all y>0 and z>0, that is, for every point E∗ located in Λ (except the positive coordinate axes and the origin). Therefore, according to Lyapunov stability theorem, the endemic equilibrium E∗ is globally asymptotically stable in the interior of Λ.


## 3.5. Existence of Two First Integrals

In this section, we will refer to the three-dimensional system (5) instead of the two-dimensional one (9). It was found that system (5) has an interesting characteristic: it has two quantities that are invariant in time along their trajectories; that is, it has two first integrals.To show this, first note that starting with (2) and using (3), (4), and (8), we obtain(36)dx+y+zdτ=0,so that (8) is in itself worth the redundancy, a relation among the variables x, y, and z that remains constant. Therefore, this is a first integral, as shown below.Theorem 5.
Dynamic system (5) has a first integral, given by(37)H1x,y,z=x+y+z.Proof.
Calculating the total derivative of (37) and using the three equations (5), we have(38)ddτH1x,y,z=∂H1∂xdxdτ+∂H1∂ydydτ+∂H1∂zdzdτ,=−R0−ρ1xy−y+R0xy+σR0yz+ρ1y2+ρ2y−σR0−ρ1yz,=yρ2+ρ1x+y+z−1,=yρ1+ρ2−1=0,since (8) is true and also ρ1+ρ2=1. The proof is complete.On the other hand, note that the first and third equations (5) are only coupled by the variable y; therefore, by forming the quotient between them, we obtain(39)dxdz=−R0−ρ1xρ2−σR0−ρ1z.This expression can be solved, in principle, by the method of separating variables, giving rise to(40)1R0−ρ1lnx=1σR0−ρ1lnρ2−σR0−ρ1z+C,where C is a integration constant. It can be shown that, after some algebraic steps, it gives rise to(41)ρ2σR0−ρ1−zx−σR0−ρ1/R0−ρ1=K,K being a constant quantity. Consequently, we have obtained a quantity, dependent on the variables x and z, which remains constant for all τ. It turns out that this quantity is also a first integral, as shown below.Theorem 6.
Dynamic system (5) has another first integral, given by(42)H2x,z=ρ2σR0−ρ1−zx−σR0−ρ1/R0−ρ1.Proof.
Calculating the total derivative of (39) and using the first and third equations (5), we have(43)ddτH2x,z=∂H2∂xdxdτ+∂H2∂zdzdτ,=σR0−ρ1R0−ρ1ρ2σR0−ρ1−zx−σR0−ρ1/R0−ρ1−1R0−ρ1xy−x−σR0−ρ1/R0−ρ1ρ2y−σR0−ρ1yz,=ρ2yx−σR0−ρ1/R0−ρ1−σR0−ρ1yzx−σR0−ρ1/R0−ρ1−ρ2yx−σR0−ρ1/R0−ρ1+σR0−ρ1yzx−σR0−ρ1/R0−ρ1=0.
The proof is complete.It should be noted that expression (42) is reduced, when α=0, to the first integral of Lemma 1 reported in [28]. In this sense, the former is a generalization of the latter.


### 3.5.1. Calculation of the Maximum ofyτ

The calculation of the first previous integrals allows us to estimate, for example, without having to know the corresponding exact solution, the maximum value that the infectious proportionyτ reaches. To do this, we first consider that we have the initial conditions x0=x0 and z0=z0. Thus, it can be shown that, from (41), we obtain(44)z=ρ2σR0−ρ1−ρ2σR0−ρ1−z0xx0σR0−ρ1/R0−ρ1.Therefore, using (8) and (44), we have the following expression for the variable y written only in terms of x:(45)y=1−x−ρ2σR0−ρ1+ρ2σR0−ρ1−z0xx0σR0−ρ1/R0−ρ1.Equating the derivative ofyτ with respect to τ to zero allows us to find its critical point; by deriving (45) in this way, it can be shown that this point is given by(46)x^=x0R0−ρ1x0ρ2−z0σR0−ρ1−R0−ρ1/1−σR0and also that it is a maximum. Therefore, substitution of (46) in (45) allows us to find the maximum of the infectious proportion yτ at any value of τ and in the presence of the parameter α, which we will denote as ymaxα; that is,(47)ymaxα=1−ρ2σR0−ρ1+R0−ρ1σR0−ρ1−1ρ2/R0−ρ1−z0σR0−ρ1/R0−ρ1x0σR0−ρ1/R0−ρ1R0−ρ1/1−σR0.It should be noted that equation (47) reduces completely, when the limit case α⟶0 is taken (that is, ρ1⟶0 and ρ2⟶1), to expression (ii) (a) of Theorem 1 in [28]. In this sense, our expression represents a generalization of that one for the case in which the parameter α of disease-induced death is taken into account.On the other hand, by means of equation (47), an estimate can be made of the difference between the maximum values of y in the presence and absence of the parameter α, that is, ymaxα−ymaxα=0. If we consider α small α≪1, we can approximate definitions (6) by ρ1≃α/δ and ρ2≃1. Therefore, for this situation, it can be shown that the said difference is given by(48)ymaxα−ymaxα=0≃1σ−1−αδσR01/R0−z0σ+z0α/δR0x0σ−α/δR01/1−σ−α/δ1−σR0−1σ−11/R0−z0σx0σ1/1−σ.A slight simplification of expressions (47) and (48) can be obtained, if we consider z0=0 in them; that is, there are no infected individuals at the beginning of the epidemic.


## 3.5.1. Calculation of the Maximum ofyτ

The calculation of the first previous integrals allows us to estimate, for example, without having to know the corresponding exact solution, the maximum value that the infectious proportionyτ reaches. To do this, we first consider that we have the initial conditions x0=x0 and z0=z0. Thus, it can be shown that, from (41), we obtain(44)z=ρ2σR0−ρ1−ρ2σR0−ρ1−z0xx0σR0−ρ1/R0−ρ1.Therefore, using (8) and (44), we have the following expression for the variable y written only in terms of x:(45)y=1−x−ρ2σR0−ρ1+ρ2σR0−ρ1−z0xx0σR0−ρ1/R0−ρ1.Equating the derivative ofyτ with respect to τ to zero allows us to find its critical point; by deriving (45) in this way, it can be shown that this point is given by(46)x^=x0R0−ρ1x0ρ2−z0σR0−ρ1−R0−ρ1/1−σR0and also that it is a maximum. Therefore, substitution of (46) in (45) allows us to find the maximum of the infectious proportion yτ at any value of τ and in the presence of the parameter α, which we will denote as ymaxα; that is,(47)ymaxα=1−ρ2σR0−ρ1+R0−ρ1σR0−ρ1−1ρ2/R0−ρ1−z0σR0−ρ1/R0−ρ1x0σR0−ρ1/R0−ρ1R0−ρ1/1−σR0.It should be noted that equation (47) reduces completely, when the limit case α⟶0 is taken (that is, ρ1⟶0 and ρ2⟶1), to expression (ii) (a) of Theorem 1 in [28]. In this sense, our expression represents a generalization of that one for the case in which the parameter α of disease-induced death is taken into account.On the other hand, by means of equation (47), an estimate can be made of the difference between the maximum values of y in the presence and absence of the parameter α, that is, ymaxα−ymaxα=0. If we consider α small α≪1, we can approximate definitions (6) by ρ1≃α/δ and ρ2≃1. Therefore, for this situation, it can be shown that the said difference is given by(48)ymaxα−ymaxα=0≃1σ−1−αδσR01/R0−z0σ+z0α/δR0x0σ−α/δR01/1−σ−α/δ1−σR0−1σ−11/R0−z0σx0σ1/1−σ.A slight simplification of expressions (47) and (48) can be obtained, if we consider z0=0 in them; that is, there are no infected individuals at the beginning of the epidemic.


## 4. Numerical Results

In this section we will present some of the most representative numerical solutions of our model that allow us to illustrate the different results obtained previously. The parameter values that we will use to calculate them correspond to estimates made for the epidemic caused by SARS-CoV-2. To carry out the numerical analysis, it is necessary to know representative values of the parameters involved in our model, namely, theβ′, δ,α, and σ rates. Regarding this disease, COVID-19 patients in Wuhan showed that the average contagious period of SARS-CoV-2 was 20 days [39]. We choose a value of 14 days to the contagious period 1/δ of an infected person so that δ≃0.071day−1. Furthermore, the disease-induced death rate has been reported in Mexico of 0.0127 and 0.0259 day−1 [40], although more studies are lacking where disease-induced death rate is estimated in populations with a high prevalence of comorbidities or in old populations. We will consider that α≃0.1; however, in our calculations, we will vary this rate and consider that 0≤α<1. It is worth mentioning that such a percentage depends on the particular characteristics of the population of each country and may be a little higher, as in the case of Mexico [8], if it presents several comorbidities [6, 7]. A meta-analysis [41] of the basic reproductive number R0 to China reported a pooled estimation of 3.32(95% confidence interval: 2.81,3.82). Knowing this quantity is important because from it, the approximate values of the other two remaining rates that are relevant in our model can be obtained. We take a value of R0≃3.11 (as in [42]). Thus, according to the parameters indicated, from relation (7), it follows that β′≃0.532day−1. Finally, it has already been mentioned that the bifurcation occurs when R0=1/σ, so a first estimate of the relapse rate, although it is for this case, is σ≃0.322; precisely in this model, 0<σ<1, and in our calculations, this quantity will vary taking values belonging to this interval. The magnitudes of the indicated rates δ, α,β′, and σ are shown in Table 1.Table 1
Estimation of the values of the main parameters used in the model, which correspond to the SARS-CoV-2 coronavirus.
ParameterMeaning (unit)Estimated valuesβ′Transmission rateday−10.532δRecovery rateday−10.071σRelative susceptibility of recovered individuals (adimensional)VaryingαDisease−induced death rateday−1VaryingSince the purpose here is to show the effects produced in the population of infected individualsyτ basically by the change of the α and σ rates, we will take the quantities indicated in Table 1 only as a reference; that is, in the different cases that we will present below, we will make some of them take values close to or different (even by several orders of magnitude) from those reported there.


### 4.1. Effects on the Trajectories ofyτ due to Variations of the α and σ Rates

We illustrate how changing the ratesα and σ affects the dynamics of the solutions of system (9), in particular those of the proportion of infected individuals yτ. For the different cases to be considered, we will indicate the values of R0 (since it varies inversely with α, as is pointed out by (7)) and their corresponding bifurcations 1/σ. Also, for each of the graphs obtained, taking into account (47), we will indicate the magnitude that reaches its peak, namely, ymaxα.Let us consider a first case whereα varies, β′=0.532, δ=0.071 (as in Table 1), and σ=0.1 will be fixed. We have chosen this last magnitude for the recovery rate, so that R0 remains below the value of the bifurcation 1/σ=10 that is constant. Figure 4 shows how the graphs of yτ change for three different values of α; specifically, for α=0.0001 (blue line), R0=7.48 and ymaxα=0.63; when α=0.1 (red line), R0=3.11 and ymaxα=0.51; and for α=0.17 (green line), R0=2.21 and ymaxα=0.41. In this figure, it can be seen that as α grows, the peak of the graphs decreases, they widen and shift to the right, and they are also less biased. Note that, in the first case α=0.0001, the rate α is so small that the effect produced is as if it were not there; in the absence of αα=0, R0 (for fixed β′ and δ) has its greatest value and the graph of y has its highest peak. Apparently, the inclusion in the model of the rate α produces a more adequate description of the evolution of infected individuals, since its absence causes overestimated values of R0 and the peak of its graph.Figure 4
Evolution ofyτ for the fixed parameters β′=0.532, δ=0.071, and σ=0.1. α=0.0001, blue line; α=0.1, red line; and α=0.17, green line.Now we consider a second case whereσ changes, β′=0.532, δ=0.071 (as in Table 1), and α=0.1 will be fixed. Note that, in this situation, for each σ, we will have different bifurcation values 1/σ so that if σ is very close to zero, the bifurcation point is much greater than one, whereas if it is very close to one, this point is much less than one. We choose the magnitudes of σ that allow to R0=3.11 (which does not change, since β′, δ, and α are constant) to be below or above the value of the bifurcation. Figure 5 illustrates how the temporal evolution of yτ changes for three different magnitudes of σ; particularly, when σ=0.1 (blue line), 1/σ=10 and ymaxα=0.51; when σ=0.3 (red line) 1/σ=3.33 and ymaxα=0.58; and for σ=0.4 (green line), 1/σ=2.5 and ymaxα=0.62. This figure shows how as σ grows, the graphs of yτ change: their peaks increase, widen, and shift slightly to the right. It should be noted that when σ=0.1 (blue line) and σ=0.3 (red line), R0 is below their corresponding bifurcation values and both trajectories tend to a disease-free equilibrium point E0 as τ tends to infinity, while when σ=0.4 (green line), R0 is above its bifurcation value and the trajectory tends to an endemic equilibrium point E∗ as τ tends to infinity.Figure 5
Evolution ofyτ for the fixed parameters β′=0.532, δ=0.071, and α=0.1. σ=0.1, blue line; σ=0.3, red line; and σ=0.4, green line.Finally, let us consider a third case, which certainly, could be considered as a combination of the previous two. Here,α varies, β′=0.532, and δ=0.071 (as in Table 1), but now, σ=0.17 will be a different fixed value. We have considered this latter magnitude in order that R0, as α increases, may exceed the value of the bifurcation 1/σ=5.88 that does not change. The different rates α selected make R0 take values below or above the bifurcation 1/σ. The effect that this has on the y curves is illustrated in Figure 6. Thus, in the figure, when α=0.0001 (blue line), R0=7.48, and ymaxα=0.65, the curve is above the value of the bifurcation and tends to an endemic equilibrium point as time tends to infinity. If α=0.1 (red line), R0=3.11 and ymaxα=0.53, whereas when α=0.17 (green line), R0=2.21 and ymaxα=0.43; in both cases, R0 is below the value of the bifurcation, so the curves tend to disease-free equilibrium points as time tends to infinity. Note that variations in the rate α can produce changes in the nature of the equilibrium points towards which the y curves tend: as α increases, the curves go from tending towards an endemic equilibrium point E∗ (when α=0.0001) to others free of disease E0 (if α=0.1 or α=0.17).Figure 6
Evolution ofyt for the fixed parameters β′=0.532, δ=0.071, and σ=0.17. α=0.0001, blue line; α=0.1, red line; and α=0.17, green line.


### 4.2. Phase Portraits before and during Forward Bifurcation

We present here some representative phase portraits of dynamic system (9), for different values of R0 between 0 and just over 1/σ, the bifurcation value, taking as reference the diagram in Figure 3. Figures 7–12 illustrate the different phase portraits that were calculated. For the realization of these graphs, the rates δ, α, and σ were considered as constants: the first with its value given in Table 1 and the other two with magnitudes equal to 0.1 and 0.17, respectively; the rate β′ was varied. The quantity R0 can increase its value, according to (7), if δ is kept constant and α decreases or β′ increases. It has been pointed out that the absence of the parameter α causes an overestimated magnitude of R0, so we will keep its value fixed and we will only vary β′ in an increasing way. For rates δ, α, and σ already mentioned, in accordance with (6), we have that ρ1≃0.5848,ρ2≃0.4152, while ρ1/σ≃3.440 and 1/σ≃5.882. For each value of R0, we will also indicate the values of their corresponding disease-free E0=0,z and endemic E∗=y∗,z∗ equilibrium points, expressions (13) and (14), respectively; since the former are infinite, we will only present one of them, namely, point E0=0,z∗.Figure 7
Phase plane of system (9) for 0≤R0<ρ1, with R0≃0.292β′=0.05; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,−0.776 is stable, and the endemic E∗≃1.776,−0.776 is an unstable node of two “straight-line” solutions; both points are indicated by a black circle.Figure 8
Phase plane of system (9), for the case where R0=ρ1≃0.585β′=0.1; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. In this case, the first vertical nullcline (15) is a line of zero slope whose intersection with the vertical axis z coincides with the first horizontal nullcline (16), that is, b=z∗ (see black line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,−0.855 is unstable, and the endemic E∗≃1.855,−0.855 is a singly degenerate equilibrium; both points are indicated by a black circle.Figure 9
Phase plane of system (9) for ρ1<R0<ρ1/σ, with R0≃2.924β′=0.5; the other rates are δ=0.071 , α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,−4.733 is unstable, and the endemic E∗≃5.733,−4.733 is a saddle point; both points are indicated by a black circle.Figure 10
Phase plane of system (9) for ρ1/σ<R0<1/σ, with R0≃4.678β′=0.8; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,1.972 is stable, and the endemic E∗≃−0.972,1.972 is an unstable node of two “straight-line” solutions; both points are indicated by a black circle.Figure 11
Phase plane of system (9), for the case where R0=1/σ≃5.882β′=1.005; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0 and the endemic E∗ coincide at the point 0,1, the latter being a doubly degenerate equilibrium; such points are indicated by a black circle.Figure 12
Phase plane of system (9) for 1/σ<R0<∞, with R0≃8.772β′=1.5; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,0.4581 is unstable, and the endemic E∗≃0.5419,0.4581 is a stable node of two “straight-line” solutions; both points are indicated by a black circle.Thus, if0≤R0<ρ1, say R0≃0.292β′=0.05, the disease-free equilibrium point E0≃0,−0.776 is stable and the endemic E∗≃1.776,−0.776 is an unstable node of two “straight-line” solutions (see Figure 7). When R0=ρ1, that is, R0≃0.585β=0.1, the disease-free equilibrium point E0≃0,−0.855 is unstable and the endemic E∗≃1.855,−0.855 is a singly degenerate equilibrium. In this case, the first vertical nullcline (15) is a line of zero slope whose intersection with the vertical axis z coincides with the first horizontal nullcline (16), that is, b=z∗; this line also consists of an infinite number of equilibrium points (see Figure 8). If ρ1<R0<ρ1/σ, say R0≃2.924β′=0.5, the disease-free equilibrium point E0≃0,−4.733 is unstable and the endemic E∗≃5.733,−4.733 is a saddle point (see Figure 9). When ρ1/σ<R0<1/σ, that is, R0≃4.678β′=0.8, the disease-free equilibrium point E0≃0,1.972 is stable and the endemic E∗≃−0.972,1.972 is an unstable node of two “straight-line” solutions (see Figure 10). If R0=1/σ, R0≃5.882β′≃1.005, the disease-free equilibrium point E0 and the endemic E∗ coincide at the point 0,1, the latter being a doubly degenerate equilibrium (see Figure 11). Finally, when 1/σ<R0<∞, say R0≃8.772β′=1.5, the disease-free equilibrium point E0≃0,0.4581 is unstable and the endemic E∗≃0.5419,0.4581 is a stable node of two “straight-line” solutions (see Figure 12).


## 4.1. Effects on the Trajectories ofyτ due to Variations of the α and σ Rates

We illustrate how changing the ratesα and σ affects the dynamics of the solutions of system (9), in particular those of the proportion of infected individuals yτ. For the different cases to be considered, we will indicate the values of R0 (since it varies inversely with α, as is pointed out by (7)) and their corresponding bifurcations 1/σ. Also, for each of the graphs obtained, taking into account (47), we will indicate the magnitude that reaches its peak, namely, ymaxα.Let us consider a first case whereα varies, β′=0.532, δ=0.071 (as in Table 1), and σ=0.1 will be fixed. We have chosen this last magnitude for the recovery rate, so that R0 remains below the value of the bifurcation 1/σ=10 that is constant. Figure 4 shows how the graphs of yτ change for three different values of α; specifically, for α=0.0001 (blue line), R0=7.48 and ymaxα=0.63; when α=0.1 (red line), R0=3.11 and ymaxα=0.51; and for α=0.17 (green line), R0=2.21 and ymaxα=0.41. In this figure, it can be seen that as α grows, the peak of the graphs decreases, they widen and shift to the right, and they are also less biased. Note that, in the first case α=0.0001, the rate α is so small that the effect produced is as if it were not there; in the absence of αα=0, R0 (for fixed β′ and δ) has its greatest value and the graph of y has its highest peak. Apparently, the inclusion in the model of the rate α produces a more adequate description of the evolution of infected individuals, since its absence causes overestimated values of R0 and the peak of its graph.Figure 4
Evolution ofyτ for the fixed parameters β′=0.532, δ=0.071, and σ=0.1. α=0.0001, blue line; α=0.1, red line; and α=0.17, green line.Now we consider a second case whereσ changes, β′=0.532, δ=0.071 (as in Table 1), and α=0.1 will be fixed. Note that, in this situation, for each σ, we will have different bifurcation values 1/σ so that if σ is very close to zero, the bifurcation point is much greater than one, whereas if it is very close to one, this point is much less than one. We choose the magnitudes of σ that allow to R0=3.11 (which does not change, since β′, δ, and α are constant) to be below or above the value of the bifurcation. Figure 5 illustrates how the temporal evolution of yτ changes for three different magnitudes of σ; particularly, when σ=0.1 (blue line), 1/σ=10 and ymaxα=0.51; when σ=0.3 (red line) 1/σ=3.33 and ymaxα=0.58; and for σ=0.4 (green line), 1/σ=2.5 and ymaxα=0.62. This figure shows how as σ grows, the graphs of yτ change: their peaks increase, widen, and shift slightly to the right. It should be noted that when σ=0.1 (blue line) and σ=0.3 (red line), R0 is below their corresponding bifurcation values and both trajectories tend to a disease-free equilibrium point E0 as τ tends to infinity, while when σ=0.4 (green line), R0 is above its bifurcation value and the trajectory tends to an endemic equilibrium point E∗ as τ tends to infinity.Figure 5
Evolution ofyτ for the fixed parameters β′=0.532, δ=0.071, and α=0.1. σ=0.1, blue line; σ=0.3, red line; and σ=0.4, green line.Finally, let us consider a third case, which certainly, could be considered as a combination of the previous two. Here,α varies, β′=0.532, and δ=0.071 (as in Table 1), but now, σ=0.17 will be a different fixed value. We have considered this latter magnitude in order that R0, as α increases, may exceed the value of the bifurcation 1/σ=5.88 that does not change. The different rates α selected make R0 take values below or above the bifurcation 1/σ. The effect that this has on the y curves is illustrated in Figure 6. Thus, in the figure, when α=0.0001 (blue line), R0=7.48, and ymaxα=0.65, the curve is above the value of the bifurcation and tends to an endemic equilibrium point as time tends to infinity. If α=0.1 (red line), R0=3.11 and ymaxα=0.53, whereas when α=0.17 (green line), R0=2.21 and ymaxα=0.43; in both cases, R0 is below the value of the bifurcation, so the curves tend to disease-free equilibrium points as time tends to infinity. Note that variations in the rate α can produce changes in the nature of the equilibrium points towards which the y curves tend: as α increases, the curves go from tending towards an endemic equilibrium point E∗ (when α=0.0001) to others free of disease E0 (if α=0.1 or α=0.17).Figure 6
Evolution ofyt for the fixed parameters β′=0.532, δ=0.071, and σ=0.17. α=0.0001, blue line; α=0.1, red line; and α=0.17, green line.


## 4.2. Phase Portraits before and during Forward Bifurcation

We present here some representative phase portraits of dynamic system (9), for different values of R0 between 0 and just over 1/σ, the bifurcation value, taking as reference the diagram in Figure 3. Figures 7–12 illustrate the different phase portraits that were calculated. For the realization of these graphs, the rates δ, α, and σ were considered as constants: the first with its value given in Table 1 and the other two with magnitudes equal to 0.1 and 0.17, respectively; the rate β′ was varied. The quantity R0 can increase its value, according to (7), if δ is kept constant and α decreases or β′ increases. It has been pointed out that the absence of the parameter α causes an overestimated magnitude of R0, so we will keep its value fixed and we will only vary β′ in an increasing way. For rates δ, α, and σ already mentioned, in accordance with (6), we have that ρ1≃0.5848,ρ2≃0.4152, while ρ1/σ≃3.440 and 1/σ≃5.882. For each value of R0, we will also indicate the values of their corresponding disease-free E0=0,z and endemic E∗=y∗,z∗ equilibrium points, expressions (13) and (14), respectively; since the former are infinite, we will only present one of them, namely, point E0=0,z∗.Figure 7
Phase plane of system (9) for 0≤R0<ρ1, with R0≃0.292β′=0.05; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,−0.776 is stable, and the endemic E∗≃1.776,−0.776 is an unstable node of two “straight-line” solutions; both points are indicated by a black circle.Figure 8
Phase plane of system (9), for the case where R0=ρ1≃0.585β′=0.1; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. In this case, the first vertical nullcline (15) is a line of zero slope whose intersection with the vertical axis z coincides with the first horizontal nullcline (16), that is, b=z∗ (see black line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,−0.855 is unstable, and the endemic E∗≃1.855,−0.855 is a singly degenerate equilibrium; both points are indicated by a black circle.Figure 9
Phase plane of system (9) for ρ1<R0<ρ1/σ, with R0≃2.924β′=0.5; the other rates are δ=0.071 , α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,−4.733 is unstable, and the endemic E∗≃5.733,−4.733 is a saddle point; both points are indicated by a black circle.Figure 10
Phase plane of system (9) for ρ1/σ<R0<1/σ, with R0≃4.678β′=0.8; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,1.972 is stable, and the endemic E∗≃−0.972,1.972 is an unstable node of two “straight-line” solutions; both points are indicated by a black circle.Figure 11
Phase plane of system (9), for the case where R0=1/σ≃5.882β′=1.005; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0 and the endemic E∗ coincide at the point 0,1, the latter being a doubly degenerate equilibrium; such points are indicated by a black circle.Figure 12
Phase plane of system (9) for 1/σ<R0<∞, with R0≃8.772β′=1.5; the other rates are δ=0.071, α=0.1, and σ=0.17. The vector field and some trajectories are indicated by lines and arrows in gray. Two nullclines are shown: one horizontal (blue line) and other vertical (red line); the horizontal and vertical remaining pair coincides with the z-axis. The disease-free equilibrium point E0≃0,0.4581 is unstable, and the endemic E∗≃0.5419,0.4581 is a stable node of two “straight-line” solutions; both points are indicated by a black circle.Thus, if0≤R0<ρ1, say R0≃0.292β′=0.05, the disease-free equilibrium point E0≃0,−0.776 is stable and the endemic E∗≃1.776,−0.776 is an unstable node of two “straight-line” solutions (see Figure 7). When R0=ρ1, that is, R0≃0.585β=0.1, the disease-free equilibrium point E0≃0,−0.855 is unstable and the endemic E∗≃1.855,−0.855 is a singly degenerate equilibrium. In this case, the first vertical nullcline (15) is a line of zero slope whose intersection with the vertical axis z coincides with the first horizontal nullcline (16), that is, b=z∗; this line also consists of an infinite number of equilibrium points (see Figure 8). If ρ1<R0<ρ1/σ, say R0≃2.924β′=0.5, the disease-free equilibrium point E0≃0,−4.733 is unstable and the endemic E∗≃5.733,−4.733 is a saddle point (see Figure 9). When ρ1/σ<R0<1/σ, that is, R0≃4.678β′=0.8, the disease-free equilibrium point E0≃0,1.972 is stable and the endemic E∗≃−0.972,1.972 is an unstable node of two “straight-line” solutions (see Figure 10). If R0=1/σ, R0≃5.882β′≃1.005, the disease-free equilibrium point E0 and the endemic E∗ coincide at the point 0,1, the latter being a doubly degenerate equilibrium (see Figure 11). Finally, when 1/σ<R0<∞, say R0≃8.772β′=1.5, the disease-free equilibrium point E0≃0,0.4581 is unstable and the endemic E∗≃0.5419,0.4581 is a stable node of two “straight-line” solutions (see Figure 12).


## 5. Discussion and Conclusions

This research work consisted of proposing a mathematical model based on a set of nonlinear ordinary differential equations to describe how a disease spreads in a population consisting of susceptible, infected, and recovered individuals. This model can be used, in principle, to study the type of epidemics generated by coronaviruses SARS-CoV, MERS-CoV, and SARS-CoV-2. In particular, two aspects of the diseases caused by these pathogens that we consider very relevant were analyzed with this approach: the first one is the possibility, given that some cases have been confirmed [19, 20], that recovered individuals have nonpermanent immunity to infection and may eventually be reinfected, and the second one is consider that infected individuals can die due to the disease and not by natural causes (as is the case of those who have some comorbidity or are older people), so they present a disease-induced death. For this purpose, it was proposed a modified SIRI epidemiological model, equation (1), that takes into account the general characteristics of these epidemics; specifically, the two aspects already mentioned: the former represented by a term with the reinfection rate σ and the latter by another with the disease-induced death rate α, respectively. In the following paragraphs, we review the most important results of our work obtained from the qualitative and numerical analysis of the aforementioned model; we discuss them and also present the conclusions that, in our opinion, emerge from them.We study a two-dimensional and dimensionless version of the aforementioned system, equation (9). We find an expression for its basic reproductive number R0, equation (7), which is in terms of two other rates: disease transmission β′ and recovery δ, and also of previous α. As a consequence of the qualitative analysis applied to system (9), it has been shown that it presents in the region of biological sense Λ an endemic equilibrium point E∗ that is globally asymptotically stable, as well as an infinite line of disease-free nonhyperbolic equilibrium points E0 which are stable and unstable; it was also shown that, in such region, it does not have periodic orbits. The analysis carried out reveals that the equilibrium points E∗ and E0 depend on R0; in particular, for the family of equilibrium points y∗ (the abscissa of E∗), equation (29), we find that R0 plays the role of a bifurcation parameter and that exhibits a forward bifurcation (see Figure 3), producing the bifurcation just when R0=1/σ. This analysis also shows that the model has two quantities that are invariant in time, that is, with two first integrals, equations (37) and (42). Using the equations of the model and these first integrals, an analytical expression can be found to quantify ymaxα, the maximum value that reaches the trajectory of the proportion of infected individuals, given by equation (47). From this quantity, an estimate can be made of the differences in peak heights between one of these graphs with α very small α≪1 and another with α=0 (see equation (48)).We find that the proposed model reduces, in the limit whenα⟶0, to the equilibrium points and the dynamics around them already known for the conventional SIRI model without vital dynamics [25, 26]. The presence of rate α in the model shows that R0, because it is in its denominator (see equation (7)), acquires a lower value than this quantity can have when α=0; it would seem then that the inclusion of α produces a more adequate value of R0, in the sense that it would not be overestimated. This effect can also be seen in the behavior of the solutions for the dimensionless proportion of infected individuals y; indeed, our numerical analysis shows that, as α grows, the peak of these trajectories decreases, they widen, and they shift to the right and are also less biased, as can be seen in Figure 4.On the other hand, the presence in the model of the rateσ produces the forward bifurcation. It should be noted that if we consider the limit case in which σ⟶0, with which the conventional SIR model without vital dynamics is recovered, the value of the bifurcation R0=1/σ tends to infinity; this is consistent because it is known that said model does not present any bifurcation [25, 26]. The numerical analysis of the solutions for the dimensionless proportion of infected individuals y shows that how has σ grows, their peaks increase, they widen, and they shift slightly to the right. It should be noted that when σ takes values such that R0 is below the bifurcation value R0=1/σ, the trajectories of y tend to a disease-free equilibrium point E0, while when σ takes values that make R0 be the abovesaid value, these trajectories tend to an endemic equilibrium point E∗, as can be seen in Figures 5 and 6.Finally, we point out that the numerical solutions of system (9) are also represented in the phase plane for different values of R0 before and after the bifurcation. The results obtained are shown in Figures 7 to 12.



---

*Source: 1570463-2021-11-26.xml*