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 number and find that the system exhibits a forward bifurcation: disease-free equilibria are stable when and unstable when , while the endemic equilibrium consist of an asymptotically stable upper branch that appears from , 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.
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 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.
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 infected , who, by interacting with the susceptible , get infected, causing the eventual spread of the disease throughout the population, and recovered , 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, , , and denote the numbers of individuals in the groups or compartments , , and , at time , respectively; the total population is considered to be of fixed size , so that . 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:
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, is the mean recovery period, and is the reduction in susceptibility due to previous infection that takes values between zero and one . It should also be noted that system (1) is in effect a modified SIRI model: the typical term represents the reinfection; that is, the recovered people, after some time, become infected again; the change comes from the disease-induced death term , 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 haveso the rate of change of the total population is not constant, it decreases proportionally with the number of infected individuals .
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 and the temporal , by the dimensionless quantities [25].
The latter form is motivated by the fact that the parameters have dimensions of (units of time)−1. From (3), we have that and ; note that the point above the quantity in question indicates its derivative with respect to time . On the other hand, from (4), using the chain rule, we obtain that . Therefore, system (1), written in terms of these new set of dimensionless variables, takes the formwhere we have also defined the dimensionless quantitieswhich 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. 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 and consequently estimate ). In (7), we have defined the quantity . It should be noted that since the parameters involved are positive and is nonnegative, then , and ; besides, we have that . Note that the new variables 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 , in the two-dimensional system
It may be shown easily that the region of biological sense isthat is, all solutions starting in remain there for all . Clearly, the set is positively invariant with respect to system (9).
2.2. Basic Reproductive Number
The procedure in which 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 is to consider that initially, with and . 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 ) is , or equivalently , wherewhere . On the other hand, 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 that we are looking for and is in complete agreement, when , 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 that satisfy this equality. As a result of this procedure, we obtain two equilibrium points. One is disease-free equilibria of the formwith ; that is, the straight line with an infinite number of points that constitute the vertical coordinate axis of the plane. The second is the endemic equilibrium
So, the dynamics develop in region , and equilibrium points must be in the first quadrant of plane , that is, and ; in other words, and . Just when , is located on the vetical -axis; that is, it is part of the points .
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 (vertical directions), given by the functionswhich consist of two straight lines: one of slope and ordered to origin and another that coincides with the vertical axis . In addition, we have the nullclines (horizontal directions), determined bywhich are again two lines: one horizontal that cuts the vertical axis at point and another vertical that also coincides with the -axis. It should be noted that, as and , the quantity is always positive; consequently, the slope of the vertical nullcline is only determined by the sign of and its ordinate to the origin by that of . As we are interested in the points located in region , it must happen that , and , as already mentioned previously; therefore, note that and . 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 to be in , ; therefore, the first horizontal nullcline (16) must be located in the upper half plane of plane .
The intersections of different nullclines coincide with the equilibrium points. Note that, at , nullclines overlap in both vertical and horizontal directions; therefore consists of an infinite line of equilibrium points, that is, the equilibria . On the other hand, another intersection between the other different nullclines gives rise to point . 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 (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.
(a)
(b)
The directions (left or right) of the horizontal nullclines are determined by substituting , 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 (which depends in turn on the signs of and ), the factor , the values of , and their order relation with . Since is satisfied that (that is, ) and ( is on the nonnegative -axis), the horizontal nullcline cuts the nonnegative vertical axis at point and the nullcline of vertical directions, first equation (15), at , which is an equilibrium point. Therefore, according to (16), if , then ; that is, in this interval, the vector field is directed to the left if we are to the right of ; when , then so that the vector field points to the right if we are to the left of and to the right of , while if , then ; that is, the direction of the field is to the left if we are to the left of . This result indicates that, on horizontal nullcline , the vector field leaves point and enters .
The abovementioned analysis can be repeated for all combinations of the quantities (which determine the sign of ) and . It can be shown that there are four possible cases. For , if , then vector field on the horizontal nullcline leaves and enters , and when , it enters and leaves ; Figure 2(a) illustrates the first case. Also, for , if , then field enters and leaves , and when , it leaves and enters ; 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
In order to determine the positive or negative sign of (18), we need to know the signs of (that depends on the signs of and ), the factor , the values of , and their relation with . Since in , we have that and , the first vertical nullcline (16) cuts horizontal nullcline at the equilibrium point and the positive vertical axis at the point and has a negative slope. Therefore, according to (18), if , then ; that is, in this interval, the direction of the vector field is up, while if , then so that, in this interval, directions are down, and when , 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 and . It can also be shown that there are four possible cases. For , if and are of opposite signs, then the direction of vector field on this vertical nullcline is down in the interval , is up in , and is down in whereas if and are of the same signs, then the vector field is up in , down in , and up in . Figure 2(a) illustrates the second case. Besides, for , if and are of opposite signs, then the vector field is down in , up in , and down in , whereas if and are of the same signs, then the vector field is up in , down in , and up in 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
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,
The Jacobian matrix (19) evaluated in the disease-free equilibria takes the simple formwhose eigenvalues are and . Note that these equilibrium points are nonhyperbolic. Due to this nature of the points 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 around their point of intersection with the vertical axis . Those results, rewritten in terms of the ideas of stability, can be written in the following way: for , if , then point is unstable, and when , it is stable, whereas for , if , then point is stable, and when , it is unstable, see Figure 2. Since we have an infinite number of points along the vertical axis , the stability along a horizontal line that cuts one of them, say , must be identical to the stability found at its intersection with the horizontal nullcline . It should be noted that the analysis based on the method of the nullclines reveals that points are stable if and unstable if , see again Figure 2.
Since and , we can establish the relations and ; besides that, we know that ; therefore, the order relation between these quantities is , which will be called segment . On the other hand, it has just been mentioned that the stability of the nonhyperbolic points depends on the signs of and , or in other words, it depends of the location that has in .
Thus, when (which determine that ) and , according to segment , the intersection of these three subintervals is . Since the previous stability result establishes that, in this situation, points are unstable, we conclude that the equilibrium points are unstable when . On other hand, when and , according to , intersection is ; since points are stable here, we have that points are unstable when . By means of an analysis similar to the previous one, it can be shown that points are unstable in and stable in . These obtained results can be summarized as follows.
Theorem 1. Disease-free equilibrium points , which constitute the vertical axis , are nonhyperbolic. These points are(i)unstable, if and (ii)stable, if and
Since when , our model is reduced to the conventional SIRI model, in this situation , so from Theorem 1, it follows that the nonhyperbolic equilibrium points are stable if and unstable if .
3.2.2. The Endemic Equilibrium,
By substituting the endemic equilibrium in system (9), the identities are obtained.which allows writing the Jacobian matrix (19), evaluated in such a point, as
The trace of (22) is given byand its determinant iswhereas it can be shown that its discriminant, defined as , iswhich is a nonnegative quantity. According to the trace-determinant plane method [36–38], the analysis of matrix (22) indicates that points , since , can be saddle points (if ) and degenerate equilibrium (if ) or nodes of two “straight-line” solutions (if ), the latter pair being stable or unstable if or , respectively. Furthermore, the eigenvalues of (22) are given bywhich can be explicitly written as
Note that eigenvalues (27) and (28) are functions of and their signs will depend on the location of this quantity in segment already indicated above. In general, they are different from zero, so the equilibrium points are hyperbolic.
Next, we will analyze the signs that and could have in . According to (27) and (28), and , if , and ; the intersection of these three subintervals occurs at , so in this region, the equilibrium points correspond to nodes of two “straight-line” solutions that are locally asymptotically stable. In this case, in which both eigenvalues are negative, points lie within region ; however, for other combinations of signs of the eigenvalues and , points lie outside of . Thus, it can be shown that and when and so that, in these regions, points correspond to unstable nodes of two “straight lines.” Additionally, and , if , so there points are saddle points; the case and does not exist because the intersection of the corresponding subintervals gives rise to the empty set. Moreover, and , when , so here, is a singly degenerate equilibrium. Lastly, , when , so 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 points are in general hyperbolic. Such points are(i)nodes (of two “straight-line” solutions) locally asymptotically stable, if (ii)a doubly degenerate equilibrium, if (in which case it is nonhyperbolic)(iii)unstable nodes (of two “straight-line” solutions), if and (iv)saddle points, if (v)a singly degenerate equilibrium, if (in which case it is nonhyperbolic)(vi)unstable nodes (of two “straight-line” solutions), if
It is important to note that, in Theorem 2, subinterval (that of condition )) is the only one that guarantees that points are in region . Again, since when , our model is reduced to the conventional SIRI model, in this situation, , so for this case, from Theorem 2, it follows that the equilibrium points are nodes (of two “straight-line” solutions) locally asymptotically stable if and unstable if .
3.3. Presence of a Forward Bifurcation
The analysis that allowed to formulate the two previous theorems shows that the coordinates of both points and are a function of the reproductive number . In particular, it is of great interest to know how the family of equilibrium points (which is the abscissa of points ) varies with and to visualize such dependence in a diagram ; this is important because could play the role of a bifurcation parameter. For our study problem, the families consist of two curves: one horizontal line of disease-free equilibrium points located on the axis and another determined byof endemic equilibrium points. By inspection, it can be seen that equation (29) has a discontinuity at , which gives rise to one left branch that cuts the vertical axis at when and moves upward indefinitely as approaches to , as well as another right branch for values greater than that comes from below cutting to the horizontal axis at and asymptotically approaches to the horizontal line as 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 according to their positions in segment defined above. Thus, for example, this figure shows how these families exhibit an abrupt change at the point ; 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 ; that is, it consists of a horizontal branch that coincides with the 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 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) asWe propose the Dulac functionIn this way, it can be easily verified thatsince as mentioned in Section 3, in , it is true that and as a consequence, in this region, . 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 point 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, , is globally asymptotically stable in the interior of .
Proof. We consider the common quadratic Lyapunov functionwhich is continuous and a positive definite function, , in a region that contains that is in turn located in . It can be verified by inspection that takes the value at so that the global minimum value of is located at that point. Calculating the total derivative of (32) and using the two-dimensional system (9), we obtainTaking into account the relations (21), the previous expression (34) takes the formNote in (35) that in region ; this is true since, as already mentioned, in such a region so that , and as , we finally obtain . Thus, since also , the negative inequality in (35) holds for all and , that is, for every point located in (except the positive coordinate axes and the origin). Therefore, according to Lyapunov stability theorem, the endemic equilibrium 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 obtainso that (8) is in itself worth the redundancy, a relation among the variables , , and that remains constant. Therefore, 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 havesince (8) is true and also . The proof is complete.
On the other hand, note that the first and third equations (5) are only coupled by the variable ; therefore, by forming the quotient between them, we obtain
This expression can be solved, in principle, by the method of separating variables, giving rise towhere is a integration constant. It can be shown that, after some algebraic steps, it gives rise to being a constant quantity. Consequently, we have obtained a quantity, dependent on the variables and , 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
Proof. Calculating the total derivative of (39) and using the first and third equations (5), we haveThe proof is complete.
It should be noted that expression (42) is reduced, when , 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 of
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 proportion reaches. To do this, we first consider that we have the initial conditions and . Thus, it can be shown that, from (41), we obtain
Therefore, using (8) and (44), we have the following expression for the variable written only in terms of :
Equating the derivative of 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 byand also that it is a maximum. Therefore, substitution of (46) in (45) allows us to find the maximum of the infectious proportion at any value of and in the presence of the parameter , which we will denote as ; that is,
It should be noted that equation (47) reduces completely, when the limit case is taken (that is, and ), 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 in the presence and absence of the parameter , that is, . If we consider small , we can approximate definitions (6) by and . Therefore, for this situation, it can be shown that the said difference is given by
A slight simplification of expressions (47) and (48) can be obtained, if we consider 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 of an infected person so that . 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 ; however, in our calculations, we will vary this rate and consider that . 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 to China reported a pooled estimation of (95% confidence interval: ). 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 (as in [42]). Thus, according to the parameters indicated, from relation (7), it follows that . Finally, it has already been mentioned that the bifurcation occurs when , so a first estimate of the relapse rate, although it is for this case, is ; precisely in this model, , 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.
Since the purpose here is to show the effects produced in the population of infected individuals 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 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 . For the different cases to be considered, we will indicate the values of (since it varies inversely with , as is pointed out by (7)) and their corresponding bifurcations . Also, for each of the graphs obtained, taking into account (47), we will indicate the magnitude that reaches its peak, namely, .
Let us consider a first case where varies, , (as in Table 1), and will be fixed. We have chosen this last magnitude for the recovery rate, so that remains below the value of the bifurcation that is constant. Figure 4 shows how the graphs of change for three different values of ; specifically, for (blue line), and ; when (red line), and ; and for (green line), and . 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 , the rate is so small that the effect produced is as if it were not there; in the absence of , (for fixed and ) has its greatest value and the graph of 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 and the peak of its graph.
Now we consider a second case where changes, , (as in Table 1), and will be fixed. Note that, in this situation, for each , we will have different bifurcation values 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 (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 changes for three different magnitudes of ; particularly, when (blue line), and ; when (red line) and ; and for (green line), and . This figure shows how as grows, the graphs of change: their peaks increase, widen, and shift slightly to the right. It should be noted that when (blue line) and (red line), is below their corresponding bifurcation values and both trajectories tend to a disease-free equilibrium point as tends to infinity, while when (green line), is above its bifurcation value and the trajectory tends to an endemic equilibrium point as tends to infinity.
Finally, let us consider a third case, which certainly, could be considered as a combination of the previous two. Here, varies, , and (as in Table 1), but now, will be a different fixed value. We have considered this latter magnitude in order that , as increases, may exceed the value of the bifurcation that does not change. The different rates selected make take values below or above the bifurcation . The effect that this has on the curves is illustrated in Figure 6. Thus, in the figure, when (blue line), , and , the curve is above the value of the bifurcation and tends to an endemic equilibrium point as time tends to infinity. If (red line), and , whereas when (green line), and ; in both cases, 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 curves tend: as increases, the curves go from tending towards an endemic equilibrium point (when ) to others free of disease (if or ).
4.2. Phase Portraits before and during Forward Bifurcation
We present here some representative phase portraits of dynamic system (9), for different values of between 0 and just over , 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 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 , 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 , while and . For each value of , we will also indicate the values of their corresponding disease-free and endemic equilibrium points, expressions (13) and (14), respectively; since the former are infinite, we will only present one of them, namely, point .
Thus, if , say , the disease-free equilibrium point is stable and the endemic is an unstable node of two “straight-line” solutions (see Figure 7). When , that is, , the disease-free equilibrium point is unstable and the endemic 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 coincides with the first horizontal nullcline (16), that is, ; this line also consists of an infinite number of equilibrium points (see Figure 8). If , say , the disease-free equilibrium point is unstable and the endemic is a saddle point (see Figure 9). When , that is, , the disease-free equilibrium point is stable and the endemic is an unstable node of two “straight-line” solutions (see Figure 10). If , , the disease-free equilibrium point and the endemic coincide at the point , the latter being a doubly degenerate equilibrium (see Figure 11). Finally, when , say , the disease-free equilibrium point is unstable and the endemic 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 , 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 that is globally asymptotically stable, as well as an infinite line of disease-free nonhyperbolic equilibrium points 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 and depend on ; in particular, for the family of equilibrium points (the abscissa of ), equation (29), we find that plays the role of a bifurcation parameter and that exhibits a forward bifurcation (see Figure 3), producing the bifurcation just when . 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 , 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 and another with (see equation (48)).
We find that the proposed model reduces, in the limit when , 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 , because it is in its denominator (see equation (7)), acquires a lower value than this quantity can have when ; it would seem then that the inclusion of produces a more adequate value of , 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 ; 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 , with which the conventional SIR model without vital dynamics is recovered, the value of the bifurcation 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 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 is below the bifurcation value , the trajectories of tend to a disease-free equilibrium point , while when takes values that make be the abovesaid value, these trajectories tend to an endemic equilibrium point , 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 before and after the bifurcation. The results obtained are shown in Figures 7 to 12.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this manuscript.