Epidemic spreading

We present an analysis of six deterministic models for epidemic spreading. The evolution of the number of individuals of each class is given by ordinary differential equations of the first order in time, which are set up by using the laws of mass action providing the rates of the several processes that define each model. The epidemic spreading is characterized by the frequency of new cases, which is the number of individuals that are becoming infected per unit time. It is also characterized by the basic reproduction number, which we show to be related to the largest eigenvalue of the stability matrix associated with the disease-free solution of the evolution equations. We also emphasize the analogy between the outbreak of an epidemic with a critical phase transition. When the density of the population reaches a critical value the spreading sets in, a result that was advanced by Kermack and McKendrick in their study of a model in which the recovered individuals acquire permanent immunization, which is one of the models analyzed here.


I. INTRODUCTION
The control and possibly the eradication of an infectious disease is not fully successful without the knowledge of the mechanisms that determine its spreading in a population. A major contribution in this direction is provided by the theoretical study of epidemic spreading through the use of deterministic and stochastic models. The theoretical study of the epidemic spreading [1][2][3][4][5][6] is not new but it did not properly start before the physical basis for the causes of the infectious disease was established in the second half of the nineteenth century [1]. By the end of the nineteenth century the mechanism of epidemic spreading, revealed by bacteriology, allowed the development of the first epidemic models [1].
Hammer in 1906, carried out an analysis of the measles epidemic considering that the infection spreads from person to person and depends on the number of the susceptible [1]. In his studies of the relation between the number of mosquitoes and the incidence of malaria, Ross starting from 1908 formulated models for the transmission of the infectious diseases [1]. In his book on the prevention of malaria of 1911 [7], he employed ordinary differential equations to describe the evolution of the number of affected and unaffected individuals by a disease, and wrote the equations in terms of rates of several types, pointing out the major role of the infection rate.
A more general theory than that of Ross, but employing similar ideas, was developed by Kermack and McKendrick in 1927 [8]. They proposed a model for an epidemic consisting of recovered individuals having permanent immunization. They were able to solve the differential equations that governs the time evolution of the number of susceptible, infected, and recovered, arriving at a remarkable theorem concerning the spreading threshold of an epidemic [1,8]. If the density of the susceptible is smaller than a certain value, the epidemic does not outbreak. In other terms, if the infection rate is smaller than a critical value the disease does not spread.
In 1929, Soper developed deterministic models for measles by using difference equations [1], and by assum-ing that the operations of transmission are analogous to the mass law action of chemistry [9]. The equations developed by Kermack and McKendrick were also in accordance with this law which would become one of the most important concepts in theoretical epidemiology [2]. The deterministic models, employing ordinary differential equations, were carried further mainly after around 1945 [1]. In 1952, Macdonald introduced a concept, which he called basic reproduction rate, concerning the threshold of spreading, with reference to the threshold theorem of Kermack and McKendrick [10,11]. However, it was not before 1975 when the concept was reintroduced by Hethcote and by Dietz, that it became widely used [12].
Around 1945, the stochastic approach for epidemic spreading, which had been originated earlier, received further development by Bartlett [1]. He transformed the Kermack and McKendrick epidemic model into a stochastic model by treating the numbers of susceptible and infective individuals as stochastic variables [13]. The corresponding master equation in these two variables was obtained by Bailey [1,14] and the stochastic version of the Kermack and McKendrick threshold theorem was advanced by Whittle [15]. The model included two processes: the recovery of infected individual, who acquires permanent immunization, and the infection of a susceptible by the contact with an infective individual. These stochastic models can be understood as discrete generic random walk in a space whose axes are the numbers of various classes of individuals [16,17].
A distinct type of stochastic models for epidemic spreading takes into account the spatial structure in which real infectious diseases spread in a population [18][19][20][21][22][23][24][25][26][27][28][29][30]. The spatial stochastic models are formulated in terms of several stochastic variables, one for each individual, and which take values corresponding to the class an individual belongs in with respect to the disease [17]. Usually, they are defined on a lattice that represents the spatial structure.
The spatial stochastic model allows us to derive the stochastic model of the generic random walk mentioned above by a suitable reduction of stochastic variables. This reduction can be made at the first level, represented by the numbers of individuals of each class, or at the second level, which takes into account the pairs of individuals [27,31]. It is also possible from the spatial stochastic models to reach the deterministic evolution equations at the first level, which are the usual evolution equations such as the ones studied here, or at the second level where the numbers of pairs of individuals are taken into account [22,28,32].
Here we only consider the deterministic approach to the epidemic spreading, which still is a relevant approach to the study of epidemic spreading of various infectious diseases [33][34][35]. We analyze models that are appropriate for the direct transmission of the disease, that is, when the transmission occurs from person to person, which is the case of most common infectious diseases such as measles, chickenpox, mumps, rubella, smallpox, common cold, and influenza, as well as the present covid-19. These models include the susceptible-infected-susceptible (SIS) [7], the susceptible-infected-recovered (SIR) [8], and the susceptible-exposed-infected-recovered (SEIR) [36]. The models that we considered here are those such that the total number of individuals is constant in time. We do not consider demographic processes, such as birth, death, and migration.
We present the main features that are used to characterize the spreading such as the reproduction number and the frequency of new cases, which is called the epidemic curve when it is plotted as a function of time. Our presentation emphasizes the analogy with the thermodynamic phase transition, the onset of spreading being understood as the critical point of a continuous thermodynamic phase transition, when the epidemic comes to an end, in which case the frequency of new cases vanishes, the size of the epidemic can be measured by the area under the epidemic curve, and can be understood the order parameter.
The spatial stochastic version of the SIS model was proposed by Harris [18], who named it the contact process. It is widely studied not only because it describes an epidemic spreading but also because of its critical behavior [20] which is distinct from models describing system in thermodynamic equilibrium such as the Ising model [17,37,38]. The SIS model describes the evolution of an infectious disease that becomes endemic. The simplest model that describes the evolution of a infectious disease that becomes extinct is the SIR model. Its spatial version was formulated by Grassberger [19], who called it general epidemic process, and its critical behavior is distinct from the Ising model and also different from the SIS model [17,37,38].
It is worth pointing out that the approach to epidemic spreading, be it deterministic or stochastic, is extended to a more general context of the population dynamics [17,39,40], as for instance in ecological problems such as the predator-prey dynamics [17,22].

II. EVOLUTION EQUATIONS
The spreading of an epidemic is understood as the time evolution of a population consisting of several classes of individuals pertaining to a certain community. The two main classes of individuals are the susceptible and the infective but other classes may also be present. As the epidemic evolves in time, an individual belonging in one class might change to another class through a transforming process. We remark that the members of the infective class are individuals that infect others. An individual that is infected by a disease but does not infect others belong in a distinct class, the exposed class.
The framework that is used here to analyze the population dynamics is borrowed from the theory of chemical kinetics. According to this theory, molecules of various species inside a vessel are subject to chemical reactions that transform molecules of one species into molecules of another species. The several classes of individuals are analogous to the chemical species, and the processes that transform an individual of one class to an individual of another class are analogous to chemical reactions.
As an example of the analogy, let us consider the process where an infective (I) individual recovers from the disease becoming a recovered (R) individual. This process is represented by which is understood as a reaction that transforms one I into one R, or the annihilation of one I and the simultaneous creation of one R. This type of reaction is called spontaneous.
Another example is the process of infection of a susceptible (S) individual who become infected by the contact with an infective (I) individual It is represented by where the symbol I over the arrow means that the reaction that transforms one S into one I needs the presence of one I. This type of reaction is called catalytic or more precisely auto-catalytic and I is the catalyst. After introducing the analogy between chemical species and chemical reactions, on one side, and classes of individuals and transforming process, in the other, we introduce a second analogy. The time evolution of the numbers of individuals of each class, determined by the transforming processes, are understood as analogous to the chemical kinetic equations. These are ordinary differential equations of the first order in time involving the concentrations or the fractions of the various chemical species and are set up by the use of the law of mass action.
The primary question one wants to answer is how many individuals are there in each class. Thus we should seek for equations in the number of individuals in each class. However, it is more convenient to express the equations in terms of the density, which is the number of individuals per unit area, or in terms of the fraction of individuals in each class. To set up the evolution equations in either one of these two types of variables, we will employ the law of mass action, presented and explained in the Appendix. Here, we choose to write down the equations in terms of the fractions of each class of individuals.
As an example of the application of the law of mass action, consider the infection process represented by equation (2), and let us denote by x the fraction of the susceptible and by y the fraction of the infective. The contribution of this process to the evolution equation of both x and y, in accordance with the law of mass action, is bxy, where b is the infection rate constant, that is, where the dots indicate the contribution coming from other processes. Notice the presence of a minus sign in the first equation, because the number of the susceptible decreases in the infection process. The evolution equations describing an epidemic consist of two or more equations that give the time variation of the fractions of the several classes of individuals. We assume that the total number of individuals is constant, so that the sum of all fractions equals one. They are solved with an initial condition that reflects what occurs in a real epidemic. At the beginning, all individuals are susceptible except a very small fraction of them that are infective.
Our presentation uses frequently the term rate, as in infection rate. It should not be confused with the variation in time, which is a derivative. One should also make a distinction between a rate and a rate constant which appears as a prefactor of a rate.

A. Epidemic curve
A relevant characterization of the evolution of an epidemic is the frequency of the number of individuals that are becoming infected, or the number of newly infected individuals per unit time, or simply the frequency of new cases. The ratio of this number and the total number of individuals we denote by f . When f is plotted as a function of the time it is called the epidemic curve. The curve initially grows exponentially until an inflexion point after which it behaves in such a way as to reach a final value which may be zero, in which case the epidemic comes to an end, or nonzero, in which case the disease becomes endemic.
To determine f , we consider all processes of the type where B represents an infected and A a non-infected individual. We then use the law of mass action to determined the rate of this process. The sum of the rates of all process of this type is f . In the case of the auto-catalytic process of the preceding section, in which an S is transformed into one I, the frequency of new cases is f = bxy.

B. Phase transition
The spreading of a disease can be viewed as a thermodynamic phase transition. To understand this point we consider that the infection rate constant, denoted by b, takes several values. For small values, there is no spreading whereas for large values the epidemic spreads. Increasing b, one passes from a non-spreading regime to a spreading regime at a critical value b c that gives the onset of spreading. This condition may seem to make no reference to the closeness of the individuals, or equivalently to the density of the population. In a real process of infection, the contact, or at least the proximity of two individuals, is an important condition for the onset of spreading. Thus this feature should in some way be included in the theory. In fact, this is the case, as explained next.
In the Appendix, it is shown that the infection rate where γ is the ratio of the total number of individuals in a community and the area of the community, or the population density. The intrinsic rate constant b * is a measure of the strength of infection or the virulence of the disease and is density-independent. Let us define the critical density γ c = b c /b * and see what happens when one increases γ = b/b * . When γ increases, b increases, and when it reaches the value γ c , the rate constant reaches the critical value b c , and the spreading outbreaks. In other words, when the density of the population reaches a critical value, the disease outbreaks. This is, in fact, the statement of the Kermack and McKendrick theorem, except that they refer to the density of the susceptible but at the beginning of spreading the population consist mostly of the susceptible. It is usual to characterize the transition from a trivial phase to the significant phase by an order parameter, which is nonzero for the significant phase and zero for the trivial one. For the cases of an epidemic spreading that comes to an end, so that the frequency of new cases f vanishes for the long term, the order parameter may be defined as the area s of the epidemic curve, where f here is understood as function of t. , recovered (R), and exposed (E). An I between parentheses indicates that the process is catalytic and I is the catalyst. The other processes are spontaneous. An exposed individual is infected but not infective. A recovered individual has permanent immunity in the SIR, SISR, and SEIR models.

C. Reproduction number
Another important characterization of the epidemic spreading is related to the number of individuals that are infected by a single infective individual. This number, called the reproduction number and denoted by R, is the ratio of two numbers and is defined more precisely as follows. During an interval of time τ , the number of new cases is N 1 = f N τ , where f is the frequency of new cases defined above and N is the total number of individuals.
The question now arises about the number N 2 of the infective that are responsible for these new infections. If the number of the infective remains the same in the interval τ , then N 2 would be equal to N 1 . But if the infective increases by a certain amount N 3 then N 2 should be smaller, and in fact equal to N 1 − N 3 , and where y is the fraction of infective we find As f is the rate of a catalytic process, in which the infective is the catalyst, f is proportional to y, that is, f = gy, which allows us to write The basic reproduction number, denoted by R 0 , is the reproduction number at the early stages of the spreading when the number of infective individuals is negligible and the whole population consists mostly of susceptible. Thus R 0 is obtained from R by setting x = 1 in the formula (9).
It is worth mentioning at this point a fundamental property that should be obeyed by the equations describing an epidemic. As the spreading of the disease needs the presence of an infected, then if y = 0 there is no evolution meaning that the other variables should not vary in time. Therefore, the state consisting of y = 0 and other variables constant should be a solution of the evolution equation. This state is the trivial solution and is called disease-free solution. In stochastic models, it is called absorbing state.
The stability analysis of the trivial solution is a way of determining the outbreak of the spreading. Linear stability reveals that the fraction of infective y is dominated by the exponential behavior, y = y 0 e αt , where α is the largest eigenvalue of the stability matrix. The onset of spreading occurs when α turns from a negative to a positive value. Thus, at the early stages of the spreading y = y 0 e αt and x near one. Replacing these results in equation (9) we reach a relationship between the basic reproduction number and the exponent α, which is where g 0 is the value of g when x = 1. As the onset of spreading occurs when α = 0, we see that the outbreak of spreading occurs at R 0 = 1, When α > 0, that is, when R 0 > 1, the epidemic spreads.

IV. MODELS OF SPREADING
The models that we will analyze consist of two or three processes involving two or more classes of individuals. They are represented by a flow diagram which gives all the possibilities of transformation of the individual from one class to another. The six models that we analyze are shown in figure 1.

A. SIS model
We start by considering the simplest model for epidemic spreading which consists of two classes of individuals: infective (I) and susceptible (S). The model comprises two transformation processes. In the first, a susceptible individual becomes infective by the contact with other infective individuals. This process is represented by and occurs with an infection rate constant b. In the second, an infective individual becomes susceptible spontaneously. More specifically, an infective individual that has recovered from the disease does not acquire immunity becoming immediately susceptible again. This process is represented by occurring with a recovery rate constant c. We denote by x and y the fraction of the susceptible and infective, respectively. As only these two classes of individuals are present, the sum of these two fractions equals one, x + y = 1. According to the rules of mass action, the equation that gives the evolution of x and y are Due to the constraint x + y = 1, these equations are not independent. Replacing x by 1 − y in equation (14), it becomes where α = b − c.
One solution of this equation is the trivial solution y = 0, and therefore x = 1, which is understood as the absence of epidemic spreading because there are no infective individuals present. To find whether this solution is stable one tries to solve the equation (15) for small values of y. The term proportional do y 2 is neglected and the equation becomes the linear equation whose solution is where y 0 is the value of y at t = 0. Therefore, if α > 0, that is, if b > c, then the number of infective increases exponentially and the spreading of the disease sets in. If α < 0, the disease does not spread. Recalling that α = b − c, it follows that b = c marks the onset of spreading. If b > c the epidemic spreads and if b < c it does not.
The equation (15) can in fact be solved in a closed form by writing it in the differential form as or yet as Integrating we find where we are considering solutions such that α ≥ by. Let us suppose that at t = 0, the number of infective is y 0 . This determines the constant k and we may write which gives y as a function of t.
We distinguish two types of solutions depending on the sign of α. Let us consider that α < 0. As we have seen above, in this case there is no spreading of the disease. The right-hand side of equation (15) is negative and y will decrease toward its asymptotic value y * which is found by taking the limit t → ∞ in equation (21), which is y * = 0. When α = 0, the exponential solution (21) is no longer valid. To find the solution in this case, one should solve the equation (15) for α = 0, which in this case reads The solution is For long times the decay towards the zero values is not exponential but algebraic, y ∼ (bt) −1 . When α > 0, the disease spreads and the fraction y of the infective as a function of time is shown in figure 2a for a small initial value of y. The fraction of the infective grows and then approaches and asymptotically the value y * . This asymptotic value is nonzero and is found by taking the limit t → ∞ in equation (15), with the result Alternatively, y * can be found as the stationary solution of equation (15) which is obtained by setting to zero the right-hand side of this equation. The frequency of new cases is the rate of the process (11) and thus given by From the solution given by equation (21), we may find f as a function of t, which is the epidemic curve shown in figure 2b. We assume that initially the number of infective individuals is small, which is the case of a real spreading of disease.
The initial increase of f is exponential, a result that reflects the initial exponential growth of y, given by the equation (17), that is, If b/c ≤ 2, the epidemic curve f increases monotonically towards its final value f * = by * (1 − y * ) or If b/c > 2, the epidemic curve increases initially, then reaches a maximum and then decreases to its final value (27). In any of these cases the final value of f * is nonzero which means that the disease becomes endemic. The order parameter s of the present model is the final fraction y * of the infective, given by (24), When b ≤ c, s = 0 and there is no spreading of the disease. When b > c, the disease spreads. A plot of s as a function of b is given in figure 2c. The reproduction number is given by equation (8). As f = bxy and dy/dt is given by equation (14) we find The basic reproduction number is found by setting x = 1,

B. SIR model
This model differs from the SIS model in an important feature. The infective individuals after healing acquires permanent immunization, meaning that they cannot be infected again. Therefore, in addition to susceptible (S) and infective (I) individuals, there are the immune individuals, which we call recovered (R). The model consists of two processes. The first is the infection, represented by occurring with an infection rate constant b, and the second is the spontaneous recovery, represented by occurring with a recovery rate constant c, which is strictly positive.
Denoting by x, y and z, the number of susceptible, infective, and recovered individuals, respectively, then according to the rule of mass action the evolution equation of these variables are dz dt = cy.
As the total number of infective, susceptible and recovered is constant, the fractions x, y, and z are related by x + y + z = 1, and the three differential equations are not all independent. The differential equations (33), (34), and (35) were introduced by Kermack and McKendrick in their study of epidemic spreading, allowing them to show the spreading threshold theorem [8], which we show next.
One solution of the evolution equations is x = 1, y = 0 and z = 0 which is the non-spreading state. To find the stability of this state, we consider small deviations from this solution. The equation for y becomes where α = b − c, whose solution is and y 0 is the value of y at t = 0. Therefore, if α > 0, that is, if b > c, the number of infective grows exponentially and the spreading of the disease occurs. If α < 0, y decreases and the disease does not spread. The onset of the spreading occurs when b = c. To solve the time evolution equations we take the ratio between equations (35) and (33) to reach the equation which can be integrated, The constant of integration was chosen by considering that at the initial times, there is no recovered individuals, z = 0, and the number of infective is very small so that the we may consider the fraction of the susceptible to be equal to one, x = 1. Replacing the result (39) into x + y + z = 1 we find the following relation between x and y, which replaced in (33) gives In the integral form it reads The integral can be solved numerically to get x as a function of t, from which we find y and z. Alternatively, these variables can be obtained by solving numerically the set of equations (33), (34), and (35). The result of y as a function of t, for a small initial value of y and z = 0, is shown in figure 3a. Initially there is an exponentially growth, as determined by the equation (37), then y reaches a maximum and then decreases with time and vanishes asymptotically. This property is a consequence of the fact that an infective individual that becomes healthy acquires immunity and cannot be infected again. For long times there will be no infective individuals but only the recovered individuals and the ones that did not have acquired the disease, the susceptible. The disease has become extinct.
The frequency of new cases is the rate of the process (31) and given by f = bxy, (43) and for the SIR model it is also f = −dx/dt. In figure 3b we show f as a function of time, the epidemic curve. It was obtained from the numerical solution of x referred to above. As initially the fraction of infective is increasing exponentially so does the quantity f . The epidemic curve increases, attains a maximum and then decreases to the zero value. The area s under the epidemic curve, which is a measure of the size of the epidemic, is the order parameter and will be determined next. When t → ∞, no infective is left, as we have already seen, because once an infective individual is healed, it acquires immunization and cannot be infected again. In other words, the infective becomes recovered and remains in this state forever. Thus when t → ∞, y → 0, and the final fraction of susceptible x * and the final fraction of recovered z * are related by x * + y * = 1. Replacing x * = 1 − z * in (39) we find We have seen above that f = −dx/dt. In we integrate this equation in time, we find the order parameter s, that is, where we have taken into account that x at the initial time is equal to one. Replacing the result z * = s in equation (44) we find the equation that determines the order parameter, In figure 3c we show s as a function of the infection rate constant b. When b ≤ c, s = 0 meaning that the size of the epidemic is zero or that there is no spreading of the disease. When b > c, the disease spreads and the size of the disease increases as the infection rate constant increases. The equation (46) is transcendental but can be solved numerically. However, a solution can be written explicitly when s is small, in which case we find 8 The reproduction number is given by (8). Using dy/dt, given by equation (34), and recalling that f = bxy, we find The basic reproduction number is found by setting x = 1,

C. SISR model
In the SIR model, whenever an infective individual becomes healed it acquires permanent immunity. Here we consider a modification of the SIR model in which an infective individual may become again susceptible [17,30]. There are three processes. The first is the infection, represented by occurring with an infection rate constant b. The second is the spontaneous recovery, represented by occurring with a rate constant c, where R represents an individual with permanent immunization. The third is the spontaneous recovery without immunization, represented by occurring with a rate constant a. An individual recovers from the disease but does not acquire immunization, becoming susceptible again. The equation for the fractions x, y, and z of susceptible, infective and recovered individuals are and again x + y + z = 1. One solution of these equations is x = 1, y = 0, and z = 0 corresponding to the absence of epidemic spreading. For small values of y and for x near one, the equation for y reads where α = b − c − a. Therefore, if b > c + a then the epidemic spreading sets in whereas if b < c + a there is no spreading of the disease. The ratio of the equations (55) and (53) gives which can be integrated, with the condition bx ≥ a. The constant of integration was chosen by considering that at the initial times, z = 0, y is negligible so that we may take y = 0 and x = 1. The frequency of new cases f is given by the rate of the process (50) and is f = bxy. The plot of f as a function of t, the epidemic curve, has a bell shape similar to that of the SIR model. To get the area under this curve, which is the order parameter s, we proceed as follows. By an appropriate combination of the evolution equations (53) and (55), we write so that Integrating f in time from zero to infinity, we get where we recall that at t = 0, x = 1 and z = 0, and x * and z * are the final values of x and z.
To determine the value z * , we observe that for long times y vanishes and x * = 1 − z * . Replacing this result in the equation (58), we find the equation for z * , Using the relation (61) and the result x * = 1 − z * , we find s = (c + a)z * /c and the size of the epidemic s is proportional to the final value of the fraction of the recovered z * . The transcendental equation (62) can be solved when z * is small, with the result The threshold of the spreading of the disease occurs when b = a + c. If b < a + c there will be no spreading and s vanishes. However, if b > a + c, the disease will spread and for long times the fraction of the recovered individuals is z * given by equation (62).
The reproduction number is given by (8). Replacing dy/dt given by equation (54) we find The basic reproduction number is found by setting x = 1, The numerical solution of the evolution equations give results for the fraction y and for the frequency of new cases f which are similar to those of the SIR model. Both these quantities grow exponentially, attain a maximum value and then decrease to their final zero values. The order parameter has also a similar behavior except that the critical point occurs when b = c + a.

D. SIRI model
Here we consider another modification of the SIR model. In the SIR model, an individual that has been infected becomes immune, which means that a recovered individual remains forever in this condition. In the present modification the recovered individual loses the immunization and may be reinfected again [17,41]. Therefore, to the two processes of the SIR model, occurring with an infection rate constant b, and occurring with a recovery rate constant c, we add the following process occurring with a reinfection rate constant a. The equation for the fractions of susceptible, infective, and recovered are and we recall that x + y + z = 1. The frequency of new cases are determined by the processes (66) and (68) and is given by One solution of the evolution equations is the trivial solution x = 1, y = 0, and z = 0, which corresponds to the absence of spreading. For small values of y and for x near one, the equation for y reads where α = b − c. Therefore, if b > c the epidemic spreads whereas if b < c there is no spreading of the disease. Dividing the equations (71) and (69) we find which can be integrated to give with the condition az ≤ c. The constant of integration was found by assuming as we did before that at the initial time, z = 0 and x = 1. For long times we have two types of solution. In one of them, the infected disappears, y * = 0, and there remain the recovered, z * = 0, and the susceptible, x * = 0. This solution occurs when a ≤ c and the asymptotic value z * is obtained by using (75) and replacing x * = 1 − z * . The result is the equation which should be solved for z * . For small values of z * the solution is The time dependent solution is similar to the solution of the SIR model. As we have mentioned, y * = 0 for this solution so that the frequency of new cases f vanishes for long times and the epidemic curve is also similar to that of the SIR model. The area s under the epidemic curve is given by the integral which can be obtained from the solution of the evolution equations. Let us consider the other solution, which is valid for a > c. The asymptotic values for this solution are x * = 0, y * = (a − c)/a, and z * = c/a, which are obtained by setting to zero the right-hand side of the evolution equations (69), (70), and (71). The time dependent solution leading to these values can be obtained by solving numerically the evolution equations. The fraction y of infected is shown in figure 4a for the case a/c = 1.1. From the solution for x, y and z we determine the frequency of new cases using (72), which is shown in figure 4b. This solution predicts a persistence of the disease because f remains finite for long times.
The reproduction number is given by (8). Replacing dy/dt given by equation (70) we find The basic reproduction number is found by setting x = 1, in which case z = 0, and

E. SIRS model
We consider here a model with three classes of individuals like the SIR model: susceptible, infective and recovered. However, the present model consists of three instead of two processes [17,26]. The first is the infection of susceptible individuals, represented by occurring with an infection rate constant b. The second is the spontaneous recovery, represented by occurring with a recovery rate constant c. These two processes are the same as the SIR model. In the present model, the recovered individuals are considered to have only partial immunity. Accordingly they may become again susceptible through an spontaneous process represented by occurring with a rate constant a. The evolution equations for x, y, and z, the fractions of susceptible, infective and recovered individuals, are Again x + y + z = 1 and the equations are not all independent. The trivial solution is x = 1, y = 0, and z = 0, corresponding to the absence of spreading. Near this solution, that is, for small values of y and for x near 1, the equation for y reads where α = b − c. Therefore, if b > c then y increases exponentially and the spreading of the disease occurs. If b < c, it does not.
To determine the asymptotic fractions x * , y * , and z * of each class of individuals at long term we set to zero the right-hand sides of the equations (84), (85), and (86). The solution with a nonzero value of y is The frequency of new cases f is determined by the process (81) and is given by (89) For long times it approaches a nonzero value given by and in this respect it is like the SIS model, that is, the model predicts a persistence of the disease for long times, or that the disease becomes endemic. However for the present model the time dependence may present damped oscillations as shown in figure 5. This behavior is shown below.
Let us linearize the evolution equations (84) and (85) around the solution (88). To this end, we define the deviations ξ = x − x * and η = y − y * and write the evolution equations in these new variables. Up to linear terms in these variables, we find The product of the roots is positive and the sum of them is negative. If the roots are real then they are both negative and the solution is stable. Now, let us take a look at the complex roots which occurs when In this case we should look at the real part of λ. Since the real part is negative, then the solution is stable. The conclusion is that the solution is stable but x, y, and z, and also f , will display damped time oscillations, as seen in figure 5, when the condition (94) is fulfilled. The values of a, b and c for this behavior is shown in figure 5c.

F. SEIR model
In some infection diseases, the individual that has been infected takes a certain time to be infective. To take into account the latent period of such a disease, we introduce a class of individual called exposed (E) who has been infected but is not infective, but eventually becomes infective. A model that includes the exposed individuals is similar to the SIR, but the directed infection is replaced by the following two processes. The first is represented by occurring with an infection rate constant b, and the second is represented by occurring with a latent rate constant h, the inverse of which is a measure of the latent period of the exposed individual. The recover of an infective is the same as the SIR model, represented by occurring with a recovery rate constant c. The evolution equations for the fractions x, u, y, and z, of susceptible, exposed, infective and recovered, respectively, are These equations are not all independent because x + u + y + z = 1. These equations were introduced by Dietz to account for the latent period of infection [36]. The frequency of new cases f is determined by the process (95) and is given by (102) Considering that bxy = −dx/dt, it follows that the area under the epidemic curve is where we have considered that at the initial times x = 1.
To find the solution of these equations we assume that u = u 0 e αt y = y 0 e αt .
Replacing this solution in the equations of evolution we find which are understood as a set of eigenvalue equations where α is understood as the eigenvalue. The possible values for α are the roots of the largest one being If α < 0, which occurs when b < c, the disease does not spread. If α > 0, which occurs when b > c, the disease spreads. In this regime, y and u grow exponentially with an exponent α given by equation (110). Let us compare this exponent with that of the SIR model, which is In both cases the onset of spreading occurs when b = c. However, when b > c, that is, in the spreading regime, α < α sir indicating that the epidemic growth rate is smaller for the SEIR model than it is for the SIR model. This is due to the necessary passage of an individual to the exposed state before reaching the infective state. The behavior of the present model concerning the epidemic curve and the time behavior of fraction of infective is qualitatively similar to those of the SIR model shown in figure 3. The behavior of the size s of the epidemic is the same as that of the SIR model and does not depend on the rate constant h as we show in the following.
We divide the equation (101) by the equation (98) which can be solved to give where the constant of integration was found by considering that at the early stages, z = 0 and x = 1. For long times the infective as well as the exposed disappear, y = 0 and u = 0. Therefore, the sum of the final values x * and z * of the fractions of the susceptible and of the recovered equals one and x * = 1 − z * . Recalling that the order parameter s = 1 − x * = z * we reach the following equation for s, which is identical to the equation (46) for the SIR model. It is worth mentioning that this equation says that s does not depend on h which means that the size of the epidemic is independent of h. In other terms, the presence of the exposed does not change the size of the epidemic. This process only slows the velocity of the spreading but not its size which is the same as that of the SIR model. The reproduction number is given by (8). Replacing dy/dt given by equation (100) we find To find the basic reproduction number we need to know the ratio u/y when x = 1, that is, in the early stages of the spreading. According to (106), this ratio equals u 0 /y 0 so that The ratio can be found from the eigenvalue equation (108). Dividing this equation by y 0 we find and R 0 acquires the form where α is the eigenvalue (110). A plot of R 0 is given in the figure 6 for various values of the latent period ℓ = c/h. Notice that, the SIR corresponds to the case when the exposed rapidly becomes infective, that is, when the latent period is zero.

V. THE ROSS THEORY OF HAPPENINGS
Ronald Ross received the Nobel Prize for Medicine in 1902 for his work on the transmission of malaria. He believed that it was not necessary to eliminate all mosquitoes to prevent the spreading of the disease but only to reduce the density of mosquitoes to a value below the critical level [11]. This conception as well as that of the threshold theorem of Kermack and McKendrick are understood as instances of the fundamental result of the statistical mechanics of interacting systems that the transition from one regime to another occurs only when one reaches a critical value of a parameter.
Ross advanced the idea of applying differential calculus to the dynamics of an infectious disease with the aim of its control and prevention. He even coined the name pathometry for what we call theoretical epidemiology. The theory introduced in his 1911 book on the prevention of malaria he called theory of happenings [7]. It is worth mentioning that the book contained the reports on the condition of malaria in several regions of the world, written by several authors. This includes the report by Oswaldo Cruz, the Brazilian epidemiologist, on the campaign on the prevention of malaria in southwest Brazil carried out in the first decade of the twentieth century.
His theory of happenings on the spreading of an infectious disease, Ross presented in an appendix to his book on the prevention of malaria of 1911 [7] and in subsequent papers [42,43], some of which with Hilda Hudson [44,45]. In a population of P individuals affected by a certain disease, the time variation dZ/dt in the number of the affected Z is hA+qZ where A = P −Z is the number of unaffected, h is the proportion of unaffected that becomes affected and q takes into account the demographic and recovery rates. For infectious diseases, Ross argues that h is proportional to the fraction x = Z/P of the affected, and writes h = cx where c is a constant which Ross calls the infection rate, and arrives at the following equation [7,42,43] The solution of this equation, given by Ross, is revealing that x increases slowly at the beginning and then very rapidly until it reaches an inflexion, and after that approaches the limiting value ℓ [42,43]. According to Ross the current proportion of new cases to the total population is [42,43] f = cx(1 − x).
The equation (119) and its solution (120) are identical to the equations (15) and (21), respectively, and the frequency of new cases (121) is identical to f given by (25). These results shows that the model considered by Ross is identified with the SIS model that we have analyzed.

VI. CONCLUSION
We have presented an analysis of deterministic models for epidemic spreading. The equations were set up by using the analogy between chemical reactions and the processes occurring in the epidemic spreading which are understood as a change of the class of an individual. The main analogy was the use of the law of mass action, which provides the rate of the several processes that define an epidemic model. We have also emphasized the analogy of the onset of an epidemic with a thermodynamic phase transition. When the infection constant changes it reaches a critical value at which the spreading takes place. The infection constant is argued to depend on the density of the population and as it increases the spreading outbreaks at a certain critical density, which is the theorem advanced by Kermack and McKendrick. We have also analyzed the quantities that characterizes the spreading of an epidemic. One of them is the frequency of new cases, which is the number of new infected individuals per unit time. When this quantity vanishes the epidemic comes to an end. In this case the area of the epidemic curve is a measure of the epidemic and may be identified as the order parameter. It may happen that the frequency of new cases does not vanishes and in this case the disease becomes endemic. We have seen that the SIRS model which is appropriate for this case, predict oscillations in the frequency of new cases.
Another way of determining the outbreak of the spreading is by means of the stability analysis of the disease-free state which is the state without any infective individuals. Any model of spreading must have this state which in stochastic models is called absorbing state. The stability analysis gives the behavior of the spreading at the beginning and shows that the growth of the epidemic is exponential. We have related the growth exponent with the basic reproduction number. When the exponent change signs, which means that the reproduction number passes from a value less than one to a value greater than one, the spreading outbreaks.
All the properties were obtained by solving the evolutions equations, which are ordinary differential equations of first order in time. This may be obtained in closed form or by numerical methods. We have finally showed that the SIS model was in fact the model originally studies by Ross in his theory of happenings.