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

Maestŕıa en Ciencias de la Complejidad, Universidad Autónoma de la Ciudad de México, Mexico City 03100, Mexico Facultad de Matemáticas, Universidad Autónoma de Guerrero, Av. Lázaro Cárdenas CU, Chilpancingo 39087, Guerrero, Mexico Sección de Estudios de Posgrado e Investigación, Escuela Superior de Medicina, Instituto Politécnico Nacional, Ciudad de México, Mexico División de Investigación, Hospital Juárez de México, Ciudad de México, Mexico


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. e 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   [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 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. e 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]. e 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]. e 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][22][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.
is 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]. e 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]. is 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. e 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][32][33]. ese 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. e 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. e rest of this paper is organized as follows. In Section 2, 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.

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 infected I, 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, S(t), I(t), and R(t) 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 S(t) + I(t) + R(t) � 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 a SIRI 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) 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 wellknown 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 so the rate of change of the total population N is not constant, it decreases proportionally with the number of infected individuals I.

Adimensionalization and Reduction of a Two-Dimensional
System. It is much more convenient to handle system (1), simplifying its writing by dimensionless. With this Discrete Dynamics in Nature and Society purpose, we propose to change the variables of number of individuals S, I, R and the temporal t, by the dimensionless quantities [25].
x � S N , e latter form is motivated by the fact that the parameters α, δ have dimensions of (units of time 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τ. erefore, system (1), written in terms of these new set of dimensionless variables, takes the form where we have also defined the dimensionless quantities 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. R 0 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 R 0 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 < R 0 ; besides, we have that ρ 1 + ρ 2 � 1. Note that the new variables x, y, z meet the condition 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 It may be shown easily that the region of biological sense is that is, all solutions starting in Λ remain there for all τ ≥ 0. Clearly, the set Λ is positively invariant with respect to system (9).

Basic Reproductive Number.
e procedure in which R 0 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 R 0 is to consider that initially, S(0) � S 0 , I(0) > 0 with I(0) ≪ S 0 and R(0) � 0. In this case, second equation (1) takes the form Hence, a necessary and sufficient condition for an initial increase in the number of infected (since where x 0 � S 0 /N. On the other hand, R 0 < 1 is the sufficient condition to have an initial decrease of infected and that the epidemic does not spread. us, equation (12) is the expression for R 0 that we are looking for and is in complete agreement, when x 0 � 1, with equation (7) found previously.

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 S I αI σβIR δI R βIS Figure 1: Flowchart of the variant of model SIRI considered. is 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. 4 Discrete Dynamics in Nature and Society with z ∈ R; that is, the straight line y � 0 with an infinite number of points that constitute the vertical coordinate axis of the yz plane. e second is the endemic equilibrium 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, σR 0 − ρ 1 > 0 and σR 0 − 1 ≥ 0. Just when σR 0 − 1 � 0, E * is located on the vetical z-axis; that is, it is part of the points E 0 .
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 which consist of two straight lines: and another that coincides with the vertical axis z. In addition, we have the nullclines dz/dτ � 0 (horizontal directions), determined by which are again two lines: one horizontal that cuts the vertical axis at point (0, ρ 2 /(σR 0 − ρ 1 )) and another vertical that also coincides with the z-axis. It should be noted that, as 0 < σ < 1 and R 0 > 0, the quantity (1 − σ)R 0 is always positive; consequently, the slope of the vertical nullcline is only determined by the sign of R 0 − ρ 1 and its ordinate to the origin by that of R 0 − 1. As we are interested in the points E * located in region Λ, it must happen that σR 0 − ρ 1 > 0, σR 0 − 1 ≥ 0, and σR 0 − ρ 1 > 0, as already mentioned previously; therefore, note that 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 Λ, σR 0 − ρ 1 > 0; therefore, the first horizontal nullcline (16) must be located in the upper half plane of plane yz. e intersections of different nullclines coincide with the equilibrium points. Note that, at y � 0, nullclines overlap in both vertical and horizontal directions; therefore y � 0 consists of an infinite line of equilibrium points, that is, the equilibria E 0 . 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. is figure also shows the different directions that the vector field has along all the nullclines. e directions (left or right) of the horizontal nullclines are determined by substituting z � z * , first equation (16), in the first equation (9), which gives rise to 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 σR 0 − 1 and σR 0 − ρ 1 ), the factor R 0 − ρ 1 , the values of y, and their order relation with y * . Since Λ is satisfied that σR 0 − ρ 1 > 0, σR 0 − 1 ≥ 0 (that is, y * ≥ 0) and R 0 − ρ 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. erefore, 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 * ).
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 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 σR 0 − 1 and σR 0 − ρ 1 ), the factor R 0 − ρ 1 , the values of y, and their relation with y * . Since in Λ, we have that (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. erefore, 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 σR 0 − 1, σR 0 − ρ 1 and Discrete Dynamics in Nature and Society R 0 − ρ 1 . It can also be shown that there are four possible cases. For y * > 0, if R 0 − ρ 1 and σR 0 − ρ 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 R 0 − ρ 1 and σR 0 − ρ 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 R 0 − ρ 1 and σR 0 − ρ 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 R 0 − ρ 1 and σR 0 − ρ 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 Λ.

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. is matrix overall is given as To perform the local stability analysis, we will consider separately the disease-free equilibria and the endemic equilibrium cases.

e Disease-Free Equilibria
e Jacobian matrix (19) evaluated in the disease-free equilibria E 0 takes the simple form whose eigenvalues are λ 1 � 0 and Note that these equilibrium points are nonhyperbolic. Due to this nature of the points E 0 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). is 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. ose results, rewritten in terms of the ideas of stability, can be written in the following way: for is stable, and when R 0 − ρ 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.
Theorem 1. Disease-free equilibrium points E 0 , which constitute the vertical axis z, are nonhyperbolic. ese points are Since when α � 0, our model is reduced to the conventional SIRI model, in this situation ρ 1 � 0, so from eorem 1, it follows that the nonhyperbolic equilibrium points E 0 are stable if 0 < R 0 < 1/σ and unstable if 1/σ < R 0 < ∞.
which allows writing the Jacobian matrix (19), evaluated in such a point, as e trace of (22) is given by and its determinant is whereas it can be shown that its discriminant, defined as which is a nonnegative quantity. According to the tracedeterminant plane method [36][37][38], the analysis of matrix (22) indicates that points E * , since Δ ≥ 0, can be saddle points (if Det J < 0 ) and degenerate equilibrium (if Det J � 0) or nodes of two "straight-line" solutions (if Det J > 0), the latter pair being stable or unstable if Tr J < 0 or Tr J > 0, respectively. Furthermore, the eigenvalues of (22) are given by which can be explicitly written as Note that eigenvalues (27) and (28) are functions of R 0 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.

Theorem 2.
e endemic equilibrium points E * are in general hyperbolic. Such points are (i) nodes (of two "straight-line" solutions) locally as- It is important to note that, in eorem 2, subinterval 1/σ < R 0 < ∞ (that of condition i)) is the only one that guarantees that points E * are in region Λ. Again, since when Discrete Dynamics in Nature and Society α � 0, our model is reduced to the conventional SIRI model, in this situation, ρ 1 � 0, so for this case, from eorem 2, it follows that the equilibrium points E * are nodes (of two "straight-line" solutions) locally asymptotically stable if 1/σ < R 0 < ∞ and unstable if 0 < R 0 < 1/σ.

Presence of a Forward Bifurcation.
e analysis that allowed to formulate the two previous theorems shows that the coordinates of both points E 0 and E * are a function of the reproductive number R 0 . 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 R 0 and to visualize such dependence in a diagram R 0 y * ; this is important because R 0 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 R 0 axis and another determined by of endemic equilibrium points. By inspection, it can be seen that equation (29) has a discontinuity at R 0 � ρ 1 /σ, which gives rise to one left branch that cuts the vertical axis at y * � 1/ρ 1 when R 0 � 0 and moves upward indefinitely as R 0 approaches to R 0 � ρ 1 /σ, as well as another right branch for values greater than R 0 � ρ 1 /σ that comes from below cutting to the horizontal axis at R 0 � 1/σ and asymptotically approaches to the horizontal line y * � 1 as R 0 tends to infinity, see Figure 3. Figure 3, 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 R 0 according to their positions in segment L defined above.
us, 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 Λ. is 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 R 0 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.

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). ese 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.

Nonexistence of Periodic Orbits.
e nonexistence of periodic orbits in the two-dimensional system (9) implies, of course, that neither system (5) has this type of trajectories. e first statement is shown in eorem 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 We propose the Dulac function In this way, it can be easily verified that since as mentioned in Section 3, in Λ, it is true that σR 0 − ρ 1 > 0 and as a consequence, in this region, R 0 − ρ 1 > 0. erefore, according to the Dulac criterion [37,38], there are no periodic orbits in the interior of Λ.

Global Stability of the Endemic Equilibrium Point.
We will establish the conditions that the endemic equilibrium point E * must meet to guarantee its global stability.
is will be carried out by following Lyapunov's second method or Lyapunov stability theorem [37,38], as is shown in eorem 4.

Theorem 4.
e unique endemic equilibrium, E * , is globally asymptotically stable in the interior of Λ.
Proof. We consider the common quadratic Lyapunov function y* 1/ρ 1 R 0 ρ 1 /σ 1/σ ρ 1 1 0 Λ Figure 3: Bifurcation diagram of system (9). Locally asymptotically stable and unstable points have been indicated, respectively, with continuous and discontinuous lines. e region of biological sense Λ is indicated with gray color. 8 Discrete Dynamics in Nature and Society which is continuous and a positive definite function, V(y, z) > 0, in a region that contains E * � (y * , z * ) that is in turn located in Λ. It can be verified by inspection that V(y, z) takes the value V(y, z) � 0 at E * so that the global minimum value of V(y, z) is located at that point. Calculating the total derivative of (32) and using the two-dimensional system (9), we obtain Taking into account the relations (21), the previous expression (34) takes the form this is true since, as already mentioned, R 0 − ρ 1 � (1 − σ)R 0 + σR 0 − ρ 1 > 0 in such a region so that R 0 − ρ 1 > (1 − σ)R 0 , and as we finally obtain R 0 − ρ 1 − (1/4)(1 − σ)R 0 > 0. us, since also (1 − σ)R 0 > 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). erefore, 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 twodimensional 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 so that (8) is in itself worth the redundancy, a relation among the variables x, y, and z that remains constant. erefore, this is a first integral, as shown below.

Theorem 5. Dynamic system (5) has a first integral, given by
Proof. Calculating the total derivative of (37) and using the three equations (5), we have since (8) is true and also ρ 1 + ρ 2 � 1. e 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 is expression can be solved, in principle, by the method of separating variables, giving rise to (40) where C is a integration constant. It can be shown that, after some algebraic steps, it gives rise to 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.
Discrete Dynamics in Nature and Society Theorem 6. Dynamic system (5) has another first integral, given by Proof. Calculating the total derivative of (39) and using the first and third equations (5), we have e 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.

Calculation of the Maximum of y(τ).
e 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 proportion y(τ) reaches. To do this, we first consider that we have the initial conditions x(0) � x 0 and z(0) � z 0 . us, it can be shown that, from (41), we obtain ) .
(44) erefore, using (8) and (44), we have the following expression for the variable y written only in terms of x: Equating the derivative of y(τ) 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 and also that it is a maximum. erefore, 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 y max α ; that is, 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 eorem 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, y max α − y max α�0 . If we consider α small (α ≪ 1), we can approximate definitions (6) by ρ 1 ≃α/δ and ρ 2 ≃1. erefore, for this situation, it can be shown that the said difference is given by .

(48)
A slight simplification of expressions (47) and (48) can be obtained, if we consider z 0 � 0 in them; that is, there are no infected individuals at the beginning of the epidemic.

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.
e 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.071 day − 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 diseaseinduced 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 (R 0 ) 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 R 0 ≃3.11 (as in [42]). us, according to the parameters indicated, from relation (7), it follows that β ′ ≃0.532 day − 1 . Finally, it has already been mentioned that the bifurcation occurs when R 0 � 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. e magnitudes of the indicated rates δ, α, β ′ , and σ are shown in Table 1.
Since the purpose here is to show the effects produced in the population of infected individuals y(τ) 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 of y(τ) 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 R 0 (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, y max α .
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 R 0 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), R 0 � 7.48 and y max α � 0.63; when α � 0.1 (red line), R 0 � 3.11 and y max α � 0.51; and for α � 0.17 (green line), R 0 � 2.21 and y max α � 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), R 0 (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 R 0 and the peak of its graph. 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 R 0 � 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 y max α � 0.51; when σ � 0.3 (red line) 1/σ � 3.33 and y max α � 0.58; and for σ � 0.4 (green line), 1/σ � 2.5 and y max α � 0.62. is 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), R 0 is below their corresponding bifurcation values and both trajectories tend to a disease-free equilibrium point E 0 as τ tends to infinity, while when σ � 0.4 (green line), R 0 is above its bifurcation value and the trajectory tends to an endemic equilibrium point E * as τ tends to infinity.

Discrete Dynamics in Nature and Society
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 R 0 , as α increases, may exceed the value of the bifurcation 1/σ � 5.88 that does not change. e different rates α selected make R 0 take values below or above the bifurcation 1/σ. e effect that this has on the y curves is illustrated in Figure 6. us, in the figure, when α � 0.0001 (blue line), R 0 � 7.48, and y max α � 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), R 0 � 3.11 and y max α � 0.53, whereas when α � 0.17 (green line), R 0 � 2.21 and y max α � 0.43; in both cases, R 0 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 E 0 (if α � 0.1 or α � 0.17).

Phase Portraits before and during Forward Bifurcation.
We present here some representative phase portraits of dynamic system (9), for different values of R 0 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. e quantity R 0 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 R 0 , 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 R 0 , we will also indicate the values of their corresponding diseasefree E 0 � (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 E 0 � (0, z * ).
e disease-free equilibrium point E 0 ≃(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 11: Phase plane of system (9), for the case where R 0 � 1/σ≃5.882(β′ � 1.005); the other rates are δ � 0.071, α � 0.1, and σ � 0.17. e 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. e disease-free equilibrium point E 0 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.

Discussion and Conclusions
is 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. is 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 R 0 , 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 E 0 which are stable and unstable; it was also shown that, in such region, it does not have periodic orbits. e analysis carried out reveals that the equilibrium points E * and E 0 depend on R 0 ; in particular, for the family of equilibrium points y * (the abscissa of E * ), equation (29), we find that R 0 plays the role of a bifurcation parameter and that exhibits a forward bifurcation (see Figure 3), producing the bifurcation just when R 0 � 1/σ. is 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 y max α , 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]. e presence of rate α in the model shows that R 0 , 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 R 0 , in the sense that it would not be overestimated. is 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 R 0 � 1/σ tends to infinity; this is consistent because it is known that said model does not present any bifurcation [25,26]. e 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 R 0 is below the bifurcation value R 0 � 1/σ, the trajectories of y tend to a disease-free equilibrium point E 0 , while when σ takes values that make R 0 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 R 0 before and after the bifurcation. e results obtained are shown in Figures 7 to 12.  Figure 12: Phase plane of system (9) for 1/σ < R 0 < ∞, with R 0 ≃8.772(β′ � 1.5); the other rates are δ � 0.071, α � 0.1, and σ � 0.17. e 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. e disease-free equilibrium point E 0 ≃(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.

Data Availability
e data used to support the findings of this study are included within the article.