An introduction to compartmental modeling for the budding infectious disease modeler

ABSTRACT Mathematical models are ubiquitous in the study of the transmission dynamics of infectious diseases, In particular, the classic ‘susceptible-infectious-recovered’ (SIR) paradigm provides a modeling framework that can be adapted to describe the core transmission dynamics of a range of human and wildlife diseases. These models provide an important tool for uncovering the mechanisms generating observed disease dynamics, evaluating potential control strategies, and predicting future outbreaks. With ongoing advances in computational tools as well as access to disease incidence data, the use of such models continues to increase. Here, we provide a basic introduction to disease modeling that is primarily intended for individuals who are new to developing SIR-type models. In particular, we highlight several common issues encountered when structuring and analyzing these models.


Introduction
Recent high profile infectious disease outbreaks around the globe -including SARS, Ebola, and Zika -have reignited a general interest in understanding the spread of disease. The utility of mathematical and computational approaches to simulate and numerically analyze infections within populations has also benefited from major advances in computational power over the past several decades. As a result, public health officials and policy makers have increasingly relied on a quantitative understanding of infectious disease dynamics.
Mathematical models of infectious disease dynamics have a rich history that dates back more than 100 years. Mathematically simple formulations that describe the transition of individuals in a population between 'compartments' that capture the infection status of individuals leads to surprisingly significant insight. Their elegance and simplicity allow the ease of expansion to more complexities through, for example, the addition of compartments. Expanding these models is often straight forward but the apparent simplicity can mask subtle, but important, model structure and parameterization choices. Furthermore, while there are many wonderful texts focused on infectious disease modeling, there exist several complexities that are rarely discussed in sufficient detail for a novice disease modeler.
Here we highlight several of these potential pitfalls, describe how they arise and, when possible, how to address them. We also point to further resources in the literature. We begin with a brief introduction to disease modeling as a means to provide basic context as well as introduce the notation that we will use throughout this paper. We often provide more mathematical detail and discussion in sections in the Appendix.
Following our introduction, we focus on analysis of the compartmental epidemiological model with an emphasis on computing the basic reproductive number, commonly known as R 0 . We focus on the following two difficulties: (i) deriving the force of infection and (ii) instances when R 0 does not represent a precise disease extinction threshold. Keep in mind, this is by no means meant to be a all-encompassing discussion of infectious disease modeling but a resource to supplement other more comprehensive texts (e.g. Anderson, May, & Anderson, 1992;Brauer, 2008;Brauer & Castillo-Chavez, 2001;Brauer & Kribs, 2015;Brauer, vanden Driessche, & Wu, 2008;Keeling & Rohani, 2008;Martcheva, 2015;Vynnycky & White, 2010).

A brief introduction to disease modeling
Mathematical models that describe the population-level dynamics of infectious diseases are typically based on the classic Susceptible -Infectious -Recovered (SIR)-type framework, a structure inspired by the seminal work of Kermack andMcKenrick in 1927 (Kermack &McKendrick, 1927). In SIR-type models, individuals are separated into compartments based on their infectious status. Thus, any individual can only be in a single compartment at a given time, although they can transition between compartments. In this section, we briefly introduce and provide an analysis of a simple SIR model. We then generalize the model analysis so that it can be applied to a large suite of compartmental disease models.

Constructing compartmental models: the SIR framework
Compartmental disease models divide a population into groups (or compartments) based on each individual's infectious status and track the corresponding population sizes through time. Here, we first consider a model with straight forward properties: susceptible individuals (S) become infected through contact with infectious individuals (I), and infectious individuals recover (R) at a fixed rate γ and confer lifelong immunity. The fixed rate γ has the interpretation that 1/γ is the average amount of time spent in the infectious class (see Appendix 1 for a detailed explanation). Although there typically is some lag between acquiring infection and becoming infectious, here we do not distinguish between infected and infectious individuals; we instead assume that individuals are infectious immediately upon infection. This assumption will be relaxed later in Section 3.1. Almost all diseases have some delay between infection and infectiousness, however, when that delay is short (hours or days) and the goal is understanding the long term dynamics (years in the future), the assumption that individuals are immediately infectious may be reasonable. Examples of some infections that are frequently modeled within the SIR framework include many common childhood diseases such as chicken pox and measles. An important common feature among these pathogens is that infection is thought to confer life-long immunity so that, from a modeling perspective, individuals move directly into (and remain in) the recovered compartment upon recovery. Now, to model the infection process we consider a system of ordinary differential equations describing the change in the population size of each compartment. Let λ be the force of infection, or the per capita rate at which susceptible individuals acquire infection. It is important to note that the force of infection is not a constant but rather a function of the size of the infectious compartments; that is, λ = λ(I).
The definition of λ(I) can indeed be tricky and depends on assumptions about the nature of interactions between individuals (see Section 4.1). However, complications in the different ways to define the force of infection most often come into play when the total population size is changing, which is not the case at the moment. Thus, we will ignore this complexity for the time being. The force of infection is comprised of a transmission rate β or, more precisely, the product of the contact rate and the probability of transmission given contact, and an interaction term with infectious individuals. In order to preserve the units of the equation, using our definition of β (which has units time −1 ), we will think about the interaction with the proportion of infectious individuals: I/N. Thus, The definition of λ(I) will change for each model, although it will always remain a function of the infectious class(es). Here, we have written λ(I) as a function of the infectious class, but frequently the dependence is not explicitly included. in model formulation.
For now, we ignore demography, i.e. natural birth and mortality, but relax this assumption later. We can now write down a system of ordinary differential equations (ODEs) of the SIR model: A schematic diagram of this model is provided in Figure 1(A). Here, the nature of our equations show us that the population is 'closed'; that is, the total population size (defined by N = S+I+R) is constant over time (dN/dt = 0). A nice property of this assumption is that the system is entirely determined by two, rather than three, of the equations. For example, we can write the recovered population as R = N−S−I so that only susceptible and infected populations need to be determined explicitly.
A final comment before we move into the model analysis is that throughout this paper S, I, and R always represent population densities (although we will refer to population density and size interchangeably throughout this paper). However, some texts define S, I, and R as proportions of the population in each class, respectively. Non-dimensionalizing the state variables in this way for disease models is a common modeling choice, and the associated derivation is demonstrated in the Appendix (see Appendix 2). Our choice to consistently write the models in terms of densities, rather than proportions, will make it more straight forward to interpret some common pitfalls in disease modeling.  (1). S, I, and R represent the total number of susceptible, infected, and recovered individuals in a population, respectively. λ(I) is the force of infection, and γ is the recovery rate. The dashed arrow indicates the necessity of interaction with the I compartment for transition from S to I. (B) Same as above, but this schematic includes demography in grey (as in Equation (2)).

Analysis of the SIR model
Many questions can be asked regarding the dynamics of an infectious disease: How quickly will the disease sweep through a population? How many people will be infected during the outbreak? Will the disease persist in the population? Compartmental ODE models such as the SIR framework introduced above can help answer these questions through analysis of the ODE systems.

Finding equilibria
A typical first step in analyzing a system of differential equations is finding the equilibria. At equilibrium we do not mean that individuals are fixed within particular compartments, but rather the rate of individuals entering and leaving a compartment is exactly balanced. This is equivalent to setting the right hand side of our equations, i.e. the rate of change of individuals into and out of each compartment, equal to zero.
In epidemiology, models generally have two important equilibria: (1) the disease-free equilibrium (DFE) and the endemic equilibrium. The DFE requires there to be no infected individuals in the population, or in the case of the SIR model I * = 0, where the star designates an equilibrium solution. In contrast, the endemic equilibrium corresponds to the state in which infected individuals persist indefinitely such that I * > 0.
In the model presented in Section 2.1, the canonical DFE is given by (S * , I * , R * ) = (N, 0, 0); in other words, the entire population is susceptible to infection. For this model, there is no endemic equilibrium as the infectious population always returns to zero once the epidemic has completed (see Appendix 3). Biologically speaking, an endemic equilibrium requires that there is a continuous supply of susceptible individuals -known as 'susceptible replacement' -in order for infection to persist endemically. This replacement process can, for example, occur through births into the susceptible population or waning of immunity. Without such replacement, the number of infected individuals returns to zero following an epidemic, as in the SIR model without demography (i.e. in Equation (1)). In this section, we therefore introduce a small modification to the SIR model of Section 2.1: demography, i.e. the inclusion of natural birth and death.
For mathematical simplicity, we make the following key biological assumptions: (i) the birth rate and death rate are equal and given by μ, (ii) all individuals are capable of reproducing and are equally subject to mortality, and (iii) all individuals are born susceptible to infection. The first two assumptions maintain a constant population size and the third assumption is made for simplicity (and is a valid assumption for many directly transmitted pathogens). It is possible to model infectious disease dynamics without these assumptions. For example, a discussion of models that include demography with variable population size can be found in Brauer (2008), Hethcote (2000), and Ledder (2017). The simple SIR model with demography is given by: and a schematic representation of this model is provided in Figure 1(B). For this model, we can find both a DFE and an endemic equilibrium, which are given by: respectively.
Finding equilibria is only the first step to understanding long-term behaviour in a system. We must also determine which of the behaviours are typically realized. This involves determining the stability of the equilibria and will show us whether we will approach the equilibria or move away from it, assuming we have started nearby. In order to determine the stability of each of these equilibria, we can apply a standard linear stability analysis (more details can be found in any standard ODEs or linear algebra textbook). In brief, we linearize our system around each equilibria using the Jacobian matrix evaluated at the chosen equilibrium of interest. We then find the associated characteristic polynomial, whose roots are the eigenvalues to the Jacobian matrix. For local stability of the equilibrium, we require that all of the eigenvalues have negative real parts; otherwise, the equilibrium is locally unstable. In some cases, the value of the eigenvalues can be calculated directly. However, under more complicated scenarios an alternative method can be used to determine the sign of the eigenvalues (e.g. with a symbolic programming language or analytically with the Routh-Hurwitz criteria, see Brauer & Castillo-Chavez, 2001).
Using linear stability analysis, we can determine the local stability of each equilibrium. However, even for simple systems, this does not always give us a clear biological interpretation on when disease will spread or die out globally. Furthermore, when using more complicated models, there may be multiple stable equilibria present simultaneously. In such cases, it is important to obtain a more global picture of the infection dynamics using, for example, bifurcation analysis to supplement the linear stability analysis (as in Section 4.2.1).
In the standard SIR model with demography, however, we can use a standard linear stability analysis to identify conditions under which the DFE and endemic equilibrium are stable as well as unstable. As it turns out, stability is closely related to an important quantity in epidemiology: R 0 , often called the basic reproductive number (or ratio). Therefore, we introduce R 0 before discussing the stability of this model.

The basic reproductive number: a biologically motivated introduction
The basic reproductive number, R 0 , is arguably the most important quantity in disease modeling and is generally defined as the 'average number of secondary cases arising from a typical primary case in an entirely susceptible population'. There is rich mathematical theory that describes how this quantity can be computed for a range of SIR-type models with varying degrees of complexity. We will introduce some of this theory beginning in Section 3. Here, however, we find R 0 for the SIR model with demography in an intuitive way motivated by one question central to disease modeling: under what conditions will a disease die out?
To begin answering this question, we assume that the population at the initial time (t = 0) has one infected individual and is otherwise susceptible. Therefore, I(0) = 1 and S(0) = N − 1. The disease will fail to produce additional infections, i.e. fail to invade, if the rate of change in the infectious compartment is decreasing at time t = 0, or dI dt t=0 < 0.
Recalling from Equation (2), the differential equation for I is given by Now, we want to evaluate this equation at t = 0 by plugging in S = N−1,I = 1 and consider when this equation is negative: We rearrange this inequality to a form which provides biological intuition, and now shifts the threshold of interest to one (rather than zero when considering whether the size in the compartment will grow or shrink). Thus, In other words, the threshold is determined by the product of the transmission rate β with the quantity 1/(γ + μ) (which is related to the average time spent in the infectious class and the average lifespan, see Appendix 1) and the initial fraction of susceptible individuals. Combined this must be 'small enough' (i.e. < 1) so that the epidemic cannot take off. As long as we are considering relatively large population sizes, we can look at the approximate behaviour in the limit N → ∞) and simplify such that implies the disease will die out. A parallel analysis that begins with the question 'under what conditions will a disease persist?' reveals that implies the disease will persist. This quantity, β/(γ + μ), is defined as the basic reproductive number for the SIR system. It is a threshold quantity where the 'critical value' is at R 0 = 1, and we are now better able to interpret our original definition: if a typical infectious individual infects less than one person on average during the time they are infected then a disease will die out. In contrast, if a typical infectious individual infects more than one person on average during the infectious period then the disease will spread. Intuitively, this makes sense: R 0 is the product of the transmission rate and -roughly -the average amount of time spent in the infectious class. Notice that this definition includes several important assumptions: (i) R 0 is the average number of onward infections produced by a single individual (and therefore excludes individual variation), (ii) the average transmissibility among individuals is used in the definition of R 0 (again ignoring individual variation), and (iii) all individuals in the population are susceptible at time t = 0. In other more complicated settings, this definition changes such as when one of these assumptions is broken by the model structure (see Section 4.1.2 on structured populations). We note that the mathematical form of R 0 will differ for each model and numerically will depend on model parameters. However, almost all methods of computing R 0 require examination of when the DFE loses stability.
As you may expect, the way we have derived the basic reproductive number here by examining whether the infectious population grows or shrinks at time zero is only possible when the model framework is simple. However, there are multiple ways to determine R 0 in a more general setting, some of which we discuss later (see Section 3).

R 0 and equilibrium stability
We are now well-suited to determine the stability of both equilibria for the SIR model with demography. Now that we have found R 0 for this model, we can actually simplify the complicated endemic equilibrium (Equation (4)) by replacing the appearance of β/(γ + μ) with R 0 : It is no coincidence that R 0 keeps appearing. In fact, this threshold is important in evaluating the stability of our equilibrium. As one might expect, as long as R 0 < 1 then the DFE is stable. More generally, we can use linear stability analysis, i.e. when the eigenvalues of the Jacobian matrix evaluated at the DFE have negative real part, to show: We emphasize that when using most standard methods for the computation of R 0 , which involve linearization around the DFE, we can really only say something definitively about the local stability of the DFE. In other words, the derivation of R 0 typically depends on the stability of the DFE. As R 0 increases past the critical threshold value of one, the DFE is destabilized -thereby often, and indeed in the case of the SIR model with demography, this leads to stability of the endemic equilibrium. We will introduce an example later where this is not the case (see Section 4.2.1).

Methods for computing R 0
In the previous section, we defined the threshold R 0 for disease invasion simply by examining the sign of the differential equation describing the change in the infectious population for the SIR model with demography (see Section 2.2.2). We then described how a linear stability analysis can be used to describe equilibrium stability in terms of R 0 of the equilibria and observed how R 0 provides information on stability properties. See Appendix 4 for a more detailed example of linear stability analysis. This is straightforward for simple models, but as models increase in complexity they often contain multiple compartments with infectious individuals, making this process difficult. In this section we introduce a commonly used method for finding R 0 : the next generation matrix method (Diekmann, Heesterbeek, & Johan, 1990;Diekmann, Heesterbeek, & Roberts, 2009). Importantly, this method is centered on considering the dynamics near the DFE.

Next generation matrix method
We introduce the next generation matrix method by demonstrating the process with a particular example where infected individuals are not immediately infectious. This is a much more biologically valid assumption for most diseases (e.g. Ebola, influenza, and pertussis), as individuals typically have some delay between infection and becoming infectious. However, the SIR model, where this delay is ignored, is often a reasonable approximation, especially when the duration of the prior to infectiousness is much shorter than the time spent infectious. Therefore, we now assume that there is a latent period prior to the onset of infectiousness in which individuals have been exposed to disease but are not yet capable of transmission. This class is commonly called the exposed, or E, class and the corresponding model is known as the SEIR model ( Figure 2).
In particular, we let σ be the rate at which individuals leave the E class (recall that implies that 1/σ is the average amount of time spent in the E class), and we can extend the model in Equation (2) to: Figure 2. Schematic representation of the SEIR model. λ(I) is the force of infection, σ is the rate at which individuals transition from the exposed class to the infectious class, γ is the recovery rate, and μ is the natural mortality rate. The dashed arrow indicates the necessity of interaction with the I compartment for transition from S to E.
where N = S+E+I+R. In this model, the force of infection λ(I) is the same as in the SIR model (λ(I) = βI/N). There are no new infections contributed through interactions with the exposed individuals, as they are not infectious. If we try to find R 0 using the same methods as we did for the SIR model, we run into trouble. More precisely, when examining the behaviour of I at t = 0, in contrast to the SIR model, we cannot factor out an I. Fortunately, the next generation matrix method allows us to generalize the process of computing R 0 . A nice feature of the next generation matrix method is that it only requires use of the DFE which is often easy to compute. In the SEIR model, the DFE is given by (S * , E * , I * , R * ) = (N, 0, 0, 0). After determining the DFE, we must create a sub-model that only considers the 'disease' compartments, a subset of the equations in the SEIR model. The disease compartments are those that include individuals that are in any stage of infection which, for the SEIR model, includes both the exposed and infectious individuals. Therefore, this sub-model will only contain the E and I equations (i.e. Equations (6)- (7)), which we write in the form d where x is a vector of the j disease compartments; in the SEIR model j = 2 because the disease compartments are E and I. The right hand sides of Equations (6)-(7) are therefore contained in the vectors F( x) and V( x).
Here, F( x) contains any terms that directly lead to new infections entering each compartment j. Notice that the second element of F is zero because no new infections enter the I compartment, rather they transition from the E compartment into the I compartment.
V( x), on the other hand, can be further broken down as where V − ( x) and V + ( x) contain all other outputs and inputs, respectively, from each disease class. This includes terms such as mortality or transition between classes. Although mathematically it is not necessary to separate V( x), we introduce this notation here as it is commonly seen in the literature on the next generation method. For the SEIR model, we can now define the vectors in Equations (9)-(10) as Notice, this is exactly equivalent to our sub-model original equations (Equations (6) and (7)). Next, we must linearize around the DFE, which involves the Jacobian, i.e. the matrix of partial derivatives, evaluated at the DFE. Recalling that the DFE is given by (S * , E * , I * , R * ) = (N, 0, 0, 0), the Jacobian matrix of the sub-model evaluated at the DFE is given by As a result, we can factor out the vector x = E I on the right hand side, leaving us with ⎡ Notice that in the F matrix the only non-zero term β arises from the derivative of the first row of F with respect to I, evaluated at S = N. To simplify, we define matrices in Equation (15) as F and V, respectively, and rewrite Equation (15) as Recall from Section 2.2.2 that the ultimate goal is to determine when the number of infectious individuals, I, is increasing or decreasing when the population is near the DFE, i.e. when a single infected individual is introduced into an otherwise susceptible population. We can now investigate a similar question: when are the compartments containing the infectious classes increasing or decreasing near the DFE? From linear stability analysis we know that stability of the DFE is equivalent to the real part of all eigenvalues of (F − V) being less than zero. From Equation (16), we can extract what is known as the next generation matrix, or Here, an analogous threshold (to the eigenvalues of the Jacobian) for when the DFE is stable is given by the spectral radius of the next generation matrix denoted by where the spectral radius (ρ) is the eigenvalue with the largest magnitude. While we omit the proof here, it can be found in other sources, e.g. van den Driessche and Watmough (2002). For the SEIR model, the next generation matrix is given by We find the eigenvalues of this matrix are βσ /(γ + μ)(σ + μ) and 0. Here, R 0 (or the spectral radius of FV −1 ) is given by .
For the SEIR model, notice that R 0 is very similar to the basic reproductive number in the SIR model with demography, (β/(γ + μ)), multiplied with one additional term, (σ/(σ + μ)). This term leads to a smaller R 0 for the SEIR model because some individuals are lost to natural mortality while in the E class and never actually reach the infectious class. While we provided a brief introduction to the next generation matrix in this section using a particular example (the SEIR model), more formal mathematical treatments of this can be found in other references, e.g. Diekmann and Heesterbeek (2000), Martcheva (2015) and van den Driessche and Watmough (2002). We also note that analytically finding the roots of the characteristic polynomial (and thereby finding R 0 ) can be challenging; it is common to use a symbolic mathematical programme (e.g. Mathematica) or to find it numerically (e.g. using MATLAB). Furthermore, a great deal of tractability, both mathematically and biologically, may be gained by first non-dimensionalizing the model. Such simplifications may reveal biologically relevant quantities that were obscured by relying solely on the mathematical computation. We refer the reader to Ledder (2017) for a more detailed description of non-dimensionalization of models.

A cautionary tale of computing R 0
The next generation matrix method is a widely used method for computing R 0 for a variety of compartmental models that have varying degrees of complexity. There are, however, instances in which the construction of a model is not amenable to this -or other -standard methods for analytically finding R 0 .
One common example in which R 0 cannot be directly computed is when there is 'importation'. Suppose a particular pathogen follows the standard SEIR structure with demography as in Equations (5)-(8). Now, assume that there is some external source of infection that does not depend on the number of infected individuals present in the population. For example, an infected individual from another population may briefly encounter a susceptible individual and transmit infection. If we assume the importation rate is given by φ, then this modifies the force of infection so that The S and E equations change to ensure that the total population size remains constant. While R 0 cannot be computed explicitly for a model containing infectious imports, some modeling studies still attempt to take advantage of the information R 0 provides in similar standard models. For example, R 0 can first be computed for the model in the absence of imports (i.e. when φ = 0). This R 0 then represents an 'intrinsic' basic reproductive number; that is, R 0 in the absence of any spatial processes. The dynamical impact of importation can then be explored numerically through simulations (e.g. Blackwood, Streicker, Altizer, & Rohani, 2013). Other models with importation focus on other aspects of the model all together, such as direct parameterizations that are independent of the precise R 0 value itself (e.g. Shrestha et al., 2015Shrestha et al., , 2013.

Other methods for finding R 0
The next generation matrix is the most commonly used method for finding R 0 in standard disease ecology studies. However, there are alternative methods such as the graph theoretic method (de Camino-Beck, Lewis, & van den Driessche, 2009). While we do not discuss them in this paper, a nice discussion can be found in Heffernan, Smith, and Wahl (2005) and Li, Blakeley, and Smith (2011).

Formulating the force of infection
Recall that the force of infection gives the per capita rate at which susceptible individuals acquire infection. Generally, the force of infection, denoted by λ(I) can be broken apart into three components: (i) the contact rate, (ii) the probability of infection given a contact, and (iii) the size of the infectious population. The product of the first two is the overall transmission rate, standardly denoted by β. Although these quantities are tightly connected we will address each of them independently, beginning with the infectious population.

Infectious population
When the total population size is constant we do not need to worry too much about how the infectious population is used within the system. However, if we expect the total population size to change due to, for example, differing birth and death rates or disease induced mortality, we need to be more careful in its formulation. Throughout our discussion thus far, the force of infection λ(I) depended on the infectious proportion of the total population: I/N.
There are two main ways to incorporate the infectious population into the force of infection λ(I). We generally assume that transmission is either 'frequency dependent' or 'density dependent'. In the former, which we have used throughout this paper, λ(I) is independent of population size and depends on the fraction of individuals who are infected in a population so that λ(I) = βI/N. In the latter, λ increases proportionally with the number of infected individuals I such that λ(I) = βI. Here, if the density of the infectious population increases, i.e. there are more infectious individuals in the same area, the force of infection also increases.
Frequency dependent transmission, also known as standard incidence, is typically used for sexually transmitted pathogens in human population (Begon et al., 2002;Keeling & Rohani, 2008). As the population size changes, the number of interactions between susceptible and infectious individuals does not change. Hence, the relevancy for sexual transmission, where the number of partners an individual has is likely independent of the number of individuals near by, at least when the population is not very low. In contrast, density dependent transmission allows for reduced interaction between susceptible and infectious individuals as the population size changes (Begon et al., 2002;Keeling & Rohani, 2008). Consider a population that is declining due to disease-induced mortality from a disease with high mortality. As fewer individuals remain, it is more difficult for individuals to come into contact, reducing the force of infection (Keeling & Rohani, 2008). Density dependent transmission is typically used to model airborne and directly transmitted infections such as measles or influenza. We remind the reader that these definitions only differ when the total population size is variable.
Although the vast majority of the literature uses one of these two functions, a variety of phenomenological non-linear formulations have also been introduced (Begon et al., 2002;McCallum, Barlow, & Hone, 2001). When altering the transmission function it is essential that careful consideration is given to the estimation of the transmission rate, β, as the units of the term must be consistent within the equation. For example, a non-linear force of infection such as the following λ(I) = βI q where β is the transmission rate, I is the number of infected individuals, and 0 < q < 1 is a constant defining a power law relationship (Hochberg, 1991;Liu, Levin, & Iwasa, 1986). The exponent on I requires that β take the units ( q √ individuals · time) −1 rather than time −1 as discussed above in Section 2.
A further complication arises when the population is subdivided by characteristics other than those that are disease-related such as risk status or age. Nothing changes when we consider density-dependent transmission because we are only concerned with the actual number of infectious individuals from a particular group. However, closer attention to the details is required when using frequency dependent transmission: the fraction of individuals infectious from a particular class refers to the number infectious of that class divided by the total number from that class rather than the fraction of the total population of all infectious individuals.
To demonstrate this, let us consider the case of an age-structured model, pictured in Figure 3. Although age changes continuously and can be more generally described using a continuous-age model with, for example, partial differential equations (Feng et al., 2016;Martcheva, 2015), we examine an age-structured model with two discrete age classes so that our population is sub-divided into children and adults. Such a model would be relevant for infections that spread much faster among children than adults, whether due to higher susceptibility or greater contact rates. We will use the subscript C to refer to children populations and the subscript A to refer to adult populations. For example, the number of susceptible children will be denoted as S C . We assume that the total number of children is given by N C = S C + I C + R C and the total number of adults is given by N A = S A + I A + R A . We also assume that children age into the adult compartment at a rate f and that only adults reproduce and all births enter the S C class. Finally, only adults are subject to natural mortality so that the total population size remains fixed. Now, A schematic representation of this model is provided in Figure 3. Notice that we distinguish between the force of infection to children λ C (I C , I A ) and the force of infection to adults λ A (I C , I A ). In the former, the rate at which susceptible children become infected will depend on the rate at which susceptible children contact both infectious adults and children as well as the fraction of the adult and child populations that are infectious in their respective populations. Additionally, in contrast to the models presented earlier, this model contains multiple infectious classes. Therefore, the forces of infection now depend on both infectious classes. We now assign two subscripts to the transmission rate so that β ij is the transmission rate from individuals in population j to those in population i, where i and j can take on values of C or A. Notice that the proportion of infectious individuals refers to the proportion of each class. Thus, the force of infection for children and adults, respectively, are

Contact rate
A key component of the force of infection is the quantity β, i.e. the transmission rate, which, as described earlier, can be broken down to the product of the probability of transmission per contact and the contacts per unit time. In many unstructured populations (for example, those that do not account for age or sex structure), β can usually be assigned a fixed value without paying direct attention to the two components that make up this quantity. However, there are also a lot of model structures that require more nuanced assumptions on the structure of β. For example, the notion of contacts per time becomes non-trivial when a population has been grouped by factors other than disease status, such as by age. We follow from the age-structured model introduced in the previous section. The different β ij values are a composition of both the probability of transmission and the contact rate. Contacts of children with other children or of adults with other adults are easiest to account for. However, contacts between adults and children must be symmetric, i.e. an adult contacting a child requires a child contacting an adult. Thus, when recording contacts between individuals as occurs in experimental attempts to estimate contact rates, the resulting data of total reported contacts will be symmetric (see Table 1 for example data). However, as long as the total number of children and adults are different, obtaining the contact rate requires that the total number of reported contacts correct for differences in population size (Table 2 computes the contact rates associated with the data in Table 1).

Conservation of biting
Models of vector-borne diseases where transmission does not occur directly between humans pose additional complications, as contacts no longer mark the interaction between two individuals. The simplest compartmental model of a vector-borne disease contains susceptible and infectious populations for both human and vector populations. In this model, the subscripts H and V denote variables and parameters that pertain to the human and vector populations, respectively. Here, we will consider standard assumptions for many mosquito-borne illnesses like malaria or dengue that new human infections only originate from bites of the infected vector population and new infected vectors only occur from biting an infected human. Furthermore, we ignore human demography to focus on infection and recovery to a susceptible state at rate γ . Vector demography, as it occurs at a much shorter time scale relative to humans, will be retained such that vector populations have a birth and death rate μ. This gives us the following system of ODEs where N V = S V + I V is the total vector population size (see schematic representation in Figure 4). Given how we defined the force of infection in previous sections, our first instinct may be to define the forces of infection for human and vector populations as follows: such that the infection of humans is dependent upon the transmission rate β HV from vectors to humans as well as fraction of the vector population that is infected; whereas the infection of vectors is dependent upon the transmission rate β VH from humans to vectors and the fraction of humans that are infected (N H = S H + I H is the total number of humans). As in the previous models, the transmission rates β ij include both the contact rate and the probability of transmission given contact. Here, contacts that are capable of generating infections occur through bites, b i , and we denote the probability of infection by α ij . In other words, we require that the transmission rate from j to i (β ij ) to be equal to the product of the contact rate b i and the infection probability α ij so that Contact through bites has the symmetrical property such that every bite received by a human is a bite taken by a vector, a concept referred to as 'conservation of biting'. We observed a similar idea in contacts between adults and children in Section 4.1.2. Consequently, the number of bites depends on the population sizes of both the human and vector populations. To demonstrate what this means, let us assume that b V is the biting rate of a single mosquito and b H is the rate at which a single human gets bitten per unit time. Then the total number of bites received by the human population and total number of bites taken by the mosquitoes should equate as follows: Solving for N V in the equation for conservation of biting (Equation (18)) and substituting in λ H (I H , I V ) into our original force of infection equations (Equation (17)) leads to As the same biting rate appears in both equations, typically the subscript is dropped. As a note of caution, the canonical model of malaria (known as the Ross-Macdonald model after Ronald Ross and George Macdonald who contributed to its formulation) often refers to biting rate by the parameter a (Anderson et al., 1992;Macdonald, 1957;Ross, 1911). Furthermore, although they use the notion of conservation of biting, the derivation is often omitted. The appearance of the ratio of mosquitoes to humans, commonly denoted by m, results from this assumption.

Disease can persist when R 0 < 1
Up to this point, we have assumed that we can define the threshold of disease spread R 0 by examining the local behaviour near the DFE. In other words, we generate conditions for which the DFE becomes unstable (i.e. when R 0 > 1) which, in turn, indicates that the endemic equilibrium exists and is locally stable. However, this scenario is not the only possibility.

Backward bifurcations
As it turns out, there exist cases where both the DFE and an endemic equilibrium stably coexist over a range of parameter values and are separated by an additional unstable equilibrium. In such cases, R 0 is no longer a strict threshold of disease decay and growth. In the infectious disease literature, these scenarios are commonly termed backward bifurcations. In these cases, R 0 < 1 remains as a requirement for the DFE to be locally stable. However, the dynamics have changed such that the condition for instability of the disease free equilibrium (that is, R 0 > 1) is no longer equivalent to the condition for the stability of a biologically relevant endemic equilibrium. Here, we describe the qualitative structure of the simplest instance where a backward bifurcation is observed. In this case, the system contains three equilibria.
When R 0 << 1 only the DFE exists in the system. As R 0 increases, there comes a point 0 < R 0 < 1 where two endemic equilibria appear. Mathematically this occurs through what is known as a saddle node bifurcation (see for example Strogatz, 2014 or most linear algebra or ODE texts). The appearance of a saddle node bifurcation often arises when the solution for the endemic equilibrium contains a radical (or square root). For 0 < R 0 << 1, the values of the endemic equilibria are imaginary, i.e. the number under the radical is negative, and are therefore ignored because biological population sizes are real-valued and non-negative. However, as R 0 increases, the term under the radical may switch from negative to positive giving rise to two real-valued, biologically relevant endemic equilibria. This type of equilibrium structure requires that the larger of the newly arisen equilibria is stable and the smaller is unstable. In the majority of backward bifurcations, as R 0 increases the smaller of the two endemic equilibria eventually intersects the DFE and creates a transcritical bifurcation (see Strogatz, 2014) at R 0 = 1 such that the DFE becomes unstable. Beyond R 0 = 1, the smaller endemic equilibrium is negative and is therefore not biologically realistic. Additionally, the positive, biologically relevant endemic equilibrium is always stable and the DFE is always unstable. Of course, there are more complex equilibrium structures that may not exactly follow this description, but the fact that R 0 = 1 is no longer a strict threshold for disease persistence always holds in backwards bifurcations.
One example of a system that exhibits this bifurcation structure is that of an Susceptible-Infectious-Susceptible (SIS) model with an imperfect vaccine (Kribs-Zaleta & Velasco-Hernandez, 2000) (see schematic in Figure 5). Assuming that the number of individuals in the vaccinated class is given by V where N = S+V +I is the total population size, λ(I) = β(I/N) is the force of infection to the susceptible class, λ V (I) = σβ(I/N) is the force of infection to the susceptible class, β is the transmission rate, μ is the natural birth/mortality rate, θ is the vaccine waning rate, σ is the ability of vaccinated individuals to produce new infections, γ is the recovery rate, and φ is the vaccination rate. When the rate of waning of immunity from β is the transmission rate, θ is the rate of immunity loss, σ represents the ability of vaccinated individuals to produce new infections, γ is the recovery rate, φ is the vaccination rate, and μ is the natural mortality rate.
the vaccine exceeds that of recovery rate, it is possible to have a stable endemic equilibrium while R 0 < 1. As found in Kribs-Zaleta and Velasco-Hernandez (2000), Varying R 0 through the variation of the vaccination rate φ yields different dynamics dependent on the value of σ , i.e. the ability of vaccinated individuals to produce new infections. Figure 6 shows an example of backward bifurcation for this system when σ = 0.02 which disappears when the vaccine is perfect, i.e. σ = 0. Backwards bifurcations have arisen in the mathematical analysis of a variety of epidemiological models (e.g. Arino, McCluskey, & van den Driessche, 2003;Dushoff, Huang, & Castillo-Chavez, 1998;Hadeler & van den Driessche, 1997;Martcheva & Thieme, 2003;van den Driessche & Watmough, 2000). However, often the constructions necessary to produce these dynamics require biologically tenuous assumptions such as in the model above where the rate of reinfection from the vaccinated state is greater than the rate of returning to the susceptible state. Another example is a model for Tuberculosis in which a backward bifurcation will occur when infection with re-exposure is more likely than infection with primary exposure (Feng, Castillo-Chavez, & Capurro, 2000;Martcheva, 2015). To our knowledge, however, there exists no real-world experimental evidence of the occurrence of backwards bifurcation in infectious disease systems.

Time varying parameters
The models we have discussed thus far have only considered constant-valued parameters. However, many of these quantities actually vary with time. Take, for example, the transmission rate, β. Some diseases have seasonal patterns, such measles and flu, which may be well represented by a periodically varying transmission rate resulting from, for example, changes in contact rates corresponding to when children are in school or on extended breaks. Over the course of a period, assuming frequency or density dependent transmission, the time varying transmission rate can be replaced, under certain conditions, by its average in the computation of the reproductive number (Wang & Zhao, 2008;Wesley & Allen, 2009). While R 0 for an entire period may be above one, there may be times when an instantaneous R 0 is actually below one. If the epidemic starts at this point in the transmission cycle with a single infectious individual, it can die out.

Conclusions
This paper provides an introduction to modeling infectious disease dynamics with a focus on several potential pitfalls. We attempt to show when these issues are likely to arise and how to proceed with models and analyses despite them. This paper provides an introduction to modeling infectious disease dynamics with a focus on several potential pitfalls. We attempt to show when these issues are likely to arise and how to proceed with models and analyses despite them. This paper is intended to serve as a supplement for the many excellent texts on infectious disease dynamics (referenced throughout) and it provides a concise reference that compiles a range of common issues in disease modeling.

Disclosure statement
No potential conflict of interest was reported by the authors. The exit rate from the I compartment, γ , can be found by measuring the average time spent infectious, d. Specifically, it turns out that γ = 1/d. To understand the origin of this result, solve Equation (A.1) with an initial condition of I(0) = I 0 to obtain the number of infectious individuals at time t Dividing both sides by I 0 allows us to rewrite this equation in terms of I(t)/I 0 , i.e. the fraction of the originally infectious individuals who are still infectious at time t. This is equivalent to the probability of remaining infectious at time t. The complement of this quantity, which we define as F(t; γ ) = 1 − exp(−γ t), is then the probability of having left the infectious class by time t. Given that the initial time is assumed to be at t = 0 (that is, there are no infections prior to the initial time), F(t; γ ) is exactly the cumulative distribution function for a probability distribution known as the exponential distribution where γ is the rate parameter. Any standard statistical text will demonstrate that the mean of the exponential distribution with parameter γ is exactly 1/γ , i.e. the average time an individual is infectious.

Appendix 2. SIR model in terms of population proportions
As discussed in the main text, the state variables in compartmental models can be written in terms of proportions of the population instead of the population density. This requires a change of variables, which we demonstrate here for the SIR model without demography (Equation (1)  A similar process applies for the I and R equations, and it can be shown that In more complex models, this change of variables can be applied in a similar way. We refer interested readers to Ledder (2017) to read more about non-dimensionalization of models in epidemiology.

Appendix 3. SIR model without demography
In Section 2.2, we introduced and analyzed the SIR model with demography. Interestingly, it turns out that if we eliminate demography by setting μ = 0 then the standard analysis of finding equilibria and R 0 no longer applies in the same way. First, recall the SIR model without demography: In this model, there is only a DFE. In fact, there are infinitely many DFEs because as long as I = 0 then the value of S and R does not matter. Suppose we try to find R 0 using, for example, the next generation matrix. In this case, the F matrix is simply a scalar, β(1/N)S and the V matrix is also a scalar given by γ (note that they are scalars because there is only one infectious class, I). Evaluated at the DFE, so that R 0 = β/γ . The magnitude of R 0 relative to one determines whether the disease initially grows (R 0 > 1) or declines (R 0 < 1). In contrast to the SIR model with demography, you always return to a disease free state. Although R 0 does not directly relate to the global asymptotic stability of the DFE, it does provide information on the shape of the epidemic curve. Namely, if R 0 < 1 then I is always decreasing. Consequently, there is no epidemic. If R 0 > 1, then I increases until S is equal to γ N/β. At that point, I decreases back towards zero. In this case, the dynamics follow a standard epidemic curveincreasing cases followed by a peak and eventually extinction of the outbreak. Therefore, although R 0 does not directly impact the global stability of the DFE it does play a substantial role in the short term epidemic dynamics.

Appendix 4. Linear stability analysis
We show how to perform a linear stability analysis for any two species model through the example of the SIR model with demography. Notice that since the total population size is constant we only need to use two equations to describe the system. Thus, our populations are represented by S and I so that with an equilibrium point denoted by (S * , I * ). Now, suppose that (x(t), y(t)) is a small perturbation away from the equilibrium such that Our goal is to determine how the perturbation (x, y) behaves as time increases, so we want to find the derivatives of x and y. Using the equalities above, this leaves us with dy dt = dI dt = F 2 (S, I) = F 2 (x + S * , y + I * ) We use the Taylor expansion to approximate these equations near the equilibrium. We now have two variables, so we need to find the Taylor expansion in two dimensions. Recall that we can approximate a function f (w, z) at a point (w * , z * ) using the following: where h.o.t. are higher order terms. We can apply this to Equations (A3) and (A4). Because the perturbations are assumed to be small, we ignore higher order terms and after some simplification we are left with the following linear approximations: dx dt = ∂F 1 (S * , I * ) ∂S x + ∂F 1 (S * , I * ) ∂I y dy dt = ∂F 2 (S * , I * ) ∂S x + ∂F 2 (S * , I * ) ∂I y We can alternatively write these approximations using slightly different notation: dx dt = ∂F 1 ∂S (S * ,I * ) x + ∂F 1 ∂I (S * ,I * ) y dy dt = ∂F 2 ∂S (S * ,I * ) x + ∂F 2 ∂I (S * ,I * ) y where the vertical bar indicates that the partial derivatives are being evaluated at the equilibrium. It is more convenient to write this in matrix form, i.e. the Jacobian. We can therefore rewrite the linear approximation as: ⎡ , n = x y Note that matrices will be indicated by bold font and vectors will have an arrow above them. We can now simply write Equation (A.5) as d n dt = J n where J is called the Jacobian matrix and here it is evaluated at the equilibrium. In the case of the SIR model the Jacobian is ⎡

S I
Note that the Jacobian can be extended to higher dimensions. As our goal is to determine the long term behaviour of the perturbation from the equilibrium, we need to understand how n will change. From calculus we know that in a one dimensional setting dn(t) dt = αn(t) =⇒ n(t) = n 0 e α The dynamics of this equation can then be determined by the sign of α. If it is positive, the perturbation grows with time. In contrast if it is negative, the perturbation shrinks with time. Now, if α represents the Jacobian evaluated around an equilibrium, these dynamics tell us about the stability of the equilibrium, i.e. when α < 0 the equilibrium is stable and when α > 0 the equilibrium is unstable.
A similar idea applies to higher order systems. Let's return to our two dimensional system Here, e λt is a non-zero constant so we can divide it out and we are left with λ n 0 = J n 0 (A8) Here, n 0 is called an eigenvector and λ is an eigenvalue of the matrix J. Details and further explanation on eigenvalues and eigenvectors can be found in any standard linear algebra or ordinary differential equation textbook. Keep in mind that our end goal is to determine how the dynamics change, which given Equation (A.6) depend on λ. To find λ, we subtract λ n 0 from both sides of Equation (A.8), leaving us with and we are left with 0 = J n 0 − λI n 0 = (J − λI) n 0 Notice that we cannot simply pull out n 0 from the above equation because J and λ have different dimensions. Therefore, we needed to introduce a 2x2 identity matrix I.
To have a non-zero solution of this equation the determinant of (J − λI) must equal zero. We further you to a standard ordinary differential equation text to find the determinant.
Thus, a general solution the two-dimensional system d n dt = J n is given by n = c 1 e λ 1 t n 1 + c 2 e λ 2 t n 2 where n 1 and n 2 are the eigenvectors corresponding to the eigenvalues λ 1 and λ 2 , respectively, and c 1 and c 2 are arbitrary constants.
In the case of the SIR model we find that the two eigenvalues of J are λ 1 = −μ and λ 2 = β − (μ + γ ). As our parameters are all positive to ensure biological realism, λ 1 < 0 always. However, λ 2 may be positive or negative dependent upon the relative sizes of β and μ + γ . The DFE equilibrium will be linearly stable when β − (μ + γ ) < 0. We can rearrange this quantity to find an equivalent comparison, i.e. β (μ + γ ) < 1 Notice, this is in fact our definition of R 0 . In fact, it can be shown that when R 0 < 1 then our eigenvalues are negative and the DFE is linearly stable and when R 0 > 1 then at least one eigenvalue is positive and the DFE is unstable.