key: cord-0119428-ljph6cao
authors: Dadlani, Aresh; Afolabi, Richard O.; Jung, Hyoyoung; Sohraby, Khosrow; Kim, Kiseon
title: Deterministic Models in Epidemiology: From Modeling to Implementation
date: 2020-04-07
journal: nan
DOI: nan
sha: d1a899fc4738bb5c248ada458e038537633bc667
doc_id: 119428
cord_uid: ljph6cao

The abrupt outbreak and transmission of biological diseases has always been a long-time concern of humankind. For long, mathematical modeling has served as a simple and yet efficient tool to investigate, predict, and control spread of communicable diseases through individuals. A myriad of works on epidemic models and their variants have been reported in the literature. For better prediction of the dynamics of a particular disease, it is important to adopt the most suitable model. In this paper, we study some of the widely-appreciated deterministic epidemic models in which the population is divided into compartments based on the health status of each individual. In particular, we provide a demographic classification of such models and study each of them in terms of mathematical formulation, near equilibrium point stability properties, and disease outbreak threshold conditions (basic reproduction ratio). Furthermore, we discuss the various influential factors that need to be considered during epidemic modeling. The main objective of this article is to provide a basic understanding of the mathematical complexity incurred in deterministic epidemic models with the aid of graphical illustrations obtained through implementation.

Throughout recorded history, human population has always been haunted by the emergence and re-emergence of infectious diseases. Several lives have been lost due to lack of knowledge on the dynamical behavior of epidemic outbreaks of contagious diseases and measures to confront them [1] , [2] . For decades, scientists have toiled to understand the transmission characteristics of such diseases so as to devise control strategies to prevent further spread of infection. The field of science that studies such epidemic diseases and in particular, the factors that influence the incidence, distribution, and control of infectious diseases in human populations is called epidemiology. In this regard, mathematical modeling has proven to serve useful in analyzing, predicting, evaluating, detecting, and implementing efficient control programs. Such analytical models accompanied by computer simulations serve as experimental tools for building, testing, and assessing theories and understanding the relationship between various parametric values involved. Practical use of epidemic models depends on how closely they realize actual biological diseases in real world. To keep the models simple and tractable, many assumptions and relaxations are taken into consideration at each level of the process. However, even such simplified models often pose significant questions regarding the underlying mechanisms of infection spread and possible control approaches. Hence, adopting the apt epidemic model for prediction of real phenomenon is of great importance.

Models that are useful in the study of infectious diseases at the population scale can be broadly classified into two types: deterministic and stochastic. Early models that were developed to study specific diseases such as tuberculosis, measles, and rubella were deterministic in nature. In deterministic models, the large population is divided into smaller groups called compartments (or classes) where each group represents a specific stage of the epidemic. Such models, often formulated in terms of a system of differential equations (in continuous time) or difference equations (in discrete time), attempt to explain what happens on the average at the population scale. A solution of a deterministic model is a function of time or space and is generally uniquely dependent on the initial data. On the other hand, a stochastic model is formulated in terms of a stochastic process which, in turn, is a set of random variables, X(t; ω) ≡ X(t), defined as {X(t; ω)|t ∈ T and ω ∈ Ω} where T and Ω represent time and a common sample space, respectively. The solution of a stochastic model is a probability distribution for each of the random variables. Such models capture the variability inherent due to demographic and environment variability and are useful under small population sizes. More specifically, they allow follow-up of each individual in the population on a chance basis [3] , [4] . Discrete-time Markov chain (DTMC), continuous-time Markov chain (CTMC), and stochastic differential equation (SDE) models are three types of stochastic modeling processes which have been deeply covered in [5] . Figure 1 shows the different classes under which epidemic models have been studied in the literature. The connecting blue lines in the figure highlight the main scope of this report.

Needless to say, both deterministic and stochastic epidemic models have their own applications. Deterministic models are used to address questions such as: what fraction of individuals would be infected in an epidemic outbreak?, what conditions should be satisfied to prevent and control an epidemic?, what happens if individuals are mixed non-homogeneously?, and so on [6] . While such models are preferable in studying a large population, stochastic epidemic models are useful for a small community and answer questions such as: how long is the disease likely to persist?, what is the probability of a major outbreak?, and the like [5] . Hence, stochastic epidemic models are generalized forms of simple deterministic counterparts. However, unlike deterministic models, stochastic mod- Figure 1 : Classification of various classes of epidemic models. els can be laborious to set up and may need many simulation runs to yield useful predictions. They can become mathematically very complex and result in misperception of the dynamics [3] . To this end, we focus on some widely-used deterministic models which are relatively easier to conceive, set up, and implement using various computer softwares at disposal. Deterministic epidemiology is believed to have started in the early twentieth century [7] . In 1906, Hamer was the first to assume that the incidence (number of new cases per unit time) is proportional to the product of the number of susceptible and infective individuals in his model for measles epidemics [8] . The exponential growth in mathematical epidemiology was boosted by the acclaimed work of Kermack and McKendrick which was published in 1927 [9] . This paper laid out a foundation for modeling infections where all members of the population are assumed to be initially equally susceptible to the disease and confer complete immunity only after recovery. After decades of neglect, the Kermack-McKendrick model was brought back to prominence by Anderson et al. in 1979 [10] . Since then several models have been developed addressing aspects such as passive and disease-acquired immunity, vaccination, quarantine, vertical transmission, disease vectors, age structure, social and sexual mixing groups, as well as chemotherapy [11, 12, 13, 14] . Improved models have also been designed for diseases such as measles, chickenpox, smallpox, whooping cough, malaria, rabies, diphtheria, filariasis, herpes, syphilis, and HIV/AIDS [14] , [15] .

The main objective of this article is to help the reader gain an insight on the basics of deterministic compartmental modeling through implementation. In this work, we formulate some well-known models and derive their steady-state solutions. Since the models under study are non-linear in nature, we investigate their qualitative behavior near their corresponding equilibria using linearization method [16] . All models discussed in this paper have been implemented using Wolfram Mathematica [17] , the codes of which are freely available.

The remainder of this article is structured as follows. Section II provides a demographic classification of deterministic models along with notations and assumptions that will be used throughout the paper. In Section III, we present the classical epidemic model known as the susceptible-infected-recovered (SIR) model which forms the basis for the extended models that follow in Sections IV to VII, accompanied by their implementation results. In Section VIII, we present some additional factors that impact the behavior of epidemic models, followed by some conclusive remarks in Section IX.

In order to analyze the structure of epidemic models as well as the relation between their structure and the resulting dynamics, it is important to classify models as clear and simple as possible. In our study, we classify and study compartmental models based upon demography or vital dynamics. Demography relates to the study of characteristics of human populations such as birth, death, incidence of disease, and so on. Epidemic models with vital dynamics consider an open population with births and deaths while models without vital dynamics have a closed and fixed population with no demographic turnover.

For better realization, we assume that the law of mass action holds for all models in this paper. This law states that if individuals in a population mix homogeneously, then the encounters between infected and susceptible individuals occur at a rate proportional to their respective numbers in the population [18] . In other words, the rate at which the susceptible population becomes infected is directly proportional to the product of the sizes of the two populations. Table 1 summarizes the notations that will be used in model derivation henceforth. Note that S , I, R, and E are used to represent the compartments in the epidemic model as well as the proportion of the corresponding compartments at any time instant t. 3 The basic SIR model

In their first paper, Kermack and McKendrick created a model in which the population is divided into three compartments: susceptible (S ), infected (I), and recovered (R) [9] as illustrated in Figure 2 . They assumed that all individuals are mutually equally susceptible to the disease and that complete immunity is obtained merely after recovery from infection. Moreover, they also assumed that the duration of the disease is same as the duration of infection with constant transmission and recovery rates. Based upon the demographic classification, the epidemic and endemic S IR models are studied below.

For a closed population of size N, we assume that the mixing of individuals is homogeneous and the law of mass action holds. Also, for large classes of communicable diseases, it is more realistic to consider a force of infection that depends on the fraction of infected population with respect to the total constant population N, rather than the absolute number of infectious subjects. Based upon this assumption, the standard disease incidence rate is defined as βS I/N and the overall rate of recovery is given as γI. In spite of the above simplifying assumptions, the resulting non-linear system does not admit a closed-form solution. Nevertheless, we shall see how significant results can be derived analytically. Figure 2 can be translated into the following set of differential equations:

Summing up (1), (2), and (3) yields zero which implies that the population is of constant size with S + I + R = N. Dividing (2) by (1) gives:

Assume that the population is susceptible up to time zero at which a relatively small number, I(0), become infected. Thus, at t = 0, S (0) = N − I(0) and R(0) = 0. As time approaches infinity, lim t→∞ I(t) = 0, lim t→∞ S (t) = S (∞), and the number of individuals that have been infected is S (0) − S (∞). Integrating (4) leads to:

where c is constant and S (∞) is the proportion of susceptibles at the end of the epidemic.

Since the initial infection is small, (5) further reduces to:

(a) R 0 = 0.7, β = 0.7, and γ = 1. where c ′ denotes some constant. Defining R 0 = β/γ as the basic reproduction ratio, we see in Figure 3 that for R 0 > 1, a small infection to the population would create an epidemic. R 0 describes the total number of secondary infections produced when one infected individual is introduced into a disease-free population. The importance of the role of R 0 can be seen by rewriting (2) as follows:

In order to avoid an epidemic, (7) should be non-positive. This is possible only if (7) is positive and thus, there will be an epidemic outbreak. Figure 3 illustrates the limiting values of the S , I, and R compartments for different values of R 0 where N is normalized to 1.

Inclusion of demographic dynamics may permit a disease to persist in a population in the long term. A disease is said to be endemic if it remains in a population for over a decade or two. Due to the long time period involved, an endemic disease model must include births as a source of new susceptibles and natural deaths in each compartment. In our study of the endemic S IR model, we consider constant birth and death rates. Using the notations in Table 1 , the scheme in Figure 4 can be expressed mathematically as: Figure 4 : Basic S IR model with vital dynamics.

Assuming that b equals µ, we can easily see that the sum of the above three equations yields zero when S + I + R = N holds in a non-varying population. Moreover, we observe that the average time of an infection is 1/(γ + µ), and since the infectious individuals infect others at rate β, R 0 is defined as β/(γ + µ).

By setting the left-hand side of the (8)-(10) to zero and solving for S , I, and R, we obtain the following two steady states (or equilibrium points) [16] :

where c 1 = µ/β and c 2 = γ/β. Points e 1 and e 2 denote the disease-free equilibrium (DFE) and endemic equilibrium (EE) points, respectively. Figure 5a depicts the system in a disease-free steady state when R 0 ≤ 1, whereas Figure 5b shows the occurrence of an endemic as the infected population reaches a limiting value of 0.225 when R 0 > 1. The stability of the system is driven by R 0 as it can be observed by linearizing the system of equations at these points.

The local stability of the model at these equilibrium points is analyzed via linearization. The Jacobian matrix for (8) and (9) is given as:

Evaluating the above matrix at e 1 and solving the characteristic equation, det(J − λI) = 0, where I is the identity matrix of size 2, results in the following pair of eigenvalues: In order to be a stable node, both eigenvalues should be negative. Therefore, e 1 is stable when β < γ + µ (or equivalently R 0 < 1) and unstable when β > γ + µ. Similarly, the Jacobian matrix evaluated at e 2 is as below:

One can easily observe that the determinant of (14) is positive as long as R 0 > 1. Hence, the endemic equilibrium is stable only when R 0 > 1.

The S-I phase portrait in Figure 6 shows how the model approaches the DFE and EE points with different initial values for S (0) and I(0). For the sake of simplicity, N has been normalized to 1. As depicted in Figure 6a , for R 0 = 0.5, the system eventually ends up at (1, 0), irrespective of the initial values of S (0) and R(0). On the other hand, an endemic occurs at (0.625, 0.225) for R 0 = 1.6 as in Figure 6b .

The phenomenon in which a parameter variation causes the stability of an equilibrium to change is known as bifurcation [16] . In continuous systems, this corresponds to the real part of an eigenvalue of an equilibrium passing through zero. With the basic reproduction number as the bifurcation parameter, Figure 7 shows a transcritical bifurcation where the equilibrium points persist through the bifurcation, but their stability properties change. Hence, we can conclude that:

A few recent studies have revealed interesting bifurcation behaviors in S IR models incorporated with factors such as varying immunity period, saturated treatment, and vaccination. We refer the interested reader to [19, 20, 21, 22] and the references therein for more details.

For viral diseases, such as measles and chickenpox, where the recovered individuals, in general, gain immunity against the virus, the S IR model is applicable. However, there exist certain bacterial diseases such as gonorrhoea and encephalitis that do not confer immunity. In such diseases, an infectious individual is allowed to recover from the infection and return unprotected to the susceptible class where he/she is prone to get infected again. Cases as such can be modeled using the susceptible-infected-susceptible (SIS) model as shown in Figure 8 , where the model variables are defined in Table 1 . 

In a fixed population, where there is no birth or death and individuals recover from the disease at the per capita rate of γ, the simplest form of the model in Figure 8 is given by:

By substituting S = N − I in (17), the system above can be reduced to:

Solving (18) analytically with I(0) = I 0 gives the solution for the complete system at time t as follows:

The behavior of the system in long-term can be inferred by looking at the possible values of (β − γ) that make (20) feasible. If (β − γ) > 0, then e −(β−γ)t → 0 as t → ∞. This can be written as:

If (β − γ) < 0, then e −(β−γ)t → ∞ as t → ∞ and thus, lim t→∞ I(t) = 0.

There exists two equilibrium points for this model which can be obtained by setting dI dt = 0 in (18) and solving for S and I. With R 0 as β/γ, we get:

where e 1 and e 2 denote the DFE and EE points, respectively. In terms of the basic reproduction number, if R 0 ≤ 1, the pathogen dies out as illustrated in Figure 9a because the infection in one individual cannot replace itself. If R 0 > 1, an existing infectious individual leads to more than one infection thus, spreading the pathogen in the population as seen in Figure 9b . 

The Jacobian matrix constructed from (16) and (17) is as follows:

Linear stability analysis for e 1 is done by solving the corresponding characteristic equation to obtain the following pair of eigenvalues:

The stability of e 1 depends on the value taken by λ 2 . The equilibrium point is a stable DFE if β < γ (or equivalently R 0 < 1) and unstable if β > γ (or R 0 > 1). Similarly, the eigenvalues of the characteristic equation for e 2 are:

In this case, the EE point is stable if β > γ and unstable if β < γ. Figure 10 depicts the vector plots for examples where e 1 and e 2 are unstable. In Figure 10a , the system converges to some state (highlighted in red) other than (1, 0) when R 0 = 2.5. Likewise, for R 0 < 1, the system converges to an invalid state as illustrated in Figure 10b . At β = γ or equivalently, R 0 = 1, a bifurcation occurs as the two equilibria collide (DFE equals EE) and exchange stability [7] . The forward bifurcation occurring at this threshold condition is similar to as seen previously in Section III.

The S IS model with varying population of constant size is as shown in Figure 11 . The corresponding system of differential equations for such a model is given below, where b and µ are assumed to be equal and S + I = N: 

To find the equilibrium points of the system, we set (26) and (27) 

where R 0 is β/(γ + µ). Figure 12 shows the system behavior for different values of R 0 . In Figure 12a , since R 0 is less than 1, the disease dies out and the system enters the diseasefree steady state. The same happens at R 0 = 1, where the two equilibria meet. For R 0 greater than 1, the disease does not die out, but instead remains in the population as an endemic with a limiting value. This can be seen in Figure 12b where an endemic occurs at (S (t), I(t)) = (0.56, 0.44) for R 0 = 1.8.

The corresponding Jacobian matrix for this model is given as: 

The eigenvalues of J in (29) for e 1 and e 2 are deduced as follows:

Linear stability analysis reveals that the disease-free equilibrium (e 1 ) is asymptotically stable if β − (γ + µ) ≤ 0 (or R 0 ≤ 1) and unstable otherwise [23] . Similarly, the endemic steady state is asymptotically stable if R 0 > 1. Figure 13 portraits the stability of e 1 and e 2 for different values of R 0 . In Figure 13a , with R 0 = 0.6, the system converges to e 1 = (1, 0) Figure 14 : The S IRS model without vital dynamics.

where it is stable. In the same manner, as in Figure 13b , the system eventually ends up at e 2 = (0.56, 0.44) for R 0 = 1.8, which is a valid endemic state. The forward bifurcation for the simple S IS model with demographic turnover and R 0 as the bifurcation parameter occurs at R 0 = 1 as studied earlier. Nevertheless, such models in presence of additional factors reveal interesting behaviors. Some recent works on S IS model that exhibit bifurcations consider factors such as non-constant contact rate having multiple stable equilibria [24] , non-linear birth rate [25] , treatment [26] , and time delay [27] .

The S IRS model is an extension of the basic S IR model in which individuals recover with immunity to the disease and become susceptible again after some time recovering. Influenza is a contagious viral disease that is usually studied using this model. In what follows, we investigate the model in both, absence and presence of demographic turnover.

The system of differential equations describing the S IRS flow diagram in Fig. 14 is as below:

With the three compartments summing up to N, we obtain the following two equilibrium points where e 1 and e 2 denote the DFE and EE, respectively. With c 1 = ν/(γ + ν), c 2 = γ/(γ + ν), and R 0 defined as β/γ, we have:

As illustrated in Fig. 15 , the infection dies out and reaches the disease-free steady state for R 0 ≤ 1 and stays as an endemic when R 0 > 1. 

The corresponding Jacobian matrix of the system is obtained by substituting R with N−S −I in (32) . Hence,

At e 1 , the matrix J yields the following two eigenvalues:

Therefore, e 1 is a stable node if λ 2 ≤ 0 or R 0 ≤ 1, and is a saddle point (unstable) if λ 2 > 0 or R 0 > 1. However, analyzing the stability of e 2 requires more care due to the structure complexity of the corresponding eigenvalues given as:

For the sake of clarity, let us denote ν(β + ν) and 2(γ + ν) as c and b, respectively. Rewriting the above eigenvalues gives:

Also, let us represent c 2 − ν(β − γ)b 2 as δ and ν(β − γ)b 2 as ζ. In order to study the stability of e 2 , we need to consider the following two cases:

• Case A: When δ ≥ 0, depending on the possible values that (β − γ) can take, we have the following two conditions :

(i) If β − γ ≤ 0 or R 0 ≤ 1, then ζ ≤ 0 which implies that the magnitude of δ is always greater than or equal to c 2 . Under this condition, the eigenvalues will always be real with opposite signs. Hence, the equilibrium point e 2 would be a saddle point as portrayed in Figure 16a , where the parametric values are the same as that in Figure 15a . (ii) If β − γ > 0 or R 0 > 1, then ζ > 0 and thus, the magnitude of δ would always be lesser than c 2 . In this case, λ 1 and λ 2 will always be real values with negative signs. As shown in Figure 16b , e 2 converges to a stable node with the same parametric values as given in Figure 15b . Here, δ = 0.019 which is lesser than c 2 = 1.519 and thus, results in a stable node at (0.333, 0.539).

• Case B: When δ < 0, β − γ should always be greater than zero and thus, the eigenvalues would be complex conjugates. Since the real parts of the eigenvalues are negative, the Figure 17 : The S IRS model with vital dynamics. equilibrium would be a stable focus point where the following condition holds:

As an example for this case, with R 0 = 2, β = 0.7, γ = 0.35, and ν = 0.25, the stable focus occurs at (0.5, 0.208) as depicted in Figure 16c .

As in Figure 17 , the S IRS model with standard incidence can be simply expressed as the following set of differential equations [23] :

It is worth mentioning that 1/γ and 1/ν can be regarded as the mean infectious period and the mean immune period, respectively. With ν = 0, the model reduces to an S IR model with no transition from class R to class S due to life-long immunity.

The system has a DFE point and a unique EE point denoted by e 1 and e 2 , respectively, as given below:

where R 0 is given by β/(γ + µ), c 1 = (ν + µ)/(γ + ν + µ), c 2 = γ/(γ + ν + µ), and b is assumed to be equal to µ. It should be noted that the basic reproduction number does not depend on the loss of immunity rate (ν). The system reaches DFE steady state for R 0 = 0.28 as shown in Figure 18a . Here, the parametric values are taken to be β = 0.04, µ = 0.043, γ = 0.1, and ν = 0.01. On the other hand, Figure 18b illustrates an example for which R 0 > 1. In this case, the system reaches the endemic state (0.625, 0.125, 0.25) for R 0 = 1.6. (41) gives the following Jacobian matrix which is obtained from (41) and (42):

Evaluating (45) at e 1 and solving its corresponding characteristic equations yields the following pair of eigenvalues:

We see that the disease-free equilibrium (e 1 ) is stable if λ 2 ≤ 0, i.e. β − (γ + µ) ≤ 0 or R 0 ≤ 1, and unstable otherwise. Similarly, on finding the eigenvalues for e 2 , we see that this unique endemic equilibrium point is stable for R 0 > 1. The instability of the equilibrium points can be seen in Figure 19 . For R 0 = 1.6, the vector plot in Figure 19a shows how the system does not converge to (S * , I * ) = (1, 0). In the same manner, Figure 19b illustrates the instability of e 2 at (S * , I * ) = (3.575, −0.892) when R 0 is less than unity. With the basic reproduction ratio as the bifurcation parameter, the system yields a transcritical forward bifurcation at R 0 = 1. More interesting behaviors have been reported when studied under factors such as stage structure [28] and non-linear incidence rates [29] , [30] .

Hitherto, we have dealt with models comprising of S , I, and R compartments. In these models, the infected individuals become infectious immediately. In the next two models, an exposed compartment in which all the individuals have been infected but are not yet infectious, is introduced. Such models take into consideration the latent period of the disease, resulting in an additional compartment denoted by E(t). The progression rate coefficient from compartment E to I is given as ε such that 1/ε is the mean latent period. Several other models with latent period such as S EIR and S EIRS have also been reported in the literature. However, such models are beyond the scope of this article. Interested readers can refer to [14] , [15] , and [31] for more on epidemic models beyond two dimensions. 6 The S EI model Unlike S IR models, the susceptible-exposed-infected (S EI) model assumes that a susceptible individual first undergoes a latent (or exposed) period before becoming infectious [32] . One example of this model is the transmission of severe acute respiratory syndrome (SARS) coronavirus.

For a fixed population of size N, the following differential equations describe the flow diagram in Fig. 20 :

The system should be analyzed asymptotically as it does not have a closed-form solution. We shall see that the population converges into a single compartment due to the straight-forward nature of the system. In what follows, we investigate the system behavior in terms of β and ε.

• Case A: When β 0 and ε 0, depending on the initial values of E(0) and I(0), the system approaches two different equilibrium points as below:

(i) If E(0) = 0 and I(0) = 0, the system remains in the following disease-free equilibrium:

This scenario is illustrated in Figure 21a where the complete population is susceptible at all times for N = 1, β = 0.8, ε = 0.5, and S (0) = 1.

(ii) If E(0) 0 or I(0) 0, then the system approaches the following equilibrium point in long-term: 

To prove this, consider the solution of (49) which is:

I(t) is a monotonically increasing function when E(t) > 0. Since all compartments are always non-negative, E(t) is greater than 0 when E(t) 0. Additionally, on solving (47), we get:

When I(t) > 0, we see that S (t) is a monotonically decreasing function. Unless when E(t) = 0 for all t, I(t) is a monotonically increasing function, and lim t→∞ S (t) = 0 as lim t→∞ t 0 I(τ)dτ = ∞. In terms of E(t), from (48), we see that the solution E(t) = ce −εt goes to 0 as t approaches infinity. In summary, if the condition E(0) 0 or I(0) 0 is satisfied, then lim t→∞ S (t) = 0. Since S = 0, lim t→∞ E(t) = lim t→∞ ce −εt = 0. Consequently, lim t→∞ I(t) = N. This case is depicted in Figure 21b .

• Case B: When β = 0 and ε 0, the reduced system yields the following solution:

Since lim t→∞ E(t) = 0, the system results in the following equilibrium solution as shown in Figure 21c , where the state of the system changes from the initial condition (S (0), E(0), I(0)) = (0.75, 0.24, 0.01) to steady state (0.75, 0, 0.25) for ε = 0.5: 

• Case D: When β 0 and ε = 0, the solution of the reduced system of differential equations is given as below: 

With b = µ in a birth-death population of total size N , the model illustrated in Fig. 22 can be written as follows:

The two set of equilibrium points obtained by setting the left-hand side of (60)-(62) to zero and solving for S , E, and I are:

with c 1 and c 2 defined as µ/(ε + µ) and ε/(ε + µ), respectively, and R 0 = βε/(µ(ε + µ)).

The DFE and EE steady-states are respectively, e 1 and e 2 . With β = 0.25, ε = 0.4, and µ = 0.2, Figure 23a depicts the disease-free steady state of the system as R 0 = 0.833. Likewise, Figure 23b shows how the system reaches the endemic equilibrium for R 0 = 1.7 when β = 0.54, ε = 0.5, and µ = 0.22. 

Considering (60) and (62), the Jacobian matrix for the system is as given below:

Evaluating the matrix at e 1 and solving the characteristic equation gives the following two eigenvalues:

For e 1 to be a stable node, both eigenvalues should be negative which means that λ 2 should be less than zero. Hence, ε(4β + ε) should be less than (ε + 2µ) 2 , which implies that R 0 < 1. In other words, e 1 is a stable point if R 0 < 1 and is a saddle point if the eigenvalues have opposite signs. Doing the same for e 2 , we get:

Similar to the above case, if R 0 > 1, then both eigenvalues are negative and thus, e 2 would be a stable node. Otherwise, e 2 would be a saddle point. In the simplest case, with R 0 as the bifurcation parameter, the system exhibits a forward transcritical bifurcation as it switches between the two equilibria. However, not much has been done in bifurcation analysis of such models in presence of other factors.

The susceptible-exposed-infected-susceptible (S EIS ) model is an extension of the S EI model such that in this model, the individual does not remain infected forever, but instead recovers and returns back to being susceptible again. Many sexually transmitted diseases (STD) and chlamydial infections are known to result in little or no acquired immunity following recovery [10] . In such cases, this model may serve as a suitable choice. 

The dynamical transfer of hosts depicted in Figure 24 can be formulated as follows, where N = S + E + I:

On solving (67)-(69) for S , E, and I, we obtain e 1 and e 2 which represent the disease-free and endemic equilibrium points, respectively:

where c 1 = γ/(ε + γ), c 2 = ε/(ε + γ), and R 0 = β/γ. As shown in Figure 25 , the system converges to e 1 when R 0 ≤ 1 and to e 2 when R 0 > 1. 

The asymptotic behavior of the model can be analyzed by studying the stability conditions of the system near its equilibrium points. The Jacobian matrix formed by using (67) and (69) is as below:

At e 1 , the matrix J yields the following eigenvalues:

Once again, in order for e 1 to be stable, both eigenvalues should be negative and since λ 1 is already negative, λ 2 should be less than zero. Hence, for λ 2 to be negative, −γ − ε + (γ − ε) 2 + 4βε should be negative, which on further simplification implies that R 0 < 1. Therefore, e 1 is a stable DFE when R 0 < 1 and is unstable when R 0 > 1. Similarly, calculating the eigenvalues of J evaluated at e 2 results in the following pair of eigenvalues:

Since λ 1 is always negative, the stability of e 2 depends on λ 2 . For λ 2 < 0, we observe that e 2 is a stable EE. On the contrary, when λ 2 > 0 (or equivalently, R > 1), the equilibrium point is unstable. This is clearly shown in Figure 26 where the stability of the system at the equilibrium points depends on R 0 . In Figure 26a is a stable endemic. As observed in previous cases, the model results in a forward bifurcation at R 0 = 1 when switching from one steady-state to another. Thus, the stability of e 1 is persistent at R 0 = 1.

In this subsection, we consider the S EIS model for a population of size N with birth and death rates that are constant. The set of equations given below refer to such a scheme illustrated in Figure 27 :

The following two equilibrium points are calculated by setting b = µ, the time-derivatives in (74)-(76) to zero, and solving for S , E, and I: e 1 : (S * , E * , I * ) = (N, 0, 0),

where c 1 , c 2 , and R 0 are defined as (γ + µ)/(γ + ε + µ), ε/(γ + ε + µ), and βε/((ε + µ)(γ + µ)), respectively. For R 0 ≤ 1, the system approaches e 1 which is disease-free, while for R 0 > 1, it ends up at e 2 .

The Jacobian matrix for this system is given as follows:

By evaluating the matrix in (78) at e 1 , we see that the equilibrium point is a stable disease-free equilibrium when both of the following eigenvalues are negative and is unstable if the eigenvalues have opposite signs. In terms of R 0 , e 1 is stable if R 0 < 1 and unstable if R 0 > 1:

Doing the same for e 2 yields a more complex pair of eigenvalues. However, on simplifying the eigenvalues, one can easily conclude that e 2 is stable when R 0 > 1 and unstable otherwise. At R 0 = 1, the system changes from e 1 to e 2 resulting in a forward bifurcation. A few works that reveal interesting bifurcation behaviors in S EIS model can be found in [33, 34, 35] .

Occurrences of certain events in nature and society may influence the behavior of the epidemic models seen so far. To guarantee that the model of choice mimics its counterpart in real-world, such influential factors should be taken into consideration. Many of these factors result in models with additional compartments which make them more complicated for mathematical analysis. In this section, a concise description on some of the most prominent factors that are considered in epidemic modeling is provided.

The time period between exposure and the onset of infectiousness is defined as the latent period. This is slightly different from the definition of incubation period which is the time interval between exposure and appearance of the first symptom of the disease in question. Thus, the latent period could be shorter or even longer than the incubation period. As seen before, S EI and S EIS are two examples of models with latent period. Models with higher dimensions such as S EIR, S EIRS , MS EIR, and MS EIRS have also been reported in the literature [36] . However, since these models cannot be reduced to planar differential equation systems due to their complexity, only a few complete analytic results have been obtained.

In absence of vaccination for an outbreak of a new disease, isolation of diagnosed infectives and quarantine of people who are suspected of having been infected are some of the few control measures available. Models such as S IQS and S IQR with a quarantined compartment, denoted by Q(t), assume that all infectives go through the quarantined compartment before recovering or becoming susceptible again [37] . However, vaccination, if available, is one of the most cost-effective methods of preventing disease spread in a population. We refer the reader to [14] , [15] , and [38] for more on models with vaccination and vaccine efficacy.

Models with time delay deal with the fact that the dynamic behavior of disease transmission at time t depends not only on the current state but also on the state of previous time [31] , [39] . Time delays are of two types namely, discrete (or fixed) delay and continuous (or distributed) delay. In the case of discrete time delay, the behavior of the model at time t depends on the state at time t − τ as well, where τ is some fixed constant. As an example, with τ being the latent period for some disease, the number of infectives at time t also depends on the number of infectives at time t − τ. On the contrary, the behavior of a model with continuous delay at time t depends on the states during the whole period prior to t as well.

Since individuals in different age groups have different infection and mortality rates, age is considered to be an important characteristic in modeling infectious diseases. Mostly, cases of sexually transmitted diseases (STDs) such as AIDS occur in younger individuals as they tend to be more active within or between populations. Likewise, malaria is responsible for nearly half of the death of infants under the age of 5 due to their weak immune system. Hence, such facts highlight the importance of age structure in epidemic modeling. Agestructured models are broadly classified into three types namely, discrete [40] , continuous [41] , and age groups or stages [42] .

In epidemiology, multi-group models describe the spread of infectious diseases in heterogeneous populations where each heterogeneous host population can be divided into several homogeneous groups in terms of geographic distributions, models of transmissions, and contact patterns. One of the pioneer models with multiple groups was investigated by Lajmanovich et al. in [43] for the transmission of gonorrhea. However, recent studies include differentiation of susceptibility to infection (DS) due to genetic variation of susceptible individuals, variation in infectiousness, and disease spread in competing populations [31] .

A central assumption in the classical models seen so far is that the rate of new infections is proportional to the mass action term. In these models, we assumed that the infected and susceptible individuals mix homogeneously. An increasingly important issue in epidemiology is how to extend these classical formulations to adequately describe the spatial heterogeneity in the distribution of susceptible and infected people and in the parameters of the spread of the infection observed in both, experimental data and computer simulations. Diffusion or migration of individuals in space are simply of two types: migration among different patches and continuous diffusion in space. In the former type, migration of individuals between patches depend on the connectivity of the patches. Models with migration between two patches [44] and n patches [45] have been reported in the literature. The latter type, on the other hand, takes into account the fact that the distribution of individuals and their interactions depend not only on the time t, but also the location in a given space.

Most of the classical epidemic models admit threshold dynamics, i.e. a DFE is stable if R 0 < 1 and an EE is stable if R 0 > 1. However, Capasso et. al [46] showed that it is likely possible for a DFE and EE to be stable simultaneously. Futhermore, periodic oscillations have been observed in the incidence of various diseases including mumps, chickenpox, influenza, and the like. The question that arises is why are classical epidemic models unable to capture these periodic phenomena? The main reason is the nature of the force of infection. Classical models frequently use mass incidence and standard incidence which imply that the contact rate and infection probability per contact are constant in time.

Nonetheless, it is more realistic (with added complexity) to consider the force of infection as a periodic function in time. As a simple example, consider the S IR model with a periodic incidence function F(I, t), and birth and death rates taken to be µ as given below [47] :

In recent years, much attention has been given to the study and analysis of chaotic behavior in epidemic models with non-linear infection forces.

Mathematical modeling of communicable diseases has received considerable attention over the last fifty years. A wide range of studies on epidemic models has been reported in the literature. However, there lacks a comprehensive study on understanding the dynamics of simple deterministic models through implementation. Aiming at filling such a gap, this work introduced some widely-appreciated epidemic models and studied each in terms of mathematical formulation, near equilibrium point stability analysis, and threshold dynamics with the aid of Mathematica. In addition, important factors that may be considered for better modeling were also presented. We believe that this article would serve as a good starting point for readers new to this research area and/or with little mathematical background.

Infectious diseases and human population history

Viruses, Plagues, and History

Deterministic modeling of infectious diseases: theory and methods

Stochastic epidemic models: a survey

An Introduction to Stochastic Epidemic Models, ser

Compartmental models for epidemics. Department of Mathematics, University of British Columbia

The mathematics of infectious diseases

Epidemic Disease in England: the Evidence of Variability and of Persistency of Type, ser

A contribution to the mathematical theory of epidemics

Population biology of infectious diseases: Part I

Epidemiological models with age structure, proportionate mixing, and cross-immunity

A Thousand and One Epidemic Models, ser

Modeling Infectious Diseases in Humans and Animals

An Introduction to Infectious Disease Modeling

Nonlinear Systems

Infectious Diseases of Humans: Dynamics and Control

Bifurcation thresholds in an SIR model with information-dependent vaccination

Bifurcation analysis in an SIR epidemic model with birth pulse and pulse vaccination

Stability and bifurcations in an epidemic model with varying immunity period

Qualitative and bifurcation analysis using an SIR model with a saturated treatment function

On the global stability of SIS, SIR and SIRS epidemic models with standard incidence

A simple SIS epidemic model with a backward bifurcation

Bifurcation analysis of an SIS epidemic model with nonlinear birth rate

Backward bifurcation in simple SIS model

Bifurcation and chaos in SIS epidemic model

Stability of Hopf bifurcation of a delayed SIRS epidemic model with stage structure

Bifurcation analysis of an SIRS epidemic model with generalized incidence

Bifurcations of an SIRS epidemic model with nonlinear incidence rate

Dynamical Modeling and Analysis of Epidemics

Population dynamics of fox rabbits in Europe

Some epidemiological models with nonlinear incidence

Complex dynamics of discrete SEIS models with simple demography

Dynamic analysis of an SEIS model with bilinear incidence rate

Mathematical Understanding of Infectious Disease Dynamics

Effects of quarantine in six endemic models for infectious diseases

Global results for an epidemic model with vaccination that exhibits backward bifurcation

Time delays in epidemic models: modeling and numerical considerations

Dynamics of a discrete age-structured SIS models

Continuous-time age-structured models in population dynamics and epidemiology, ser

An age-structured model for pertussis transmission

A deterministic model for gonorrhea in a nonhomogeneous population

Qualitative analysis for communicable disease models

An epidemic model in a patchy environment

Analysis of a reaction-diffusion system modeling man-environment-man epidemics

Melnikov analysis of chaos in an epidemiological model with almost periodic incidence rates