Dynamics of infectious diseases: A review of the main biological aspects and their mathematical translation


 In this paper, the main biological aspects of infectious diseases and their mathematical translation for modeling their transmission dynamics are revised. In particular, some heterogeneity factors which could influence the fitting of the model to reality are pointed out. Mathematical tools and methods needed to qualitatively analyze deterministic continuous-time models, formulated by ordinary differential equations, are also introduced, while its discrete-time counterparts are properly referenced. In addition, some simulation techniques to validate a mathematical model and to estimate the model parameters are shown. Finally, we present some control strategies usually considered to prevent epidemic outbreaks and their implementation in the model.


Introduction
Infectious or communicable diseases are caused by a biological pathogen such as virus, bacteria, protozoa, toxins, etc. Different transmission mechanisms affect the spread of the disease as, for instance, direct physical contact, aerosol droplets of an infected individual, passive vectors (water, food, etc.), active vectors (mosquitoes, ticks, etc.). In addition to these ways of horizontal transmission, there could be vertical transmission, i.e. from a mother to her child before or at the moment of the birth.
A mathematical model of an infectious disease is a (mathematical) formula which describes the transmission process of the disease. To formulate a mathematical model, the understanding of the main biological features of the disease (states and parameters involved) and the relation among them is essential. This is usually informally reflected by a graphical representation called the transfer diagram of the model. From the transfer diagram, the model is often formalized by a system of equations.
The mathematical analysis of the model provides a qualitative long-term view of the transmission process, while simulations help to understand quantitatively the short-term behavior of the disease. The model parameters can be calibrated by comparing with empirical data. Subsequent validation of a model also allows us to test its underlying hypotheses. This task may be difficult or even impossible to be done, depending on the quality of the available data. Understanding the transmission process characteristics of a specific disease can help to decide what control measures can be taken in order to prevent its outbreaks.
The purpose of this review is to illustrate the key biological features that should be considered to formulate and analyze infectious processes. While reviewing the existing literature, we provide the main biological aspects and their mathematical translation, although this may be someway incomplete. We will mainly focus on deterministic models in continuous-time, while discrete-time counterparts will be properly referenced.
The review is organized as follows: Section 2 presents the most relevant biological features, which lead to different compartmental approaches and their corresponding transfer diagrams. In Section 3, we revise the most relevant classes of mathematical models, according to the main biological aspects affecting the disease dynamics. In particular, we deal with the inclusion of heterogeneity in infectious disease models. That is, special factors that could be relevant for the modeling of the phenomena, although they make the model more complex than when homogeneously considered. Section 4 briefly describes the mathematical tools and techniques that can be used for a qualitative study of an infectious disease model formulated by a system of ordinary differential equations, while referencing their counterparts for a system of difference equations. Section 5 is devoted to shown some simulation techniques to validate the mathematical model and to estimate the model parameters. In Section 6 we present some control strategies usually considered to prevent epidemic outbreaks and their implementation within the model through new variables or parameters. Finally, we present the main conclusions from this review.

Main biological features of infectious disease models
In this section, we review the most important (biological) features that should be taken into account to reflect the disease behavior. More specifically, we describe the main (possible) states of a disease and the forms of transition among them in relation to fundamental parameters. This allows us to establish the class of the disease, according to the states and transitions initially observed. Likewise, other relevant features as demography or other heterogeneity factors are introduced. Some transfer diagrams are shown, as examples of schematic representation. A transfer diagram should be always obtained as a conclusion of the observed biological features.
R. It is usually expressed as γI, where 1/γ represents the average residence time in the I compartment, i.e. the time for the disease to be discovered in the infected individual. When the disease is fatal and every infected host finally died, the compartment R is interpreted as the class of death or removed individuals.
Of course, other rates or coefficients different from the incidence rate and the recovery rate could be considered as model parameters, if we observe transference between other compartments or states of the diseases under study.
Usually, a compartmental approach leads to a schematic representation by means of a transfer diagram. Each compartment in a transfer diagram is represented by a box labeled with the initial capital of the corresponding population class. Then, arrows indicate the movement of individuals from one compartment to another, depending on the above indicated transference rate, as shown in figure 1 [91]. The transfer diagram in  [81][82][83]. It is a simple epidemic model of SIR type, predicting the evolution of the number of cases of an infectious disease, as it spreads through a population. It represents the first compartmental model properly studied [26,91]. For a historical review of the development of models for epidemics, there exists a wide variety of papers and books describing the different views, such as [16] or [56].

Demography and heterogeneity factors
The transmission of the infection is usually firstly described by a homogeneous approach. That is, as a first step, it is usually assumed: total constant population, homogeneous mixing, with all the hosts distributed uniformly, and having homogeneous (horizontal) infectiveness, with identical rates leading to new infections at any time.
In particular, variations of the total population are referred as demography. To add demographic effects in compartmental models, assumptions are made about the migration, birth and death rates. If demographic effects are not considered in a compartmental model, this model will only be appropriate to simulate epidemic outbreaks of relatively short duration (under a year). However, in many cases, in spite of demographic effects, it is assumed that the total population remains constant, and homogeneously mixed. The assumption of homogeneous mixing is a simplification that renders the analysis more manageable. However, sometimes, it does not fit reality as much as desired [7,20,26,65]. Nevertheless, several homogeneous approaches have proved to have a strong predictive power [20].
In view of the homogeneous mixing and homogeneous infectiveness oversimplifications, which do not capture other more subtle features of reality, in the last years, other approaches have been developed based on non-homogeneous or heterogenous assumptions. That is, spatial heterogeneity, population heterogeneity, transmissional heterogeneity and seasonal heterogeneity have been considered as relevant biological features to be reflected in the transfer diagram describing the disease. Spatial heterogeneity refers to non-homogeneous mixing, i.e. non-uniform distribution of the host within the habitat. Population heterogeneity concerns the nonhomogeneous infectivity with respect to different population groups (age, sexual or social groups). Transmissional heterogeneity is related to other routes of transmission different from the horizontal one, as vertical (from mother to child) or vectorial (by means of animals/insects or fomites). Seasonal heterogeneity is associated with differential infectivity depending on periods of time.
Likewise, some infectious agents, such as bacteria, can replicate outside their hosts in environmental reservoirs (e.g. water, food, etc.). Humans acquire the infection indirectly from the pathogen-contaminated environmental reservoir, as is the case of cholera (see [30,57] and the references therein). In fact, even common airborne infections can spread indirectly by a virus transferred through an intermediate object, for instance, contaminated hands or fomites [25]. To reflect this feature, a box is added in the transfer diagram, for the environmental pathogen concentration, which is denoted by P. The incidence rate then becomes λ (P)S [122]. Figure 2 shows the transfer diagram of an approach with an intermediate infectious vector (e.g. mosquitoes, ticks, fleas, rodents, etc.). The indirect transmission is represented by the terms T h S h I M /N M and T M S M I h /N h [59]. An example, with several compartmental approaches for the malaria transmission dynamics, with different levels of complexity of host-vector-parasitic interactions, is found in [90]. In case of dengue transmission in humans, a review of compartmental approach is given in [9]. The Zika virus instead can be transmitted both indirectly and directly through sexual contact at rate β S h I h [28]. In addition, other models where, in the transmission of the disease, the virus is eliminated from the infected individual and is acquired by the susceptible [20,25,28,98], models with a saturation effect in the transmission rate, for which a higher density of infected individuals decreases their per capita infectiousness, and situations where multiple exposures to an infected individual are required for an effective transmission to occur [84], have been considered.

Mathematical formalization of infectious diseases.
This section presents the different formulations to (mathematically) model the main biological aspects mentioned in the previous section. Likewise, this leads to a mathematical model classification, so differentiating deterministic versus stochastic, continuous versus discrete, etc. In addition, we show how the different heterogeneous aspects, which could influence their evolution (e.g. seasonality, age structure, spatial mobility, etc.) can be implemented using various mathematical tools.

Deterministic and Stochastic Models
Deterministic models are those whose parameters have fixed (proportional) values and whose state variables evolve due to fixed rules, which are continuous functions of time. They are usually established by continuous or discrete-time formulations, i.e. differential or difference equations (see, for instance, [7,17,27,29,32,44,65,79,91]). They are appropriate if the epidemics occurs in large populations, but they describe the global behavior of complex systems consisting of multiple elements, without taking into account the local interactions between individuals.
Stochastic models instead are more effective for small populations [79]. To formulate such models, one needs to specify a probability law for the time and type transition. In contrast with determinist ones, they allow us to capture the stochastic nature of one-to-one transmission. They make use primarily of Markov chains and probability distributions, with continuous or discrete states. They better account for individual interactions: most of such models are based on cellular automata [118], individual agents [101] and networks [98,100]. A detailed introduction to stochastic epidemic models and their statistical treatment can be found in [3] and [8], respectively.
On the other hand, in [5], one can find a good comparison between deterministic and stochastic SIS and SIR models in discrete-time.
In the rest of the paper, we focus on deterministic models and mainly in those given by differential equations, although we reference properly the counterparts for those given by difference equations.

Continuous-time Models
In continuous-time models, their variables can take an infinite number of values within a given range. The vast majority of compartmental deterministic models for studying the spread of diseases are based on the use of non-linear differential equations. For instance, the SIR Kermack and McKendrick model [81], with the assumption that the population size is constant, i.e. S(t) + I(t) + R(t) = N, by replacing the compartments in Figure 1 with their mathematical descriptions, becomes: However, due to the constraint on the total population size, the evolution of the number of individuals in the infected class I can be studied just using the first two differential equations [91].
The behavior of this model depends fundamentally on the threshold parameter R 0 , which determines the stability of the disease-free and the endemic equilibria. Another important asymptotic property in epidemic models is the size of the epidemic, defined as the total number of individuals infected at the end of an epidemic, namely I(∞) [91]. If R 0 < 1 the infection will be reduced, the number of infected individuals will decrease until disease eradication, and the system will attain an equilibrium state with no infected. If, on the contrary, R 0 > 1, then the infection will spread, i.e. the number of infected individuals will grow, to reach a stable endemic equilibrium, where the disease persists. The number R 0 depends on the parameters involved in the system: for the model (1), it is R 0 = β N/γ. This can be interpreted as the amount of the population that has effective contacts, β N, in the time unit with an infectious individual within a totally susceptible population of size ≈ S 0 , during its average infectious period, 1/γ. If R 0 ≤ 1, it follows β S ≤ β N ≤ γ. Thus, from the second equation in (1), I(t) decreases monotonically to zero as t → ∞, so that the disease gets extinguished. Instead of that, if R 0 > 1, from β N −γ > 0, initially I(t) starts growing leading to an epidemic outbreak, reaching its peak and then decreasing as t → ∞. In view of this result and the settings of (1), S(t) and R(t) are, respectively, a monotonically decreasing and a monotonically increasing function. The maximum number of infected individuals at any time is obtained from I = 0, i.e. when S m = γβ −1 . This maximum is given by For this model, we have σ = β /γ. The number of replacements R at the initial time t = 0 is the product of the number of contacts multiplied by the number of susceptible individuals at the initial time, i.e. R = σ S 0 = (β /γ)S 0 = R 0 . As claimed earlier, at the initial time the basic reproductive number and the number of replacements are equal. Based on the model of Kermack and McKendrick, many other types of models have been developed, incorporating more features and details to describe the transmission dynamics of the infectious diseases [44,65]. A more complex model than (1) may include demographic aspects in the population, i.e. birth and death rates assumed to be equal to each other, µ. This reformulation of (1) reads: In such a case, R 0 = β (γ + µ) −1 . Thus, the value of R 0 will depend on the epidemiological characteristics of the disease and the demographic features of the population, µ.
The previous two models consider only direct horizontal transmission, i.e. contacts among susceptible and infected individuals. But for certain diseases such as AIDS, toxoplasmosis or the Zika virus the infected mother transfers the infection to her offsprings, i.e. vertical transmission occurs [32]. Model (2) can be modified as follows to account for this feature. A fraction q of descendants of infected individuals are born infected, and correspondingly a p = 1 − q fraction are born susceptible. Birth and mortality rates are still constant in the class of susceptible and recovered, b > 0, whileb > 0 denotes the birth and death rate for the class of infected. The disease does not induce mortality, and once again the total population size is fixed, N = S + I + R. The model takes the following form: Here note that from the second equation in (3), the basic reproduction number can be written in two different ways: The first formulation shows that the infection disappears, namelyR 0 < 1, when the number of new infections, obtained either by horizontal β N or verticalbq transmission, per average length of the infectious period, 1/(γ + b), is less than one. The second one indicates that the vertical transmission has the effect of reducing R 0 , namelỹ Indeed, if vertical transmission is allowed, there will be less horizontal contacts needed to spread the disease. The first model incorporating a nonlinear transmission rate appears in [34]. Formulated to understand the rapid curbing of the cholera epidemics of 1973 in Bari, Italy, it contains a contact rate that decreases as the infected population grows, such as a gamma function In this way the individuals response to the spread of epidemics is modeled, for which suitable measures are taken to avoid possible contacts with infected individuals.

Discrete-time models
A good reason for working with epidemic models in discrete-time is that data are collected and/or reported in units of discrete-time and the models can be adjusted to the collected data [2]. In addition, numerical scanning of discrete-time epidemic models is quite simple and can, therefore, be easily implemented also by nonmathematical experts [27]. The advantage of discrete-time over continuous-time modeling is that the former one allows an easier comparison of their results with empirical data, which are given discretely. Further, these types of models are best adapted for populations that have several generations a year. For instance, insects reproduce several times in a summer, so that N k would represent their population at generation k. Transfers between compartments or states can depend on the size of the compartments in the past as well as at the instant of the transfer. Numerical simulations of continuous-time models are discretized versions of these models and therefore related to discrete systems. A simple relationship of this sort is given by the Euler scheme [12,70,71,111,119] or other non-standard difference schemes [13,75,130].
A discrete model to describe the evolution of West Nile encephalitis with weekly data from New York City in 1999 is proposed in [120]. It allows us to determine the specific amounts of periodic fumigation which are needed for virus eradication. A comparison of the similarities between single-outbreaks in classical continuoustime epidemic models and their discrete-time counterparts, as well as the final size of an epidemic are studied in [29]. A discrete model for the transmission of Babesiosis between cattle and ticks shows that similar conclusions can be obtained for both discrete and continuous models. However, the former one needs some parametric restrictions that are not present in the continuous version. But the parameter constraints are natural and reflect the epidemiological robustness of the model [12]. The model is described by a system (4) of non-linear difference equations and it is epidemiologically meaningful if and only if the following parametric restrictions are satisfied: The transmission and control of the SARS infection in China was studied in [133] via a discrete mathematical model whose numerical simulations fitted the real data well. An age-structured model for the spread of tuberculosis in China, where the discrete-time intervals correspond to the reproductive periods or the average duration of a particular stage, can be found in [33].

Networks models
While global models study the behaviour of complex systems composed of multiple elements, disregarding the local interactions between individuals, individual-based models, on the contrary, take these individuals interactions into account. Network modeling approaches can be useful to incorporate a heterogeneous contact structure. They are most useful when each individual is in contact with a small fraction of the population [20,78].
The first network-based epidemic approaches were developed to study the spread of sexually transmitted diseases, the so-called couple formation models [63]. Extending these models to include long-term simultaneous interactions, led to contacts network models, described by a graph with nodes representing individuals and links representing their contacts (see Figure 4). The different structural properties of the network can be exploited to study the epidemic propagation speed [78,98,100]. In the so-called small world networks, most contacts between individuals are local. But it is enough that some long-distance contacts exist, to ensure a fast global spread of an epidemic, which is a threat in an increasingly globalized world, where individuals move for economic, social, cultural and political reasons. These network models have been applied to study the spread of respiratory diseases such as mycoplasma pneumoniae [95], influenza A(H1N1) [60], as well as vector-borne diseases [132] and cholera [131].

Spatial heterogeneity
The mobility of the population affects to a large degree the speed of transmission of the disease. To study the relation between the mobility of the population and the propagation of the disease, space must be included in the mathematical description [14,97]. In a compartmental model, this can be reflected by partitioning the population into spatially different sites (e.g. cities, countries, towns), called patches, giving rise to metapopulation models [10,14,91,96]. Transfer diagram for an SIR metapopulation compartmental model (see [91]).
Metapopulation models were first introduced in ecology in 1997 [64]. The patches are connected through transference of individuals. This approach has been effectively extended also to modeling infectious diseases [14,79]. For epidemic metapopulation models, an adequate modified version of the SIR epidemic model is given where An SEIR model that combines several of the above discussed features, specifically metapopulations with integrodifferential or differential equations of fractional order, has been formulated for the spread of mumps infection in four New Zealand cities in [49].
Space can be incorporated as an independent variable on a continuous or discrete basis. Continuous spatial models are based mainly on partial differential equations, concretely, on reaction-diffusion equations [22,97,105]. However, more recently, in reaction-diffusion models, the nearest neighbor contact, as the main infection transmission mechanism, has been considered. In this approach, the infection is considered in a population of immobile individuals, and the propagation occurs by a local contact mechanism between individuals and the part of the population nearest to them [103].
To describe the spatial structure in the SIR model (1) discussed above, we can consider that the individuals move randomly in a limited one-dimensional habitat 0 ≤ x ≤ L. For this purpose, we can assume that the densities of susceptible, infectious and removed depend on space as well as time, S(t, x), I(t, x) and R(t, x). The dynamics of (1) is then described by the following system (6) of partial differential equations with corresponding initial conditions where D is the diffusion coefficient related to spatial dispersion. The sub-indices indicate the partial derivatives taken with respect to t and x, respectively. In addition, boundary conditions need to be specified. They can be of Dirichlet, Neumann, or Robin type [91]. System (6) was successfully used to describe the geographic spread of the bubonic plague in Europe in the years 1347-1350, as well as the wolf rabies in England in 1985 [91,97]), the spatial spread of the West Nile virus that reached New York in 1999 and spread to the West Coast of North America in 2004 [88], or the bovine babesiosis in [1]. Spatial heterogeneity leads to consider the possibility of heterogeneous contacts Recently, different approaches based on epidemiological networks have been developed to explicitly model the heterogeneity of disease contacts independently of host contact rates. Each individual has a finite number of potential interactions by which the infection can be transmitted. This contact structure profoundly affects the dynamics of the infection [20,78,100]. Diseases affecting interacting populations have been reviewed in [129].

Population heterogeneity
Several characteristics could contribute to population heterogeneity. A disease can present a differential degree of infectivity depending on the social, ethnic, sexual or age group. In particular, age usually gives a source of heterogeneity and, due to that, age structure models are commonly studied. One of the most important characteristics in the modeling of infectious diseases is the age, because individuals respond differently to diseases according to their age. Age, of course, usually influences the single individual reproductive and survival capacities. But this occurs also for diseases. In general, older individuals are indeed weaker than young ones to fight against an infection. The same could be said for newborns. As diseases may have different infection rates and mortality rates for different age groups, the risk of acquiring and transmitting a contagious disease through multiple pathways changes with the biological age [4,79,91]. We must not confuse biological age of individuals, with the "age" of the disease, i.e. the time elapsed from the moment in which the individual has become infected. A model that considers both biological age as well as the age of the disease, i.e. the stage of the disease is presented in [128]. An example of an SIR model structured by age groups. The population is partitioned into two age groups, juvenile, and adult, and in each of these age groups the subpopulation is again partitioned into susceptible (S), infected (I) and recovered (R). Each of these age groups has a different level of interaction within their own age group and with people from the other groups.
Models incorporating population heterogeneity may exhibit fundamentally different behaviors compared to homogeneous infectivity models. These models can be used to compare, for instance, different alternative agespecific immunization programs (see Figure 5). A broad spectrum of age-structured compartmental models can be found in [4,79].
Age structure can be described as a continuous variable. In fact, this description already appears in the first epidemic model [93]. This approach leads to hyperbolic partial differential equations models [7,44,90]. Thus, these equations are analytically more difficult to be handled than ordinary differential equations. However, a numerical solution can be obtained by finite difference schemes [91].
Let us introduce the population density in age and time N(t, a), with t > 0 and a ∈ [0, A), where A < ∞ represents the maximum age of individuals. Further, where S 0 (a), I 0 (a) and R 0 (a) are suitable nonnegative known functions, and µ(a) and γ(a) represent the natural mortality and recovery rates, respectively. The sub-indices indicate once again the partial derivatives, now with respect to t and a, respectively. In this situation, the problem is an initial-boundary value problem defined in (t, a) ∈ [0, ∞) × [0, A], as the populations densities need to be specified at age 0. The boundary conditions defined at age 0 for each subpopulation indicate their respective newborns. Therefore, we define a maternity function b(a) ≥ 0, giving the number of expected offsprings for individuals of age a. Assuming that every individual in the population, even if infected, can reproduce and all newborns are healthy, i.e. in the absence of vertical disease transmission, we set Assuming intracohort mixing [79], the infection rate is given by λ (t, a), expressing the fact that the infection can be transmitted only between infected individuals of the same age. This may be appropriate for some childhood illnesses. The strength of the infection is given by With intercohort mixing, the transmission occurs between all age classes, since contacts can occur between an infected of age α with a susceptible aged a. In such a case, the strength of the infection is defined by Assuming the separable form β (a, α) = β 1 (a)β 2 (α), it can be rewritten as For further details, see [4,31,74,91].
Compartmental models can be adapted to reflect age groups [15,79,117]. Figure 5 shows a compartmental model with just two age stages; juvenile and adult. This age-structured deterministic SIR model is represented by the following system (8) of differential equations: A broad spectrum of models using this approach can be found for instance in [4,79,[108][109][110].

Transmissional heterogeneity
Although in traditional literature it does not appear as an heterogeneity factor, it is obvious that vertical transmission and indirect transmission by means of vectors or fomites, leads to a non-homogeneous infectivity. This is the case of the majority of vector-host diseases where, although homogeneously distributed, infected hosts cannot infects other hosts, but only vectors (see for instance [11,92,106]. Nevertheless, the modeling of these kinds of heterogeneity can be usually set out by means of new state variables and incidence rates when transmission is possible. For instance, the following system (9) of ordinary differential equations was the first established one for modeling the Babesiosis on bovine and ticks populations [11].
In this model, the variables with subindex B represent the bovine population, while the variables with subindex T represent the ticks population. Concerning the parameters, µ B is the birth and mortality rates for bovines, which are assumed identic. The same can be said for the parameter µ T , but in the case of ticks. In addition, β B and β T denote, respectively, the transmission coefficients from bobine to ticks and from ticks to bovine populations. Likewise, λ B is related to proportion of recovered bovines, and α B to the proportion of recovered bovines that go back to the susceptible class since, as occurs in many vector-host diseases, the model includes temporal immunity for the interesting species [91]. On the other hand, p is associated with the vertical transmission in ticks population. That is, p denotes the probability of having a susceptible tick from an infected one.

Seasonal heterogeneity
The transmission of childhood infectious diseases has been proven to vary seasonally, peaking at the beginning of the school year and decreasing significantly in the summer months [17]. Human influenza exhibits seasonality in temperate zones and tropical climates, showing seasonal patterns of infection. Influenza and measles epidemics tend to occur during the winter, and in tropical regions. Hence, these diseases have a seasonal variation that seems to have a strong relationship with the local factor. Differences in the efficiency of influenza transmission at different levels of absolute humidity have been postulated as a possible reason for seasonal variation in influenza transmission. Possible drivers of seasonality are often location-specific. Diseases transmitted by vectors (malaria or dengue) are often highly seasonal [50]. In [79], a detailed literature review of this compartmental modeling approach is provided. It is therefore generally assumed that climatic conditions or seasonal variations lead to seasonal variation in incidence.
The seasonality of propagation is modeled by incorporating periodicity, between susceptible and infected individuals, by a variable contact rate β (t), the seasonally forced function [76]. Models with different infectious periods would require the use of several infectious compartments. However, a more flexible formulation can be achieved by memory models [79], which explicitly take into account the time elapsed from infection, thereby tracking the history of the infection. Most of them are formulated via delay differential equations or integrodifferential equations (see Chapter 4 of [32]). More recently, also fractional order derivatives have been used in epidemiology to account for memory effects. For the inclusion of this kind of derivative in the simple SIR and SIS models see [52,61], while in [36,58,99,112] more specific diseases are investigated with the same tools. Let denote the derivative of fractional order θ of the function U in the sense of Caputo [35], where Γ is the Gamma function.
The reformulation of the SIR model (1) via a system of fractional order differential equations becomes as shown in (10).
While, in general, the equilibria remain the same as for the ordinary system, their stability is usually affected by this reformulation.

Mathematical study of deterministic models
This section reviews some mathematical tools and methods used to carry out a qualitative study of a (deterministic) epidemic model, mainly of those given by a system of differential equations.
For the analysis of epidemic models described by differential equations, where f (x) is an autonomous nonlinear vector field of the form the theory of differential equations and continuous dynamical systems can be applied [27,68,85,91,102]. The analysis of a model usually begins with the study of the feasible region Ω in which we could situate and find the stationary or equilibrium solutions, that are constant solutions of the model. For the model (1), of dimension three, the variables S, I, R are the components of the vector x and the components of the function f are the explicit expressions in terms of the variables S, I, R on the right hand side of (1). In epidemic models, there is usually a disease-free equilibrium and one or more endemic equilibria, with persistence of infected individuals. Their stability is related to the function f (x). Linearization near each equilibrium leads to the assessment of the eigenvalues of the Jacobian matrix of f at each stationary point. If all of them have negative real part, the equilibrium is stable. In the discrete case, this should be replaced by the condition of having modulus less than 1 at the equilibrium or fixed point. Generally, for epidemic models, if the parameter R 0 representing the basic reproduction number does not exceed the unity, R 0 < 1, the only equilibrium is the disease-free equilibrium and this equilibrium is asymptotically stable. On the contrary, if R 0 > 1, the usual situation is that the diseasefree equilibrium is unstable and there is a unique endemic equilibrium that is asymptotically stable. This is commonly known as the threshold phenomenon for infectious diseases. This approach can also be applied to models with vertical transmission [32], or to vector-host diseases [27,91]. We present a brief review of the main mathematical tools and methods needed in epidemiology, focusing just on the models described by a equation like in (11). For further information, see for instance [68,102] and the references therein. We will deal with the following topics: positively invariant regions, equilibria and their local stability, global stability and finally bifurcation and existence of periodic solutions.
This defines a vector at each point (x 1 , x 2 , . . . , 0, . . . , x n ) of this face. A path of the solution must have this direction, and if the vector points inward to R n + , the path can not leave R n + . Therefore, it is sufficient to show that the vector field points towards the interior of the boundary of R n + , and, by the existence and uniqueness theorem of solutions [68], trajectories cannot intersect each other.
More generally, let Ω be a region in R n + and let x be a point on the boundary of Ω. Let n denote the outer normal to the boundary of Ω. Then, if the direction of the vector field v = f (x) and n have non-positive inner product, v · n ≤ 0, Ω is positively invariant for the system (11).

Equilibria and local stability analysis
When analyzing a system of the form (11), the first step is usually to determine its equilibria. An equilibrium is a (constant) point x * satisfying the equation. Since it does not vary in time, it can be found by solving the nonlinear (algebraic) system f (x * ) = 0, obtained by annihilating the derivatives on the left hand side of (11).
Secondly, one must assess the stability of the equilibria. An equilibrium x * of (11) is stable if any solution starting at t = t 0 sufficiently near this equilibrium point remains in its neighborhood (stability) or tends back to the equilibrium (asymptotic stability) as t → ∞. Assessing stability is important in order to know the behavior of the system. To assess the local stability, the linearization method is illustrated below.
Note that exogenous disturbances have important effects on equilibria, but they are absent in the case of autonomous systems, as it is the case here. In general, time-dependent forcing functions alter the system behavior and determine whether or not the steady state will persist under the inevitable disturbances that occur in any real situation being modeled.
Linearization method: Let x * be an equilibrium of the system (11) and let J(x * ) = D f x * denote the Jacobian matrix of f evaluated at x * . The linearized system of (11) near x * is An equilibrium is said to be hyperbolic when none of the eigenvalues of J has zero real part. In the discrete case, this condition corresponds to eigenvalues with modulus different from 1. Hartman-Grobman's theorem proves that the solutions of an n × n autonomous system of ordinary differential equations in a neighborhood of a stable hyperbolic equilibrium can be qualitatively seen as the solutions of its linearized system near that equilibrium point [102]. Thus, the local behavior of the solutions of (11) are determined by the eigenvalues λ of J, i.e. the roots of the algebraic equation Specifically, the condition for local asymptotic stability is that all eigenvalues of J have a negative real part. It should be noted that local asymptotic stability of an equilibrium does not mean that all solutions tend to this point when t → ∞, but only those solutions which at some time are close to the equilibrium. A given system may have several equilibria, some stable and others unstable. In addition, it could behave also in other different ways, such as exhibiting periodic or unbounded solutions. Some of these possibilities are explored in the following subsections. Knowledge of the local stability of all equilibria gives an information on the dynamic behavior of the system. If the regions of attraction of the stable equilibria are small, the local information may be of little help to explain the whole behavior of the system. A global stability analysis would then be required. This should be underlined because most of the stability analysis in epidemic models are restricted to local stability. For the assessment of the nature of the roots of (13) the Routh-Hurwitz criterion can be used [72,91,107,125]. They give a set of necessary and sufficient conditions for stability. In particular, for the case n = 2, all eigenvalues of the matrix J(x * ) of 2 × 2 have real negative parts if and only if tr(J(x * )) < 0, det(J(x * )) > 0; (14) for the case n = 3, they become tr(J(x * )) > 0, det(J(x * )) < 0, tr(J(x * ))M 2 > det(J(x * )), where M 2 denotes the sum of all the principal minors of order 2 of J(x * ). In the discrete-time case, the Jury criterion offers similar conditions to determine the modulus of the eigenvalues at an equilibrium or fixed point [27,77,91].
In compartmental approaches of diseases, three important disease threshold parameters can appear in relation to the parameters in the description of the system: the basic reproduction number R 0 , the number of contacts σ and the replacement number R. The concept of the basic reproduction number was firstly introduced by Ronald Ross in a study of the transmission of malaria [90]. R 0 represents the most important threshold for many epidemic models, being defined as the average number of secondary infections produced by an infected individual in a totally susceptible population. This parameter often represents the threshold determining whether a disease can invade and persist in a population [127]. More specifically, if R 0 > 1 an epidemic outbreak will occur, while if R 0 < 1 the disease will become extinguished. If R 0 = 1, each individual will replace himself and there will be no epidemic outbreak [48]. After the spread of the disease, the fraction of susceptible individuals is less than 1. Thus, not all suitable contacts turn out to be a new case. When populations include infected individuals of different types, a sophisticated approach defines R 0 as the maximum eigenvalue of the Next Generation Matrix (NGM) (see [6,46,47,126] and Chapter 5 in [91]). Note that, at the beginning of the spreading of the disease, these three threshold quantities coincide, if all the individuals are susceptible. R 0 is calculated at the beginning of the spreading of the disease, while σ , and R can be assessed at any time. The replacement number R is always smaller than the number of contacts after the spreading, i.e. R 0 ≥ σ ≥ R [48,65].

Global Stability Analysis
Given a continuous system of the form (11) in the plane, i.e. n = 2, there can be two behaviors of great relevance for the global behavior of the trajectories: equilibria (or fixed) points and periodic (or closed) orbits.
We now present some criteria to prove non-existence of periodic or closed orbits. This situation implies that if there is only one equilibrium point in Ω, it turns out to be globally asymptotically stable [32].
A periodic solution is a function such that x(t + p) = x(t) for all t ≥ 0 where p is a fixed nonvanishing constant called the solution period. In the space R n , this solution is represented by a closed curve because the path starting at x(t 0 ) returns to this point after a time p. The closed curve is also called a limit cycle if it attracts all nearby paths containing some points not lying on the limit cycle. Consider a bidimensional system. Let x(t) = (x 1 (t), x 2 (t)) be a solution, t n a sequence of real numbers such that t n → ∞ when n → ∞, and consider the point P n = x(t n ). P is called an ω-limit point of the curve x(t) if P n → P, while if P n → P for t n → −∞, it is an α-limit point. The sets of all these points are respectively called the ω-limit and α-limit sets of x(t). The importance of the ω-limit set lies in the property that it is an attractor for the system trajectories. All trajectories in a bounded region of the plane will indeed spiral toward the ω-limit set.
If a curve x(t) lies in a bounded set for all t ≥ 0, the ω-limit set is nonempty; further, it is closed, bounded and connected. The qualitative theory of differential equations aims to classify all limited sets of a given system. The results constitute the Poincaré-Bendixson theory [102]. Essentially, the Poincaré-Bendixson theorem tells us that ω(x 0 ) can contain a fixed point or closed orbits of the flow.
This theorem can be used to show the global asymptotic stability of an equilibrium point if periodic orbits can be ruled out in a positively invariant set. A technique to achieve this result is the Bendixon-Dulac criterion [68].
Dulac-Bendixson Criterion: For the system (11), assume g(x) to be a continuously differentiable real valued function, defined in a simply connected region D ⊂ R 2 . Then, D does not contain periodic orbits if As mentioned, the Poincaré-Bendixson theorem holds only in dimension 2. For higher order nonlinear systems, it is not valid. The system can exhibit many more complicated dynamic behaviors. Although in the last years some researchers have obtain generalization for systems with higher dimensions [54,89,116], the global analysis of the trajectories is much more difficult to assess.
An efficient method that can sometimes be applied is the Lyapunov method. This method consists in constructing a suitable function that forces the system trajectories toward the equilibrium point x = 0.
Lyapunov's direct method: Consider the system (11) and assume that there is a continuously differentiable function V : D → R such that Then, Let G * the maximum invariant subset of the set {x ∈ G;V (x) = 0}. As V (x(t)) decreases along a solution, x(t, x 0 ) decreases along a system solution of (11). LaSalle Invariance Principle: For the system (11), assume that a compact and positively invariant set D and a function V : D → R, such that V ≤ 0 along the trajectories of the system exist. Let G * = {x ∈ D : V (x) = 0} and let G be the maximum invariant set contained in D.
This result is important in the stability analysis of dynamical systems. It provides a tool to prove global asymptotic stability when V ≤ 0. If x * is locally stable and G is positively invariant, i.e. all solutions starting in G remain in G, then x * is globally asymptotically stable in the region G.
This method applies when the equilibrium point is located at the origin, but if x * = 0, the change of variableŝ x = x − x * can be made. In this way, the system (11) becomesx = g(x), where g(x) = f (x + x * ), and the equilibrium x * shifts to the origin, i.e. g(0) = f (x * ) = 0. At this point, the method could be applied.
The Lyapunov functions method can be used to establish global stability results for epidemic models [53]. However, it is often difficult to be applied, since there is no general technique to build such Lyapunov functions. One of the general forms of Lyapunov functions used in the literature is where the coefficients c i must be determined in such a way that the derivative of V along the system solutions is nonnegative, a rather challenging task in higher dimensional models. Similar techniques have been recently found out to deal with discrete-time systems. One can check these results in detail in [73]. If the disease model possesses a network structure, in particular if it can be considered as a coupled system in which each subsystem has a Lyapunov function, then the network theory method can be used to construct a Lyapunov function for the whole system. This method has been successfully applied to establish the overall stability of the endemic equilibrium point for several complex disease models [115]. However, the same problem faced above for higher order systems subsists here as well, since a systematic approach is lacking to guide in the construction of Lyapunov functions for each subsystem.

Bifurcation and existence of periodic solutions.
In this section, we address the question of how the qualitative behavior of (11) changes as a function of parameters associated with the description of the vector field f . Thus, from now on, we will denote (11) by x = f (x, µ), where µ = (µ 1 , . . . , µ m ) ∈ R m represents such parameters. If the qualitative behavior remains the same for all nearby vector fields, the system (family) given by x = f (x, µ) or the vector field (family) f (x, µ), with µ ∈ R m , is said to be structurally stable [102].
If a particular vector field f (·, µ) ∈ C (D) is not structurally stable, it belongs to the bifurcation set in C (D). The qualitative structure of the solution set or global phase portrait of (11) changes as the vector field f crosses a critical point in the bifurcation set. In particular, we study bifurcations at non-hyperbolic equilibrium points (see [18,19,85]). Simple bifurcation types occur when there is a variation in the number and stability of the equilibria as µ changes.
In epidemic models, bifurcation analysis is used to establish the feasible epidemic evolutions. The most common situation is a transcritical bifurcation at the parameter value R 0 = 1, with the appearance of an asymptotically stable endemic equilibrium with an infected individuals amount that continuously depends on R 0 [38,55,69,91]. Here, a change of stability occurs between the disease-free equilibrium, which is asymptotically stable for R 0 < 1 and unstable for R 0 > 1, and the endemic equilibrium, which, in contrast, shows the opposite stability behavior, although it is biologically significant only for R 0 > 1.
Another simple way in which non-linearity can lead to periodic behavior is the Hopf bifurcation [85,124]. The technique consists in looking for parameter thresholds below which an equilibrium is asymptotically stable, and above which it becomes unstable and a periodic solution emerges. The real part of an eigenvalue of the Jacobian must vanish, and originate a pair of purely imaginary eigenvalues. If all the remaining eigenvalues are negative or have negative real parts, the periodic solution that arises is stable, while, conversely, it will be a repelling limit cycle.
The bifurcation behavior can be described graphically. It contains the size of the dependent variable at equilibrium as a function of one of the bifurcation parameters. For epidemic models, the bifurcation curve is the graph that plots the values of the infectious population I at the equilibria in terms of the parameter given by the basic reproductive number R 0 .
In some epidemic models with multiple sets and asymmetry between sets or mechanisms of multiple interactions, it is possible to have different bifurcation behaviors at R 0 = 1 [32]. Also, there can be several positive endemic equilibria for values of R 0 < 1 [80] and a subcritical (also called backward) bifurcation at R 0 = 1 [24].
The qualitative behavior of a system with a subcritical (backward) bifurcation differs from the one of a system with a supercritical (forward) bifurcation [62]. The epidemiological consequence of a subcritical bifurcation lies in the fact that the condition R 0 < 1 is necessary, but not any longer sufficient, to ensure disease eradication. Therefore, the presence of this phenomenon in the disease makes the control of propagation difficult [62]. Figure  6 shows the bifurcation diagrams for supercritical and subcritical bifurcations, respectively [62]. An epidemic model has periodic oscillations if and only if the model is cyclical and the disease has temporary immunity by which the return to susceptibility of an individual can be significantly delayed [62]. Several epidemiological mechanisms may lead to periodic solutions. The most direct way in which the periodicity arises is through the contact rate, but periodicity can also emerge autonomously [76]. Infectious disease models with non-linear incidences of certain general forms may have periodic solutions. Some models with variable population size and disease-related deaths have periodic solutions. Most of them are vector-host models, where the vector's life is much shorter than the one of the host [50].
Recently, periodic solutions have numerically been found in age-structured models with cross-immunity between two viral strains [66]. For more details, see also [51,91].

Parameter estimation and simulation of an epidemic model
In general, after a qualitative study of an epidemic model a quantitative study is needed as well, to confront the theoretical results with the real data and to validate the model. To this end, numerical simulations are performed to observe the short-term behavior of the disease spread in the host population. Classical adaptive numerical methods for solving ordinary differential equations can be used, such as the Runge-Kutta-Fehlberg method, or other routines implemented for instance in MATLAB, e.g. ode45 or ode15s for stiff problems [11,12,59,104,121].
To run simulations, parameter values which have to be estimated directly from epidemiological data are used. In a mathematical model, it is important to investigate the relevance of each parameter for the system outcome. The objective of sensitivity analysis is to assess the most influential parameters on the dynamical evolution of the system [91]. In epidemiological models, this is especially relevant for the disease transmission dynamics. The parameters of an epidemic model, such as the transmission coefficient in the incidence term and the basic reproduction number R 0 , can be estimated from data points collected at successive moments during an outbreak. If reliable annual prevalence data are available for a specific disease, an estimation of these parameters can be obtained by finding the best fit of the model outcome to the available data, using standard optimization techniques. They usually minimize the sum of the square errors (SSE); for this purpose several Matlab or Mathematica routines can be used [91]. An alternative approach to solve the least-squares problem is to use a genetic algorithm [59].
For outbreaks of an emerging disease, the main disadvantage of estimating parameters is the long enough time that must be elapsed for collecting data. Also, stochastic effects are important, and need to be taken into account particularly at the start of an outbreak [45,91].
6 Implementation of prevention and control strategies in compartmental models.
Once the dynamics of a model is known, one can search for control strategies that are relevant for each specific epidemic.
One of the most effective ways to prevent and control a disease is to reduce contacts. However, the emergence and re-emergence of infectious diseases, climate change, population growth, unplanned urbanization and the ability of humans to migrate on a large scale have increased. Therefore, it is not easy to reduce interactions among individuals. Vaccines and medical treatments are the two widely used prevention tools that can potentially reduce transmission of a disease and control its spread. Prevention and control measures are used to keep the number of infected individuals as low as possible. Other measures may include isolation or quarantine, education, and biological control. They lead to a reduction in the number of medical treatments, hospitalizations and absences from work due to illness, so contributing to a reduction of total costs, which is the main goal of adopting alternative measures.
Mathematical models can help to quantify the benefits of possible prevention scenarios. Actually, it results as necessary to accurately assess the benefits of these programs. Different mathematical approaches that include vaccination are discussed for different infections in [21]. For the human papilloma virus, vaccination is shown to be effective in reducing the number of infected and the number of fatalities [117].
Another form of intervention is through education-induced behavioral changes. They can alter the dynamics of an outbreak of a disease. For instance, for the Ebola virus in Sudan, a deterministic model for a community partitioned into two types of individuals, those that are educated about Ebola and take precautions to avoid contracting the disease, and other ones not taking precautions, shows the importance of education as a preventive measure [87]. A transfer diagram for the inclusion of two compartments representing the intervention by education in an SIR model is shown in Figure (7). The effects of public health education on the evolution of a disease have also been studied for HIV and the Hepatitis C virus transmission [87]. Isolation and quarantine were proved as possible strategies to stop the spread of SARS in Toronto [37]. Early treatments to significantly reduce the transmission of new HIV infections are discussed in [104]. A mathematical model for the prevention of the avian influenza pandemic with three controls (namely, isolation of infected individuals, antiviral treatment of infected individuals, and elimination of infected birds) is presented in [123]. As we mentioned in the previous sections, R 0 allows us to determine the severity of a disease and to study the qualitative behavior of the epidemic model. It can be quite easy to estimate the average duration of the infection and the initial number of infectious individuals in order to produce an estimation of R 0 .
Different prevention and control strategies can be implemented, e.g. vaccination at the beginning of an outbreak. If a perfect vaccine is available for a specific disease and a fraction p of the population is vaccinated, then, in an SIR model like (2), the term pµN must be included in the R class and, consequently, it disappears from the susceptibles. That is, an SIR model with vaccination is as given in (16): In this context, the disease will not spread if (1 − p)R 0 < 1. This indicates that if a fraction p > 1 − R −1 0 of the population is vaccinated, herd immunity is attained. This simple model predicts that diseases with large values of R 0 require a large proportion of the population successfully vaccinated in order to be eradicated. Smallpox is the only human disease that currently has been eradicated, due to vaccination campaigns [91].
The different mathematical approaches useful for capturing a diversity of early epidemic growth profiles, obtained from empirical data, ranging from sub-exponential to exponential growth dynamics, are presented in [40]. Early epidemic forecasts consisting of the likely short-term evolution of an ongoing outbreak can help guide the type and intensity of interventions, including diagnostic needs of health infrastructure, isolation of infected individuals, and contact tracing activities.
Recently, optimal control theory has been employed to find the best combination of the different diseasecurbing alternatives in such a way that the control, the economic and the social costs are minimized. Generally, the optimal control problems for real life applications found in the literature can only be solved numerically [91]. Using the Pontryagin Maximum Principle, the necessary optimality conditions for the existence are usually derived (see [113]).
The prevention of the avian influenza pandemic has been studied by adjusting three control functions in the human-to-human transmission model to minimize the impact of the disease. They consist of isolation of infected individuals, antiviral treatment of infected individuals, and elimination of infected birds [123]. A new deterministic model to assess the population-level impact by quarantining individuals suspected of being exposed to the disease during the ebola outbreaks in 2014-2015 is presented in [43]. Early treatments to significantly reduce the transmission of new HIV infections are introduced in [104]. The optimal policy for distributing chlorine tablets for water purification during a cholera outbreak to assess the fraction of susceptible individuals who should have access to them in order to minimize the total number of new infections is investigated in [86], with an application to real data of the cholera outbreak in Yemen in 2017-2018.
A preventive control for ebola in the form of an education campaign and two treatment controls applied to an infected and advanced stage of the human population are shown to lead to a reduction of the ebola infection, giving a good framework for planning and designing cost-effective strategies for good interventions in the treatment of this disease [23]. A mathematical model with optimal strategies for dengue reduction and prevention in Cali, Colombia is presented in [114]. A new version of an optimal control problem for a vaccination strategy with two dengue serotypes appears in [41].
Another type of control strategy that can be implemented is the feedback control [39,42]. The results show that, when appropriately choosing the feedback control variables values, the endemic disease may be eradicated. Alternatively, it may remain endemic, but the level of endemicity can be kept at a low value.
All these measures, including population-level immunity as well, taken jointly or separately, can provide useful information for the proper implementation of effective treatment programs for a pandemic.

Conclusions
In this brief review paper, different compartmental paradigms that are generally used in mathematical epidemiology were discussed. Biological aspects used to make a qualitative analysis of the model were also described to support the relevance of this features for understanding the dynamics of diseases modeled by ordinary differential equations.
The vast majority of the contributions in the literature for the spread of infectious diseases that have been discussed are of a deterministic nature, continuous in time, and modeled by ordinary non-linear differential equations or by difference equations. Furthermore, these models do not simulate the individuals dynamics, for which individual-based or network (stochastic) models would be more appropriate. Nonetheless, the latter are more complex to be analyzed. Stochastic approach is of importance for modeling the probability of contact with the infected as well as other relevant biological features such as the immune system of each host, and so on.
Mathematical models may effectively help in assessing various biological features to implement or suppress heterogeneities, control measures and environment conditions, that are key links in the dynamics of disease transmission models. Consolidating our understanding of the evolution of different modeling approaches helps in the development of new models, by incorporating different biological aspects, and finally contributes to improve predictions and to take more realistic control measures. They will require the development of innovative and manageable mathematical and statistical tools and techniques for their overcoming. Contact networks, spatial heterogeneity and early epidemic growth profiles represent some of these shortcomings at the present state of knowledge.
Collection and use of large empirical data sets for parameter estimation is another relevant current challenge. Efforts to improve the modeling of disease transmission patterns will provide a better understanding and, as a consequence, better prevention strategies, having a positive impact on the calibration of the intervention and contingency plans, in order to face the threats of the emergence of new infectious diseases, and the reemergence of old ones.