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

In compartmental approaches, the population (N) is partitioned into several mutually exclusive subsets or classes, which exhibit the same properties with respect to the disease (disease states of the population), e.g. all individuals that can contract the disease, called Susceptibles (S); those which have been exposed to the pathogen but are not yet contagious, named Exposed (E); the ones which are able to transmit the disease, denominated Infectious (I); those which have overcome the disease, called Recovered (R); the ones which show symptoms and are therefore isolated, named Quarantined (Q); those which have received a vaccine, called Vaccinated (V); and so on. These subsets represent the various compartments, and the transference of individuals from one compartment to another one occurs at transition rates that are possibly estimated from empirical data. In general, they depend on the size of the compartment from which the migration occurs. The state variables which translate the compartments description represent the number of individuals in each compartment. A state of the system at any time is the vector whose components are all these variables. On the other hand, the transition rates are translated as function of parameters.

Concerning the class E, the period of time from the moment at which the host acquires the infection until the moment when the host becomes infectious, i.e. when it is able to transmit the pathogenic agent to another individual is called latency. On the other hand, the interval of time between the moment of the exposure to an infectious agent and the time at which the disease symptoms begin is said to be the incubation period. That is, the incubation period is the time required for the infectious agent to multiply sufficiently to produce symptoms or evidence of the infection. The incubation and latency periods do not necessarily match. For example, for the influenza, people become infectious approximately one day before the symptoms appear.

The most common approaches frequently found in the literature are as follows: SI (Susceptible–Infected), SIS (Susceptible–Infected–Susceptible), SIR (Susceptible–Infected–Recovered), SIRS (Susceptible–Infected–Recovered–Susceptible), SEIR (Susceptible–Exposed–Infected–Recovered), SEIRS (Susceptible–Exposed–Infected–Recovered–Susceptible), together with other variants, such as SEIRV, SEIQR [91], etc. For life lasting infections, which describe viral diseases like HIV/SIDA, the usual approach is of type SI. The disease for which a recovered individual acquires lifelong immunity, like measles, rubella, mumps or smallpox are treated via an SIR approach. Diseases caused by bacterial agents (meningococcal meningitis, pest, sexually transmitted diseases, etc.) and protozoa (malaria, zika, dengue, etc.) or even a simple cold, for which no immunity after the infection exists, are dealt with an SIS approach [4, 65, 91].

A disease outbreak arises when the number of infected individuals in any part of the population exceeds the average one in a short period of time. An epidemic is a contagious disease spreading rapidly in (all) the population or community, producing a large number of infected individuals during a short period of time. An epidemic can collapse a health system, when the propagation of the infection is uncontrolled, as it happened for instance with the outbreak of ebola in West Africa in the year 2014 (WHO). If the disease remains long time in the population, it is said to be endemic. If the outbreak affects large geographic areas, such as continents, then the disease becomes a pandemic; as in the case of HIV or the current Covid-19. The prevalence of the disease represents the number of infected individuals at time t whereas the incidence represents the rate at which new infections occur. Generally, it is assumed that the new cases of infected individuals are generated by direct contacts through homogeneous mixing, meaning that all the individuals in a set have the same probability to be infected [20]. This is known as mass action law (or frequency-dependent) incidence. Mathematically, it is expressed as βIS, where β is the transmission coefficient. This rate can be applied to diseases like covid19, influenza or tuberculosis, for which the direct (horizontal) transmission occurs through the droplets released in the air, mainly when someone coughs [65], or also to diseases (vertically) transmitted by the placenta of a mother to her child before or at the moment of the birth such as HIV, hepatitis B, and syphilis [32]. For other diseases, for instance those sexually transmitted, the contact number cannot vary with the density of the population. This is modeled by the standard incidence (or density-dependent) rate βIS/N [65, 79]. For further information, refer to [4, 79, 84, 91]).

The recovery rate is the proportion at which the diseased individuals are transferred into the compartment 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 figure 1 represents one of the first compartmental models of direct transmission. It was proposed by Kermack-MacKendrick [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].

Fig. 1

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 non-homogeneous 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.

For instance, one can find in the recent literature models with spatial heterogeneity [14, 22, 49, 88, 97, 103, 105], age-structured models, [4, 7, 15, 31, 33, 67, 74, 91, 117, 119], or seasonal models [17, 76, 79].

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_hS_hI_M/N_M and T_MS_MI_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_hI_h [28].

Fig. 2

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.

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.

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: (1) {S′(t)=−βS(t)I(t), S(0)=S0≥0,I′(t)=βS(t)I(t)−γI(t), I(0)=I0≥0,R′(t)=γI(t), R(0)=R0≥0, \left\{ {\matrix{ {S'(t) = - \beta S(t)I(t),\quad S(0) = {S_0} \ge 0,} \hfill \cr {I'(t) = \beta S(t)I(t) - \gamma I(t),\quad I(0) = {I_0} \ge 0,} \hfill \cr {R'(t) = \gamma I(t),\quad R(0) = {R_0} \ge 0,} \hfill \cr } } \right.

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 ℛ₀, 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 ℛ₀ < 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, ℛ₀ > 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 ℛ₀ depends on the parameters involved in the system: for the model (1), it is ℛ₀ = β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₀, during its average infectious period, 1/γ. If ℛ₀ ≤ 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 ℛ₀ > 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 = γβ⁻¹. This maximum is given by Imax=γβ[lnγβS(0)−1]+I(0)+S(0). {I_{max}} = {\gamma \over \beta }\left[ {\ln {\gamma \over {\beta S(0)}} - 1} \right] + I(0) + S(0). 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₀ = (β/γ)S₀ = ℛ₀. 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: (2) {S′(t)=μN−βS(t)I(t)−μS(t),I′(t)=βS(t)I(t)−γI(t)−μI(t),R′(t)=γI(t)−μR(t), \left\{ {\matrix{ {S'(t) = \mu N - \beta S(t)I(t) - \mu S(t),} \hfill \cr {I'(t) = \beta S(t)I(t) - \gamma I(t) - \mu I(t),} \hfill \cr {R'(t) = \gamma I(t) - \mu R(t),} \hfill \cr } } \right. In such a case, ℛ₀ = β (γ + μ)⁻¹. Thus, the value of ℛ₀ 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, while b˜>0 \tilde b > 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: (3) {S′(t)=b(S(t)+R(t))+b˜pI(t)−βS(t)I(t)−bS(t),I′(t)=βS(t)I(t)+b˜(1−p)I(t)−γI(t)−b˜I(t),R′(t)=γI(t)−bR(t), \left\{ {\matrix{ {S'(t) = b(S(t) + R(t)) + \tilde bpI(t) - \beta S(t)I(t) - bS(t),} \hfill \cr {I'(t) = \beta S(t)I(t) + \tilde b(1 - p)I(t) - \gamma I(t) - \tilde bI(t),} \hfill \cr {R'(t) = \gamma I(t) - bR(t),} \hfill \cr } } \right. Here note that from the second equation in (3), the basic reproduction number can be written in two different ways: ℛ˜0=βN+b˜qγ+b˜=βNγ+b˜p. {\tilde {\cal R}_0} = {{\beta N + \tilde bq} \over {\gamma + \tilde b}} = {{\beta N} \over {\gamma + \tilde bp}}. The first formulation shows that the infection disappears, namely ℛ˜0<1 {\tilde {\cal R}_0} < 1 , when the number of new infections, obtained either by horizontal βN or vertical b˜q \tilde bq transmission, per average length of the infectious period, 1/(γ+b˜) 1/(\gamma + \tilde b) , is less than one. The second one indicates that the vertical transmission has the effect of reducing ℛ₀, namely ℛ˜0=ℛ0γγ+b˜p≤ℛ0. {\tilde {\cal R}_0} = {{\cal R}_0}{\gamma \over {\gamma + \tilde bp}} \le {{\cal R}_0}. 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 β=β(I)=Ime−αI. \beta = \beta (I) = {I^m}{e^{ - \alpha I}}.

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.

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 continuous-time 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 (4) {S¯B(t+1)=S¯B(t)+(μB+αB)C¯B(t)−βBS¯B(t)I¯T(t)NT(t),I¯B(t+1)=I¯B(t)+βBS¯B(t)I¯T(t)NT(t)−λBI¯B(t),C¯B(t+1)=C¯B(t)+λBI¯B(t)−(μB+αB)C¯B(t),S¯T(t+1)=S¯T(t)+μTpI¯T−βTS¯T(t)I¯B(t)NB(t),I¯T(t+1)=I¯T(t)+βTS¯T(t)I¯B(t)NB(t)−μTpI¯T(t), \left\{ {\matrix{ {{{\bar S}_B}(t + 1) = {{\bar S}_B}(t) + \left( {{\mu _B} + {\alpha _B}} \right){{\bar C}_B}(t) - {\beta _B}{{\bar S}_B}(t){{{{\bar I}_T}(t)} \over {{N_T}(t)}},} \hfill \cr {{{\bar I}_B}(t + 1) = {{\bar I}_B}(t) + {\beta _B}{{\bar S}_B}(t){{{{\bar I}_T}(t)} \over {{N_T}(t)}} - {\lambda _B}{{\bar I}_B}(t),} \hfill \cr {{{\bar C}_B}(t + 1) = {{\bar C}_B}(t) + {\lambda _B}{{\bar I}_B}(t) - \left( {{\mu _B} + {\alpha _B}} \right){{\bar C}_B}(t),} \hfill \cr {{{\bar S}_T}(t + 1) = {{\bar S}_T}(t) + {\mu _T}p{{\bar I}_T} - {\beta _T}{{\bar S}_T}(t){{{{\bar I}_B}(t)} \over {{N_B}(t)}},} \hfill \cr {{{\bar I}_T}(t + 1) = {{\bar I}_T}(t) + {\beta _T}{{\bar S}_T}(t){{{{\bar I}_B}(t)} \over {{N_B}(t)}} - {\mu _T}p{{\bar I}_T}(t),} \hfill \cr } } \right. and it is epidemiologically meaningful if and only if the following parametric restrictions are satisfied: 1−(μB+αB)≥0, 1−βB≥0, 1−λB≥0, 1−βT≥0, 1−μTp≥0. 1 - ({\mu _B} + {\alpha _B}) \ge 0,\quad 1 - {\beta _B} \ge 0,\quad 1 - {\lambda _B} \ge 0,\quad 1 - {\beta _T} \ge 0,\quad 1 - {\mu _T}p \ge 0.

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].

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].

Fig. 3

Fig. 4

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]. (9) {SB′(t)=μB(SB(t)+RB)+αBRB(t)−μBSB(t)−βBSB(t)IT(t)NT(t),IB′(t)=μBIB+βBSB(t)IT(t)NT(t)−μBIB−λBIB(t),RB′(t)=λBIB(t)−(μB+αB)RB(t),ST′(t)=μT(ST(t)+pIT(t))−βTST(t)IB(t)NB(t)−μTST(t),IT′(t)=βTST(t)IB(t)NB(t)+(1−p)μTIT(t)−μTIT(t), \left\{ {\matrix{ {{S_B^\prime}(t) = {\mu _B}({S_B}(t) + {R_B}) + {\alpha _B}{R_B}(t) - {\mu _B}{S_B}(t) - {\beta _B}{S_B}(t){{{I_T}(t)} \over {{N_T}(t)}},} \hfill \hfill \cr {{I_B^\prime}(t) = {\mu _B}{I_B} + {\beta _B}{S_B}(t){{{I_T}(t)} \over {{N_T}(t)}} - {\mu _B}{I_B} - {\lambda _B}{I_B}(t),} \hfill \hfill \cr {{R_B^\prime}(t) = {\lambda _B}{I_B}(t) - ({\mu _B} + {\alpha _B}){R_B}(t),} \hfill \hfill \cr {{S_T^\prime}(t) = {\mu _T}({S_T}(t) + p{I_T}(t)) - {\beta _T}{S_T}(t){{{I_B}(t)} \over {{N_B}(t)}} - {\mu _T}{S_T}(t),} \hfill \hfill \cr {{I_T^\prime}(t) = {\beta _T}{S_T}(t){{{I_B}(t)} \over {{N_B}(t)}} + (1 - p){\mu _T}{I_T}(t) - {\mu _T}{I_T}(t),} \hfill \hfill \cr} } \right.

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.

A solution for the system (11) is a function ϕ = ϕ(t) of class 𝒞 defined in a range J ⊂ ℝ, ϕ : J → ℝⁿ such that ϕ′(t) = f (ϕ(t)). Since in epidemic models the dependent variables S(t), E(t), I(t), R(t) represent populations, they must be nonnegative for all t ≥ 0, given any initial condition S(t₀),E(t₀),I(t₀),R(t₀) ≥ 0 for the system (11). Let Ω=ℝ+n \Omega = {\mathbb R}_ + ^n be the first orthant given by Ω={(x1,x2,…,xn):x1≥0,x2≥0,…,xi≥0,…,…,xn≥0}. \Omega = \left\{ {({x_1},{x_2}, \ldots ,{x_n}):{x_1} \ge 0,{x_2} \ge 0, \ldots ,{x_i} \ge 0, \ldots , \ldots ,{x_n} \ge 0} \right\}. Then, any solution ϕ(t) starting in Ω should remain in Ω for all t ≥ 0.

The positive invariance of the non-negative orthant ℝ+n {\mathbb R}_ + ^n for the system (11) ensures that no trajectory escapes from ℝ+n {\mathbb R}_ + ^n by crossing one of its faces. Assume that, on the face x_i = 0, the vector field f in (11) has components (f1(x1,x2,…,0,…,xn),f2(x1,x2,…,0,…,xn),…,fn(x1,x2,…,0,…,xn)). ({f_1}({x_1},{x_2}, \ldots ,0, \ldots ,{x_n}),{f_2}({x_1},{x_2}, \ldots ,0, \ldots ,{x_n}), \ldots ,{f_n}({x_1},{x_2}, \ldots ,0, \ldots ,{x_n})). This defines a vector at each point (x₁,x₂,...,0,...,x_n) of this face. A path of the solution must have this direction, and if the vector points inward to ℝ+n {\mathbb R}_ + ^n , the path can not leave ℝ+n {\mathbb R}_ + ^n . Therefore, it is sufficient to show that the vector field points towards the interior of the boundary of ℝ+n {\mathbb R}_ + ^n , and, by the existence and uniqueness theorem of solutions [68], trajectories cannot intersect each other.

More generally, let Ω be a region in ℝ+n {\mathbb 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).

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₀ 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 (12) y′=J(x*)y. y' = J({x^*})y. 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 (13) det[J(x*)−λI]=0. \det [J({x^*}) - \lambda I] = 0. 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 (14) tr(J(x*))<0, det(J(x*))>0; {\rm tr}(J({x^*})) < 0,\quad \det (J({x^*})) > 0; for the case n = 3, they become (15) tr(J(x*))>0, det(J(x*))<0, tr(J(x*))M2>det(J(x*)), {\rm tr}(J({x^*})) > 0,\quad \det (J({x^*})) < 0,\quad {\rm tr}(J({x^*})){M_2} > \det (J({x^*})), where M₂ 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 ℛ₀, 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]. ℛ₀ 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 ℛ₀ > 1 an epidemic outbreak will occur, while if ℛ₀ < 1 the disease will become extinguished. If ℛ₀ = 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 ℛ₀ 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. ℛ₀ 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 [48, 65].

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 ℝⁿ, this solution is represented by a closed curve because the path starting at x(t₀) 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₁(t),x₂(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₀) 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 ⊂ ℝ². Then, D does not contain periodic orbits if ∇⋅h(x)<0 (or>0), x∈D, h(x)=(gf)(x). \nabla \cdot h(x) < 0\;(or > 0),\quad x \in D,\quad h(x) = (gf)(x).

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 → ℝ such that V(0)=0, V(x)>0, x∈D, x≠0. V(0) = 0,\quad V(x) > 0,\quad x \in D,\quad x \ne 0. Then, a)

If V′(x) ≤ 0,∀x ∈ D, x = 0 is stable;

Let G ⊂ ℝⁿ an open set. A function V (x) is said to be a Lyapunov function with respect to G if V′(x)≤0, ∀x∈G. V'(x) \le 0,\quad \forall x \in G. 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₀) 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 → ℝ, 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. If x₀ ∈ D, then x(t) → G as t → ∞.

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 variables x^=x−x* \hat x = x - {x^*} can be made. In this way, the system (11) becomes x^′=g(x^), \hat x' = g(\hat x), , where g(x^)=f(x^+x*) g(\hat x) = f(\hat 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 V=∑inci[ xi−xi*−xi*ln(xixi*) ], V = \sum\limits_i^n {c_i}\left[ {{x_i} - x_i^* - x_i^*\ln \left( {{{{x_i}} \over {x_i^*}}} \right)} \right], 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.

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 μ = (μ₁,..., μ_m) ∈ ℝ^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 μ ∈ ℝ^m, is said to be structurally stable [102].

If a particular vector field f (·, μ) ∈ 𝒞 (D) is not structurally stable, it belongs to the bifurcation set in 𝒞 (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 ℛ₀ = 1, with the appearance of an asymptotically stable endemic equilibrium with an infected individuals amount that continuously depends on ℛ₀ [38, 55, 69, 91]. Here, a change of stability occurs between the disease-free equilibrium, which is asymptotically stable for ℛ₀ < 1 and unstable for ℛ₀ > 1, and the endemic equilibrium, which, in contrast, shows the opposite stability behavior, although it is biologically significant only for ℛ₀ > 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 ℛ₀.

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 ℛ₀ = 1 [32]. Also, there can be several positive endemic equilibria for values of ℛ₀ < 1 [80] and a subcritical (also called backward) bifurcation at ℛ₀ = 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 ℛ₀ < 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].

Fig. 6

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].