Using asymptotics for efficient stability determination in epidemiological models

1.   Introduction

Local stability analysis is an important part of the investigation of any endemic disease model. When possible, it is ideal to identify a small set of inequalities that divide the parameter space into regions that define the set of equilibria that are locally asymptotically stable. Often the stability structure is simple, with a unique stable equilibrium at each point in the parameter space and relatively simple defining inequalities as well. A typical case is for there to be two regions, one with a globally stable disease-free equilibrium (DFE) and the other with a globally stable endemic disease equilibrium (EDE), with the boundary occurring where the basic reproduction number ${\cal{R}}_0$ is unity. More complicated disease models can produce more complicated behavior. For example, some models exhibit backward bifurcation, in which there is a region in the parameter space with both a stable EDE and a stable DFE ^[1,2,3], while other models have parameter space regions where there are no stable equilibria (see, for example, ^[4,5]).

Local asymptotic stability for a given set of parameter values can be determined by numerical evaluation of the eigenvalues, but general results require analytical methods. The direct method is to compute the characteristic polynomial as $P(\lambda) = \det (J-\lambda I)$, where $J$ is the Jacobian matrix and $I$ the identity matrix, followed by calculation of eigenvalues. This method requires complicated calculations; hence, it is seldom used for systems with more than two components. Systems of three components can be done without too much difficulty by finding the characteristic polynomial in the usual way and then employing the Routh-Hurwitz conditions ^[3,6,7]. In principle, the same method can be used for higher-order models, but this has rarely been done and no examples are given in the principal monographs of mathematical epidemiology; e.g., ^[1,3].

Our ability to determine stability requirements analytically depends on the efficiency of the methods we use. If the Jacobian is reducible, as is often the case for a DFE, then the eigenvalue problem decouples into a set of lower-dimensional eigenvalue problems. Otherwise, the Routh-Hurwitz conditions are more efficient than full calculation of eigenvalues. An example is in ^[3], page 101, where a 5th order eigenvalue problem is reduced to a 3rd order one, with the Routh-Hurwitz conditions used for the latter. There are also more efficient methods for calculating the characteristic polynomial (see Section 3). With these tools, it is frequently feasible to do analytical stability calculation for systems with four components, as in ^[8]. However, the additional tool of asymptotic approximation can extend the range of problems for which analytical stability analysis is possible to systems of five components, as shown below, and occasionally even six components, as in ^[9]. This gain comes at the minor cost of a degree of error that is insignificant when compared to error caused by uncertainty of parameter values or by neglecting minor features of disease natural history. In the development that follows, we provide detailed guidelines for the method and two illustrative examples. In so doing, we note the specific requirements for a model to be amenable to the method.

We begin in Section 2 with a general discussion of the characteristics of epidemiological models with two timescales. This is followed by a description in Section 3 of the three mathematical tools that form the basis for our method and a set of guidelines for employing the method. Section 4 consists of an example that goes through the process in some detail, culminating in a proof of stability subject to a simple, previously unknown, restriction on one of the parameter values. Section 5 consists of an example in which there is a very small region in the parameter space for which there are no stable equilibria. We conclude with a discussion of the benefits of the method.

2.   Two timescales in epidemiological models

Figure 1 shows a typical simulation for the infectious disease model,

where $I_0$ is very small. This is the basic SEIR epidemic model with demographic terms for death rates proportional to the population and a constant birth rate just large enough to maintain the population. The initial conditions are chosen so that the state variables represent population fractions. The scenario is a typical example of an infectious disease model with disease duration far less than host lifespan and a predominantly susceptible initial population. The plots help us to understand the importance of timescales in understanding epidemiological models.

The left panel shows the epidemic phase of the simulation. We see that

1) Infected population fractions ($E$ and $I$) are significant.

2) There is no indication of the endemic disease behavior that follows: we cannot tell whether there is a stable EDE.

In contrast, the plot on the right shows the endemic behavior. Here,

1) Infected population fractions are very small.

2) There is no indication of the epidemic behavior that preceded the endemic phase. Everything that happened in the first 100 days is located at the left edge of the plot.

3) Recurrences at approximately 48 and 72 years, and at shorter intervals thereafter, show rapid changes on the epidemic timescale.

The difference in model behavior on the epidemic and endemic timescales is a feature that can be utilized to simplify a given task. For investigations focused on the epidemic phase, anything that takes a long time to change the result, such as births and natural deaths, can be ignored. For investigations focused on the endemic phase, the low populations of the infected classes can be exploited through asymptotic approximation. This requires that the variables used to represent uninfected classes are scaled as population fractions, while the variables used to represent infected classes are rescaled under the assumption that they are vanishingly small in some asymptotic limit.

The parameters used for the simulation correspond to a mean lifespan of 78.3 years, a mean infectious period of 10 days, a mean incubation period of 5 days, and a basic reproduction number of 4, all typical values for infectious diseases. The ratio of disease duration to mean lifespan is 0.00035. This timescale ratio is responsible for the qualitative differences between the two panels in Figure 1. Even if the mean infectious period is as large as 30 days, this ratio is still only 0.001, which is more than small enough for the qualitative differences to be apparent.

3.   Mathematical tools

The method we employ for equilibrium analysis is based on a combination of three mathematical tools that we describe here: an efficient scheme for calculating the characteristic polynomial, asymptotic approximation to eliminate terms with negligible impact, and the Routh array for constructing asymptotic approximations to the Routh-Hurwitz conditions.

The Characteristic Polynomial Theorem  The usual way to compute the characteristic polynomial for a matrix $J$ is to find the determinant of $\lambda I-J$; ¹ however, there is a more efficient method based on a theorem from matrix theory (see ^[10], for example):

¹$\lambda I-J$ is preferable to the more usual $J-\lambda I$ because its leading order coefficient is always positive.

Theorem 1. For an $n \times n$ matrix $J$, let $S$ be the set of all nonempty subsets $K$ of the integers $1, 2, \ldots, n$. For each possible index set $K$, let $J_K$ be the determinant of the sub-matrix of $J$ that contains the entries in the rows and columns indicated by $K$. Then the characteristic polynomial of $J$ is

where

While this formula technically requires computation of all the sub-determinants of J, this is still easier than computing the determinant of $\lambda I-J$. The next tool allows us to neglect some of those sub-determinants.

Asymptotic Approximation  Asymptotic approximation is simply a matter of Taylor series expansion about $\epsilon \to 0$, or $\Gamma = 1/\epsilon \to \infty$ (see Chapter 1 of ^[11], for example). Each coefficient in the characteristic polynomial has an asymptotic expansion of the form

to determine the stability requirements to leading order, we need only to compute the $k_m$ and $p_m$. For systems that have dynamic equations on two different timescales, there are some rows in the Jacobian with a common factor of $\Gamma$ and some without. The presence of factors of $\Gamma$ in some rows of the Jacobian means that some of the sub-determinants in the calculations of $c_1$ through $c_{n-1}$ are dominated by others. For example, suppose row 1 in the Jacobian has a factor of $\Gamma$ while rows $2, 3, \ldots, n$ do not. Then at least one sub-determinant of the form $J_{1j}$ is $O(\Gamma)$, while all sub-determinants $J_{ij}$ with $i > 1$ are $O(1)$.

While one can retain terms beyond the leading order, if $\Gamma$ is sufficiently large then only terms that contribute to the leading order expansions of the characteristic polynomial coefficients need be retained to get a good approximation of stability requirements. This generally means that only dominant terms need to be kept in all calculations; however, all Jacobian terms of $O(1)$ should be kept long enough to calculate the characteristic polynomial coefficients, even if they are dominated in their differential equation by a term of $O(\Gamma)$. The necessity of this precaution is seen in Section 4.

The Routh Array  The Routh-Hurwitz (RH) conditions use the algebraic signs of combinations of the coefficients of the characteristic polynomial to determine if all of its roots are in the left half of the complex plane. The RH conditions for degrees 2 and 3 are well-known; for higher degrees it is convenient to construct them from the Routh array, which is a two-dimensional array with $n+1$ rows and $\lceil n/2 \rceil$ columns ^[12].² Let $R_{i, j}$ be the element in row $i$, column $j$ of the array. The array is populated from a characteristic polynomial

²Other computational schemes could also be used, such as the Hurwitz matrix; however, the Routh array is particularly convenient when combined with asymptotic approximation.

according to the following rules:

1) The first row consists of the leading 1 from the characteristic polynomial, followed by $c_2$, $c_4$, and so on.

2) The second row consists of the coefficients $c_1$, $c_3$, and so on, with a trailing 0 if $n$ is even.

3) For $i > 2$, the entries $R_{i, j}$ are given by

where formula terms requiring a column outside the dimensions of the array are taken as 0.

As an example, consider the Routh array for a fourth degree polynomial $P(\lambda) = \lambda^4 + c_1 \lambda^3 + c_2 \lambda^2 + c_3 \lambda + c_4$, which takes the simple form³

³It is common to leave 0 entries as blanks.

where

Expressing some of the entries as rational quantities, with a previously defined quantity in the denominator and a newly defined quantity in the numerator, as has been done here, is a convenient modification to facilitate the next step.

The Routh theorem picks the RH stability conditions out of the Routh array.

Theorem 2 (Routh). The critical point with characteristic polynomial $P(\lambda)$ is locally asymptotically stable if and only if all entries in column 1 of the Routh array are positive.

For the fourth degree case, the conditions are

Asymptotic approximation sometimes reduces the formulas that arise in the Routh array to simpler ones. For example, suppose $c_1 = k_1\Gamma$, $c_2 = k_2\Gamma^2$, and $c_3 = k_3\Gamma^2$, all as $\Gamma \to \infty$ with $k_j = O(1)$. Then the product $c_1c_2$ dominates $c_3$, and therefore $q_1 = k_1k_2\Gamma^3+O(\Gamma^2)$. Since $k_1 > 0$ is already required, we can replace $q_1 > 0$ with the simpler $k_2 > 0$. Simplifications such as this emerge naturally when the Routh array is constructed from the characteristic polynomial approximation

and only the leading order elements of the array are given.

3.1.   Guidelines for the asymptotic approach

The theorems for calculation of the characteristic polynomial and the RH conditions can be used on any problem. Obtaining asymptotic approximations to the stability conditions for an equilibrium requires that some guidelines be followed so that the method can be properly employed and so that algebraic complexity is kept to a minimum. We present a set of guidelines here; the starting point is the dimensional dynamical system model. See ^[6] for more detail on steps 1–4, which comprise the scaling.

1) Scale populations so that the DFE population total is 1. This can often be done during the initial model construction by choosing a birth rate coefficient to match the natural death rate.

2) Identify one principal fast timescale and one principal slow timescale.⁴ Use the principal slow timescale to scale the time.⁵

⁴For most infectious diseases in human populations, fast processes have timescales on the order of days or weeks, while slow processes have timescales on the order of years, or perhaps months. Usually the principal fast timescale is the amount of time spent in the principal infectious state and the principal slow timescale is the lifespan of the disease host.

⁵The fast timescale is preferable for simulations, but the slow timescale is preferred when the primary interest is in long-term behavior.

3) Define $\epsilon$ as the ratio of the principal slow rate to the principal fast rate. Make all other parameters dimensionless by reference to the appropriate principal scale; that is, rate parameters in fast processes, such as infection, should be referred to the principal fast rate, while rate parameters in slow processes, such as loss of immunity, should be referred to the principal slow rate.⁶

⁶Very few, if any, authors who write about scaling pay any attention to the definition of dimensionless parameters. This is unfortunate. The full set of parameters should be chosen to at least make sure there are not any parameters whose value is close to 1, such as $\gamma/(\gamma+\mu)$ and that only one parameter is "small" or "large" in the sense of being suitable for an asymptotic series.

4) Rescale "small" populations by a factor of $\epsilon$. Small populations are generally those for which the dynamics of increase and decrease are dominated by fast processes; for example, in the SEIR model, the dynamics of the latent ("exposed") class is due to the fast processes of infection and incubation, and the dynamics of the infectious class is due primarily to the fast processes of incubation and recovery.⁷ A population can be seen to be small if it appears only in terms with factors of $\epsilon^{-1}$ relative to the other terms in the equation. These terms look large, but only because the corresponding state variable has an inherent, albeit invisible, factor of $\epsilon$ when viewed on the long timescale. The rescaling serves to make all variables $O(1)$, so that the importance of terms is determined entirely by the visible factors of $\epsilon$. If no populations meet this criterion, then asymptotics has limited value, if any.

⁷The infectious class in the SI model is not small in this sense because there is no recovery process to drain away the infectious population.

5) It is often helpful to choose an equation order for the Jacobian that groups the fast equations at the beginning of the list. These are equations for which the derivative side has a factor of $\epsilon$ compared to the terms on the right side.

6) When equilibria cannot be determined by linear equations alone, use the equilibrium relations to try to write all but one of the dependent variables in terms of a principal dependent variable $z$. Derive the resulting polynomial equation $G(z) = 0$ for that dependent variable, but do not solve it. Parameter range requirements for solution existence can be determined in terms of the function $G$ evaluated at critical values when direct comparison of the critical value with the equilibrium values is difficult.

7) Define auxiliary quantities, as needed, to simplify formulas. Examples include quantities that represent weighted averages of infectious or susceptible populations or quantities that combine together in one entry of the Jacobian matrix. Because of the scaling, there are often terms that include a variable plus 1; in these cases, a more compact notation is facilitated by the convention of using a bar over the top of any quantity to indicate that quantity plus 1. A general rule of thumb is that having more entries in the table of symbols is preferable to having more terms in an equation or factors in a term.

8) Look for opportunities to simplify formulas in the Jacobian by replacing dependent variables with versions that include an extra factor. See (21) in Section 4 below.

9) During the stability analysis, it is important to keep checking whether combinations of auxiliary quantities can be simplified or whether it would be helpful to introduce additional combinations.

10) Unless the eigenvalues are immediately apparent from the Jacobian, use the appropriate RH conditions to determine stability. In doing so, use Theorem 1 to compute the coefficients of the characteristic polynomial to leading order. For anything over three variables, use the Routh array (Theorem 2) to determine the RH conditions to leading order.

4.   An SEIR model with two risk groups

As an example of the method outlined in Section 3.1, we consider a model previously studied by Gumel ^[2]. The analysis in that paper identifies a complicated inequality that is both necessary and sufficient for backward bifurcation to occur, but does not attempt to find a simpler necessary condition or to show that the inequality is satisfied in a meaningful portion of the parameter space. Here we show that 1) backward bifurcation requires a mortality probability of greater than 75%, an extreme condition that calls some of the model assumptions into question, and 2) when mortality is not that high, there is always a unique EDE that is locally asymptotically stable.

The model considers a disease scenario having two risk groups of changing composition, along with standard incidence and disease-induced mortality. The schematic appears in Figure 2. We euphemistically describe the subgroups of susceptibles as "protected" and "unprotected, " to be taken in the sense that "protected" people are less susceptible than "unprotected" people, for whatever reason. Individuals are born/recruited into either of these categories, with some fixed proportion of each. In the general case, they can also move between these categories through behavior changes. By choosing parameters so that all new individuals are initially unprotected and protected individuals never switch to unprotected, we obtain the special case of a disease with a vaccine that reduces the risk of infection for individual encounters without conferring immunity. The protected are less susceptible to the disease by a factor of $\sigma < 1$.

As needed during the analysis, we define new quantities to reduce the number of symbols appearing in formulas (see guidelines 7 and 8). For example, the total rate of infection can be written as $\beta QI/N$ where $Q = U+\sigma P = (1-\sigma) U+\sigma S$. Note that this definition implies $U \le Q \le S$. Quantities that appear in the original form of the model are listed in Table 1, while quantities that appear in the scaled model are listed in Table 2. Additional notation defined merely for algebraic convenience is introduced in the analysis and listed in Table 3 for convenience.

4.1.   Problem definition

From Figure 2, the model is

with

and with $T$ used for dimensional time in order to reserve $t$ for dimensionless time. The constant birth rate coefficient has been assigned the same symbol ($\mu$) as that for the natural death rate, thereby automatically scaling the model so that the disease-free population total is $N = 1$ whenever the initial conditions sum to 1.

The total susceptible population $S$ and the total population $N$ are more important quantities for understanding the model than $P$ and $R$; hence, we rewrite the model using $S$ and $N$ in place of these quantities:

where we have introduced the additional variable

to simplify the notation. We can think of $Q$ as a "susceptibility fraction" that indicates the size of an unprotected population that would have the same total susceptibility in the absence of a protected subclass as the actual overall population.

Guidelines 1–5 in Section 3.1 direct a recasting of the system to render it optimal for analysis. Our original system satisfies the first guideline. Next we identify the various rate parameters and use them to choose a timescale for nondimensionalization and an appropriate set of dimensionless parameters.

Natural demographic processes operate on a long timescale of years, giving us a natural long timescale $1/\mu$. Following guideline 2, we choose this as the timescale; hence, the dimensionless time variable is $t = \mu T$; this variable change is implemented as

The infection, incubation, recovery, and death processes operate on a short timescale of days or weeks, giving us fast rates $\beta$, $\eta$, $\gamma$, and $\delta$, respectively. The most important fast rate is $\gamma+\delta+\mu$, which is the rate of departure from the infectious class. Following guideline 3, we define the small parameter $\epsilon$ and dimensionless infection, incubation, and death parameters $b$, $\rho$, and $m$ as

Note that dimensionless versions of parameters often have a significant meaning of their own; here, $m$ is the probability that an infected individual dies and $b$ is the basic reproduction number for the special case where all susceptibles are high risk.⁸ These parameters are more significant to the model than the original parameters $\delta$ and $\beta$ that they represent. While $\rho$ does not have a special interpretation, we can see that it is the ratio of expected duration of infectiousness to expected incubation time for those individuals who do not die during the incubation period.

⁸The basic reproduction number if all susceptibles are low risk is $\sigma b$. The actual basic reproduction number $\sigma b < {\cal{R}}_0 < b$ must be determined in the analysis.

It is not immediately clear how to classify the status-change rate parameters $\Psi$ and $\Omega$, as it seems likely that the expected amount of time required to change risk class would be somewhere between the short time of infection duration and the long time of demographic change. With the benefit of hindsight (see below), we tentatively take them to be slow rates and justify this choice later; hence, we define dimensionless counterparts

where the combined rate constant $\Omega+\Psi$ is used rather than $\Omega$ because $\Omega$ occurs in the model (10) only as $\Omega+\Psi$. Note that $\Sigma$ can be interpreted as the expected number of times an individual changes their risk classification, while $\psi$ is the expected fraction of lifespan required for the transition from low risk to high risk.

With all the definitions given here, the scaled model is

where the prime symbol indicates the derivative with respect to dimensionless time $t$.

The asymptotic parameter $\epsilon$ appears on the left side of the differential equations for $E$ and $I$ as a timescale parameter, done in anticipation of these being identified as fast equations, that is, equations for variables that change rapidly in the epidemic phase, as seen in Figure 1.

Asymptotics is based on the idea that terms of $O(\epsilon)$ are less important than terms of $O(1)$ in the limit $\epsilon \to 0$. This only works if no quantities other than $\epsilon$ are arbitrarily small.⁹ If a problem is properly scaled, the reduced problem obtained from $\epsilon = 0$ should make sense on the long timescale. Here it does not. In the $N$ equation, for example, there is a term $-\epsilon^{-1}mI$ that appears to dominate the others if quantities other than $\epsilon$ are $O(1)$. It therefore appears that the leading order approximation to this equation is $I = 0$, which contradicts the expectation that the term $\epsilon^{-1}mI$ is of $O(\epsilon^{-1})$. On closer inspection, we see that every term with $I$ appears to be more important than it actually should be. Terms with $E$ are similarly mis-scaled.

⁹See ^[6].

The solution to this dilemma is to realize that the full population size is not the correct scale for the infected classes $E$ and $I$, as was apparent in Figure 1. Everything makes sense if we assume these quantities are $O(\epsilon)$ rather than $O(1)$. Thus, we rescale the infected class variables using the substitutions

With these new variables, we obtain the correctly scaled model,

where $\bar{\Sigma} = \Sigma+1$ and the system has been reordered to put the fast equations first, as per guideline 5.

At this stage, we can see why it was best to think of the transitions between low and high risk susceptibles, $\Omega U$ and $\Psi P$, as slow processes. Had we assumed they were fast processes and defined $\tilde{\psi} = \Psi/(\gamma+\delta+\mu)$ instead of the dimensionless parameter $\psi$ that we actually used, this would have given us the term $\epsilon^{-1} \tilde{\psi} S$ instead of $\psi S$. Since we cannot assume $S = O(\epsilon)$ without messing up the $S$ equation, we would have had to fix the problem by redefining $\psi = \epsilon^{-1}\tilde{\psi} = O(1)$, taking us to where we ended up with the original choice of scale for $\Psi$. This is not something we can expect to anticipate. Regardless of whether we make the right guess for the timescale of risk class change, we end up in the right place if we are careful to check that all terms have sizes given by the number of factors of $\epsilon$.

Note that the total population equation, after rescaling, is

To leading order, this is simply $R = N-S-O(\epsilon)$. Hence, $S \le N$ is required, but no required upper bounds for $X$ and $Y$ are apparent at this stage.

4.2.   The DFE and the basic reproduction number

For algebraic convenience (guideline 7), we define

where $w$ is any generic quantity.

The DFE has $X = Y = 0$, $S = N = 1$, $U = (\psi+f)/ \bar{\Sigma}$. The basic reproduction number, which would have been $bS$ if all susceptibles were high risk, is instead¹⁰

¹⁰See Appendix A.

It is convenient to define a parameter $c$ to be the numerator of the fraction representing ${\cal{R}}_0-1$; thus,

Here we omit the routine demonstration via eigenvalues of the Jacobian that the DFE is stable if $c < 0$ and unstable if $c > 0$. We should expect a unique EDE when $c > 0$, but this needs to be demonstrated. Similarly, there are usually no endemic disease equilibria when $c < 0$, but this property has to be checked.

4.3.   EDE existence and uniqueness

The equilibrium equations for $Y' = 0$, $X' = 0$, $N' = 0$, and $S' = 0$ yield

where the asymptotic operator $\sim$ indicates equality to leading order. We define new quantities

these quantities add new symbols to the lexicon of the analysis, but their use reduces the number of symbols that appear in formulas used in the analysis. In terms of these quantities, we have

to leading order. From the relations $q = \sigma s+(1-\sigma) u$ and $q = s+(1-\sigma) p$, we obtain additional useful formulas

To complete the specification of the EDE, ¹¹ we substitute the formulas for $1/N$, $s$, and $u$ into the equation $U' = 0$ to get the leading order equation,

¹¹See Appendix B.

or

where

Gumel found a complicated formula, involving all the parameters of the original model, for the condition required for a backward bifurcation ^[2]. The complexity of that condition makes interpretation difficult. Here we show that backward bifurcation requires unrealistically extreme parameter values.¹² Analysis of the leading order model is sufficient for this purpose, as the corresponding result for the original model would be merely an $O(\epsilon)$ correction to the leading order result.

¹²This is not a pejorative statement, nor is it a value judgment, but is based on the assumptions used to build the model. At the end of this subsection, we explain why the model we are using would be inappropriate with such a large disease-related death rate.

Proposition 1. In the asymptotic limit $\epsilon \to 0$, endemic disease equilibria for the model (17) with ${\cal{R}}_0 \le 1$ can occur only with the combination $m > \kappa$ and $m > 0.75$.

Proof. Guideline 6 is the key to establishing Proposition 1. We first note that solutions of the equation for endemic disease equilibria depend continuously on the model parameters, with the consequence that positive solutions with $c < 0$ (i.e., ${\cal{R}}_0 < 1$) occur if and only if the second solution for the case $c = 0$ (the other being $z = 0$) is positive. As the leading coefficient is positive, the non-zero solution is positive if the middle coefficient is negative; hence, we assume that the middle coefficient of $G(z)$ is negative while $c = 0$ and identify the implications of this assumption. Multiplying through by $b$ and using $\kappa = \sigma b$, this means that we require

while

Multiplying the inequality by $1-\kappa$ and substituting from the equality yields the inequality

which we can rearrange to get

We see that this inequality requires $m > \kappa$. With the restrictions $\psi \ge 0$ and $f \le 1$ that are part of the model assumptions, we obtain

as a necessary condition for EDE existence. Rearranging yields

The left-hand side of this inequality is already known to be positive, so the right-hand side must also be positive. Thus, we must have

To complete the description of existence and uniqueness properties, we consider the possibilities for solutions of (25) when ${\cal{R}}_0 > 1$.

Proposition 2. In the asymptotic limit $\epsilon \to 0$, there is a unique EDE for the model (17) that satisfies all restrictions on the signs of the state variables whenever ${\cal{R}}_0 > 1$.

Proof. It is immediately clear that the equation (25) has a unique positive solution when $c > 0$, but it still needs to be shown that the solution is always biologically meaningful. The critical requirement is $p \ge 0$, which through (24) and (23) corresponds to

where $z^*$ is the positive root of the function $G$ defined in (26). Given that $G$ is a quadratic function that is negative at $z = 0$, $0$ at $z = z^*$, and becoming infinite as $z \to \infty$, the requirement $z^* \le \hat{z}$ is equivalent to

Replacing products $(1-m)\hat{z}$ by $b-1$, we have

which simplifies to

Given the parameter restriction $f \le 1$, along with (27), we have

as required. □

The results of Propositions 1 and 2 are very important from a modeling standpoint. While backward bifurcation and nonuniqueness of equilibria are important features of the mathematical problem (8), they are not important features of the epidemiological setting from which the mathematical problem was derived. These interesting mathematical results are only possible if $m > 0.75$. This means that more than 75% of people infected with the disease die of it; of known epidemic diseases, only the pneumonic plague might possibly have had a mortality fraction this high.¹³ Worse yet, the model assumptions of recruitment and encounter rates that do not depend on the total population are inappropriate for a disease with such high mortality. In reality, a population drastically reduced by disease mortality can be expected to have a much lower birth rate and a much lower average number of encounters per day than one not so reduced. Different assumptions about these processes are required for such a disease, leading to a different model than the one considered here.

¹³According to the World Health Organization, plague killed 30–60% of the population in affected cities in Europe. People who were not infected are included in this statistic, so the actual mortality fraction among the infected would have been a little bit higher, even possibly higher than 75% ^[13].

4.4.   EDE stability

Given that we now know that there is a unique EDE if and only if ${\cal{R}}_0 > 1$, it seems likely that this EDE is always stable. To confirm this conjecture, we begin by identifying a set of three stability requirements, given as inequalities in terms of the EDE solution $y$ and the other parameters. These follow from the use of the three tools described in Section 3 in accordance with guideline 10.

Proposition 3. An EDE for the model (17) is locally asymptotically stable in the asymptotic limit $\epsilon \to 0$ if and only if $A, B, C > 0$, where

Proof. The first step in the determination of these stability criteria is to obtain the characteristic polynomial, which is (see Section C for a detailed derivation)

to leading order, where

Neglecting all lower-order terms, the corresponding Routh array takes the simplified form

where

The quantity $k_1$ is already known to be positive. Given the requirements $q_1 > 0$ and $k_4 > 0$, the requirement $q_2 > 0$ shows that $k_3 > 0$ is also required. Since the formula for $k_3$ is simpler than that for $q_1$, we need only to check $k_3 > 0$ and $q_2 > 0$ to confirm $q_1 > 0$, and of course we have to check $k_5 > 0$. The requirement $k_4 > 0$ is automatically satisfied when $k_3 > 0$ and $k_5 > 0$. Using the simplifying notation introduced earlier, we are left with three conditions: $A > 0$, $B > 0$, and $C > 0$, where

Substituting in the formulas for the $k$'s yields the formula given in (29). □

The three stability criteria given in Proposition 3 can be used to determine stability of the equilibrium solution for any specific set of parameters. However, with careful algebra we can establish a strong result (see Section D for the proof):

Proposition 4. If $m \le \kappa$ or $m \le 0.75$, any EDE for the model (17) is locally asymptotically stable in the asymptotic limit $\epsilon \to 0$.

The combination of this result with Propositions 1 and 2 establishes that the result of the model analysis is straightforward in the realistic case $m < 0.75$, with a unique stable EDE if ${\cal{R}}_0 > 1$ and a stable DFE if ${\cal{R}}_0 < 1$.

5.   Variant competition

We now consider a model of variant competition, inspired by, although not specifically designed for, the possibility that the first omicron variant of COVID-19 might have been less infectious to fully susceptible individuals than the delta variant that it replaced. The model is a simplified version of one designed to illustrate cross-immunity ^[14].

For this simplest model of variant competition, suppose there is a resident variant $I$ and an invading variant $J$. Individuals who recover from the resident variant move to the partially-recovered class $P$, where they remain susceptible to the invading variant. Individuals who recover from the invading variant move to the fully-recovered class $F$, regardless of whether they have previously been infected by the resident variant. Relative to the infectiousness of the resident variant toward the fully susceptible, the infectiousness of the invading variant is a factor $0 < \kappa < 1$ toward susceptibles and a factor $0 < \delta \le \kappa$ toward partially-recovered individuals. Fully-recovered individuals can lose immunity to the invading variant, moving them into the partially-recovered class, but nobody can lose immunity to the resident variant. Finally, we assume that there is no mortality due to either variant¹⁴ and that the mean recovery time from both variants is the same.

¹⁴The model is inspired by COVID-19, but it is obviously not appropriate for this disease.

More features could be added to the model if we had a particular disease in mind and detailed data on a particular pair of disease variants or if we wanted to focus on the effects of differences in mortality or recovery rate. In this investigation, we want to see if the limited features we are including are sufficient to yield an outcome that favors the variant with the lower basic reproduction number. The invading variant has the advantage of being infectious toward some recovered individuals to compensate for the disadvantage of being less infectious toward fully susceptible individuals. With these assumptions, it seems plausible that different portions of the parameter space could yield a variety of outcomes, including stable equilibria that include only one of the two variants, a stable coexistence equilibrium, and possibly even an unstable coexistence equilibrium with a stable limit cycle.

The model is illustrated in the schematic diagram of Figure 3.

5.1.   Problem definition

From Figure 3, the model is

Note that the constant birth rate has been assigned the same symbol as the natural death rate. This makes the population constant; without loss of generality, we can think of the populations having been scaled so that the total population is $N = 1$.

There are four dependent variables after using the algebraic equation to eliminate $F$ from the system. A full stability analysis for systems of this size can be done using the RH criteria without asymptotics, but the algebra can be obscure and close to unmanageable. Instead, we can achieve approximate stability boundaries in a relatively simple form by making use of asymptotic approximation. In preparation for this, we need to scale the time and rescale the two infected class populations.

We consider the loss of immunity and natural demographic processes to operate on a long timescale of years, or at least months, while the infection and recovery processes operate on a short timescale of weeks. Thus, we have slow rates $\mu$ and $\theta$ and fast rates $\gamma$ and $\beta$, along with rates that combine multiple processes, such as the fast rate $\gamma+\mu$ that represents the rate by which individuals leave either infectious class. We choose $\mu$ and $\gamma+\mu$ as the best rates to represent slow and fast dynamics and choose $1/\mu$ as the timescale for the dimensionless version of the problem. In the nondimensionalization, we refer other rates to the appropriate principal rate according to whether the rate in question is fast or slow.

Our choices dictate parameter definitions and variable change

As in the Two Risk-Group model, the parameter $\epsilon$ is the ratio of the principal slow rate $\mu$, which represents mean lifespan, to the principal fast rate $\gamma+\mu$, which represents mean time in the infectious class. In addition, $b$ proves to be the basic reproduction number for the resident variant and $h$ is the ratio of lifespan to mean time for loss of immunity. We expect $\epsilon$ to be less than 0.001 (which, with a life expectancy of 80 years, would require an infections period lasting a whole month) and $b > 1$, while $h$ could vary greatly from one disease to another. To simplify algebraic notation, we adopt the convention of putting a bar over the top to indicate adding 1 to any other quantity; for example,

At this stage, the scaled model is

As with the Two Risk-Group model, stability analysis of the Variant Competition model requires a rescaling of the infected populations, which we do using the substitutions

With these new variables, we obtain the correctly scaled model,

where

have been introduced because it is algebraically better to have more symbols in a model than more terms in a formula.

As before, $\epsilon$ appears on the left side of the first two differential equations as a timescale parameter, where it signifies a fast variable, and on the right side of the last equations as a regular parameter, where it represents a small perturbation. In the analysis that follows, we define $\Gamma = 1/\epsilon$ and do asymptotic approximation in the limit $\Gamma \to \infty$. Removing terms of $O(\epsilon)$ from (44) yields the leading order problem

5.2.   Equilibria

The model (46) has one DFE and as many as three endemic disease equilibria, which we call EDE-Y, EDE-Z, and EDE-YZ, with the letters after the dash indicating which of the infectious populations are non-zero. In calculating these equilibria, we must be careful to check that $S$, $P$, $F$ and $Y$ are all non-negative, with the non-negativity of $Z$ following from that for $F$.

It is convenient to define quantities

as the equilibrium equations for the endemic disease equilibria are a little simpler with these modified variables.

5.2.1.   DFE and EDE-Y

The DFE is easily determined by setting $Y = Z = 0$, thereby obtaining also $F = P = 0$ and $S = 1$.

With $Z = 0$ and $Y \neq 0$, we obtain the solution

for EDE-Y. The only requirement for existence is $b > 1$. Thus, the basic reproduction number for the resident variant is $R_{0r} = b$.

5.2.2.   EDE-Z

With $Y = 0$ and $Z \neq 0$, we quickly obtain results

The equations $\kappa s+\delta p = 1$ and $s+p+f = b$ then allow us to write $p$ and $f$ in terms of $s$:

We then have

Substituting this result into the remaining equation, $sx+s = b$, and multiplying by $\delta$ to clear any fractions, yields a quadratic equation for $s$:

This equation has a unique positive solution $s^*$, but we still need to check that $p$ and $f$ are non-negative.

The requirement $p \ge 0$ is equivalent to $\kappa s^* \le 1$. Now note that $G < 0$ for $s < s^*$ and $G > 0$ for $s > s^*$; hence, the requirement $\kappa s^* \le 1$ is equivalent to $G(1/\kappa) \ge 0$. Now

hence, $p \ge 0$ iff $\kappa b \ge 1$.

Similarly, $f \ge 0$ is equivalent to $(\kappa-\delta) s^* \ge 1-\delta b$, or

A similar calculation as before shows that this inequality is always satisfied when $\kappa b \ge 1$. Hence, the existence of the EDE-Z equilibrium requires $\kappa b \ge 1$, which confirms that the basic reproduction number for the invading variant is $R_{0i} = \kappa b$.

5.2.3.   EDE-YZ

If neither $y$ nor $z$ is 0, then we must have $s = 1$ and $q = 1$. The equations work out to the solution

For this equilibrium to exist, we must have both $f > 0$ and $y > 0$. These requirements lead to bounds on $x$: $x > p$ and $x > \kappa z$. The first of these comes down to $\delta(b-1) > 1-\kappa$, while the second is $(\bar{h} \kappa-1)\delta x < \bar{h} \kappa \delta p$. This could be satisfied either by $\bar{h} \kappa < 1$ or by

In the latter case, we have two upper bounds on $\delta$, the other being $\delta \le \kappa$. By identifying the point where the two upper bounds match, we obtain the following:

1) If $\frac{1}{b} \le \kappa \le \frac{h+b}{ \bar{h} b}$, then existence of EDE-YZ requires $\frac{1-\kappa}{b-1} < \delta \le \kappa$.

2) If $\frac{h+b}{ \bar{h} b} \le \kappa \le 1$, then existence of EDE-YZ requires $\frac{1-\kappa}{b-1} < \delta < \frac{ \bar{h} \kappa (1-\kappa)}{(\bar{h} \kappa-1)(b-1)}$.

Figure 4 shows the existence regions for the endemic equilibria with $b = 3$, $h = 0$ (left), $b = 3$, $h = 2$, (center), and $b = 6$, $h = 4$ (right) in the $\kappa \delta$ plane. Larger values of $b$ move the boundary between regions 1 and 2 to the left, while larger values of $h$ move the boundary between regions 3 and 4 to the left. Thus, shorter times for immunity loss tend to suppress the coexistence equilibrium in favor of the invader-only equilibrium, although it is not yet clear what that does to the stability of the various endemic disease equilibria. Note that the region $\delta > \kappa$ is not appropriate for the biological problem, as that would mean that prior infection by the resident variant increases susceptibility to the invading variant.

5.3.   Asymptotic stability analysis

For stability, we can use the equation

to eliminate $F$ from (46), leaving state variables $Y$, $Z$, $S$, and $P$. The Jacobian for this system is

where

5.3.1.   DFE stability

The Jacobian for the DFE is lower triangular, with diagonal entries $(-1+b)\Gamma$, $(-1+\kappa b)\Gamma$, $-1$, and $- \bar{h}$. Given the problem requirement $\kappa < 1$, the binding stability requirement is $b < 1$; of course $b$ is the basic reproduction number for variant 1 and $\kappa b$ is that for variant 2 when variant 1 is absent.

5.3.2.   EDE-Y stability

For the EDE-Y equilibrium, the Jacobian is

The eigenvalue problem for this system decouples into three smaller problems, a degree 1 problem with eigenvalue $(-1+q)\Gamma$, a degree 1 problem with eigenvalue $- \bar{h}$ and a degree 2 problem using the sub-Jacobian that omits rows 2 and 4 and columns 2 and 4,

This matrix has negative trace and positive determinant, so the corresponding eigenvalues are negative. The sole stability criterion is $q < 1$, or

This equilibrium is therefore stable in regions 1 and 2, but not in regions 3 and 4.

Note that the instability of the EDE-Y equilibrium when (55) is not satisfied means that the inequality

is the criterion for successful invasion of the invading variant into a population with the resident variant. One might have expected a greater rate of immunity loss ($h$) to benefit the invader; it may still lead to a larger invader population, but it does not enhance the ability of the invading variant to invade.

5.3.3.   EDE-Z stability

EDE-Z has Jacobian

where

and $s$ comes from the quadratic equation $G(s) = 0$ (49). As with EDE-Y, the system is decoupled, this time with eigenvalue $(-1+s)\Gamma$ and three eigenvalues coming from

Given the characteristic polynomial formula

we have

The requirements $c_1 > 0$, $c_2 > 0$, and $c_1c_2 > c_3$ are immediate, and $c_3 > 0$ follows from $\kappa \xi > \delta h$ This leaves the criterion $s < 1$, which is equivalent to $G(1) > 0$, where

This comes down to

Thus, EDE-Z is stable in region 4 and unstable in regions 2 and 3. Note that region 4 disappears if $h = 0$; hence, with no loss of immunity to the invading variant, the resident strain is always persistent.

5.3.4.   EDE-YZ stability

The Jacobian for EDE-YZ is

where $p$, $z$, and $y$ are given in (50). Instead of computing the characteristic polynomial for $J$ as the determinant of $\lambda I-J$, we use the characteristic polynomial theorem ^[6] along with asymptotic approximation to find just the leading order characteristic polynomial coefficients.

There are only a small number of non-zero terms in the matrix, and some combine together in the same way for more than one coefficient. One such example is the product of the 4, 2 and 2, 4 elements, which appear in $J_{24}$, $J_{234}$, and $\det(J)$. It is helpful to look for such combinations and give them separate symbols; here we define

The coefficient $c_1$ is the negative of the trace:

For $c_2$, we have four non-zero $2 \times 2$ sub-determinants, yielding

Continuing, we obtain

With leading order terms for each coefficient, the characteristic polynomial is

where

The RH conditions are

where

The first two are clearly satisfied, and

This leaves only the condition $q_2 > 0$, or

We now know that there is a unique stable equilibrium (EDE-Y) in regions 1 and 2 of Figure 4, a unique stable equilibrium (EDE-Z) in region 4, and a possible stable equilibrium in region 3 (EDE-YZ). For the EDE-YZ equilibrium, there remains one inequality that needs to be checked. The form of the inequality as written is simple; in practice, we have to use numerical methods since the criterion is nonlinear. The parameter spaces obtained by asymptotic approximation for $b = 5$, $h = 0$ and $b = 5$, $h = 1$ are shown in the left panels of Figure 5. Each has a small region of instability (red with black border), bounded in part by the line $\delta = \kappa$ and lying toward the midpoint of the region where EDE-YZ exists. Sample points in the unstable region that are used for illustrative simulations are marked by black dots. In the $h = 1$ case, the region of instability is so small that it is entirely covered by the black dot.

The center and right panels of Figure 5 examine the solution at the points of instability marked on the region plots. The initial conditions were chosen to hit the point in the limit cycle where the resident population ($Y$) is largest. The fully susceptible population is low at that point, relative to its average, while the partly recovered population is high. These levels cause a decrease in the resident population and an increase in the invader population. At roughly 7 years, the resident population reaches its low point while the invader population is still decreasing; however, the reduced competition from the decreasing invader is enough to reverse the decline of the resident. The resident population continues to increase while the invader population reaches its low point and begins to recover. This instability only happens for a narrow range of $\kappa$ and $\delta$ very close to $\kappa$. Otherwise, the EDE-YZ solution is stable in its region of existence.

Note that loss of immunity for fully-recovered individuals ($h > 0$) has no effect on persistence of the invader, serving instead to increase the equilibrium size of the group infected by the invader and help the invader drive the resident to extinction. As $h$ increases beyond 1 (which corresponds to a very slow loss of immunity), the large region where both variants persist quickly diminishes in favor of a large region where the invader alone persists. By $h = 10$, corresponding to loss of immunity with mean time of 8 years, the invader becomes dominant in most of the region where it is persistent.

To get a picture of just how changes in $b$ and $h$ affect the region of instability, we plot instability regions in the $b\kappa$ plane with various values of $h$ and with $\delta = \kappa$; these appear in Figure 6. The top group is two pairs of curves, a blue pair that shows the instability region for $h = 0$ with asymptotics and a black pair that shows the slightly smaller instability region for $h = 0$ with $\epsilon = 0.001$. As $h$ increases, the region of instability quickly shrinks.

6.   Discussion

We have now seen two examples of local stability analysis facilitated by asymptotic approximation in conjunction with the simplest formula for calculating the characteristic polynomial of a matrix and a mechanism for constructing RH conditions for polynomials of any degree. Using the characteristic polynomial coefficient theorem and the Routh array come at no cost in accuracy, while the use of asymptotic approximation comes with the cost of using an approximate result in place of an exact result. We now discuss in more detail the advantages and disadvantages of approximation.

6.1.   Advantages and disadvantages for qualitative results

The main advantage of asymptotic approximation is that it can simplify complex formulas to the point where intractable ones become tractable, as in the stability demonstration for the Two Risk-Group model. While this demonstration required some messy algebra, it is instructive to compare it to the alternative computation performed without benefit of asymptotics. Using the generic form for the characteristic polynomial that includes the $c$ coefficients rather than products of the $k$'s with factors of $\Gamma$, the entry in row 5, column 1 of the Routh array is given by $q_4 > 0$, where $q_4$ is calculated from the successive formulas

If these are combined into a single formula with all parentheses removed, the result is

This is sufficiently complicated to deter anyone from attempting analytical stability analysis for systems with five components.

Now suppose we apply asymptotics from this starting point. When we replace the $c$'s with factors of $k$'s and $\Gamma$, we discover that the first, fourth, and fifth terms in $q_4$ are $O \left(\Gamma^6 \right)$, while the other four are $O \left(\Gamma^5 \right)$ or smaller. Removing the lower order terms and simplifying reduces the formula to $q_4 = c_1c_4q_2$. Other RH conditions have already required $c_1$ and $q_2$ to be positive, so the condition reduces further to the single condition $c_4 > 0$, corresponding to the condition $k_4 > 0$ that appeared in the same position of the Routh array with asymptotics. In our development, we showed that this condition is never binding because $k_4 = k_3+k_5$, a result that would have been impossible to obtain from (69).

It would have been far more difficult to obtain results in the general case similar to what we obtained using asymptotics. There remains the question of whether results that require an asymptotic limit have general value. In problems such as the Two Risk-Group model, where stability is guaranteed for any reasonable parameter values ($m \le 0.75$), asymptotic stability as $\epsilon \to 0$ must at least hold for $\epsilon$ values below some finite threshold. In epidemiological models where the slow and fast timescales are as different as the duration of a typical human disease and the typical human life span, actual values of $\epsilon$ are quite small, so there is a very high likelihood that the results hold for any realistic values of $\epsilon$.

There is also the question of whether retaining generality in parameter values is sufficiently important to justify the complicated calculations required for a general analysis. In the Variant Competition model, the analysis found a very small portion of the parameter space where there are no stable EDEs. Had we instead looked for such a region by viewing example simulations, not knowing what parameter ranges to examine, we would have had to view quite a lot of examples in order to be lucky enough to find one with instability.

6.2.   Advantages and disadvantages for quantitative results

Even in the case of quantitative results, such as in the Variant Competition model of Section 5, the error caused by the approximation $\epsilon \to 0$ is insignificant compared to the natural error in results caused by inexact knowledge of parameter values and neglect of minor epidemiological features that were left out of the model. This is seen in Figure 6, where the black region corresponding to $\epsilon = 0.001$ is almost visually indistinguishable from the corresponding blue region obtained by asymptotics. The value $\epsilon = 0.001$ is not unreasonably small; with an 80-year human lifespan, this corresponds to an infectious period of duration 0.08 years, corresponding to almost a month. Most diseases have a shorter infectious period, which would yield a smaller value of $\epsilon$.

It is also worth noting that the most convenient method for obtaining the boundary for $\epsilon > 0$ relies on prior knowledge of the asymptotic stability result. Points on the boundary curves in Figure 6 are found by converting the defining inequality into a function $q = f(b, h, \kappa, \delta)$ whose value is positive for stability and negative for instability and then numerically solving the equation for $\kappa$ using the fixed values of $b$ and $h$ and $\delta = \kappa$. The asymptotic work identified which of the RH conditions fails at instability, allowing the function $f$ to be obtained from a single entry of the Routh array, calculated numerically. Absent knowledge of which specific condition fails, the RH conditions could not have been used to create the function $f$ and stability would instead have had to be determined by numerical computation of the dominant eigenvalue.

6.3.   The method itself

Two minor observations on the guidelines are warranted. First, one could avoid the rescaling by assuming a regular perturbation structure and accepting the non-standard possibility of a leading order term of 0, that is, using $I = I_0+\epsilon I_1$ rather than $I = \epsilon Y$. This only works if we begin the analysis by determining the result $I_0 = 0$ and then retaining the "correction" term $I_1$ while not using a regular perturbation expansion for $S$, $U$, and $N$. The result is identical to what we have done here, except with the symbol $I_1$ in place of the symbol $Y$. The two choices are thus mathematically equivalent, but the rescaling choice has the advantage of being conceptually more clear. Second, the introduction of so much additional notation adds a degree of inconvenience, as one must keep a list of symbol definitions handy. However, the advantage of guidelines 7 and 8 becomes clear if the reader replaces the extra symbols in the formula for the Jacobian with their definitions and then tries to obtain the characteristic polynomial coefficients. It is far better to have more symbols in the table of notation than to have more symbols in the formulas needed for analysis.

6.4.   Broader application

There is no reason why the method presented here should be restricted to epidemiological models. Any model that has processes with timescales that are significantly different could be approached with the same method. In a food web model, for example, there could be organisms with a very short life span, such as microorganisms, along with organisms with much longer life spans. It is common for predator species to have a longer lifespan than their prey, but we should be wary of using the asymptotic approach if the ratio of the lifespans is on the order of 10 rather than 1000, as we generally have in epidemiology models. But even in these less clear cases, the asymptotic approach can give useful approximations, as long as they are not taken to be quantitatively accurate.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

Glenn Ledder is a guest editor for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. The author declares there is no conflict of interest.

Appendix

A.   Basic reproduction number for the two risk-group model

The most straightforward way to calculate the basic reproduction number is to use its definition as the expected number of secondary infections per primary infective in an otherwise fully-susceptible population. In the model (17), new infections are created at rate $\epsilon^{-1}bQY/N$, or $\epsilon^{-1}bQ/N$ per infectious individual. The expected duration of time in the infectious compartment is $\epsilon^{-1}$; hence, the expected number of secondary infections per infective is $bQ/N$, leading to the result

where the state variables $Q$ and $N$ are evaluated at the DFE. We have $N = 1$ at the DFE, along with

Multiplying by $b$ and using (18) yields (19).

B.   Derivation of the EDE equation (26) for the two risk-group model

Beginning with the $U' = 0$ equation, we divide by $N$ and substitute from (21) and (22) to get

Multiplying by $h = (1-\sigma)b$ then yields

From (23), we have

Substituting these formulas and (24) into (71) yields

after which gathering like terms yields (26).

C.   Proposition 3

The Jacobian, given in terms of the additional symbol definitions of (18), (21), and (22), takes the compact form

where $\Gamma = 1/\epsilon$. The leading order characteristic polynomial is obtained from this matrix by a combination of the characteristic polynomial theorem (Theorem 1) and asymptotic approximation.

The coefficient $c_1$ is the negative of the trace. Some terms are $O(\Gamma)$, while others are $O(1)$. To leading order, we need only to consider the former; hence,

At first glance, $c_{12}$ appears to be dominated by an $O(\Gamma^2)$ term coming from $J_{12}$; however, the terms of that order cancel, giving $J_{12} = \Gamma$. This means that the other $O(\Gamma)$ terms, coming from sub-determinants with indices 1 or 2 paired with 3, 4, or 5, also contribute to the leading order approximation. The result is

Note that we need the second term in the 11 entry even though it is subdominant to the first term. That second term is not needed for any other characteristic polynomial coefficients, but it was necessary to retain it in the Jacobian.

The coefficient $c_3$ is dominated by terms of $O(\Gamma^2)$ that come from the three sub-determinants of form $J_{12j}$. These conform to a simple pattern:

where

Thus,

where $A$ was defined in (28).

The full determinant calculation yields

with $B$ given in (28).

The coefficient $c_4$ is the sum of the $4 \times 4$ sub-determinants of $J$. The three (of the total of five) that include both rows 1 and 2 contribute terms of $O(\Gamma^2)$. Each of the three required sub-determinants is described by the template

Cofactor expansion on the first column yields

Careful accounting of the terms in this formula from each of the three relevant sub-matrices yields the convenient result

To summarize, we have found the leading order characteristic polynomial

where

Note that asymptotic simplification only applies to terms within each individual coefficient. We cannot compare the magnitudes of the different terms because there are eigenvalues of different orders. As an example, suppose we assume $\lambda = O(1)$. The leading order characteristic polynomial is then quadratic and has two roots. Hence, only two of the eigenvalues are $O(1)$. If instead we assume $\lambda = O(\Gamma)$, then the dominant terms are the first and second, leading to one eigenvalue, $-k_1 \Gamma$. The remaining two eigenvalues are of an intermediate order. We need not be concerned with details, as we do not need to find the eigenvalues to determine stability.

D.   Proposition 4

The first step in establishing Proposition 4 is to obtain a tight lower bound for $C$ in terms of a formula with form similar to that of $A$ and $B$, that is, $C > (\kappa-m)C_a+C_b$, where $C_a$ and $C_b$ are non-negative. We begin by multiplying (70) by $b$, rewriting in terms of $w$, and using (24) and $p+u = s$ to get

Rewriting (70) as

allows us to complete the calculation as

leading to a useful computational result:

From

we obtain

Applying (80) gives the result

Next, we write

Then $u \le q$ and $bq = 1$ combine to give $bu \le 1$, so

Finally, we substitute (82) and (83) into (29) to obtain the desired lower bound

From (28) and (84), we see immediately that stability is guaranteed if $m \le \kappa$. It remains to be shown that $m \le 0.75$ is sufficient to guarantee stability when $m > \kappa$. The first step in this demonstration is to rewrite (84) as

where

Given $h < b$, $hu+\sigma = (1-\sigma)bu+\sigma \le (1-\sigma)+\sigma = 1$, and $m > \kappa$, we have $A \ge C_1$ and $C_2 \ge C_1$, leaving $C_1 > 0$, $C_3 > 0$, and $B > 0$ to be confirmed.

From the formula for $C_1$ and $hu \sim 1-\kappa+(1-m)\sigma z \ge 1-\kappa$, we have

For $C_3$, we first note that $s \ge 0$ requires $y \le 1/(1-m) \le 4$, and we also have $bhu > hu \ge 1-\kappa$. If $m > bhu$, then

and

If $m \le bhu$, then $1+y(m-bhu) \le 1$ and

To show $B \ge 0$, we first note the result

This allows us to eliminate the $z$ terms from (28), leaving

Multiplying by $1-\kappa$ and using $b > h$, this becomes

With $m \le 0.75$, there is no backward bifurcation; hence, solutions occur only if $c > 0$; along with the modeling restriction $f \le 1$, we then have

Substituting into the previous inequality yields

This completes the proof of Proposition 4.