The effect of the incidence function on the existence of backward bifurcation

In modelling, the dynamics of infectious disease, the choice of the specific mathematical formulation of disease transmission (i.e. the incidence function) is one of the initial assumptions to be made. While inconsequential in many situations, we show that the incidence function can have an effect on the existence of backward bifurcation (the phenomenon where a disease can persist even when the basic reproductive number is less than 1). More specifically, we compare mass action (MA) and standard incidence (SI) (the most common incidence functions) versions of two hallmark models in the backward bifurcation literature and an original combination model. Our findings indicate that the SI formation of disease transmission is more conducive to backward bifurcation than MA, a trend seen in all the models analysed.


Introduction
The basic reproductive number, denoted R 0 , is one of the most important quantities in the study of infectious disease. Defined as the average number of secondary infections caused by a single infection introduced into a completely susceptible population, R 0 is a key measure of a disease's potential to spread in a given population. In mathematical epidemiology, the basic reproductive number determines the stability of the disease-free equilibrium (DFE). Specifically, if R 0 > 1, the infection replaces itself and contributes additional infection leading to larger numbers of infections and an unstable DFE. If R 0 < 1, the infection does not replace itself and the DFE is stable.
In most situations, it is safe to overlook the 'introduced into a completely susceptible population' phrase in R 0 's definition to conclude that an average infection in general does (or does not) replace itself and consequently the disease will persist (or die out). This typical situation is realized through a forward, or supercritical, bifurcation at R 0 = 1 and is illustrated in Figure 1  However, having R 0 < 1 is not always enough to guarantee eradication of a given disease. Through a phenomenon known as backward, or subcritical, bifurcation depicted in Figure 1(b), it is possible for a disease to persist despite R 0 < 1. Here, the asymptotic stability of the DFE is identical to that of the forward bifurcation (i.e. asymptotically stable for R 0 < 1 and unstable for R 0 > 1). The difference is that the endemic equilibrium that appears at R 0 = 1 is unstable and exists for values of R 0 < 1. When backward bifurcation occurs, this unstable endemic equilibrium is typically accompanied by a larger asymptotically stable endemic equilibrium. In this situation, the initial amount of infection in the population determines whether the disease dies out or persists. Specifically, the disease persists if the amount of initial infection is above the unstable endemic equilibrium and dies out if it is below the unstable equilibrium.
While the existence of backward bifurcation is interesting from a purely mathematical perspective, the phenomenon has important practical consequences for public health as well. Summarizing points made by Brauer (2004) and Dushoff, Huang, and Castillo-Chavez (1998), these include: • For an endemic setting in which control measures are reducing the reproductive number, the condition that R 0 < 1 is not sufficient to eradicate the disease. Rather, the reproductive number must be reduced beyond an eradication threshold, which we denote R * 0 (see Figure 1(b)), in order to eliminate the disease. • For a setting in which R 0 increases across the threshold of R 0 = 1, the equilibrium prevalence is a discontinuous function of R 0 . As the reproductive number crosses 1 from below, the introduction of an arbitrarily small number of infected individuals results in a prevalence of P * > 0 (see Figure 1(b)).
While most of the literature on backward bifurcation has focused on simply proving the existence of the counterintuitive behaviour, this paper takes a deeper look into backward bifurcation by examining its sensitivity to the choice of incidence function. Our current work is motivated by that of Sharomi, Podder, Gumel, Elbasha, and Watmough (2007) which provides examples of HIV models with imperfect vaccination that demonstrate backward bifurcation when using a standard incidence (SI) formulation but do not exhibit the behaviour with an equivalent mass action (MA) formulation. In this paper, we test the generality of their observation that SI favours the existence of backward bifurcation over MA incidence by establishing analytic threshold conditions for both SI and MA formulations of the models in two hallmark backward bifurcation papers and comparing the resulting conditions. While many types of incidence functions exist McCallum, Barlow, and Hone (2001), we focus our attention on SI (also referred to as frequency-dependent transmission) and MA (also referred to as density dependent) as they are the two common most. In a SI formulation, the force of infection is given by βS(t)I(t) N(t) , where S(t) and I(t) denote the number of susceptible and infectious individuals, respectively, N(t) is the total population size and β is a transmission coefficient that includes the average number of contacts that an infectious individual makes per unit time and the proportion of those contacts that result in disease transmission. Biologically, this reflects a situation where each infectious individual would infect β people per unit time but only a proportion S(t) N(t) of their contacts will be with individuals that are currently susceptible to the disease. Hence, the rate of infection depends on the proportion of the population that is currently susceptible. In the MA formulation, the force of infection is given by σ S(t)I(t) (where σ is the transmission coefficient) and each infectious individual transmits the disease to σ S(t) people per unit time. Consequently, the rate of transmission depends on the number of susceptible individuals in the population rather than the proportion. Notably, the SI and MA formulations are mathematically equivalent if the size of the total population is constant (i.e. N(t) = N * ) as the constant population size can be absorbed into the transmission coefficient (i.e. σ = β N * ). The structure of our paper is as follows: in Section 2, we examine the effect of the incidence function in the backward bifurcation in the model of Feng, Castillo-Chavez, and Capurro (2000); in Section 3, we do the same for the model of Kribs-Zaleta and Velasco-Hernández (2000); and in Section 4, we analyse a model that combines aspects of Feng et al. (2000) and Kribs-Zaleta and Velasco-Hernández (2000). Throughout this work, we rely on the general theorem for the existence of backward bifurcation of Castillo- Chavez and Song (2004) which is closely related to that of van den Driessche and Watmough (2002). The theorem proves the existence of backward bifurcation using the centre manifold theory to determine the local behaviour of the system at the DFE when R 0 = 1. A slightly modified form of their result is given in Theorem 1 of the Supplementary Material (SM) for direct applicability to our systems.

Backward bifurcation caused by exogenous reinfection
In Feng et al. (2000), the authors present a model for the epidemiology of tuberculosis that includes exogenous reinfection (the characteristic of TB where a latently infected individual can acquire a new infection from another infectious individual) and showed that exogenous reinfection can produce backward bifurcation. We revisit this hallmark work in the backward bifurcation literature to examine the effect of the incidence function on the existence of backward bifurcation.

The Feng et al. model
A slightly generalized form of the model proposed in Feng et al. (2000) is given by   Feng et al. (2000), Blower et al. (1995) k Rate of slow progression .005 .001-.008 Feng et al. (2000), Blower et al. (1995)  where N = S + E + I + T is the total population size. The model divides the population into four classes depending on disease state; S (susceptible), E (latent stage of disease), I (infectious stage of disease), and T (treated, but still susceptible). The model utilizes a constant recruitment rate, , and a per-capita disease induced rate of mortality, d. The rate of progression from latency to infectivity is represented as k, the per-capita natural mortality rate of the population as μ, the per-capita treatment rate as r, the factor of susceptibility after treatment as σ and the factor of susceptibility during latent infection as p. As p determines the level of exogenous reinfection, Feng et al. use p as a threshold parameter for the existence of backward bifurcation (i.e. How much exogenous reinfection is required to produce backward bifurcation?). In Feng et al. (2000), the authors use a SI force of infection given by λ(X, I, N) = cβ XI N for X ∈ {S, E, T}, where c is the average number of contacts per unit time and β is the average transmission probability of a single contact with a susceptible individual. For simplicity, we use the single term β as the overall transmission coefficient without loss of generality to formulate the force of infection as:

mass action version)
A flow diagram for the Feng et al. model is given in Figure 2 and a summary of model parameters is given in Table 1.

Backward bifurcation threshold; SI Feng et al. model
We begin by establishing the condition for the existence of backward bifurcation for the Feng et al. model as proposed in Feng et al. (2000) (i.e. the SI formulation) which is given by System (1) with λ(X, I, N) = β XI N , for X ∈ {S, E, T}. We use the general theorem for the existence of backward bifurcation of Castillo-Chavez and Song (2004) which was restated in Theorem 1 in the SM in a modified form for more direct applicability to our systems. The basic reproductive number, R 0 , of the model is given by the product of the transmission rate and average duration of infection, 1/(μ + r + d), multiplied by the probability that a new infection progresses from latency to active infection before death k/(μ + k); resulting in As backward bifurcation describes the local behaviour of the DFE at the reproductive number threshold, we begin our analysis by finding the DFE and β * , the transmission rate for which R 0 = 1 to get Linearizing around the DFE, we calculate the Jacobian of the system at the DFE to get As we are concerned with the system's behaviour near R 0 = 1, we evaluate the Jacobian at β = β * and calculate its eigenvalues to get λ 1,2,3, The left and right eigenvectors, v and w, respectively, corresponding to the eigenvalue λ = 0, are then found to be In accordance with the conditions of Theorem 1 in the SM, we have that v · w = 1. With this condition met, we rewrite our system in the form x = f ( x) for a clearer application of Theorem 1 in the SM by renaming the state variables as S = x 1 , E = x 2 , I = x 3 , and T = x 4 . Our system becomes From this reformulated system, we may directly calculate a and b as described in the Theorem 1 in the SM as As the parameters involved in the expression for b are all positive when assigned biologically feasible values, we have that b > 0. Therefore, the existence of backward bifurcation depends entirely on the sign of a. Noting that the necessary condition for the existence of backward bifurcation from Theorem 1 in the SM is that a > 0, we see from the form of (3) that backward bifurcation cannot occur without exogenous reinfection (i.e. p = 0). More specifically, we can use (3) to establish the level of exogenous reinfection necessary to produce backward bifurcation for the SI formulation of the Feng et al. model: Notably in Feng et al. (2000), explicit conditions for the existence of backward bifurcation are only found for a simplified model with no disease-induced mortality and where reinfection after recovery occurs at the same rate as for susceptible (i.e. d = 0, σ = 1). With these simplifications, the threshold p * SI = k(μ+k+r) μ(μ+r) is identical to that stated in Feng et al. (2000).

Backward bifurcation threshold; MA Feng et al. model
In order to examine the effect of MA vs SI on the existence of backward bifurcation for the Feng et al. model, we perform the same mathematical analysis on the MA version of the model which is given by System (1) with λ(X, I, N) = βXI, for X ∈ {S, E, T}.
By considering the case where the infectious population is zero and solving for each of the state variables, we find the DFE to be The basic reproductive number of the MA version of the Feng et al. model is found by multiplying that of the SI version by the total population size at the DFE (as N does not appear in the denominator of the MA term). Consequently, the basic reproductive number and β * (the transmission rate for which R 0 = 1) are given by To linearize the system around the DFE, we evaluate the the Jacobian at the DFE to get The eigenvalues of the Jacobian evaluated at β = β * are then given by λ 1,2,3, The left and right eigenvectors, v and w, respectively, corresponding to λ = 0, are then given by Noting that v · w = 1, we rename the state variables as S = x 1 , E = x 2 , I = x 3 and T = x 4 to reformulate our system as Calculating a and b from our differential system, according to the Theorem 1 in the SM, we have that As was the case in the model formulated for SI, b > 0 for all biologically reasonable parameters, and thus the existence of backward bifurcation hinges solely upon a > 0. Using this condition, we establish the following threshold level of exogenous reinfection necessary for the MA version of the Feng et al. model to exhibit backward bifurcation:

Comparison of backward bifurcation thresholds
From the SI and MA thresholds in (5) and (10), we immediately see that exogenous reinfection is the mechanism causing backward bifurcation. Without exogenous  Table 1. The threshold rates of exogenous reinfection for the existence of backward bifurcation evaluate to p * SI = .3224687 and p * MA = .3270838 for the standard incidence and mass action versions, respectively. Notes: Solid (dashed) curves indicate stable (unstable) equilibria. Green (red) curves show endemic equilibria for the standard incidence (mass action) version. Behaviour of the disease-free equilibrium (shown in black) is the same for standard incidence and mass action versions. reinfection (i.e. p = 0), backward bifurcation is not exhibited by either the SI or MA system. While the results of Feng et al. (2000) showed that backward bifurcation can occur for biologically feasible parameter values in a model for tuberculosis, the question of whether the required rates of exogenous reinfection are biologically realistic has been debated (see Lipsitch & Murray, 2003;Wang, Feng, Aparicio, & Castillo-Chavez, 2014). To some degree, its unintuitive nature has led to questions about the realism of model assumptions/parameter values required to produce backward bifurcation. Here, we note that a model's formulation of disease transmission (i.e. choice of incidence function) affects the threshold conditions for backward bifurcation and consequently makes the behaviour more (or less) biologically realistic.
A numerical illustration of the backward bifurcation phenomenon is shown in Figure 3. Evaluating (5) and (10) using the baseline parameter values in Table 1, we have that the critical levels of exogenous reinfection necessary to produce backward bifurcation are p * SI = .3224687 and p * MA = .3270838 in the SI and MA versions, respectively. In Figure 3(a), we see that neither version of the model exhibits backward bifurcation when the exogenous reinfection rate, p, is such that p < p * SI < p * MA . When p * SI < p < p * MA as in Figure 3(b), we see that the standard version of the model (green curve) exhibits backward bifurcation but the MA version (red curve) does not. When p * SI < p * MA < p as in Figure  3(c), both versions exhibit the phenomenon.
An initial comparison of backward bifurcation thresholds for the Feng et al. model, shows that the MA formulation always requires a higher rate of exogenous reinfection to produce backward bifurcation than a corresponding SI formulation. As the SI formulation relaxes the necessary condition on the required rates of exogenous reinfection (the central question regarding the realism of backward bifurcation), we say that the SI formulation 'favors the existence of backward bifurcation' or 'makes backward  (12)). Notes: In each plot, a single parameter is varied over its range in Table 1. All others are fixed at baseline values. bifurcation more likely to occur' and we reiterate that this does not mean that SI causes backward bifurcation.
From (11), we see that thresholds, p * MA and p * SI , are equivalent if either k = 0 or d = 0. The situation of k = 0 is not interesting because without slow progression from latent to active TB, the infection is guaranteed to go extinct (note from (2) that k = 0 =⇒ R 0 = 0, because the initial infection introduced into a susceptible population cannot progress from latent to active TB). When there is no disease-induced mortality (i.e. d = 0), the total population is constant. As was noted in Section 1, the MA and SI models are equivalent because constant population size in the SI can be absorbed into the transmission coefficient β. Therefore, it follows that the threshold for backward bifurcation must be identical. From (11), we see that amount by which SI favours backward bifurcation is an increasing function of the rate of slow progression (k) and a decreasing function of the natural mortality rate (μ) and the recovery rate (r). Differentiating (11) with respect to d to get k 2 μ+r μ 2 μ+r+d have that the difference in thresholds also increases as a function of the disease-induced mortality rate.
To gain quantitative insight into the amount by which SI favours backward bifurcation for the Feng et al. model, we consider the ratio of the backward bifurcation thresholds: Evaluating at the baseline model parameters from Table 1, we have that p * MA /p * SI = 1.014 showing that the threshold level of exogenous reinfection is only about 1.4% higher for the MA model. In Figure 4, we see that ratio of backward bifurcation thresholds is significantly more sensitive to the disease-induced mortality rate (d) and the recovery rate (r) than other variables. For biologically reasonable parameters, we see that MA can require up to ≈ 6% faster rates of exogenous reinfection than the corresponding SI formulation.

Backward bifurcation caused by imperfect vaccination
In the previous section, we showed that the SI formulation of the Feng et al. model is more conducive to backward bifurcation (i.e. requires a lower level of exogenous reinfection) than the equivalent MA formulation. To see if this relationship holds when backward bifurcation is caused by a different mechanism (other than exogenous reinfection as in Feng et al. (2000)), we now investigate the model proposed by Kribs-Zaleta and Velasco-Hernández (2000). In their model, which we will refer to as the Kribs-Zaleta et al. model, backward bifurcation is caused by a partially protective vaccine.

The Kribs-Zaleta et al. model
The model proposed in Kribs-Zaleta and Velasco-Hernández (2000) is given by where N = S + I + V is the total population size. The model divides the population into three classes; S, I, and V , representing susceptible, infectious, and vaccinated individuals, respectively. Model parameters include the rate of vaccination, γ , the factor of protection conferred by vaccination, ψ, the waning rate of the vaccine's protection, θ, and the natural rate of mortality, μ. The model analysed in Kribs-Zaleta and Velasco-Hernández (2000) assumed no disease-induced mortality and recruitment into the population at a rate identical to that of mortality (i.e. μN) which produces a constant population size. In the last paragraph of Section 3 of Kribs-Zaleta and Velasco-Hernández (2000), the authors pose a modified model with a constant recruitment rate (which we refer to as for consistency with the Feng et al. model) and disease-induced mortality at a rate of d. As the SI and MA formulations are equivalent when the total population size is constant, it is this   Figure 5 and a summary of the parameter values we use is given in

Comparison of thresholds; Kribs-Zaleta et al. model
Again, we rely on the general theorem of Castillo-Chavez and Song (2004) Table 2. The threshold rates of recovery for the existence of backward bifurcation evaluate to r * SI = .4819918 and r * MA = 4.219441 for the standard incidence and mass action versions, respectively. Notes: Solid (dashed) curves indicate stable (unstable) equilibria. Green (red) curves show endemic equilibria for the standard incidence (mass action) version. Behaviour of the disease-free equilibrium (shown in black) is the same for standard incidence and mass action versions.
Feng et al. model, we move directly to the resulting thresholds. The full calculations are provided in Section 2 of the SM. The resulting thresholds are given by and A numerical illustration of the backward bifurcation phenomenon is shown in Figure 6. Evaluating (14) and (15) using the baseline parameter values in Table 2, we have that the critical rates of recovery necessary to produce backward bifurcation are r * SI = .4819918 and r * MA = 4.219441 in the SI and MA versions, respectively. In Figure 6(a), we see that neither version of the model exhibits backward bifurcation when the exogenous reinfection rate, r, is such that r < r * SI < r * MA . When r * SI < r < r * MA as in Figure 6(b), we see that the standard version of the model (green curve) exhibits backward bifurcation but the MA version (red curve) does not. When r * SI < r * MA < r as in Figure 6(c), both versions exhibit the phenomenon.
In both thresholds, we see that backward bifurcation requires partial protection conferred by vaccination (i.e. ψ = 0, 1) and a sufficiently fast recovery rate (r). Given that these are the mechanisms causing the behaviour, we proceed to find which form of the incidence function is more conducive to backward bifurcation (i.e. which incidence function produces a less restrictive threshold on r.) Considering the difference between the two thresholds, we see that the MA formulation is greater than the SI formulation, requiring a greater recovery rate in order to produce backward bifurcation than the SI formulation of the model. As with the Feng et al. model, the MA formulation can be said to be more  (17)). Notes: In each plot, a single parameter is varied over its range in Table 2. All others are fixed at baseline values.
restrictive towards backwards bifurcation, and that the SI formulation favours the existence of backward bifurcation. From (16), we see that the thresholds are identical when d = 0. However, as stated before, the total population size is constant without disease-induced mortality, rendering the SI and MA formulations equivalent.
For a more quantitative analysis, we consider the ratio of the thresholds given by Evaluating at the baseline parameter values from Table 2, we have that r * MA /r * SI = 7.924, showing that backward bifurcation for the MA formulation requires recovery rates nearly eight times as large as the SI version. Hence, a SI formulation of the Kribs-Zaleta et al. model is much more likely to exhibit backward bifurcation than an equivalent MA version. Figure 7 shows the ratio of thresholds when individual parameters are varied through their range in Table 2. We see that ratio of backward bifurcation thresholds is by far most sensitive to disease-induced mortality, d, and that the ratio is a decreasing function of both the natural morality rate, μ, and the waning rate of the vaccine's protection, θ. Interestingly, the ratio of the thresholds is non-monotonic as a function of both the vaccine rate, γ , or the level of vaccine protection, ψ, and exhibits a maximum value within the biologically feasible range.

Combining exogenous reinfection and imperfect vaccination
To this point, we have shown that the SI formulation of disease transmission is more conducive to backward bifurcation than the MA formulation for two hallmark models in the backward bifurcation literature. While the simple vaccination model of Kribs-Zaleta et al. is quite general and not designed for any particular disease, it is worth noting that its main assumptions (a vaccine that offers imperfect protection that wanes with time and an infectious disease that does not confer immunity after recover) are completely applicable in the case of tuberculosis. Hence, combining aspects of the Feng et al. and Kribs-Zaleta et al. models produces a model that remains relevant to the dynamics of TB infection. We now proceed to examine the effect of the choice of incidence function on the existence of backward bifurcation for this new model that contains two characteristics that are associated with backward bifurcation (exogenous reinfection and imperfect vaccination).

The combination model
Incorporating aspects of the Feng et al. and Kribs-Zaleta et al. models, our combination model for the dynamics of TB is given by where N = S + V + E + I + T is the size of the total population. The model divides the population into five classes; S, V , E, I and T, representing susceptible, vaccinated, latently infected, infectious and treated (but still susceptible) individuals, respectively. Model parameters include a constant recruitment rate, , the natural mortality rate, μ, per-capita disease-induced rate of mortality, d, the rate of vaccination, γ , the factor of protection conferred by vaccination, ψ, the waning rate of the vaccine's protection, θ, rate of progression from latency to infectivity, k, the factor of susceptibility during latent infection as p, per-capita treatment rate, r, and factor of susceptibility after treatment as σ .
The force of infection is formulated as: where β is the transmission coefficient. A flow diagram for the combination model is given in Figure 8 and a summary of model parameters is given in Table 3.

Comparison of thresholds; combination model
In Section 3 of the SM, we use Theorem 1 in the SM to establish conditions for the existence of backward bifurcation for the SI and MA versions of the combination model. To compare the thresholds of the combination model to those of the previous models, we examine both the exogenous reinfection rate, p, (as in Section 2.4 for Feng et al. model) and the rate of recovery, r (as in Section 3.2 for the Kribs-Zaleta et al. model). As the calculations for the   Feng et al. (2000), Blower et al. (1995) k Rate of slow progression .005 .001-.008 Feng et al. (2000), Blower et al. (1995) where c 1 , . . . , c 4 are positive terms (see (14) and (19) in Sections 3.1 and 3.2 of the SM for details). As the c 2 r term is positive, we see that a high rate of recovery r is no longer sufficient (as was the case for the Kribs-Zaleta et al. model) to produce backward bifurcation for the combination model. Rather, it is the rates of reinfection, both after recovery (σ ) and  Table 3. The threshold rates of exogenous reinfection for the existence of backward bifurcation evaluate to p * SI = .1690434 and p * MA = .1719923 for the standard incidence and mass action versions, respectively. Notes: Solid (dashed) curves indicate stable (unstable) equilibria. Green (red) curves show endemic equilibria for the standard incidence (mass action) version. Behaviour of the disease-free equilibrium (shown in black) is the same for standard incidence and mass action versions. exogenous (p), that drive backward bifurcation. For this reason, we focus our attention on the threshold for exogenous reinfection, p, necessary to produce the behaviour.
A numerical illustration of the backward bifurcation phenomenon is shown in Figure 9. Evaluating (15) and (20) from the Supplemental Material using the baseline parameter values in Table 3, we have that the critical levels of exogenous reinfection necessary to produce backward bifurcation are p * SI = .1690434 and p * MA = .1719923 for the SI and MA versions, respectively. In Figure 9(a), we see that neither version of the model exhibits backward bifurcation when the exogenous reinfection rate, p, is such that p < p * SI < p * MA . When p * SI < p < p * MA as in Figure 9(b), we see that the standard version of the model (green curve) exhibits a very small backward bifurcation and the MA version (red curve) does not. When p * SI < p * MA < p as in Figure 9(c), both versions exhibit the phenomenon.
While the exact expressions for the threshold rates of exogenous reinfection are unwieldy (see (15) and (20) in the SM), the difference between the MA and SI thresholds simplifies to Therefore, we have that the MA formulation always requires a higher rate of exogenous reinfection to produce backward bifurcation than the corresponding SI formulation. To quantify the amount by which SI favours backward bifurcation for the combination model, we consider the ratio of the backward bifurcation thresholds, p * MA /p * SI , where p * MA and p * SI are (15) and (20) in the SM, respectively. Evaluating at the baseline model parameters from Table 3, we have that p * MA /p * SI = 1.017. In Figure 10, we illustrate the effect of varying each individual model parameters over its range in Table 3. As with the Feng et al. model, we see that disease-induced mortality rate, d and the recover rate, r, have the largest impact on the ratio of the MA and SI backward bifurcation thresholds. As disease-induced mortality increases, we have that MA threshold increases from being the same as the SI threshold to about higher than the SI threshold. As the recovery rate Figure 10. The ratio of the levels of exogenous reinfection required to produce backward bifurcation in the mass action and standard incidence formulation of the combination model (solid curves). Notes: Corresponding ratios of thresholds for the Feng et al. model (dashed curves) are included for comparison when available. In each plot, a single parameter is varied over its range in Table 3. All others are fixed at baseline values.
increases, the MA threshold decreases from about 6% higher down to about 1.5% larger than the SI threshold.
Lastly, we compare the results of the combination and Feng et al. models (solid and dashed curves in Figure 10, respectively) to examine the effect of adding imperfect vaccination into the Feng et al. model. Overall, the general behaviour of the two models is very similar but we do see that the ratios are indeed larger for the combination model. Therefore, we see that the relative restrictiveness of the MA thresholds vs the SI thresholds is increased by incorporating the additional model complexity.

Discussion
In all of the models we have considered, we have found that the SI formulation of disease transmission favours backward bifurcation as opposed to the MA formulation. More precisely, we found that the existence conditions for backward bifurcation are always less restrictive for SI than for MA. For both the Feng et al. and combination models, the MA formulation was only modestly more restrictive, requiring exogenous reinfection rates to be no more than 8% larger than the SI formulation. For the Kribs-Zaleta et al. model, the effect is huge as the MA formulation requires recovery rates up to an order of magnitude larger than those of the SI formulation. For all of the models considered, we found that the ratio of the MA and SI thresholds was most sensitive to the disease-induced mortality rate, d, and that the ratio was an increasing function of disease-induced death. For the Feng et al. and combination models, the threshold ratio decreases significantly as the recovery rate, r, increases.
It must be stated, however, that our results only hold when one considers the case of a dynamic population. If the total population is constant, N * , the SI formulation β 1 SI N can be seen to be equivalent to the MA formulation β 2 SI by setting β 2 = β 1 N * (i.e. by absorbing the constant population size N * into the transmission coefficient). In the models considered in this work, the size of the total population changes in response to disease dynamics because we have considered only fatal diseases (i.e. those with disease-induced mortality). Hence, the sensitivity of our results to disease-induced mortality seems appropriate.
In this work, we have found that the SI formulation of disease transmission favours the existence of backward bifurcation compared to a corresponding MA formulation. Like much of the backward bifurcation literature, we have focused on the mathematical details of the phenomenon rather than biological explanation. To conclude our work and offer some intuition regarding our findings, we propose thinking of backward bifurcation in terms of the infectious potential of an average infected individual, which we refer to as Patient A. In the standard case of forward bifurcation, Patient A's infectious potential is maximized when the rest of the population is uninfected (i.e. Patient A's maximum infectious potential is equal to R 0 ). If other individuals besides Patient A are infected, the pool of uninfected is depleted and Patient A's infectious potential is reduced. In the situation of backward bifurcation, the presence of other infected individuals in the population can actually increase the infectious potential of Patient A. This occurs if the depletion of the uninfected pool is somehow counteracted by a change in the structure of the population that increases Patient A's infectious potential (Dushoff et al., 1998).
Mathematically, we can make sense of the main observation of this work by examining the structure of the SI and MA forms of the incidence function. When using the MA formation of incidence, Patient A's rate of infecting others is βS. If individuals other than Patient A are infected, the number of susceptible individuals is reduced, thus reducing the infectious potential of Patient A. If SI were used instead, Patient A's rate of infecting others is β S N . When additional individuals are infected, the number of susceptibles is reduced, but the disease-induced mortality of those infections also reduces the size of the total population. This reduction in the denominator of the incidence term, means that having additional infection in the population has a smaller effect (possibly only slightly) on the infectious potential of Patient A. Thus, the conditions for the depletion of the uninfected pool to be counteracted by a change in population structure will be less restrictive (possibly only slightly). This would confirm our result that for a dynamic population, backward bifurcation appears more readily in a model utilizing the SI form of the incidence function, as opposed to the MA form.
Despite being the foundation of any mathematical model of infectious disease, the choice of incidence function is often made with insufficient consideration and justification. In many circumstances, the behaviour of interest is insensitive to this decision. In this work, we have shown if one is interested in studying backward bifurcation, it is very important to carefully choose and justify the formulation of disease transmission that most realistically represents the situation being modelled.