Bistability in deterministic and stochastic SLIAR-type models with imperfect and waning vaccine protection.

Various vaccines have been approved for use to combat COVID-19 that offer imperfect immunity and could furthermore wane over time. We analyze the effect of vaccination in an SLIARS model with demography by adding a compartment for vaccinated individuals and considering disease-induced death, imperfect and waning vaccination protection as well as waning infections-acquired immunity. When analyzed as systems of ordinary differential equations, the model is proven to admit a backward bifurcation. A continuous time Markov chain (CTMC) version of the model is simulated numerically and compared to the results of branching process approximations. While the CTMC model detects the presence of the backward bifurcation, the branching process approximation does not. The special case of an SVIRS model is shown to have the same properties.

Introduction

According to the World Health Organization World Health Organization (2021b), at the time of writing, the total number of confirmed cases of Coronavirus Disease (COVID-19) has exceeded 170 million, about a year and a half after the disease and the virus that causes it (SARS-CoV-2) were officially named. Among the confirmed cases are more than 3.7 million deaths World Health Organization (2021b).

Of particular importance in the fight against COVID-19 was the development of vaccines. It is a testament to the amount of effort that went into this endeavour that, barely one year after the start of the crisis, at least 7 different vaccines have been administered World Health Organization (2021a). Three vaccines, mRNA-based vaccines from Pfizer-Biontech and Moderna along with a viral vector vaccine from Janssen, have been granted emergency use status in the United States United States Centers for Disease Control and Prevention (2021). Worldwide, several other vaccines are available at the time of writing: Astra-Zeneca, Sputnik, Sinovacc, Covishield, according to World Health Organization (2021a). All these vaccines have different characteristics, which are best summarized using three features: the number of doses required to generate immunity, the level and type of immunity provided and the duration of the protection.

While the number of doses is an important feature, for COVID-19, it has this far generated specific issues mostly at the time of vaccine roll-out. The second characteristic of a vaccine is its efficacy: a vaccine does not always confer full immunity to the disease, in which case is is called imperfect. Such vaccine imperfections have different sources, but in with COVID-19, one of the main reason is the emergence of SARS-CoV-2 variants less sensitive to some of the vaccines. Finally, it is possible that the immunity provided by the vaccine does not last: the vaccine is said to wane. Clearly, the three characteristics are intimately linked: vaccines against some diseases require booster shots because the protection they afford is known to wane; efficacy may diminish because a full dose regimen is not followed, etc.

Vaccination for COVID-19 is presenting heretofore unseen challenges. It is the first time in the history of vaccination that the cycle from development to global worldwide roll-out occurs over such a short time period. Vaccination against poliomyelitis and smallpox were on the same scale, but much evidence had been collected about efficacy and waning over the years. Vaccination against the 2009 pandemic H1N1 influenza strain became possible very soon but a lot was already known about influenza vaccines because of annual vaccination campaigns.

Another facet of the fight against COVID-19 has been the unprecedented reliance on mathematical models to study critical scientific questions about the dynamics of spread of the disease and thereby help guide public policy. Even in 2020, there were quite a lot of mathematical models published to study COVID-19. See, e.g., Mohamadou et al. (2020), Ogden et al. (2020) and the references therein. This has only increased since; see, e.g., the extensive review in Cao and Liu (2021), Xiang et al. (2021) for models regarding public health interventions or Arino (2022) for models related to the spatial spread of the disease. It is interesting to note that many of the mathematical models used to study SARS-CoV-2 and COVID-19 have their roots in the classical Susceptible-Infected-Removed (SIR) model of Kermack and McKendrick (1927). Some authors have adapted the classical SIR structure; see, e.g., Cooper et al. (2020), Nguemdjo et al. (2020). Others have extended the model to include a latent or exposed compartment, such as Chen et al. (2020), Yang and Wang (2020), while others have gone further to include both latent/exposed and asymptomatic compartments, as in Arino et al. (2020), Arino and Portet (2020), Basnarkov (2021), Li et al. (May 2020), Rădulescu et al. (2020), Tsay et al. (2020). We refer to the general class of models that extend the classical SIR model structure to include latent/exposed (L) and asymptomatic (A) compartments as SLIAR-type models. SLIAR-type models were first developed to study H5N1 influenza by Arino et al. (2006), Longini et al. (2005) combining characteristics of models for SARS-CoV suggested by Brauer (2006). Of importance to justify our work here is that most models have made the (completely justified) assumption that considering the vital dynamics (demography) of the population was not required. However, as the crisis drags on and there is more and more talk about vaccination becoming recurrent, the validity of this hypothesis is coming under question.

This is all the more important that, because of uncertainties about COVID-19 vaccine characteristics, the use of mathematical modeling to investigate the potential effects of vaccines is a necessity. However, it has been known for some time that in models with demography, the compounding of imperfection and waning of the protection impacts the dynamics of the spread of the disease, sometimes leading to so-called backward bifurcation situations as in Hadeler and van den Driessche (1997), which were first exhibited as a result of vaccination in an SIS model in Kribs-Zaleta and Velasco-Hernández (2000) and in an SIRS model in Arino et al. (2003). In a backward bifurcation scenario, there can be multiple endemic equilibria, often occurring in subcritical situations, i.e., when the reproduction number is less than unity. This greatly complicates the mathematical analysis; Arino et al. (2003) is the only work to our knowledge where global conclusions are drawn in a backward bifurcation region. The phenomenon of backward bifurcations has been shown for a variety of epidemic models, see Villavicencio-Pulido et al. (2015). The existence of a backward bifurcation also has important implications more generally for the dynamics of disease spread and, ultimately, control, since in a backward bifurcation region, solutions are initial-condition dependent. Therefore, the first question we investigate is the existence of a backward bifurcation in a deterministic SLIARS model with demography and vaccination, i.e., a deterministic SVLIARS model.

Our second question concerns the link between the bistable behavior of a deterministic SVLIARS and its related continuous time Markov chain analogue. Indeed, another important teaching of the COVID-19 crisis is the importance of case importations (or introductions in the ecological lingo) in the global spread of infectious diseases. Highlighted in Arino et al. (2020), Arino et al. (2021), this phenomenon is best described using a stochastic model. While the effect of importations was of course known before, the amount of scrutiny and the wide availability of data had made studying this aspect possible. In the context of COVID-19 and more generally, the emergence or re-emergence of a disease, case importations and the initial spread in communities involve a relatively small number of infectious individuals, leading to disease dynamics driven by randomness. Randomness is also present when there is a degree of variability associated with the mechanisms driving the disease dynamics such as transmission or recovery or if environmental randomness affects disease outcomes. Stochastic epidemic models date back to Bernoulli (1760) and have been studied in earnest at least since Kermack and McKendrick (1927). An important feature of these models is that they allow to compute the probability of a minor epidemic and its complement, the probability of a major epidemic. These probabilities were calculated for continuous time Markov chain (CTMC) SI, SIR models by Whittle (1955). CTMC and stochastic differential equation (SDE) models related to deterministic models exhibiting backward bifurcation have been shown to detect bistability Allen and van den Driessche (2006). Relying on knowledge of the bifurcation structure of the corresponding ODE model, we will investigate the capacity to detect bistability in the SVLIARS CTMC.

When direct analysis of CTMC models is intractable, it is necessary to approximate the probabilities of minor and major epidemics. One common approximation technique is to construct and analyze a branching process model. Branching process approximation of the probability of a minor epidemic has been used to great effect Allen and Bokil (2012), Allen and van den Driessche (2013), Edholm et al. (2018) and has been related to the basic reproduction number of a corresponding deterministic model Allen and van den Driessche (2013). However, this technique is a linearization and has been shown to be less successful for highly nonlinear models and away from the disease-free equilibrium in Milliken (2017). Another way to approximate the probability of a minor epidemic in a CTMC model is to make calculations from an ensemble of sample paths simulated via Gillespie algorithm. This brings us to the third question considered in this work, namely the consideration of whether the CTMC and the BP approximation are equally able to detect a backward bifurcation.

To summarize, in this work we investigate the dynamical behavior of CTMC models in the presence of a backward bifurcation and whether this behavior can be detected using branching process approximation techniques. To this end, we formulate an ODE SLIAR-type model with vital dynamics that includes imperfect vaccination and waning immunity and establish conditions for the existence of a backward bifurcation in a biologically relevant parameter regime. A related continuous time Markov chain (CTMC) model is then presented, and branching process approximation and simulation via the Gillespie algorithm are used to approximate the probability of a minor epidemic. By way of example, we show that branching process approximation does not detect the change in dynamics in the CTMC model associated to the backward bifurcation. In an appendix, we consider the same properties in the special case of an SVIRS model.

The deterministic SVLIARS model

Model formulation

The model is an SLIAR-type model to which both demography and vaccination are added, i.e., essentially a cross between the SVIRS model of Arino et al. (2003) and the SLIARS model of Arino et al. (2006). We interpret the R compartment here as containing individuals who are immune (potentially for a limited time) to the disease because of disease-acquired immunity, rather than as removed individuals as is often done in SLIAR-type models. See Table 1 for a list of state variables used in the model and Table 2 for a list of parameters. The flow diagram of the system takes the form shown Fig. 1.

Table 1

Description of model variables

Variable	Description
S	Number of susceptible individuals
V	Number of individuals with vaccine induced partial immunity
L	Number of latently infected individuals
I	Number of symptomatic infected individuals
A	Number of asymptomatic infected individuals
R	Number of individuals with disease induced immunity

Table 2

Description of model parameters

Parameter	Description
Demography
B	Natural birth rate (set to dN_0)
d	Natural mortality rate
Disease characteristics
β	Rate of transmission
η	Modification of transmission for asymptomatics
π	Proportion of asymptomatic cases
ε	Rate of transition from latency to infectious stage
γ	Rate of recovery
ω_r	Rate of waning of disease induced immunity
μ	Disease induced mortality rate
Vaccination
p	Proportion of newborns vaccinated
e	Rate of vaccination of adults
σ	Vaccine efficacy
ω_v	Rate of waning of vaccine induced immunity

Fig. 1

Flow diagram of the SVLIARS model (1). Demography flows are shown using dashed lines

The model has flow diagram shown in Fig. 1 and incorporates imperfect and waning vaccine immunity as well as waning disease-induced immunity. While it lasts, disease-induced immunity is perfect, contrary to the vaccine-induced one. Disease induced mortality is assumed to only occur for symptomatic infected individuals.

System (1) is considered with nonnegative initial conditions. To avoid trivial solutions, it is generally assumed that .

Preliminary analysis of the SVLIARS model

Before proceeding further, let us briefly consider properties of model (1) without vaccination, i.e., (1) without equation (1b) and with . This system has disease-free equilibrium () and using the same method as used later, one finds the basic reproduction numberThis quantity is useful to understand the effect of vaccination. We now turn to the system with vaccination (1).

Proposition 1

The disease-free equilibrium of system (1) iswhere and are given by

Proof

Suppose there is no disease, i.e., . Then from (1f), . So at a disease free equilibrium, . Now consider the positively invariant subsystem of (1) without disease given byIt follows that the is given byi.e., the expression (3).

The disease-free equilibrium allows us to define a quantity useful in practical applications, vaccine coverage. The situation of COVID-19 is of course different, since there was no pre-existing vaccination coverage, but having access to this quantity is nonetheless helpful for simulations as it can be used to set parameters. Vaccine coverage is the fraction (typically expressed as a percentage) of the population vaccinated, i.e., . Vaccine coverage can be evaluated at any point in state space, for instance at endemic equilibria, although a closed-form solution is rarely obtained. When the system is at the disease-free equilibrium, though, the following expression holds, in which the dependence on is indicated:This expression also provides a useful form of the nonzero components of the DFE asTo continue the analysis, it is useful to define the combination parameterUsing , (2) takes the form

Theorem 1

The basic vaccination reproduction number for system (1) is given bywhere and are given by (3). If , then (3) of (1) is locally asymptotically stable, while it is unstable if .

Let us use the method of van den Driessche and Watmough (2002). Infected compartments are L, I and A (presented in that order) and the matrices F and V obtained by taking the Fréchet derivatives of the subsystem in the infected variables x written as and evaluating at take the formandComputing the dominant eigenvalue of gives (8). It is straightforward to show that conditions for applying van den Driessche and Watmough (2002, Theorem 2) are satisfied, giving local asymptotic stability of when and instability when .

Remark that , with equality when , i.e., when vaccine efficacy is zero. This implies that , the inequality being strict unless . Thus, vaccination always lowers the reproduction number. Actually, the following relationship links as given by (2) and for the full model (1):

Proposition 2

System (1) admits endemic equilibria according to the following cases:

A unique positive endemic equilibrium if or , and ;

Two positive endemic equilibria if , and ;

Zero positive equilibria otherwise, and in particular when .

Write (1) as , where . Solving the nonlinear system of equations , we find that, at a positive equilibrium, each of the components of can be written as functions of the equilibrium value of L, itself the root of the quadratic polynomialwhere Before considering endemic equilibria of the full system, let us focus on roots of (9). Since ,Since , it follows thatTherefore, (unless when , in which case ). Clearly,Note that is equivalent towhere is given in (2) and the right equality only holds when ; recall that when . Thus, when , , since . Since , this implies in turn that when , so by Descartes’ Rule of Signs, there are no positive equilibria when and . The case is discussed at the end of Section 2.2 and shown to preclude existence of positive equilibria; thus, there are no positive equilibria whenever .

Our next conclusion follows readily from the Rule of Signs of Descartes: regardless of the sign of , if , i.e., if , then there is a single positive root of (9). So now assume , i.e., . Using the Rule of Signs of Descartes again, if then (9) has no positive root.

So there remains the case where and to consider, where the Rule of Signs provides the unsatisfactory answer that there are either two (potentially with multiplicity 2) or zero positive roots to (9). Using the Rule of Signs on indicates that there are no negative roots of (9) when and . As a consequence, the situation is fully determined by the discriminant : if it is positive, there are two distinct positive roots of (9); if it is zero, there is a single (double) root; and if it is negative, there are no real roots.

Putting all this together gives the conditions in the result. Now return to the equilibrium for the whole system (1), expressed as a function of L as

The components I(L), A(L) and R(L) are positive scalar multiples of L. Furthermore, for all and when the equilibrium value of L is positive, implies that for all , implying in turn that for all . V(L) is a monotone decreasing function of L, so given a value of L, V(L) and S(L) take on unique positive values. Thus, if L is a positive root of (9), then it is the L component of a positive endemic equilibrium with other components given by (14).

When the root of (9) is unique, we denote it , whereas when there are two distinct positive roots, we denote and these values, with the convention that . The resulting endemic equilibria obtained by using (14) are denoted and . From the properties of (14), all components of are larger than those of , except .

Remark that the second condition in Proposition 2(i) defines the value of at which the saddle node bifurcation seen in the left part of Fig. 2 takes place; see the discussion in Section 2.3.

Fig. 2

(a) Subcritical () and (b) supercritical () transcritical bifurcations at . varies as a function of . N is the population size at the disease-free equilibrium. All other parameters can be found in Table 7 Appendix A

As a function of the transmission rate parameter , the basic vaccine reproduction number is monotone and there exists unique such that . Solving for yieldsThis is useful in simulations, as values of are easier to grasp than those of .

At this point, Proposition 2 tells us of situations where multiple endemic equilibria exist, so-called backward bifurcation situations, but says nothing of the stability properties of these equilibria. This is investigated (close to ) in Section 2.3.

Before doing so, though, it is worth briefly investigating particular cases, to check that the model behaves in a way that is similar to existing models exhibiting backward bifurcations. Of particular interest is the situation that prevails when the vaccine is fully efficacious (), or not efficacious at all (). These situations should preclude the existence of multiple endemic equilibria. Suppose . Then in (9), so P(L) has at most one positive root and there cannot exist multiple endemic equilibria. Suppose now that . Then , so that , implying there can be no endemic equilibria when .

Existence of a backward bifurcation

The existence of multiple non-trivial equilibria when given in Proposition 2 suggests that the transcritical bifurcation at the disease-free equilibrium can be subcritical. This phenomenon is more commonly referred to as a backward bifurcation Hadeler and van den Driessche (1997) and has important implications for control and ultimate eradication of the disease.

Recovered individuals develop a degree of protection from reinfection in the form of (possibly temporary) natural immunity as a result of their bodies’ own immune response Lumley et al. (2021). One goal of vaccination campaigns is to further increase the proportion of protected individuals in the population by generating a vaccine-induced immune response in individuals with no immunity or weak immunity. This provides vaccinated individuals with protection against infection but also can have a population-level effect, the so-called herd immunity. The latter occurs when a sufficient proportion of the population becomes immune, implying that infected individuals make most of their contacts with immune individuals, thereby reducing the risk of disease spread to the point that the remaining susceptible individuals are indirectly protected Desai and Majumder (2020). In an ideal setting, once herd immunity is achieved, the number of cases declines to extinction or near extinction. In reality, the possibility of reemergence of the disease remains an ever-present threat.

In mathematical models, it is often the case that the transcritical bifurcation at (or, in the case of models with vaccination, at is supercritical. In these models, the disease-free equilibrium can typically be shown to be globally asymptotically stable when the basic (vaccine) reproduction number is less than one. Therefore, vaccination campaigns are often designed with the reproduction number in mind, with the herd immunity threshold given by Aschwanden (2020)The effective reproduction number is the average number of new cases generated by a single infectious case in a population consisting of a combination of susceptible and non-susceptible individuals. The effective reproduction number is related to the basic reproduction number, but evaluated away from the disease-free equilibrium. The presence of a backward bifurcation implies that the effective reproduction number has to be reduced further, below a constant . In this section, we give an analytical proof of the existence of a backward bifurcation and illustrate it numerically for biologically reasonable parameters.

For convenience, write the state variable vector as . Then system (1) can be written . At each pair , we can calculate the Jacobian matrix . Writing and , we have and the Jacobian matrix evaluated at the disease-free equilibrium can be written . The following Lemma is the result of direct calculation.

Lemma 1

At , the Jacobian matrix admits a simple zero eigenvalue with left eigenvector and right eigenvector given bywhere is free.

Remark 1

Note that for . This follows from the fact that

For the purpose of the proof of the backward bifurcation, define the following hypothesis associated to the existence of two positive equilibria:where is given by (36), is the disease-free equilibrium and is given by (5).

We can now proceed with establishing the existence of the backward bifurcation analytically.

Theorem 2

System (1) undergoes a backward bifurcation at and whenever condition () holds.

Choose the bifurcation parameter . Then can be viewed as a function of with . By Lemma 1, at , admits a simple zero eigenvalue with associated eigenvectors and given in equations (17) and (18), respectively. Following Castillo-Chavez and Song (2004, Theorem 4.1), we calculate the bifurcation parametersLet the free component of , . From the Lemma, we have for . Thus, we can restrict ourselves to the mixed partials with respect to . Upon inspection we find that the only non-zero mixed partials areApplying this to (20), along with the fact that and due to Clairaut’s Theorem,Since andNow consider b. Among , only exhibits a nonzero derivative with respect to . The non-zero partials areHence,As noted in Remark 1, and when . Thus, if and only ifwhich is equivalent to hypothesis ().

Note that van den Driessche and Watmough (2002, Theorem 4) could also have been used to obtain this result.

Remark 2

By (3), the disease-free equilibrium number of susceptible and vaccinated individuals is a function of e, the vaccination rate of adults. It follows from (5) that the vaccination coverage is also a function of e. Since the disease-free equilibrium, and therefore the vaccination coverage, are independent of , the vaccination coverage is constant for each graph depicted in Fig. 2. The vaccine coverage, , is monotone increasing in e. However, the right hand side of hypothesis () is a non-monotone function of e.

Assuming that parameters other than e and are as in Table 7 in Appendix A, we find that the backward bifurcation exists (() is true) for . Therefore, the backward bifurcation exists for a wide range of vaccine coverage from approximately 35.4% to approximately 99.5%.

Table 7

Parameters for numerical simulations in Appendix 3.2

Parameter	Value
N	1×10^4
d	3.65×10^-5
B	3.65
γ	0.048
ω_r	0.032
μ	1.0e-4
p	0
σ	0.9
ω_v	5.5×10^-4

Remark 3

The existence of a backward bifurcation corresponds to a region of bistability for , i.e., existence of two positive endemic equilibria described in Proposition 2 (ii). As decreases, the nullclines corresponding to these equilibria come together at the vertex of the parabola (9). As seen in the graph on the left of Fig. 2, the vertex is also the point at which attractive and repulsive branches of equilibria collide, suggesting the point is a saddle-node bifurcation. We refer to the value of that corresponds to the vertex (or, equivalently, to this bifurcation point) as .

The stochastic SVLIARS model

We now consider a stochastic model in the form of a continuous time Markov chain which is related to the deterministic SVLIARS model (1). Informed by that model, we construct the CTMC, , with the infinitesimal transition probability to state j from state i,where is the transition rate associated with transition to j from i and can be found in Table 3.

Table 3

State transitions and rates for the CTMC SVLIARS model

Description	Transition i↦j	Rate ρ(i,j)
Birth of S	(S,V,L,I,A,R)↦(S+1,V,L,I,A,R)	B(1-p)
Waning of V	(S,V,L,I,A,R)↦(S+1,V-1,L,I,A,R)	ω_vV
Waning of R	(S,V,L,I,A,R)↦(S+1,V,L,I,A,R-1)	ω_rR
Vaccination of S	(S,V,L,I,A,R)↦(S-1,V+1,L,I,A,R)	eS
I infects S	(S,V,L,I,A,R)↦(S-1,V,L+1,I,A,R)	βSI
A infects S	(S,V,L,I,A,R)↦(S-1,V,L+1,I,A,R)	βηSA
Natural death S	(S,V,L,I,A,R)↦(S-1,V,L,I,A,R)	dS
Birth of V	(S,V,L,I,A,R)↦(S,V+1,L,I,A,R)	Bp
I infects V	(S,V,L,I,A,R)↦(S,V-1,L+1,I,A,R)	(1-σ)βVI
A infects V	(S,V,L,I,A,R)↦(S,V-1,L+1,I,A,R)	(1-σ)βηVA
Natural death V	(S,V,L,I,A,R)↦(S,V-1,L,I,A,R)	dV
Progress to I	(S,V,L,I,A,R)↦(S,V,L-1,I+1,A,R)	(1-π)εL
Progress to A	(S,V,L,I,A,R)↦(S,V,L-1,I,A+1,R)	πεL
Natural death L	(S,V,L,I,A,R)↦(S,V,L-1,I,A,R)	dL
Recovery of I	(S,V,L,I,A,R)↦(S,V,L,I-1,A,R+1)	γI
Disease death I	(S,V,L,I,A,R)↦(S,V,L,I-1,A,R)	μI
Natural death I	(S,V,L,I,A,R)↦(S,V,L,I-1,A,R)	dI
Recovery of A	(S,V,L,I,A,R)↦(S,V,L,I,A-1,R+1)	γA
Natural death A	(S,V,L,I,A,R)↦(S,V,L,I,A-1,R)	dA
Natural death R	(S,V,L,I,A,R)↦(S,V,L,I,A,R-1)	dR

The CTMC model characterized by (26) captures dynamics not captured by the deterministic model (1). In particular, model (26) is well suited to questions related to emergence and reemergence of a disease, when a small population of infectious individuals results in a high degree of randomness driving interactions. It also tracks integer numbers of individuals, making transitions into a state where the disease is extinct possible, rather than approached as a limit in ODE models.

The difference between these models can be seen in relation to the backward bifurcation, when it exists. In the case of model (1), the nullclines which are represented in Fig. 2 divide the solution space into basins of attraction which are forward invariant. However, in model (26), realizations can move between analogous regions with positive probability.

Consider the deterministic setting of model (1) and suppose that the model is at endemic equilibrium in the presence of a backward bifurcation and that is slightly greater than one. As a vaccination strategy is deployed, it reduces the reproduction number. However, the population will asymptotically return to endemic equilibria along the nullcline of Fig. 2 until decreases below , which is less than one.

In the remainder of this section, we explore two important natural questions related to the backward bifurcation and model (26). Does model (26) detect the backward bifurcation shown above for (1)? In other words, do numerical simulations of the probability of a minor epidemic in (26) change quantitatively and qualitatively in the parameter region for which the backward bifurcation exists for (1)? We also investigate whether the backward bifurcation can be detected in a branching process model commonly used to approximate the probability of a minor epidemic.

Multitype Branching Process approximation

As noted in the Introduction, it is common to formulate a branching process to approximate the probability of extinction in a Markov chain model. One of our primary objectives is to investigate whether CTMC models and the branching processes commonly used to approximate CTMCs detect a change in the dynamics associated to the presence of a backward bifurcation in the mean field system (1). For epidemic models with a single infectious individual type, a Bienaymé-Galton-Watson branching process (BGWbp) is used, and a multitype branching process (MTbp) is used for multiple infectious individual types.

Branching process approximation is a linearization technique. In the context of epidemic modeling, the branching process used to approximate a CTMC model captures birth and death probabilities of each infectious type at the disease-free equilibrium. As a result, branching process approximations are not appropriate to approximate behavior far away from the disease-free equilibrium and do not capture behavior driven by non-linear dynamics Milliken (2017).

By assuming the number of susceptible and vaccinated individuals are fixed at disease-free levels, we can construct the MTbp with the offspring probability generating function (p.g.f.)where is the offspring p.g.f. for , respectively. Then has the formwhere is the probability that an individual of type i gives birth to individuals of type j, for . Specifically, and

The probability of a minor epidemic relative to the CTMC with infinitesimal transition probabilities (26) can be approximated by the extinction probability of the multitype branching process with generating function (27). Suppose that individuals of types L, I, A, respectively, are introduced at the disease free quasi-stationary distribution. Suppose also that the multitype branching process satisfies the conditions of Harris (1963, Theorem 7.1). Then the probability of a minor epidemic iswhere is the fixed point ofwhere is the standard partial order in .

The threshold theorem of Allen and van den Driessche (2013) associates the criticality of the branching process to the reproduction number. Let be as above. Then the following theorem is a corollary to Harris (1963, Theorem 7.1) and the Threshold Theorem of Allen and van den Driessche (2013).

Theorem 3

The probability of extinction in the multitype branching process with generating function (27) is given by (30). If , then . If , then is the unique vector such that

First considerIt follows that, if with , thenThis together with the fact , implies that the multitype branching process is not singular. Furthermore, the matrix of first moments is primitive, since is positive Berman and Plemmons (1979). Therefore, the branching process is positive regular. The result follows from application of the Threshold Theorem of Allen and van den Driessche (2013) together with Harris (1963, Theorem 7.1).

The previous result illustrates an important implication of the Threshold Theorem of Allen and van den Driessche (2013) relative to branching process approximation of major epidemics in the presence of a backward bifurcation in the mean field system (1). The threshold theorem implies that the disease will go extinct almost surely in the approximating branching process. Since the mean field system is a kind of limit of the CTMC, the presence of a backward bifurcation there suggests we should be able to detect a backward bifurcation in the CTMC. However, the branching process approximation, according to the previous result, provides the same approximation of almost sure extinction throughout the region of bistability. In the next section, we use numerical simulations of sample paths of the CTMC to determine, by way of example, whether the backward bifurcation is detected in the CTMC or not.

Numerical results

In order to determine whether the CTMC detects the dynamical change suggested by the backward bifurcation in system (1), we estimate the probability of a minor epidemic as the frequency of minor epidemics in an ensemble of one thousand sample paths. The results are compared to the branching process approximation for the same choice of parameters. These experiments are performed for the parameters found in Table 7, except for and e, which are allowed to vary. The parameter e is allowed to vary to investigate both the super-critical (no backward bifurcation, ) and sub-critical (backward bifurcation, ) paradigms. In the absence of the backward bifurcation, varies to investigate the two cases and . In the presence of the backward bifurcation, varies to investigate the three cases , , and . The case represents the region of bistability in the presence of a backward bifurcation. There is no difference between this region and in the absence of the bifurcation. The case always corresponds to no positive equilibria and the case always corresponds to a unique positive equilibrium.

Numerical experiments are initialized at the disease-free quasi-stationary distribution following the introduction of a single individual of infectious type L, I or A.

Table 4 contains the results of branching process approximation and numerical simulation in the absence of a backward bifurcation. These results correspond to the graph on the right in Fig. 2. We see that the branching process approximation captures both the character and accurately quantifies the probability of a minor epidemic for all three sets of initial conditions and for both values of .

Table 4

Numerical results in the absence of a backward bifurcation. , varies with and the remaining parameter values are presented in Table 7. The branching process approximation of the probability of a minor epidemic appears in the column denoted BP Approx and the approximation of that probability based on the frequency of outcomes in one thousand sample paths appears in the column denoted Gillespie

(L_0,I_0,A_0)	No. of Pos. Eq.	R_v	BP Approx	Gillespie
(1,0,0)	0	0.99	1.000	1.000
(1,0,0)	1	1.10	0.990	0.990
(0,1,0)	0	0.99	1.000	1.000
(0,1,0)	1	1.10	0.990	0.991
(0,0,1)	0	0.99	1.000	1.000
(0,0,1)	1	1.10	0.990	0.991

Table 5 contains the results of the branching process approximation and numerical simulation in the presence of the backward bifurcation and correspond the graph on the left of Fig. 2. They tell a different story. As noted in the discussion following Theorem 3, the branching process provides the same prediction across the entire region , namely, that a major epidemic occurs with probability zero. Therefore, the branching process approximation does not detect any dynamical changes due to the presence of the backward bifurcation. Numerical simulation of the CTMC by the Gillespie algorithm does detect such changes. This indicates that a major epidemic occurs in one out of one thousand sample paths in the region where the deterministic model (1) exhibits bistability. We conclude that the CTMC has different dynamical behavior in this region that is not captured by the branching process approximation.

Table 5

Numerical results in the presence of a backward bifurcation. e=0.05, varies with and the remaining parameter values are presented in Table 7.The branching process approximation of the probability of a minor epidemic appears in the column denoted BP Approx and the approximation of that probability based on the frequency of outcomes in one thousand sample paths appears in the column denoted Gillespie

(L_0,I_0,A_0)	No. of Pos. Eq.	R_v	BP Approx	Gillespie
(1,0,0)	0	0.92	1.000	1.000
(1,0,0)	2	0.99	1.000	0.999
(1,0,0)	1	1.10	0.990	0.989
(0,1,0)	0	0.92	1.000	1.000
(0,1,0)	2	0.99	1.000	0.999
(0,1,0)	1	1.10	0.990	0.982
(0,0,1)	0	0.92	1.000	1.000
(0,0,1)	2	0.99	1.000	0.999
(0,0,1)	1	1.10	0.990	0.989

Results for the region provide further insight. In this region, the branching process approximation captures the qualitative behavior of the CTMC model – that major epidemics occur with positive probability . However, numerical simulations suggest that the branching process may not accurately capture the behavior quantitatively. The evidence for this can be seen in the last two rows of Table 5 in the rows corresponding to a unique positive equilibrium. Here we see that the probability approximated by simulation (column Gillespie) is consistently less than the branching process approximation (column BP Approx). This signal in the data can be considered weak since the sample mean has standard deviation on the order of the reciprocal of the square root of the number of simulations. In this case, only 1,000 simulations were performed. The weakness of the signal may also be due to the complexity of the system.

A similar analysis is performed for the SVIRS model studied in Appendix B. In that case, the signal that numerical simulation predicts a lower probability of minor epidemic that branching process approximation is stronger. For that model, the branching process is a single-type BGWbp and the approximation corresponds to Whittle’s approximation for an SIR model Whittle (1955),This relationship between and the probability of a minor epidemic is more representative of a model similar to Fig. 2(a) than to Fig. 2(b). That is, it is more similar to the case of the absence of a backward bifurcation.

Discussion

The phenomenon of bistability due to existence of a backward bifurcation has long been known to occur in SIR-type models that include vaccination, so it is not too surprising that an ODE SVLIARS model including imperfect waning vaccination, waning natural immunity, and disease-induced death, would exhibit the bistable behavior we proved here. The analysis of the SVLIARS model (1) is similar to that of similar SIR-type models, but serves as an important reminder that this dynamical behavior is present in many SIR-type epidemic models when they include demography. While vaccination always lowers the reproduction number, its imperfection and waning is also the source of the backward bifurcation in the deterministic model (1), with the important implication that the reproduction number must be lowered below the critical value to guarantee the elimination of the disease.

There are many situations in which a stochastic model provides a very valuable complement or alternative to a deterministic ordinary differential equations system, including when considering the emergence or re-emergence of a pathogen or variant. Indeed, they provide trajectories that are more realistic and have a nontrivial variance, and provide access to quantities such as mean first passage times that cannot be evaluated with ODEs. One class of stochastic models used in this context is continuous time Markov chains. Indeed, ODEs and CTMCs are very closely related and have many features in common. For instance, equilibria in a deterministic model are related to quasi-stationary distributions in CTMC models. However, there are also fundamental differences. For one, there is no well-defined separatrix or basins of attraction in CTMC models, so while the mean over a large number of realizations of a CTMC often tracks close to the solution of the corresponding ODE, in a case like the present one where the ODE exhibits bistability, the situation is more complicated.

The epidemiological consequences of a backward bifurcation to public health strategies and outcomes have been previously discussed in the literature; see, e.g., Dushoff et al. (1998); Greenhalgh and Griffiths (2009); Safan et al. (2006). Here we will just review the implications of this dynamical regime in two situations. In the case of an already occurring epidemic, the presence of a backward bifurcation means that the effective reproduction number must be reduced below a threshold rather than simply below the threshold at one in order to guarantee the elimination of the disease. Therefore, the presence of the backward bifurcation indicates that additional effort in the way of control measures are needed to eradicate the disease. In the case of emergence or re-emergence of the disease, if the basic reproduction number is below one but above , then the introduction of a few cases may not result in a major outbreak, but there is a critical number of cases, determined by the parameters, which would lead to a major outbreak. This unexpected behavior is even more likely due to the destabilizing effects of randomness.

For example, according to the deterministic model, the introduction of a few infected individuals generally results in a return to the disease-free equilibrium. However, our numerical experiments show that, in the case of the CTMC model, there is a positive probability of a major outbreak, consistent with approaching the quasi-stationary distribution associated to the stable positive equilibrium. This is true despite the fact that the reproduction number is less than unity. This phenomenon occurs in the presence of the backward bifurcation. Simulation of the CTMC in the absence of the backward bifurcation shows the disease almost surely goes extinct following the introduction of a small number of cases when the reproduction number is less that one. The behavior of the CTMC model depends on whether or not there is a backward bifurcation in the related ODE model, particularly in the region . We say the CTMC model detects the bifurcation.

However, as models become more complex, analysis and simulation of CTMC models becomes more difficult, often requiring the use of approximation methods. One such method, relatively simple and easy to implement, is referred to as branching process approximation. Herein lain the third broad question in this work: given a deterministic system exhibiting a bistable situation, would a CTMC and its branching process approximation detect the feature in the same way. The answer to this is negative: while the CTMC detects the backward bifurcation, the relevant branching process approximation does not. In both Sect. 3 and Appendix B, we have shown that, as proven by the combination of the Threshold Theorem of Allen and van den Driessche (2013) and the Criticality Theorem for branching processes Harris (1963, Theorem 7.1), the branching process indicates almost sure extinction of the disease when the reproduction number is below unity, regardless of the dynamical paradigm. In addition to the fact that the branching process approximation fails to detect the backward bifurcation when the reproduction number is less than one, it is also qualitatively inaccurate in the presence of the backward bifurcation when the reproduction number is greater than one. This is true for both the SVLIARS and the SVIRS models, but is more evident in the SVIRS example. This approximation error most likely stems from the fact that the branching process assumes a paradigm in which the probability of a minor epidemic is related to Whittle’s approximation (32). Near in this paradigm, the probability of a new infection (birth) and the probability of a recovery (death) are nearly equal. However, in the presence of a backward bifurcation (see Fig. 2(a)), the probability of a new infection is, in reality, much higher than the probability of a recovery. This results in greater directed motion away from the disease-free boundary and a higher probability of a major outbreak (lower probability of a minor outbreak). Thus, the major conclusion about bistability exhibited by both the ODE and the CTMC is obfuscated when the branching process approximation is used.

Table 6

Parameters for numerical simulations in Sect. 3.2

Parameter	Value
N	1×10^5
d	3.65×10^-5
B	0.365
η	1
π	0.4
ε	0.14
γ	0.071
ω_r	0.032
μ	1.0e-4
p	0
σ	0.9
ω_v	5.5×10^-4

Table 8

Description of model variables in SVIRS model

Variable	Description
S	Number of susceptible individuals
V	Number of individuals with vaccine induced partial immunity
I	Number of symptomatic infected individuals
R	Number of individuals with disease induced immunity

Table 9

State transitions and rates for the CTMC SVIRS model

Description	Transition i↦j	Rate ρ(i,j)
Birth of S	(S,V,I,R)↦(S+1,V,I,R)	B(1-p)
Waning of V	(S,V,I,R)↦(S+1,V-1,I,R)	ω_vV
Waning of R	(S,V,I,R)↦(S+1,V,I,R-1)	ω_rR
Vaccination of S	(S,V,I,R)↦(S-1,V+1,I,R)	eS
I Infects S	(S,V,I,R)↦(S-1,V,I+1,R)	βSI
Natural death S	(S,V,I,R)↦(S-1,V,I,R)	dS
Birth of V	(S,V,I,R)↦(S,V+1,I,R)	Bp
I Infects V	(S,V,I,R)↦(S,V-1,I+1,R)	(1-σ)βVI
Natural death V	(S,V,I,R)↦(S,V-1,I,R)	dV
Recovery of I	(S,V,I,R)↦(S,V,I-1,R+1)	γI
Disease death I	(S,V,I,R)↦(S,V,I-1,R)	μI
Natural death I	(S,V,I,R)↦(S,V,I-1,R)	dI
Natural death R	(S,V,I,R)↦(S,V,I,R-1)	dR

Table 10

Numerical results in the absence of a backward bifurcation. , varies with and the remaining parameter values are presented in Table 6

I_0	No. of Pos. Eq.	R_v	BP Approx	Gillespie
1	0	0.99	1.000	1.000
1	1	1.01	0.990	0.991

Table 11

Numerical results in the presence of a backward bifurcation. e=0.05, varies with and the remaining parameter values are presented in Table 6

I_0	No. of Pos. Eq.	R_v	BP Approx	Gillespie
1	0	0.92	1.000	1.000
1	2	0.99	1.000	0.992
1	1	1.01	0.990	0.977