Persistence, Coexistence and Boundary Transcritical Relays for Multi-Strain Epidemic Models

1. Introduction

Definition 1

(Positive / non-negative ODEs [Pos]). A dynamical system is calledpositive[1] ornon-negative[2] if the non-negative orthant

$${\mathbb{R}}_{\ge 0}^n:=\{x\in {\mathbb{R}}^n:x_i\ge 0,i=1,\dots ,n\}$$

is forward invariant under the flow.

Our research started as an investigation of the escape paths from the boundary of some multi-strain models, a mathematical epidemiology (ME) topic which may be traced back at least to [3,4,5]. Multi-strain models generalize both ME single-strain ODEs and Lotka–Volterra ODEs, and render their common underlying structure visible; a reasonable first step towards making that structure precise is the unifying CRN language adopted below.

1.1. Why CRN Language?

An important characteristic of the models in mathematical epidemiology (ME), ecology, and immuno-virology, etc is that important phenomena occur when trajectories reach the boundary of the positive orthant, and hence that the notion of forward–invariant coordinate faces is crucial (these appear in ME as disease–free or partially disease–free manifolds). Now, while boundary faces have been touched upon also in classical ODE theory, the systematic study of phenomena arising for positive ODEs was initiated in CRN works, see notably [6,7], who showed that all positive polynomial models admit (non-unique) mass-action representations.

The next crucial development for the CRN-ME interface was the discovery in [8] that forward–invariant coordinate faces correspond to the combinatorial object of siphons in the theory of Petri nets, and that boundary $\omega -$ limit points with interior starting points must lie in siphon faces. Note also that we observed in [9] that minimal siphons are close conceptually to “infection strains", with the DFE face being the intersection of all minimal siphon faces, that on all siphon faces the Jacobian has a triangular structure [10], and that a transversal Jacobian block on a siphon face is Metzler [11].

1.2. Main Results

Beyond results for specific multi-strain epidemic models, the notions introduced below (reproduction functions attached to siphon faces, relay graphs, boundary transcritical relays) are meant as a common language for invasion, persistence, coexistence, and replacement, applicable beyond mathematical epidemiology; see the abstract for the broader motivation.

Our first new result is a persistence theorem for two-strain models, proved as part of a general acyclic-boundary-Morse-decomposition framework (Theorem 1 and Corollary 1 of §Section 3) that shows, under explicit boundary assumptions, that mutual invasibility implies uniform persistence. The criterion is expressed entirely in terms of siphon-face invasion quantities and provides a practical persistence test based on boundary dynamics.

A second group of results concerns coexistence. For the permanent-immunity Feng–Gavish model, Theorem 3 (Coexistence criterion for the permanent-immunity Feng–Gavish model) shows that invasibility of the less fit strain implies coexistence; the converse is conjectured (the permanent-immunity instance of Conjecture 1) but not proved in general. More generally, Theorem 4 (Coexistence criterion for two-strain models with essentially rank-one siphon matrices and one input species) yields an explicit coexistence criterion for a broad class of models in terms of invasibility numbers and a reduced scalar polynomial.

Our third contribution concerns boundary transcritical relay (BTR) phenomena. Such transitions are frequently observed in explicit models: a resident equilibrium loses stability at the same parameter value at which a new equilibrium appears on a larger invariant face. We show that this phenomenon is a structural consequence of three features shared by a large class of positive ODEs: (a) face invariance induced by siphons; (b) factorization of the vector field in transversal variables; (c) Metzler structure of the corresponding transversal Jacobian blocks.

Theorem 7 (Local boundary transcritical law) formalizes this mechanism by proving that, under natural structural assumptions, crossing an invasion threshold simultaneously destroys transversal stability of a resident equilibrium and creates a successor equilibrium on an adjacent inhabited siphon face.

These results lead naturally to relay graphs and relay tables, which organize equilibria according to the lattice of inhabited siphon faces and the associated invasion inequalities. For the Feng–Gavish model, Theorem 10 (LAS–CEP for the Feng–Gavish model) yields a complete local-stability and exclusion partition of the boundary equilibria.

These results separate into two logically different layers. Reproduction numbers, local invasion (loss of transversal hyperbolicity), and the boundary transcritical relay theorem are local: read off directly from the transversal Jacobian at a single boundary equilibrium, with no global argument. Persistence and coexistence are global: persistence needs, on top of local invasion at every boundary equilibrium, a global acyclicity hypothesis on the boundary Morse decomposition (Theorem 1); coexistence needs, on top of local invasion at the relevant boundary equilibria, a global algebraic existence argument (Theorem 4). Local invasion is necessary for both but sufficient for neither. This local/global split reflects a broader shift of viewpoint, from isolated equilibria to the partially ordered family of inhabited siphon faces and the invasion relations between them, that leads to algorithmic procedures for enumerating equilibria, building relay graphs and tables, and predicting replacement and coexistence as parameters vary.

1.3. Structure of the Paper

Section 2 presents some background results from CRN and ME.

Section 3 presents the general persistence framework and its two-strain specialization, with a worked example in §Section 3.1.

Section 4 introduces a two-strain model with cross-immunity and ADE of Gavish, chosen to illustrate the relay and persistence results.

Section 5 is dedicated to the question of whether mutual invasibility implies the existence of coexistence equilibria.

Section 6 provides a boundary transcritical relay result, and a relay table for the Feng-Chung-Gavish model.

Section 7 provides some conclusions.

2. Mathematical Background

The ODEs we work on belong to the following class:

Definition 2

(Stoichiometric representation and chemical ODEs [Sto]). A stoichiometric representation of f is a pair $(\mathrm{\Gamma },r)$ such that

$$f\left(x\right)=\mathrm{\Gamma }r\left(x\right),r\left(x\right)\ge 0,$$

where $\mathrm{\Gamma }\in {\mathbb{R}}^{n\times n_R}$ is constant and $r:{\mathbb{R}}_{\ge 0}^n\to {\mathbb{R}}_{\ge 0}^{n_R}$ is locally Lipschitz.

An ODE is called chemical if it admits a stoichiometric representation such that

$${\mathrm{\Gamma }}_{i\rho }<0\Rightarrow r_{\rho }\left(x\right)=x_i{\tilde{r}}_{\rho }\left(x\right),{\tilde{r}}_{\rho }\left(x\right)\ge 0,$$

and, moreover, $r_{\rho }$ depends only on variables in $supp\left(r_{\rho }\right)\subseteq {\mathrm{\Gamma }}_i^{-}$ whenever ${\mathrm{\Gamma }}_{i\rho }<0$. Furthermore, we assume that each rate $r_{\rho }\left(x\right)$ is monotone non-decreasing with respect to each variable $x_k$ for $k\in supp\left(r_{\rho }\right)$.

The class of chemical ODEs is very general and is in fact a classic favorite in ME (more than mass-action, which is best viewed in our context as a polynomial approximation of chemical ODEs). Note that for computational implementations, we restricted always below to rational models, which we call “perturbed" mass-action, and define informally as models which become mass-action after a subset of parameters is set to 0. This restriction provides the possibility of investigating first the computationally easier mass-action case.

2.1. Siphons and Persistence

Angeli, de Leenheer and Sontag [8], Anderson [12] and Shiu and Sturmfels [13] proved that a the boundary face $F_{\Sigma }$ associated to a nonempty set $\Sigma$ of variables being 0 is forward-invariant for a chemical ODE iff it is a siphon in the sense of definition 3. This property may be taken as the definition of siphons.

For completeness, we provide also the “combinatorial siphon definition" encountered in the CRN literature.

Definition 3

(siphon, minimal siphon, total siphon/DFE [sip])).

• Asiphon set$\Sigma \subseteq \mathcal{S}$ is a nonempty subset of species such that whenever a species in Σ appears in a “product complex" in the RHS of a reaction, at least one species in Σ must appear also in the LHS of the reaction (the “reactant complex").
• A siphon isminimalwhen it contains no other siphon included within.
• The union of all minimal siphon will be calledtotal siphon/DFE .

The total siphon/disease free equilibrium set (DFE) is a fundamental object of study in ME. However, its rigorous definition is only available via CRNT.

The siphon property is fundamental for understanding persistence.

Definition 4

(persistence and uniform persistence). Let $\dot{x}=f\left(x\right)$ be a positive ODE defined on a compact, forward-invariant set $E\subset {\mathbb{R}}_{+}^n$, with semiflow ${\Phi }_t$. The ODE ispersistentif for every component i and every trajectory starting in $\mathrm{int}\left(E\right)$,

$$\underset{t\to \infty }{lim\; inf}{x}_i\left(t\right)>0.$$

The ODE isuniformly persistentif there exists $\epsilon >0$ such that

$$\underset{t\to \infty }{lim\; inf}{x}_i\left(t\right)\ge \epsilon$$

uniformly for all trajectories starting in $\mathrm{int}\left(E\right)$.

[8] showed that if the $\omega$-limit set does not intersect siphon boundaries except at equilibria, and if all trajectories starting on non-semilocking boundaries eventually leave those boundaries, then the system exhibits persistence.

The fundamental problems of persistence, permanence, and extinction have been studied extensively, in particular in the context of generalized Lotka–Volterra systems [14,15,16,17,18], for more general mathematical ecology-epidemiology models [19,20,21] and also for particular CRNs [8,22,23,24,25]. Intriguingly, these authors used different methods, and their results have little intersection.

The potential of siphons in ME is illustrated by a recent proof of the folklore result that the Jacobian has a triangular block form on siphon faces [10], and [11] further noted that transversal Jacobians on siphons have the fundamental Metzler/cooperativity property. This property further ensures the existence of right and left Perron eigenvectors which provide escape directions from siphon faces, and Lyapunov functions at each resident fixed point, respectively. The result of [10] is reviewed in the next subsection.

2.2. The Next-Generation Matrix Theorem and Beyond: From Siphons to Right and Left Perron Eigenvectors of Metzler Transversal Jacobians

There is no doubt that the next-generation matrix (NGM) theorem is a fundamental law of mathematical epidemiology. This classical result [26,27,28], probably the most cited in the field of positive ODEs, is emblematic of the fragmentation of positive ODE sciences, by being not enough known outside ME.

It provides a striking invasion criterion of the form $R_0>1$ at the DFE ($R_0$ is called reproduction number), which is motivated by the theory of population models [29] and their branching process approximations [30]. This result implicitly relies on three distinct mechanisms that are not enough disentangled:

(i) The fact that forward invariance of a boundary face enforces a triangular block structure of the Jacobian evaluated on it [10], with the obvious ensuing reduction of stability to that of the lower–dimensional tangential and transversal blocks – see for example [31].

(ii) The Metzler structure of the transversal Jacobian block on a siphon face, recently proved in [11].

(iii) Regular splitting. A convenient feature of Metzler matrices, and hence of the models studied in ME, is that the classic stability criteria for the transversal Jacobian may be replaced by a criterion specific to Metzler matrices which admit regular splittings $M=F-V$, with $F,V^{-1}$ non-negative componentwise. For such matrices, the equivalence between the sign of the spectral abscissa of M and that of $R_0-1$, where $R_0$ is the spectral abscissas of $F{V}^{-1}$, goes back to [27,32,33]. At first, this equivalence suggests preferring spectral abscissa, like in ecology, which are unique, while regular splittings and reproduction numbers aren’t. However, in this paper we use also the latter, both due to empirical observations that sparseness renders reproduction numbers often simpler, and to possible probabilistic interpretations which might turn out useful in future works.

Definition 5

(ME-type siphon, next generation matrix [ME-]). Let $\mathbf{x}$ denote the variables which are 0 on a siphon face $L_{\sigma }$, let $\mathbf{y}$ denote the positive (resident) variables, and recall that the Jacobian $J_x\left(\mathbf{y}\right)$ on $L_{\sigma }$ is Metzler (equivalently, the ODE is locally increasing in $\mathbf{y}$ on $L_{\sigma }$) [10,11]. We will call ME-type siphon a siphon which admits some regular splitting $J_x=F-V$ in the sense of [27,32] (F is non-negative, an $V^{-1}$ exists and is non-negative).

The matrix $K=F{V}^{-1}$ is called next generation matrix (NGM).

Remark 1

(Three aspects of NGM theory are relevant at every siphon face). With NGM decomposed as above, we realize immediately that its three parts identified above are relevant at every siphon face: the first two apply always, and the third may be occasionally useful, when regular splitting to the transversal Jacobian on that face may be applied.

Definition 6

(transversal Jacobian blocks [tra]). For a chemical dynamical system, put the (Metzler) Jacobian of the DFE equations with respect to the DFE variables in its Perron-Frobenius normal form. The irreducible diagonal blocks of this decomposition will be called transversal Jacobian blocks.

2.3. What is a Multi-Strain Model?

We have now the building bricks to define multi-strain systems.

The next three definitions separate notions that are often bundled together: a purely combinatorial siphon decomposition, a multiplicative factorization of the resulting blocks, and the further ME-admissibility (regular splitting) of those blocks. Each is strictly weaker than the next, and results below are stated at the weakest level they actually require.

Definition 7

(Strain decomposition). Ak-strain decompositionof a chemical ODE with variables $(s,x^1,\dots ,x^k)\in {\mathbb{R}}_{\ge 0}^m\times {\prod }_{j=1}^k{\mathbb{R}}_{\ge 0}^{d_j}$ is a collection of k pairwise disjoint minimal siphons ${\Sigma }_1,\dots ,{\Sigma }_k$, with $x^j$ the variables of ${\Sigma }_j$ and s the remaining (susceptible/uninfected) variables.

Definition 8

(Multiplicative strain representation [Mul]). A k-strain decomposition ismultiplicativeif the dynamics factor as

$$\dot{s}=f_s(s,x^1,\dots ,x^k),{\dot{x}}^j=M_j(s,x^1,\dots ,x^k)x^j,j=1,\dots ,k,$$

(1)

with $f_s,M_j$ locally Lipschitz.

Definition 9

(k–strain ME system, and reproduction functions [k]). A multiplicative k-strain decomposition is ak-strain ME systemif each block is ME-admissible: $M_j=F_j-V_j$ is a regular splitting (regular splitting) of the Metzler transversal Jacobian $J_{{\Sigma }_j}^{\perp }$, i.e. $F_j,V_j\ge 0$ and $V_j$ is a nonsingular M-matrix.

For each strain j, define the reproduction function

$$R_j\left({\mathbf{y}}_j\right):=\rho \left(K_j\left({\mathbf{y}}_j\right)\right),$$

where the variables ${\mathbf{y}}_j$ which are not zero on block j are left free, will be calledR-reproduction functions(associated with the block j), and where $K_j\left({\mathbf{y}}_j\right)=F_j\left({\mathbf{y}}_j\right)V_j{\left({\mathbf{y}}_j\right)}^{-1}$ is the next–generation matrix associated with the Metzler transversal Jacobian of strain j.

This gives the hierarchy

$$k-\mathrm{strain}\mathrm{decomposition}\supseteq \mathrm{multiplicative}k-\mathrm{strain}\mathrm{decomposition}\supseteq k-\mathrm{strain}\mathrm{ME}\mathrm{system},$$

so that a property proved from Definition 7 or 8 alone applies more broadly than one requiring the full ME-admissible structure of Definition 9.

Definition 10

(reproduction numbers and invasion numbers [rep]).

1.

The evaluations

$$R_i=R_i\left(E_0\right),i=1,...,k,$$

(2)

where $E_0$ is the DFE, and $R_i\left({\mathbf{y}}_i\right)$ are the reproduction functions, will be calledbasic reproduction numbers .

2.

In the case when each minimal siphon has a unique resident fixed point $E_j,j=1,...,k,$

$$R_i^{\tilde{j}}:=R_i\left(E_j\right)$$

(3)

will be called theinvasion number of invading block i on (resident) block j .

2.4. Rank one ME Models with Essentially Rank One Transversal Jacobians

The regular splitting $J_{\Sigma }^{\perp }=F_{\Sigma }-V_{\Sigma }$ ($F_{\Sigma }\ge 0$, $V_{\Sigma }$ a nonsingular M-matrix) used throughout is the standard van den Driessche–Watmough splitting, already required by hypothesis (RS) above; we do not redefine it here. Because $J_{\Sigma }^{\perp }$ is Metzler, the trivial regular splitting $J_{\Sigma }^{\perp }=(J_{\Sigma }^{\perp }-D_{\Sigma })-(-D_{\Sigma })$, with $D_{\Sigma }$ the diagonal part, always exists whenever the diagonal entries are strictly negative; this is elementary and only shows regular splittings always exist. The interesting case is when the resulting infection matrix $F_{\Sigma }=J_{\Sigma }^{\perp }-D_{\Sigma }$ has rank one:

Definition 11

(essentially rank one irreducible Metzler matrix [ess]). An irreducible Metzler matrix will be called essentially rank one if its off-diagonal part is of rank one, i.e. $J_{\Sigma }^{\perp }=w_{\Sigma }{\pi }_{\Sigma }^{\top }-D_{\Sigma }$ for some $w_{\Sigma },{\pi }_{\Sigma }\gg 0$ and diagonal $D_{\Sigma }$.

Remark 2

(Which matrix each Perron vector belongs to). Throughout, $w_{\Sigma }$ and ${\pi }_{\Sigma }$ are the right and left Perron vectors of thenew-infectionsmatrix $F_{\Sigma }=w_{\Sigma }{\pi }_{\Sigma }^{\top }$, with ${\pi }_{\Sigma }$ normalized as a probability vector; this is the convention of the shared notation file, where the siphon index is written σ or k rather than Σ. Every other Perron vector occurring below is obtained from these two by at most one application of $V_{\Sigma }^{-1}$:

Here ${\mathcal{R}}_{\Sigma }={\pi }_{\Sigma }^{\top }{V}_{\Sigma }^{-1}{w}_{\Sigma }$ is the per-susceptible reproduction coefficient, and

$${\tilde{\pi }}_{\Sigma }^{\top }:={\displaystyle \frac{{\pi }_{\Sigma }^{\top }{V}_{\Sigma }^{-1}}{{\mathcal{R}}_{\Sigma }}},{\tilde{\pi }}_{\Sigma }^{\top }{w}_{\Sigma }=1,$$

(4)

is the normalized left Perron eigenvector of $K_{\Sigma }$ – equivalently the normalized left nullvector of $s_{\Sigma }{F}_{\Sigma }-V_{\Sigma }$ at $s_{\Sigma }=1/{\mathcal{R}}_{\Sigma }$ – which is the vector that carries the Lyapunov weights of §Section 3. The last row of the table restates Lemma 5.

Two identifications must therefore be resisted: the right Perron vector of $K_{\Sigma }$ is $w_{\Sigma }$, whereas the escape direction $ker{J}_{\Sigma }^{\perp }$ is $V_{\Sigma }^{-1}{w}_{\Sigma }$; and the left Perron vector of $F_{\Sigma }$ is ${\pi }_{\Sigma }^{\top }$, whereas that of $K_{\Sigma }$ is ${\pi }_{\Sigma }^{\top }{V}_{\Sigma }^{-1}$. Both distinctions collapse exactly when $V_{\Sigma }=v_{\Sigma }I$ is scalar, which is the case in every worked example below (each strain having the single removal rate $v_k={\gamma }_k+\mu$), so no computation of this paper is affected; they must be kept for blocks carrying several distinct removal rates.

Remark 3

(Why rank-one matters). Multi-strain models generalize Lotka–Volterra models along two axes: not every equation need have multiplicative structure, and the scalar multiplications of Lotka–Volterra (scalar siphons) become matrix multiplications in (1) for vector siphons. Rank-one blocks are the natural boundary case of this generalization, and they matter because all the relevant Perron eigenvectors are then explicit, being obtained from the two factors $w_{\Sigma },{\pi }_{\Sigma }$ by a single application of $V_{\Sigma }^{-1}$ (Remark 2). The escape direction from the boundary is $V_{\Sigma }^{-1}{w}_{\Sigma }$, which spans $ker{J}_{\Sigma }^{\perp }$ at criticality by Lemma 5 and appears directly in endemic-point coordinates in rank-one models [11,34,35,36,37]; the normalized left Perron eigenvector ${\tilde{\pi }}_{\Sigma }^{\top }={\pi }_{\Sigma }^{\top }{V}_{\Sigma }^{-1}/{\mathcal{R}}_{\Sigma }$ of (4) builds the local repelling Lyapunov function used throughout §Section 3 (see [38] for a unification: this property holds for systems with disjoint scalar minimal siphons, with at most one exception, which must be rank one). This is the notion used by Lemma 5 and the coexistence results built on it (§Section 5.3).

Definition 12

(simple ME models [sim]). An essentially rank one ME model for which the NGM of the DFE has only zero eigenvalues and unconditionally positive eigenvalues, one for each minimal siphon block, which were called invasion functions in [9], will be called simple.

2.5. Is the Persistence Criterion Explicit for ME Models?

A fundamental result of [39] – see also [40] – based on Morse decomposition theory, implies that when the maximal invariant sets on the boundary are fixed points, uniform persistence follows from two ingredients: transversal repulsivity of every boundary equilibrium,

$$\underset{(\Sigma ,E)}{min}{\alpha }_{\Sigma }\left(E\right)>0,{\alpha }_{\Sigma }\left(E\right):=\alpha \left(J_{\Sigma }^{\perp }\left(E\right)\right),$$

(5)

equivalently every relevant ME invasion number exceeding one; and acyclicity of the boundary connection graph, which prevents tangential motion along siphon faces from compensating for transversal repulsion. Because the transversal Jacobians are Metzler, positivity of the invasion numbers yields more than instability: Perron–Frobenius theory supplies the positive left eigenvectors used to build the local repelling Lyapunov functions of §Section 3. In the one- and two-siphon cases acyclicity is essentially automatic (Corollary 1); for general n-strain models it becomes a genuinely nontrivial combinatorial problem.

2.6. Coexistence, the Relay Property, and the Complete Exclusion Partition Principle

Beyond persistence, another fundamental question in ME is coexistence, i.e. existence of a strictly positive endemic equilibrium (EE), which typically appears only after certain boundary equilibria lose stability. In all explicit examples known to us, persistence and EE existence coincide, but whether this holds beyond low-dimensional examples remains an open problem.

Along the “escape paths from the boundary” generated by positive invasion numbers (a viewpoint already emphasized in ecology by Schuster and Sigmund [41,42]), one repeatedly observes a relay principle: the same invasion inequality that destroys transversal stability of a boundary equilibrium E also guarantees existence of a successor equilibrium on the adjacent larger invariant face, both triggered by the same threshold crossing ${\alpha }_{\Sigma }\left(E\right)=0$ (equivalently $R_{\Sigma }\left(E\right)=1$). We call this a boundary transcritical relay (BTR); unlike classical transcritical bifurcation (e.g. [43]), the crossing occurs along an entire siphon block rather than a single coordinate, governed by the spectral abscissa of a Metzler transversal block. The precise statement, with its nondegeneracy hypotheses, is Theorem 7 below; the resulting relay tables (Table 3) partition parameter space into regions where exactly one boundary or interior equilibrium is locally asymptotically stable, which we call a complete exclusion partition.

The relay graphs used below are related to the invasion graphs of Hofbauer and collaborators [44,45], which encode which resident communities can be invaded through positive invasion rates, but differ in two respects: nodes are equilibria indexed by inhabited siphon faces, partially ordered by the minimal-siphon lattice, and edges are labeled by explicit invasion numbers, so that adjacency corresponds to distance-one moves in that lattice. This makes the relay graph algorithmic: once the minimal siphons and their reproduction functions are known, one can in principle enumerate candidate equilibria, determine their stability, and predict relay transitions between them.

Figure 1. The relay graph for the two-strain Feng–Gavish model. Nodes correspond to equilibria indexed by inhabited siphon faces: $E_0$ (DFE), one-strain equilibria $E_1\in \mathrm{relint}{F}_1$ and $E_2\in \mathrm{relint}{F}_2$, and the coexistence equilibrium $E_{12}\in \mathrm{relint}\mathrm{\Gamma }$. Dashed arrows are boundary transcritical relays (BTRs): the invasion inequality on each arrow simultaneously destroys transversal stability of the source equilibrium and creates a successor equilibrium on the adjacent siphon face (Thm. 7). LAS conditions at $E_1$ and $E_2$ include the strain-ordering requirement (${\mathcal{R}}_1>{\mathcal{R}}_2$ or ${\mathcal{R}}_2>{\mathcal{R}}_1$, respectively) from Table 3. Dotted arrows lead to two downstream results. Persistence Cor. 1: all four invasion inequalities together imply uniform persistence; in this model the first two (${\mathcal{R}}_k>1$) are implied by the last two, but that implication is not established in general (open problem). Coexistence Thm. 3: the ordered invasion condition guarantees existence of $E_{12}\in {\mathbb{R}}_{>0}^n$; the equivalence persistence ⇔ coexistence holds in all known explicit cases but remains open in general.

Figure 1. The relay graph for the two-strain Feng–Gavish model. Nodes correspond to equilibria indexed by inhabited siphon faces: $E_0$ (DFE), one-strain equilibria $E_1\in \mathrm{relint}{F}_1$ and $E_2\in \mathrm{relint}{F}_2$, and the coexistence equilibrium $E_{12}\in \mathrm{relint}\mathrm{\Gamma }$. Dashed arrows are boundary transcritical relays (BTRs): the invasion inequality on each arrow simultaneously destroys transversal stability of the source equilibrium and creates a successor equilibrium on the adjacent siphon face (Thm. 7). LAS conditions at $E_1$ and $E_2$ include the strain-ordering requirement (${\mathcal{R}}_1>{\mathcal{R}}_2$ or ${\mathcal{R}}_2>{\mathcal{R}}_1$, respectively) from Table 3. Dotted arrows lead to two downstream results. Persistence Cor. 1: all four invasion inequalities together imply uniform persistence; in this model the first two (${\mathcal{R}}_k>1$) are implied by the last two, but that implication is not established in general (open problem). Coexistence Thm. 3: the ordered invasion condition guarantees existence of $E_{12}\in {\mathbb{R}}_{>0}^n$; the equivalence persistence ⇔ coexistence holds in all known explicit cases but remains open in general.

3. Persistence for Eco-Epidemic Models

This section develops general results that generalize results obtained originally for specific examples, but apply to broad classes of positive ODEs. Taken together, they illustrate three distinct ways in which CRNT concepts enter the analysis. Persistence uses siphons and invariant faces; relay phenomena use siphons together with Metzler transversal Jacobians; competitive exclusion partitions use reproduction functions attached to resident faces.

The first result below is a persistence theorem for multi-strain models, stated for two minimal siphons. The underlying persistence mechanism is more general: it depends on the boundary Morse decomposition, on the existence of at least one invading missing siphon at each boundary resident equilibrium, and on acyclicity of the boundary Morse graph.

Theorem 1

(Persistence from acyclic boundary Morse decompositions).

Let Ω be a compact positively invariant set.

Assume that the boundary semiflow admits a finite Morse decomposition whose Morse sets are resident equilibria

$$E_1,\dots ,E_N.$$

Assume further:

(RS)

For each resident equilibrium $E_i$ and each minimal siphon σ absent at $E_i$, the transversal Jacobian admits a regular splitting

$$J_{\sigma }^{\perp }\left(E_i\right)=F_{\sigma }\left(E_i\right)-V_{\sigma }\left(E_i\right),F_{\sigma }\left(E_i\right)\ge 0,V_{\sigma }\left(E_i\right)\mathrm{a}\mathrm{nonsingular}M-\mathrm{matrix}.$$

Put

$$K_{\sigma }\left(E_i\right)=F_{\sigma }\left(E_i\right)V_{\sigma }{\left(E_i\right)}^{-1}.$$

(IC)

For every boundary resident equilibrium $E_i$, there exists at least one absent minimal siphon σ such that

$$\rho \left(K_{\sigma }\left(E_i\right)\right)>1.$$

(A)

The boundary Morse graph is acyclic.

Let ${\partial }_{\Sigma }\Omega$ denote the union of the siphon faces on which the Morse sets $E_1,\dots ,E_N$ lie. Then the system is uniformly persistent with respect to

$${\rho }_{\partial }\left(x\right)=d(x,{\partial }_{\Sigma }\Omega ),$$

that is, there exists $\epsilon >0$ such that ${lim\; inf}_{t\to \infty }{\rho }_{\partial }\left(x\left(t\right)\right)\ge \epsilon$ for every interior solution. Moreover, since ${\partial }_{\Sigma }\Omega$ is the union of the individual missing-siphon faces $\{x_{\sigma }=0\}$ and $d(x,A\cup B)=min\left\{d\right(x,A),d(x,B\left)\right\}$, this is equivalent to ${lim\; inf}_{t\to \infty }\left|x_{\sigma }\left(t\right)\right|\ge {\epsilon }^{\prime }$ (for a possibly smaller ${\epsilon }^{\prime }>0$) separately for every siphon σ that is missing at some Morse set: persistence of the aggregate boundary distance automatically decomposes into persistence of each individual siphon.

Proof. 

Fix a boundary resident equilibrium $E_i$.

By assumption (IC), there exists an absent minimal siphon $\sigma$ such that

$$\rho \left(K_{\sigma }\left(E_i\right)\right)>1.$$

By assumption (RS) and the next-generation-matrix sign theorem,

$$\rho \left(K_{\sigma }\left(E_i\right)\right)>1⟺s\left(J_{\sigma }^{\perp }\left(E_i\right)\right)>0.$$

Hence the equilibrium $E_i$ possesses an unstable transversal direction corresponding to a missing minimal siphon.

Since $J_{\sigma }^{\perp }\left(E_i\right)$ is Metzler [11], there exists a positive left Perron eigenvector ${\pi }_i\gg 0$ satisfying

$${\pi }_i^{\top }{J}_{\sigma }^{\perp }\left(E_i\right)={\lambda }_i{\pi }_i^{\top },{\lambda }_i>0.$$

Let $x_{\sigma }$ denote the variables of the missing siphon. In local coordinates near $E_i$,

$$x_{\sigma }^{\prime }=J_{\sigma }^{\perp }\left(E_i\right)x_{\sigma }+o\left(\right|x_{\sigma }\left|\right).$$

Therefore the linear functional

$$L_i\left(x\right):={\pi }_i^{\top }{x}_{\sigma }$$

satisfies

$${\dot{L}}_i\left(x\right)>0$$

for all interior points sufficiently close to $E_i$.

Consequently no trajectory starting in the interior can remain in a sufficiently small neighborhood of $E_i$ for all future time. Thus every boundary Morse set is isolated and weakly repelling with respect to the interior semiflow.

Collecting the ingredients: $\Omega$ is compact and positively invariant (standing assumption); the boundary semiflow has a finite Morse decomposition $\left\{E_1\right\},\dots ,\left\{E_N\right\}$, each Morse set a single equilibrium and, by the definition of Morse decomposition, compact and isolated as an invariant set; each Morse set is weakly repelling for the interior semiflow, by the Lyapunov argument just given; and the Morse graph is acyclic by (A). These are exactly the hypotheses of the Freedman–Ruan acyclic persistence theorem [39], and no more. It follows that the boundary

$$\partial \Omega$$

is a uniform repeller for interior trajectories.

Hence there exists $\epsilon >0$ such that every interior solution satisfies

$$\underset{t\to \infty }{lim\; inf}d(x\left(t\right),\partial \Omega )\ge \epsilon .$$

Therefore the system is uniformly persistent.

   □

Corollary 1

(Two-siphon persistence theorem).

Assume that there exist exactly two minimal siphons

$${\sigma }_1,{\sigma }_2,{\sigma }_0={\sigma }_1\cup {\sigma }_2,$$

with associated invariant faces [8]

$$F_1=\{x_{{\sigma }_2}=0\},F_2=\{x_{{\sigma }_1}=0\},F_0=F_1\cap F_2.$$

Assume that there exist resident equilibria

$$E_0\in F_0,E_1\in relint\left(F_1\right),E_2\in relint\left(F_2\right),$$

which form a finite Morse decomposition

$$\left\{E_0\right\},\left\{E_1\right\},\left\{E_2\right\}$$

of the boundary semiflow on $(F_0\cup F_1\cup F_2)\cap \Omega$: that is, $E_0,E_1,E_2$ are the only chain-recurrent points of the boundary flow, so that every boundary orbit $x(·)\subset (F_0\cup F_1\cup F_2)\cap \Omega$ has $\omega \left(x\right)\subseteq \{E_0,E_1,E_2\}$, with $\omega \left(x\right)\subseteq \{E_0,E_1\}$ for orbits in $F_1\cap \Omega$ and $\omega \left(x\right)\subseteq \{E_0,E_2\}$ for orbits in $F_2\cap \Omega$ (equivalently, $E_0$ is the unique equilibrium on $F_0$, and $E_1,E_2$ are each the unique equilibrium in the relative interior of $F_1,F_2$ and globally attracting there).

Assume also that hypothesis(RS)of Theorem 1 holds and that

$$R_1\left(E_0\right)=\rho \left(K_{{\sigma }_1}\left(E_0\right)\right)>1,R_2\left(E_0\right)=\rho \left(K_{{\sigma }_2}\left(E_0\right)\right)>1,$$

$$R_1^{E_2}=\rho \left(K_{{\sigma }_1}\left(E_2\right)\right)>1,$$

$$R_2^{E_1}=\rho \left(K_{{\sigma }_2}\left(E_1\right)\right)>1.$$

Assume furthermore that the boundary Morse graph associated with the boundary dynamics is acyclic.

Then the system is uniformly persistent with respect to ${\rho }_{\partial }\left(x\right)=min\left\{\right|x_{{\sigma }_1}|,|x_{{\sigma }_2}\left|\right\}$: there exists $\epsilon >0$ such that ${lim\; inf}_{t\to \infty }\left|x_{{\sigma }_1}\left(t\right)\right|\ge \epsilon$ and ${lim\; inf}_{t\to \infty }\left|x_{{\sigma }_2}\left(t\right)\right|\ge \epsilon$ for every interior solution, by Theorem 1’s decomposition of ${\rho }_{\partial }$ into the two individual siphon faces.

Proof. 

The hypothesis directly supplies the finite boundary Morse decomposition

$$E_0,E_1,E_2$$

required by Theorem 1.

The four invasion inequalities imply condition (IC) of Theorem 1: at $E_0$ both missing siphons are individually repelling, and at $E_1$, $E_2$ the unique missing siphon is repelling.

Hypothesis (A) is supplied by the assumed acyclicity of the boundary Morse graph.

Therefore all assumptions of Theorem 1 hold, and uniform persistence follows.

   □

Remark 4

(Role of the Morse decomposition).

The persistence mechanism depends on three ingredients:

1.

a finite boundary Morse decomposition;

2.

transversal instability of every boundary Morse set in at least one absent-siphon direction;

3.

absence of directed cycles in the boundary Morse graph.

The invasion inequalities

$$\rho \left(K_{\sigma }\left(E\right)\right)>1$$

provide the local repulsion mechanism, while the acyclicity condition controls the global organization of the boundary dynamics.

The theorem does not use relay graphs. The relevant graph is the boundary Morse graph arising from the Morse decomposition of the boundary semiflow. Establishing acyclicity of that graph is generally a separate dynamical problem.

Remark 5

(Why Corollary 1 asks for both DFE invasion inequalities separately). Condition(IC)of Theorem 1 only requires, at each resident equilibrium, thatsomemissing siphon be repelling. At $E_0$ the two missing siphons are ${\sigma }_1$ and ${\sigma }_2$, so(IC)at $E_0$ would already follow from the weaker disjunction $\rho \left(K_{{\sigma }_1}\left(E_0\right)\right)>1$or$\rho \left(K_{{\sigma }_2}\left(E_0\right)\right)>1$, equivalently $\rho \left(K_{{\sigma }_0}\left(E_0\right)\right)>1$ whenever $K_{{\sigma }_0}\left(E_0\right)$ is block diagonal in the two siphon blocks. The corollary states the stronger conjunction $R_1\left(E_0\right)>1$and$R_2\left(E_0\right)>1$ instead, for two reasons unrelated to the abstract theorem’s logic: it is the form directly verifiable from the two one-strain reproduction numbers without first checking block-diagonality of $K_{{\sigma }_0}\left(E_0\right)$, and it is automatically satisfied whenever $E_1,E_2$ are genuine one-strain endemic equilibria of a model without single-strain backward bifurcation (existence of a positive $E_k\in relint\left(F_k\right)$ then already forces $R_k\left(E_0\right)>1$). Models with single-strain backward bifurcation, where $E_1$ or $E_2$ can persist with $R_1\left(E_0\right)<1$ or $R_2\left(E_0\right)<1$, are exactly the case where the two hypotheses part ways and the conjunction is the one that should be checked directly.

3.1. A Worked Example: Persistence in a Two-Strain SI2V Vaccination Model

We consider epidemic systems

$$x^{\prime }=f\left(x\right),x\in \mathrm{\Gamma },$$

where $\mathrm{\Gamma }$ is a compact forward-invariant set (often a simplex induced by $N^{\prime }\le \mathrm{\Lambda }-\mu N$), where $\mathrm{\Lambda }$ and $\mu$ denote the birth and death rates, respectively.

The general persistence machinery needed here — an acyclic boundary Morse decomposition theorem (Theorem 1) and its two-siphon equilibrium specialization (Corollary 1), both proved via the Freedman–Ruan theorem [39] by constructing a local repelling Lyapunov function from the Perron eigenvector of each transversal Metzler block $J_{\Sigma }^{\perp }\left(E\right)$ at every boundary equilibrium, plus an acyclicity argument for the boundary connection graph — is stated and proved once, for general multi-strain ME models, in §Section 3. In particular, the remark following Corollary 1 there explains why the DFE invasion condition is stated as the conjunction $R_1\left(E_0\right)>1$and$R_2\left(E_0\right)>1$ of the two one-strain reproduction numbers, rather than as a single combined spectral radius condition $\rho \left(K_{{\Sigma }_0}\left(E_0\right)\right)>1$ — exactly how it is verified below in the Feng–Gavish application (Theorem ).

Example 1

(Two-strain SI2V vaccination model with scalar blocks [Two]). We consider the system, on $\mathrm{\Gamma }=\{(s,i_1,i_2,v)\ge 0:s+i_1+i_2+v\le \mathrm{\Lambda }/\mu \}$,

$$\left\{\begin{array}{c}{s}^{\prime }=\mathrm{\Lambda }-(\mu +\nu )s-{\beta }_{1}s{i}_1-{\beta }_{2}s{i}_2,\hfill \\ i_1^{\prime }={\beta }_{1}s{i}_1-\mu i_1,\hfill \\ i_2^{\prime }={\beta }_{2}s{i}_2-\mu i_2,\hfill \\ v^{\prime }=\nu s-\mu v,\hfill \end{array}\right.$$

with $s,i_1,i_2,v$ susceptible, strain-1-infected, strain-2-infected, and vaccinated individuals, $\mathrm{\Lambda },\mu ,\nu ,{\beta }_1,{\beta }_2>0$ the recruitment, death, vaccination, and transmission rates. The two minimal siphons are ${\Sigma }_1=\left\{i_1\right\}$, ${\Sigma }_2=\left\{i_2\right\}$, with faces $F_1=\{i_2=0\}$, $F_2=\{i_1=0\}$, $F_0=\{i_1=i_2=0\}$.

Since each infection block is scalar, $K_{{\Sigma }_j}\left(y\right)={\mathcal{R}}_{j}s$ (${\mathcal{R}}_j={\beta }_j/\mu$), so $R_j\left(y\right)={\mathcal{R}}_{j}s$; at the DFE $E_0=(s_0,0,0,v_0)$, $s_0=\mathrm{\Lambda }/(\mu +\nu )$, this gives $R_j={\mathcal{R}}_j{s}_0$, and at the one-strain equilibrium $E_k=(s_k,i_k^{*},0,v_k^{*})$ ($s_k=1/{\mathcal{R}}_k$, when $R_k>1$) the invasion numbers are $R_2^{P_1}=R_2\left(E_1\right)={\mathcal{R}}_2{s}_1$, $R_1^{P_2}=R_1\left(E_2\right)={\mathcal{R}}_1{s}_2$.

Checking Corollary 1. On $F_k=\{i_{3-k}=0\}$ the system reduces to a standard SIS-type model in $(s,i_k,v)$, whose force of infection is linear in s; a standard Lyapunov argument gives a unique equilibrium ($E_0$ if $R_k\le 1$, else $E_k$) that is globally asymptotically stable on $F_k$, so $\left\{E_0\right\},\left\{E_1\right\},\left\{E_2\right\}$ form the required Morse decomposition, with (RS) and acyclicity immediate since all transversal blocks are scalar. Persistence therefore reduces exactly to $R_1>1,R_2>1,R_2^{P_1}>1,R_1^{P_2}>1$, the usual two-strain invasion conditions.

4. A Two Strain Model with ADE and Immunity Waning, That May Exhibit Hopf Bifurcations: [3,4,5,46,47,48,49,50];GavScan.wl

4.1. Background

The two strain model with ADE may be traced back to [3,4], and is appropriate for modelling simultaneous epidemics with different pathogens, like for example Dengue and Zika. Subsequently, two-strain models which add further compartments allowing for temporary cross-immunity have been developed in the works of Aguiar, Stollenwerk and Kooi [46,51,52,53,54], and [40,55] examined the effects of single-strain vaccination on the dynamics of an epidemic multi-strain Dengue model (see also [56] for a first public notebook). This model has also been used for several strains of pathogens (without immunity-effectors compartments), and in ecology [57,58].

The Feng-Chung-Gavish model illustrates the relay and persistence results in a simple, but non-trivial setting.

Figure 2. Schematic diagram of disease dynamics for two co-circulating strains. The model includes: autocatalytic primary infections (${\beta }_i$), green arrows: primary infections catalyzed by secondary infections (${\beta }_i{\eta }_i$), amplified secondary infections catalyzed by primary infections (${\beta }_j{\sigma }_j$), orange arrows: amplified auto-catalyzed secondary infections (${\beta }_j{\sigma }_j{\eta }_j$), recoveries (${\gamma }_i$), and waning immunity (${\theta }_i$).

Figure 2. Schematic diagram of disease dynamics for two co-circulating strains. The model includes: autocatalytic primary infections (${\beta }_i$), green arrows: primary infections catalyzed by secondary infections (${\beta }_i{\eta }_i$), amplified secondary infections catalyzed by primary infections (${\beta }_j{\sigma }_j$), orange arrows: amplified auto-catalyzed secondary infections (${\beta }_j{\sigma }_j{\eta }_j$), recoveries (${\gamma }_i$), and waning immunity (${\theta }_i$).

4.2. Model and Siphons

Examining the reaction network (23 reactions) reveals that the minimal siphons are precisely the two infections $i_1,i_{21}$ and $i_2,i_{12}$. This is emphasized by writing the ODE with $X=(s,i_1,i_{21},i_2,i_{12},r_1,r_2,r_{12})$ ordered so that each minimal siphon occupies a contiguous block:

$$\left\{\begin{array}{c}{s}^{\prime }=\mathrm{\Lambda }-\mu s-{\displaystyle \sum _{k=1}^2}{\beta }_k{J}_{k}s+{\displaystyle \sum _{k=1}^3}{\theta }_k{r}_k,\hfill \\ \left(\begin{array}{c}{i}_1^{\prime }\\ i_{21}^{\prime }\end{array}\right)=\underset{M_1(s,r_2)}{\underbrace{\left[{\beta }_1\left(\begin{array}{c}s\\ {\sigma }_1{r}_2\end{array}\right)\left(\begin{array}{cc}1& {\eta }_1\end{array}\right)-v_1\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\right]}}\left(\begin{array}{c}{i}_1\\ i_{21}\end{array}\right),\hfill \\ \left(\begin{array}{c}{i}_2^{\prime }\\ i_{12}^{\prime }\end{array}\right)=\underset{M_2(s,r_1)}{\underbrace{\left[{\beta }_2\left(\begin{array}{c}s\\ {\sigma }_2{r}_1\end{array}\right)\left(\begin{array}{cc}1& {\eta }_2\end{array}\right)-v_2\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\right]}}\left(\begin{array}{c}{i}_2\\ i_{12}\end{array}\right),\hfill \\ r_1^{\prime }={\gamma }_1{i}_1-\left({\beta }_2{\sigma }_2{J}_2+\mu +{\theta }_1\right)r_1,\hfill \\ r_2^{\prime }={\gamma }_2{i}_2-\left({\beta }_1{\sigma }_1{J}_1+\mu +{\theta }_2\right)r_2,\hfill \\ r_{12}^{\prime }={\gamma }_1{i}_{21}+{\gamma }_2{i}_{12}-(\mu +{\theta }_{12})r_{12},\hfill \end{array}\right.$$

(6)

with $v_k:={\gamma }_k+\mu$, and forces of infection

$$J_1=i_1+{\eta }_1{i}_{21},J_2=i_2+{\eta }_2{i}_{12}.$$

.Definition 13 (Antibody-Dependent Enhancement (ADE) [Ant]).For a multi-strain model, we say a secondary infection parameter η has Antibody-Dependent Enhancement (ADE) if it is bigger than 1 [3,4].

Since the system is positive and

$$N^{\prime }=\mathrm{\Lambda }-\mu N,$$

where N is the sum of all compartments, it follows that

$$N\left(t\right)\le max\left\{N\left(0\right),{\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\right\},\mathrm{and}\underset{t\to \infty }{lim}N\left(t\right)={\displaystyle \frac{\mathrm{\Lambda }}{\mu }}.$$

Hence, the system is dissipative on the positively invariant compact set

$$\mathrm{\Gamma }=\left\{x\in {\mathbb{R}}_{\ge 0}^8:N\le {\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\right\}.$$

Denote the two minimal siphons and their union by

$${\Sigma }_1=\{i_1,i_{21}\},{\Sigma }_2=\{i_2,i_{12}\},{\Sigma }_0={\Sigma }_1\cup {\Sigma }_2,$$

with invariant faces

$$F_1=\{x_{{\Sigma }_2}=0\},F_2=\{x_{{\Sigma }_1}=0\},F_0=\{x_{{\Sigma }_0}=0\}.$$

4.3. Reproduction Functions

The next-generation matrices are rank one on each block, yielding the reproduction functions

$$R_1\left(y\right)={\mathcal{R}}_1\left(s+{\eta }_1{\sigma }_1{r}_2\right),R_2\left(y\right)={\mathcal{R}}_2\left(s+{\eta }_2{\sigma }_2{r}_1\right),{\mathcal{R}}_i={\displaystyle \frac{{\beta }_i}{v_i}},v_i={\gamma }_i+\mu ,i=1,2,$$

which are independent of the immunity parameters $\theta$.

At the DFE

$$E_0=(s_0,0,\dots ,0),s_0={\displaystyle \frac{\mathrm{\Lambda }}{\mu }},$$

$$R_1:=R_1\left(E_0\right)=s_0{\mathcal{R}}_1,R_2:=R_2\left(E_0\right)=s_0{\mathcal{R}}_2.$$

4.4. Boundary Equilibria and Invasion Numbers

If $R_1\left(E_0\right)>1$, there exists a unique equilibrium $E_1\in relint{F}_1$:

$$s_1={\displaystyle \frac{1}{{\mathcal{R}}_1}},i_1^{*}=(1-{\tilde{\gamma }}_1)(s_0-s_1),r_1^{*}={\tilde{\gamma }}_1(s_0-s_1),{\tilde{\gamma }}_1={\displaystyle \frac{{\gamma }_1}{{\gamma }_1+{\theta }_1+\mu }}.$$

(7)

Remark 6

(Explicit invasion numbers).

It is conjectured and proved below in some examples that persistence and the relay table are completely determined by $R_1\left(E_0\right),R_2\left(E_0\right),$ and the two quantities

$$\begin{array}{cc}\hfill R_2^{E_1}:=R_2\left(E_1\right)={\mathcal{R}}_2\left(s_1+{\eta }_2{\sigma }_2{r}_1^{*}\right)& ={\displaystyle \frac{R_2}{R_1}}\left(1+{\eta }_2{\sigma }_2{\tilde{\gamma }}_1\left({\displaystyle \frac{s_0}{s_1}}-1\right)\right),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill R_1^{E_2}:=R_1\left(E_2\right)={\mathcal{R}}_1\left(s_2+{\eta }_1{\sigma }_1{r}_2^{*}\right)& ={\displaystyle \frac{R_1}{R_2}}\left(1+{\eta }_1{\sigma }_1{\tilde{\gamma }}_2\left({\displaystyle \frac{s_0}{s_2}}-1\right)\right).\hfill \end{array}$$

(9)

which we will call invasion numbers, following the ME literature (and diverging from the ecology literature, which reserves this term for spectral abscissas of Jacobians).

Both are increasing in the ADE parameters ${\eta }_j$ in the relevant ranges, and when ${\eta }_j=0$ reduce to $R_j/R_i$, a classical formula.

Remark 7

(Explicit ADE invasion thresholds). Since $R_1\left(E_0\right)>1$ forces $s_1<s_0$ (susceptible depletion at $E_1$) and ${\tilde{\gamma }}_1,{\sigma }_2>0$, the prefactor ${\sigma }_2{\tilde{\gamma }}_1(s_0/s_1-1)$ in (8) is strictly positive whenever $E_1$ exists. Solving $R_2^{E_1}=1$ for ${\eta }_2$ therefore gives an explicit invasion threshold:

$${\eta }_2>{\eta }_{2,\mathrm{inv}}⟺R_2^{E_1}>1,{\eta }_{2,\mathrm{inv}}:={\displaystyle \frac{R_1/R_2-1}{{\sigma }_2{\tilde{\gamma }}_1(s_0/s_1-1)}}.$$

Symmetrically, whenever $E_2$ exists ($R_2\left(E_0\right)>1$, so $s_2<s_0$),

$${\eta }_1>{\eta }_{1,\mathrm{inv}}⟺R_1^{E_2}>1,{\eta }_{1,\mathrm{inv}}:={\displaystyle \frac{R_2/R_1-1}{{\sigma }_1{\tilde{\gamma }}_2(s_0/s_2-1)}}.$$

When $R_1>R_2$ (strain 1 fitter), ${\eta }_{2,\mathrm{inv}}>0$: this quantifies precisely how much antibody-dependent enhancement of the less fit strain 2 is required to invade the fitter resident’s equilibrium $E_1$, as an explicit function of the two one-strain reproduction numbers, the cross-enhancement susceptibility ${\sigma }_2$, and the strain-1 recovery fraction ${\tilde{\gamma }}_1$. In the same regime ${\eta }_{1,\mathrm{inv}}<0$, so $R_1^{E_2}>1$ automatically for every ${\eta }_1\ge 0$: the fitter strain needs no enhancement to invade the less fit resident’s equilibrium $E_2$. The roles reverse symmetrically when $R_2>R_1$.

4.5. Persistence for the Feng–Gavish Model

Lemma 1

(Boundary-face geometry for two-strain models). Let the two minimal strain siphons be disjoint: ${\Sigma }_1\cap {\Sigma }_2=⌀$. Let $F_1=\{x_{{\Sigma }_2}=0\}$, $F_2=\{x_{{\Sigma }_1}=0\}$, $F_0=F_1\cap F_2$. Then $F_0$ is the disease-free face. Since $F_1$ and $F_2$ are forward invariant, any boundary orbit connecting the relative interior of $F_1$ to the relative interior of $F_2$ must pass through $F_0$.

Argument:

Forward invariance of each face prevents an orbit starting in $F_1$ from acquiring the missing ${\Sigma }_2$-block. Thus a boundary connection between the two one-strain faces can only occur through $F_0$.

Lemma 2

(Global attractivity of $E_1$ on $F_1$, and symmetrically of $E_2$ on $F_2$). If $R_1\left(E_0\right)>1$, every orbit on $F_1\cap \mathrm{\Gamma }$ converges to $E_0$ or to $E_1$; the symmetric statement holds on $F_2$ for $E_2$ when $R_2\left(E_0\right)>1$. Consequently $\left\{E_0\right\},\left\{E_1\right\},\left\{E_2\right\}$ form the boundary Morse decomposition required by Corollary 1.

Proof. 

On $F_1$ ($i_2=i_{12}=0$, hence $J_2=0$), the equations for the auxiliary variables $r_2,i_{21},r_{12}$ decouple from any strain-2 source term:

$$r_2^{\prime }=-({\beta }_1{\sigma }_1{J}_1+\mu +{\theta }_2)r_2\le -(\mu +{\theta }_2)r_2,$$

so $r_2\left(t\right)\to 0$ exponentially, for every orbit on $F_1\cap \mathrm{\Gamma }$, regardless of $J_1$. Since $\mathrm{\Gamma }$ is compact, $i_1$ is bounded, so the forcing term ${\beta }_1{\sigma }_1{r}_2{i}_1$ in $i_{21}^{\prime }={\beta }_1{\sigma }_1{r}_2{i}_1-v_1{i}_{21}$ decays exponentially to 0; since $v_1>0$, a standard comparison argument gives $i_{21}\left(t\right)\to 0$. The same argument applied to $r_{12}^{\prime }={\gamma }_1{i}_{21}-(\mu +{\theta }_{12})r_{12}$ then gives $r_{12}\left(t\right)\to 0$.

Consequently the $(s,i_1,r_1)$-subsystem on $F_1$ is asymptotically autonomous with limit equation the classical SIRS-with-waning-immunity system

$${\overline{s}}^{\prime }=\mathrm{\Lambda }-\mu \overline{s}-{\beta }_1\overline{s}{\overline{ı}}_1+{\theta }_1{\overline{r}}_1,{\overline{ı}}_1^{\prime }={\beta }_1\overline{s}{\overline{ı}}_1-v_1{\overline{ı}}_1,{\overline{r}}_1^{\prime }={\gamma }_1{\overline{ı}}_1-(\mu +{\theta }_1){\overline{r}}_1,$$

whose unique equilibrium is the disease-free state when $R_1\left(E_0\right)={\beta }_1{s}_0/v_1\le 1$ and the endemic state $E_1$ of (7) when $R_1\left(E_0\right)>1$, each globally asymptotically stable on the limit system by a standard Lyapunov argument for SIRS models [59]. By Thieme’s convergence theorem for asymptotically autonomous semiflows with a globally stable limit equilibrium [60], every orbit of the full $(s,i_1,i_{21},r_1,r_2,r_{12})$-system on $F_1\cap \mathrm{\Gamma }$ converges to the same equilibrium, extended by $r_2=i_{21}=r_{12}=0$: to $E_0$ if $R_1\left(E_0\right)\le 1$, to $E_1$ if $R_1\left(E_0\right)>1$. The argument on $F_2$ is symmetric under $1\leftrightarrow 2$.    □

Theorem 2

(Persistence for the Feng–Gavish model [Gav-pers]). [][Gav-pers] If

$$R_1\left(E_0\right)>1,R_2\left(E_0\right)>1,R_2^{E_1}>1,R_1^{E_2}>1,$$

then the system is uniformly persistent:

$$\exists \epsilon >0:\underset{t\to \infty }{lim\; inf}\left|x_{{\Sigma }_i}\left(t\right)\right|\ge \epsilon ,i=1,2.$$

Argument:

By Lemma 2, $\left\{E_0\right\},\left\{E_1\right\},\left\{E_2\right\}$ form the boundary Morse decomposition required by Corollary 1. By Lemma 1, the boundary connection graph is acyclic. The four inequalities imply every boundary equilibrium repels, so no interior orbit has its omega-limit set on the boundary. Uniform persistence follows from Corollary 1 of §Section 3 (the model falls under its hypotheses with ${\Sigma }_1=\{i_1,i_{21}\}$, ${\Sigma }_2=\{i_2,i_{12}\}$).

5. Is Existence of Coexistence Equilibria Equivalent to the Invasibility of the Fittest Strain, or to Mutual Invasibility?

5.1. Invasibility of the Less Fit Strain Implies Coexistence When Immunity Is Permanent

We investigate first the permanent-immunity case

$${\theta }_i=0,i=1,2,3$$

(10)

because computations are shorter in this setting and already exhibit the mechanism underlying the general boundary-factorization criterion of §Section 5.3.

The kernel structure of an essentially rank-one Metzler siphon block, used repeatedly below, is recorded once in general form as Lemma 5 of §Section 5.3.

Theorem 3

(Coexistence criterion for the permanent-immunity Feng-Gavish model [Coe]). Consider the permanent-immunity model

$${\theta }_i=0,i=1,2,3,$$

with

$${\beta }_i,{\sigma }_i,{\gamma }_i,\mu >0,i=1,2.$$

Define

$$v_i={\gamma }_i+\mu ,{\mathcal{R}}_i={\displaystyle \frac{{\beta }_i}{v_i}},s_i={\displaystyle \frac{1}{{\mathcal{R}}_i}}.$$

Let

$$R_1\left(y\right)={\mathcal{R}}_1(s+{\eta }_1{\sigma }_1{r}_2),R_2\left(y\right)={\mathcal{R}}_2(s+{\eta }_2{\sigma }_2{r}_1)$$

denote the invasibility numbers.

Then

$$\left\{\begin{array}{cc}{R}_2^{E_1}>1,\hfill & {\mathcal{R}}_1>{\mathcal{R}}_2,\hfill \\ R_1^{E_2}>1,\hfill & {\mathcal{R}}_2>{\mathcal{R}}_1\hfill \end{array}\right.⟹E_{12}\in {\mathcal{R}}_{>0}^n.$$

The converse implication is proved below under the same hypotheses, but the proof of that direction contains a step (comparing the coexistence balance to the one-strain balance) that has not been fully justified in general; see the remark following the proof and Conjecture 1 of §Section 5.3.

Proof. 

Recall that

$$s_j={\displaystyle \frac{1}{{\mathcal{R}}_j}},j=1,2,$$

by the rank-one degeneracy equations applied to the siphons ${\sigma }_j$ – see (7).

Lemma 3 (Ordering of the susceptible coordinate at coexistence)If s is the susceptible coordinate of a coexistence equilibrium, then

$$0<s<\overline{s}=min\{s_1,s_2\}.$$

Proof: At coexistence, the rank-one degeneracy equations yield

$${\beta }_1(s+{\eta }_1{\sigma }_1{r}_2)=v_1⟹r_2\left(s\right)={\displaystyle \frac{v_1-{\beta }_{1}s}{{\beta }_1{\eta }_1{\sigma }_1}},{\beta }_2(s+{\eta }_2{\sigma }_2{r}_1)=v_2⟹r_1\left(s\right)={\displaystyle \frac{v_2-{\beta }_{2}s}{{\beta }_2{\eta }_2{\sigma }_2}}.$$

(11)

Therefore

$$0<s<min\{s_1,s_2\}.$$

□

Using the fixed-point equations, the system reduces to the scalar RUR polynomial

$$P\left(s\right)=a_0+a_{1}s+a_2{s}^2+a_{12}{s}^3.$$

Since the coefficients are rather long, we display them only in the special case

$${\eta }_1={\eta }_2=1,\mathrm{\Lambda }=\mu .$$

$$\begin{array}{cc}\hfill a_0& ={\sigma }_1{\sigma }_2({\gamma }_1+\mu )2({\gamma }_2+\mu )2,\hfill \\ \hfill a_1& =({\gamma }_1+\mu )({\gamma }_2+\mu )({\sigma }_1({\gamma }_1+\mu )({\gamma }_2+\mu )\hfill \\ \hfill & -{\sigma }_2\left({\sigma }_1({\gamma }_2+\mu )({\beta }_1+{\gamma }_1+\mu )+{\beta }_2{\sigma }_1({\gamma }_1+\mu )-({\gamma }_1+\mu )({\gamma }_2+\mu )\right)),\hfill \\ \hfill a_2& =(-{\beta }_2){\sigma }_1({\gamma }_2+\mu )({\gamma }_1+\mu )2\hfill \\ \hfill & +{\beta }_1({\gamma }_2+\mu )({\gamma }_1+\mu )\left(({\sigma }_1-1){\sigma }_2({\gamma }_2+\mu )-\mu {\sigma }_1\right)\hfill \\ \hfill & +{\beta }_2{\sigma }_2({\sigma }_1\left({\gamma }_1\mu ({\beta }_1+2({\gamma }_2+\mu ))+\mu ({\beta }_1+\mu )({\gamma }_2+\mu )+{\gamma }_1^2({\gamma }_2+\mu )\right)\hfill \\ \hfill & -\mu ({\gamma }_1+\mu )({\gamma }_2+\mu )),\hfill \\ \hfill a_{12}& ={\beta }_1{\beta }_2\mu \left({\sigma }_1({\gamma }_1+\mu )+{\sigma }_2\left(-{\sigma }_1({\gamma }_1+{\gamma }_2+\mu )+{\gamma }_2+\mu \right)\right).\hfill \end{array}$$

(12)

Lemma 4 (Boundary factorization of the scalar RUR polynomial)The scalar RUR polynomial satisfies

$$P\left(0\right)>0,$$

and

$$P\left(s_1\right)=-{\displaystyle \frac{{\sigma }_1{v}_1^3{v}_2^2}{{\beta }_1^2}}\left(R_2^{E_1}-1\right).$$

(13)

Similarly,

$$P\left(s_2\right)=-{\displaystyle \frac{{\sigma }_2{v}_2^3{v}_1^2}{{\beta }_2^2}}\left(R_1^{E_2}-1\right).$$

(14)

Hence

$$R_2^{E_1}>1⟺P\left(s_1\right)<0,$$

and

$$R_1^{E_2}>1⟺P\left(s_2\right)<0.$$

Proof: The constant term satisfies

$$P\left(0\right)={\sigma }_1{\sigma }_2{v}_1^2{v}_2^2>0.$$

Direct substitution of

$$s_1={\displaystyle \frac{v_1}{{\beta }_1}}$$

into (12) gives

$$P\left(s_1\right)=-{\displaystyle \frac{{\sigma }_1{v}_1^3{v}_2^2}{{\beta }_1^2}}\left(R_2^{E_1}-1\right).$$

Similarly, substituting

$$s_2={\displaystyle \frac{v_2}{{\beta }_2}}$$

gives

$$P\left(s_2\right)=-{\displaystyle \frac{{\sigma }_2{v}_2^3{v}_1^2}{{\beta }_2^2}}\left(R_1^{E_2}-1\right).$$

Since all prefactors are positive, the sign equivalences follow. □

Assume first that

$${\mathcal{R}}_1>{\mathcal{R}}_2,$$

equivalently

$$s_1<s_2.$$

Invasion implies coexistence.

Assume

$$R_2^{E_1}>1.$$

By (13),

$$P\left(s_1\right)<0.$$

Since

$$P\left(0\right)>0,$$

continuity gives

$$s^{*}\in (0,s_1)$$

such that

$$P\left(s^{*}\right)=0.$$

Because

$$s_1<s_2,$$

we obtain

$$s^{*}<s_1<s_2.$$

Hence the reconstruction formulas (11) give

$$r_2\left(s^{*}\right)>0,r_1\left(s^{*}\right)>0.$$

The corresponding singular strain blocks therefore have positive kernel generators

$$\left({\displaystyle \genfrac{}{}{0pt}{}{s^{*}}{{\sigma }_1{r}_2\left(s^{*}\right)}}\right),\left({\displaystyle \genfrac{}{}{0pt}{}{s^{*}}{{\sigma }_2{r}_1\left(s^{*}\right)}}\right).$$

By Lemma 5, all infected coordinates are positive. The remaining coordinates are positive by their linear balance equations. Therefore

$$E_{12}\in {\mathcal{R}}_{>0}^n.$$

On the converse: does coexistence imply the ordered invasion condition?

Suppose

$$E_{12}\in {\mathcal{R}}_{>0}^n,$$

and let $s^{*}$ denote its susceptible coordinate; by the ordering lemma, $0<s^{*}<s_1<s_2$. The following computation was offered as a proof of $R_2^{E_1}>1$, but it contains an unjustified step, which we now isolate explicitly rather than suppress.

At coexistence,

$$0={\gamma }_1{i}_1^{*}-\mu r_1^{*}-{\beta }_2{\eta }_2{\sigma }_2{r}_1^{*}(i_2^{*}+i_{12}^{*}),$$

while at the one-strain equilibrium $E_1$,

$$0={\gamma }_1{i}_1^{E_1}-\mu r_1^{E_1}.$$

These are balance equations at two different equilibria ($s^{*}\ne s_1$ in general, so $i_1^{*}\ne i_1^{E_1}$ in general); comparing them directly does not, by itself, yield a relation between $r_1^{E_1}$ and $r_1\left(s_1\right)=(s_2-s_1)/\left({\eta }_2{\sigma }_2\right)$ — the missing ingredient is control of how $i_1$ (hence the production term ${\gamma }_1{i}_1$) changes between the two equilibria, which the depletion term alone does not supply. This is precisely the permanent-immunity instance of the general obstruction recorded as Conjecture Section 5.3 of §Section 5.3, rather than an assertion proved here.

The case

$${\mathcal{R}}_2>{\mathcal{R}}_1$$

is symmetric (forward implication only).    □

This is exactly the permanent-immunity instance of Conjecture of §Section 5.3: if $E_{12}\in {\mathcal{R}}_{>0}^n$ and ${\mathcal{R}}_1>{\mathcal{R}}_2$, then conjecturally $R_2^{E_1}>1$ (symmetrically for ${\mathcal{R}}_2>{\mathcal{R}}_1$), making the implication of Theorem 3 an equivalence. Numerically (see GavC.wl, GavScan.wl), we have not found a counterexample among the sampled permanent-immunity parameter sets; the four disjoint boundary regimes of Table 1 are consistent with the conjectured equivalence, but this is not a proof.

5.2. The General Boundary-Factorization criterion

The permanent-immunity computation above is the concrete instance, for the Feng–Gavish model, of a general algebraic mechanism: Theorem 4 (Boundary-factorization criterion for coexistence), Lemma 5 (kernel of an essentially rank-one Metzler block), Conjecture 1 (the converse gap), and Question 1 (structural origin of the scalar reduction) are stated and proved once, for general two-strain models with essentially rank-one siphon invasion matrices and one input species, in §Section 5.3. The candidate hypotheses for Question 1 listed there include sign-regularity of the Cramer matrix that produces $\mathrm{\Delta }\left(s\right)$ in Lemma 6 below, where the Feng–Gavish model supplies the worked example.

5.3. A Boundary-Factorization Criterion for Coexistence in Rank-One Two-Strain Models

This section isolates the algebraic mechanism common to coexistence proofs for two-strain models whose siphon blocks are essentially rank-one: a scalar polynomial obtained by elimination, whose boundary values factor through the invasibility deficits, and whose sign change on a feasible interval reconstructs a positive coexistence equilibrium by the intermediate value theorem.

Lemma 5

(Singularity and kernel of an essentially rank-one Metzler matrix). Let

$$M=w{\pi }^T-Diag\left(v\right),w,\pi ,v\in {\mathcal{R}}_{>0}^m.$$

Then

$$M\sin \mathrm{gular}⟺{\pi }^{T}Diag{\left(v\right)}^{-1}w=1.$$

If this condition holds, then

$$rankM=m-1,kerM=〈Diag{\left(v\right)}^{-1}w〉.$$

Consequently, if

$$M{I}_d=0,I_d>0,$$

then

$$I_d=cDiag{\left(v\right)}^{-1}w,c>0.$$

Proof. 

Put $D=Diag\left(v\right)$. Then

$$M=w{\pi }^T-D=-D(I_d-D^{-1}w{\pi }^T).$$

Since D is invertible, M is singular iff $I_d-D^{-1}w{\pi }^T$ is singular. The rank-one matrix $D^{-1}w{\pi }^T$ has unique nonzero eigenvalue

$${\pi }^T{D}^{-1}w.$$

Hence singularity occurs iff

$${\pi }^T{D}^{-1}w=1.$$

Assume this condition holds. If $Mx=0$, then

$$Dx=w\left({\pi }^{T}x\right),$$

so

$$x=\left({\pi }^{T}x\right)D^{-1}w.$$

Thus every kernel vector is proportional to $D^{-1}w$. Conversely,

$$M\left(D^{-1}w\right)=w{\pi }^T{D}^{-1}w-w=0.$$

Therefore

$$kerM=〈D^{-1}w〉,rankM=m-1.$$

If $M{I}_d=0$ and $I_d>0$, then necessarily

$$I_d=c{D}^{-1}w,c>0.$$

   □

The theorem below is, admittedly, close to an abstract intermediate-value lemma: its hypotheses (RANK ONE)–(POSITIVE RECONSTRUCTION) already contain most of the work needed for the conclusion, and we have not shown how to derive them from more primitive structural (network) conditions. We therefore name it for what it is — a criterion, not a theorem deriving coexistence from siphon/CRN structure alone — and record the derivation of its hypotheses from finer structure as Question Section 5.3 below.

Theorem 4

(Boundary-factorization criterion for coexistence, for two-strain models with essentially rank-one siphon invasion matrices and one input species). Consider a two-strain model with one susceptible/input coordinate s, two minimal siphons ${\Sigma }_1,{\Sigma }_2$, and siphon variables

$$I_k\in {\mathcal{R}}_{\ge 0}^{m_k},k=1,2.$$

Let $E_1,E_2$ be the one-strain boundary equilibria, with susceptible coordinates $s_1,s_2$, and let

$$R_2^{E_1},R_1^{E_2}$$

denote the invasibility numbers.

Assume:

• (A1) (RANK ONE)The strain equations may be written as

  $$I_k^{\prime }=M_k(s,z_k)I_k,M_k(s,z_k)=w_k(s,z_k){\pi }_k^T-Diag\left(v_k\right),$$
  with

  $$w_k(s,z_k)>0,{\pi }_k>0,v_k>0.$$
• (A2) (SCALAR RUR POLYNOMIAL)After imposing the rank-one degeneracy equations

  $${\pi }_k^{T}Diag{\left(v_k\right)}^{-1}{w}_k=1,k=1,2,$$
  all variables outside

  $$\left\{s\right\}\cup {\Sigma }_1\cup {\Sigma }_2$$
  are algebraically eliminated, and the remaining fixed-point equations reduce to a scalar RUR polynomial

  $$P\left(s\right)=0.$$
• (A3) (NORMALIZATION)

  $$P\left(0\right)>0.$$
• (A4) (POSITIVE BOUNDARY FACTORIZATION)Direct evaluation at $s=s_1,s_2$ gives

  $$P\left(s_1\right)=-K_1\left(R_2^{E_1}-1\right),$$

  $$P\left(s_2\right)=-K_2\left(R_1^{E_2}-1\right),$$
  with

  $$K_1,K_2>0.$$
• (A5) (POSITIVE RECONSTRUCTION)Every root

  $$s^{*}\in (0,min\{s_1,s_2\})$$
  of $P\left(s\right)=0$ yields positive values for all eliminated variables and positive amplitudes

  $$c_1^{*},c_2^{*}>0$$
  in the two siphon blocks.

Then

$$\left\{\begin{array}{cc}{s}_1<s_2,R_2^{E_1}>1\hfill & ⟹\exists E_{12}\in {\mathcal{R}}_{>0}^n,\hfill \\ s_2<s_1,R_1^{E_2}>1\hfill & ⟹\exists E_{12}\in {\mathcal{R}}_{>0}^n.\hfill \end{array}\right.$$

Remark 8

(role of elimination theory in coexistence results). This coexistence result requires not only CRNT assumptions, but also certain intriguing elimination conditions, which require further study.

Proof. 

Assume first

$$s_1<s_2.$$

If

$$R_2^{E_1}>1,$$

then by (A4),

$$P\left(s_1\right)<0.$$

Since

$$P\left(0\right)>0,$$

the intermediate value theorem gives

$$s^{*}\in (0,s_1)$$

such that

$$P\left(s^{*}\right)=0.$$

Because $s_1<s_2$,

$$s^{*}<s_1<s_2,$$

hence

$$s^{*}\in (0,min\{s_1,s_2\}).$$

By (A5), all eliminated variables are positive and the amplitudes satisfy

$$c_1^{*},c_2^{*}>0.$$

At $s=s^{*}$, the rank-one degeneracy equations make each $M_k(s^{*},z_k\left(s^{*}\right))$ singular. By Lemma 5,

$$ker{M}_k(s^{*},z_k\left(s^{*}\right))=〈Diag{\left(v_k\right)}^{-1}{w}_k(s^{*},z_k\left(s^{*}\right))〉.$$

Hence

$$I_k^{*}=c_k^{*}Diag{\left(v_k\right)}^{-1}{w}_k(s^{*},z_k\left(s^{*}\right)),c_k^{*}>0.$$

Since

$$v_k>0,w_k(s^{*},z_k\left(s^{*}\right))>0,$$

we obtain

$$I_k^{*}>0,k=1,2.$$

Thus the reconstructed fixed point belongs to ${\mathcal{R}}_{>0}^n$.

The case $s_2<s_1$ is symmetric, using $R_1^{E_2}>1$ and the boundary value $P\left(s_2\right)$.    □

Conjecture 1

(Converse of the boundary-factorization coexistence criterion). Under hypotheses i–v of Theorem 4, if $s_1<s_2$ and there exists $E_{12}\in {\mathcal{R}}_{>0}^n$, then $R_2^{E_1}>1$ (symmetrically for $s_2<s_1$); equivalently, the implication in Theorem 4 is in fact an equivalence. The obstruction to proving this in general is comparing the coexistence value $I_1^{*}$ (at $s=s^{*}$) to the one-strain boundary value $I_1^{E_1}$ (at $s=s_1\ne s^{*}$): the two are balance-equation solutions at different points of s, and no general monotonicity principle controlling this comparison is currently known. This exact gap has been isolated explicitly in worked model instances (e.g. the permanent-immunity two-strain case), where no counterexample has been found numerically.    □

Question 1

(Structural origin of scalar coexistence reduction). Find verifiable network (CRN) conditions on a two-strain model implying:

(i)

the coexistence ideal has elimination dimension one;

(ii)

the eliminated equation may be chosen in the common input variable s;

(iii)

its boundary values factor through the invasion deficits $R_j^{E_i}-1$;

(iv)

the reconstruction map is positive on a specified feasible interval.

Candidate hypotheses to investigate: one common input species; rank-one strain kernels; linear balance equations for all non-siphon variables; a tree or acyclic graph among resident recovery variables; sign-regularity of the resulting Cramer matrix.    □

5.4. Existence of Coexistence Equilibria for the Feng–Gavish Model

We now revisit the full Feng–Gavish model. The structure is the same as in the permanent-immunity case:

• the coexistence equations reduce to a scalar RUR polynomial,
• the boundary evaluations factor through the invasion deficits,
• a sign change on $(0,\overline{s})$ produces a coexistence root.

The only new feature is that positivity of the reconstructed equilibrium requires further work to show that the Cramer determinant $\mathrm{\Delta }\left(s\right)$ (15) of the equations $r_1^{\prime }=0,r_2^{\prime }=0$ does not change sign – see next Lemma.

Lemma 6

(Univariate reduction and sign structure for the Gavish model). Consider the Gavish model (6). Put

$$v_i={\gamma }_i+\mu ,d_i=\mu +{\theta }_i,d_{12}=\mu +{\theta }_{12},$$

and

$$s_i={\displaystyle \frac{v_i}{{\beta }_i}},\overline{s}=min\{s_1,s_2\}.$$

At a coexistence fixed point, the two rank-one degeneracy equations are

$${\beta }_1(s+{\eta }_1{\sigma }_1{r}_2)=v_1,{\beta }_2(s+{\eta }_2{\sigma }_2{r}_1)=v_2.$$

Hence

$$r_2\left(s\right)={\displaystyle \frac{s_1-s}{{\eta }_1{\sigma }_1}},r_1\left(s\right)={\displaystyle \frac{s_2-s}{{\eta }_2{\sigma }_2}}.$$

Consequently every positive coexistence point satisfies

$$0<s<\overline{s}.$$

Moreover,

$$\left({\displaystyle \genfrac{}{}{0pt}{}{i_1}{i_{21}}}\right)=c_1\left({\displaystyle \genfrac{}{}{0pt}{}{s}{{\sigma }_1{r}_2\left(s\right)}}\right),\left({\displaystyle \genfrac{}{}{0pt}{}{i_2}{i_{12}}}\right)=c_2\left({\displaystyle \genfrac{}{}{0pt}{}{s}{{\sigma }_2{r}_1\left(s\right)}}\right).$$

The equations $r_1^{\prime }=0,r_2^{\prime }=0$ reduce to

$${\gamma }_{1}s{c}_1-{\sigma }_2{v}_2{r}_1\left(s\right)c_2=d_1{r}_1\left(s\right),$$

$$-{\sigma }_1{v}_1{r}_2\left(s\right)c_1+{\gamma }_{2}s{c}_2=d_2{r}_2\left(s\right).$$

Let

$$\mathrm{\Delta }\left(s\right)={\gamma }_1{\gamma }_2{s}^2-{\sigma }_1{\sigma }_2{v}_1{v}_2{r}_1\left(s\right)r_2\left(s\right).$$

(15)

Then

$$c_1\left(s\right)={\displaystyle \frac{r_1\left(s\right)\left({\gamma }_{2}s{d}_1+{\sigma }_2{v}_2{d}_2{r}_2\left(s\right)\right)}{\mathrm{\Delta }\left(s\right)}},$$

$$c_2\left(s\right)={\displaystyle \frac{r_2\left(s\right)\left({\gamma }_{1}s{d}_2+{\sigma }_1{v}_1{d}_1{r}_1\left(s\right)\right)}{\mathrm{\Delta }\left(s\right)}}.$$

Finally,

$$r_{12}\left(s\right)={\displaystyle \frac{{\gamma }_1{\sigma }_1{r}_2\left(s\right)c_1\left(s\right)+{\gamma }_2{\sigma }_2{r}_1\left(s\right)c_2\left(s\right)}{d_{12}}}.$$

Substitution into $s^{\prime }=0$ gives a univariate equation

$$F_G\left(s\right)=0.$$

With the sign convention

$$P_G\left(s\right):=-\mathrm{\Delta }\left(s\right)F_G\left(s\right),$$

the coexistence equation is equivalently

$$P_G\left(s\right)=0.$$

Moreover, putting

$$\kappa ={\displaystyle \frac{{\gamma }_1{\gamma }_2{\eta }_1{\eta }_2}{v_1{v}_2}},$$

one has

$$\mathrm{\Delta }\left(s\right)={\displaystyle \frac{v_1{v}_2}{{\eta }_1{\eta }_2}}\left[\kappa s^2-(s_1-s)(s_2-s)\right].$$

There exists a unique number

$$s_{\mathrm{\Delta }}\in (0,\overline{s})$$

such that

$$\mathrm{\Delta }\left(s_{\mathrm{\Delta }}\right)=0,$$

namely

$$s_{\mathrm{\Delta }}={\displaystyle \frac{2s_1{s}_2}{s_1+s_2+\sqrt{{(s_1+s_2)}^2+4(\kappa -1)s_1{s}_2}}}.$$

Furthermore,

$$\mathrm{\Delta }\left(s\right)<0⟺0<s<s_{\mathrm{\Delta }},$$

and

$$\mathrm{\Delta }\left(s\right)>0⟺s_{\mathrm{\Delta }}<s<\overline{s}.$$

Consequently,

$$c_1\left(s\right),c_2\left(s\right)>0⟺s>s_{\mathrm{\Delta }}.$$

Proof. 

The infected blocks have the rank-one forms

$$M_1={\beta }_1\left({\displaystyle \genfrac{}{}{0pt}{}{s}{{\sigma }_1{r}_2}}\right)\left(\begin{array}{cc}1& {\eta }_1\end{array}\right)-v_1{I}_2,$$

$$M_2={\beta }_2\left({\displaystyle \genfrac{}{}{0pt}{}{s}{{\sigma }_2{r}_1}}\right)\left(\begin{array}{cc}1& {\eta }_2\end{array}\right)-v_2{I}_2.$$

At coexistence $M_k{I}_k=0$ with $I_k>0$. By Lemma 5,

$${\beta }_1(s+{\eta }_1{\sigma }_1{r}_2)=v_1,{\beta }_2(s+{\eta }_2{\sigma }_2{r}_1)=v_2.$$

This gives

$$r_2\left(s\right)={\displaystyle \frac{s_1-s}{{\eta }_1{\sigma }_1}},r_1\left(s\right)={\displaystyle \frac{s_2-s}{{\eta }_2{\sigma }_2}}.$$

Positivity of $r_1,r_2$ immediately implies

$$0<s<\overline{s}.$$

The same lemma gives

$$\left({\displaystyle \genfrac{}{}{0pt}{}{i_1}{i_{21}}}\right)=c_1\left({\displaystyle \genfrac{}{}{0pt}{}{s}{{\sigma }_1{r}_2\left(s\right)}}\right),\left({\displaystyle \genfrac{}{}{0pt}{}{i_2}{i_{12}}}\right)=c_2\left({\displaystyle \genfrac{}{}{0pt}{}{s}{{\sigma }_2{r}_1\left(s\right)}}\right).$$

Substitution into $r_1^{\prime }=0,r_2^{\prime }=0$ gives the displayed linear system for $c_1,c_2$. Its determinant is $\mathrm{\Delta }\left(s\right)$, and Cramer’s rule gives the formulas for $c_1\left(s\right),c_2\left(s\right)$.

The formula for $r_{12}\left(s\right)$ follows from $r_{12}^{\prime }=0$, and substitution into $s^{\prime }=0$ gives $F_G\left(s\right)=0$. Multiplication by $-\mathrm{\Delta }\left(s\right)$ gives the polynomial $P_G\left(s\right)$.

Using the formulas for $r_1\left(s\right),r_2\left(s\right)$,

$$\begin{array}{cc}\hfill \mathrm{\Delta }\left(s\right)& ={\gamma }_1{\gamma }_2{s}^2-{\sigma }_1{\sigma }_2{v}_1{v}_2{r}_1\left(s\right)r_2\left(s\right)\hfill \\ & ={\gamma }_1{\gamma }_2{s}^2-{\displaystyle \frac{v_1{v}_2}{{\eta }_1{\eta }_2}}(s_1-s)(s_2-s),\hfill \end{array}$$

which yields the displayed factorization.

Hence $\mathrm{\Delta }\left(s\right)=0$ is equivalent to

$$(\kappa -1)s^2+(s_1+s_2)s-s_1{s}_2=0.$$

The positive solution is

$$s_{\mathrm{\Delta }}={\displaystyle \frac{2s_1{s}_2}{s_1+s_2+\sqrt{{(s_1+s_2)}^2+4(\kappa -1)s_1{s}_2}}}.$$

Since

$$\mathrm{\Delta }\left(0\right)=-{\displaystyle \frac{v_1{v}_2}{{\eta }_1{\eta }_2}}{s}_1{s}_2<0,$$

while

$$\mathrm{\Delta }\left(\overline{s}\right)={\gamma }_1{\gamma }_2{\overline{s}}^2>0,$$

the stated sign pattern follows.

Finally, on $(0,\overline{s})$,

$$r_1\left(s\right)>0,r_2\left(s\right)>0,$$

and the numerators in the formulas for $c_1\left(s\right),c_2\left(s\right)$ are strictly positive. Therefore

$$c_1\left(s\right),c_2\left(s\right)>0⟺\mathrm{\Delta }\left(s\right)>0⟺s>s_{\mathrm{\Delta }}.$$

   □

Lemma 7

(Boundary factorization for the Gavish RUR polynomial). Let

$$s_1={\displaystyle \frac{v_1}{{\beta }_1}},s_2={\displaystyle \frac{v_2}{{\beta }_2}}.$$

The one-strain equilibria satisfy

$$r_1^{E_1}={\displaystyle \frac{{\gamma }_1(\mathrm{\Lambda }-\mu s_1)}{\mu (v_1+{\theta }_1)}},r_2^{E_2}={\displaystyle \frac{{\gamma }_2(\mathrm{\Lambda }-\mu s_2)}{\mu (v_2+{\theta }_2)}}.$$

The invasion numbers are

$$R_2^{E_1}={\displaystyle \frac{{\beta }_2}{v_2}}\left(s_1+{\eta }_2{\sigma }_2{r}_1^{E_1}\right),$$

$$R_1^{E_2}={\displaystyle \frac{{\beta }_1}{v_1}}\left(s_2+{\eta }_1{\sigma }_1{r}_2^{E_2}\right).$$

With the convention $P_G=-\mathrm{\Delta }{F}_G$,

$$P_G\left(0\right)>0,$$

and

$$P_G\left(s_1\right)=-{\displaystyle \frac{{\gamma }_2{s}_1^2\mu (v_1+{\theta }_1)v_2}{{\beta }_2{\eta }_2{\sigma }_2}}\left(R_2^{E_1}-1\right),$$

$$P_G\left(s_2\right)=-{\displaystyle \frac{{\gamma }_1{s}_2^2\mu (v_2+{\theta }_2)v_1}{{\beta }_1{\eta }_1{\sigma }_1}}\left(R_1^{E_2}-1\right).$$

Thus

$$R_2^{E_1}>1⟺P_G\left(s_1\right)<0,R_1^{E_2}>1⟺P_G\left(s_2\right)<0.$$

Proof. 

At $s=s_1$,

$$r_2\left(s_1\right)=0,\mathrm{\Delta }\left(s_1\right)={\gamma }_1{\gamma }_2{s}_1^2,$$

and

$$r_1\left(s_1\right)={\displaystyle \frac{s_2-s_1}{{\eta }_2{\sigma }_2}}.$$

Using $P_G=-\mathrm{\Delta }{F}_G$, the boundary value reduces to

$$P_G\left(s_1\right)=-{\gamma }_1{\gamma }_2{s}_1^2\left(\mathrm{\Lambda }-\mu s_1-{\displaystyle \frac{\mu (v_1+{\theta }_1)}{{\gamma }_1}}{r}_1\left(s_1\right)\right).$$

Since

$$r_1^{E_1}={\displaystyle \frac{{\gamma }_1(\mathrm{\Lambda }-\mu s_1)}{\mu (v_1+{\theta }_1)}},$$

we get

$$P_G\left(s_1\right)=-{\gamma }_2{s}_1^2\mu (v_1+{\theta }_1)\left(r_1^{E_1}-r_1\left(s_1\right)\right).$$

Moreover,

$$R_2^{E_1}-1={\displaystyle \frac{{\beta }_2}{v_2}}{\eta }_2{\sigma }_2\left(r_1^{E_1}-r_1\left(s_1\right)\right).$$

Therefore

$$P_G\left(s_1\right)=-{\displaystyle \frac{{\gamma }_2{s}_1^2\mu (v_1+{\theta }_1)v_2}{{\beta }_2{\eta }_2{\sigma }_2}}\left(R_2^{E_1}-1\right).$$

The calculation at $s=s_2$ is symmetric. The sign of the displayed prefactors is positive, so the stated sign equivalences follow. Finally, $P_G\left(0\right)>0$ follows by direct evaluation with the same sign convention.    □

Lemma 8

(Positive reconstruction). Let

$$s^{*}\in (s_{\mathrm{\Delta }},\overline{s})$$

satisfy

$$P_G\left(s^{*}\right)=0.$$

Then the reconstructed equilibrium belongs to

$${\mathcal{R}}_{>0}^n.$$

Proof. 

Since

$$s^{*}\in (0,\overline{s}),$$

Lemma 6 gives

$$r_1\left(s^{*}\right)>0,r_2\left(s^{*}\right)>0.$$

Moreover, since

$$s^{*}>s_{\mathrm{\Delta }},$$

the sign statement in Lemma 6 gives

$$\mathrm{\Delta }\left(s^{*}\right)>0.$$

The numerators in the Cramer formulas

$$c_1\left(s\right)={\displaystyle \frac{r_1\left(s\right)\left({\gamma }_{2}s{d}_1+{\sigma }_2{v}_2{d}_2{r}_2\left(s\right)\right)}{\mathrm{\Delta }\left(s\right)}}$$

and

$$c_2\left(s\right)={\displaystyle \frac{r_2\left(s\right)\left({\gamma }_{1}s{d}_2+{\sigma }_1{v}_1{d}_1{r}_1\left(s\right)\right)}{\mathrm{\Delta }\left(s\right)}}$$

are strictly positive on $(0,\overline{s})$. Hence

$$c_1\left(s^{*}\right)>0,c_2\left(s^{*}\right)>0.$$

Therefore

$$\left({\displaystyle \genfrac{}{}{0pt}{}{i_1^{*}}{i_{21}^{*}}}\right)=c_1\left(s^{*}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{s^{*}}{{\sigma }_1{r}_2\left(s^{*}\right)}}\right)>0,$$

and

$$\left({\displaystyle \genfrac{}{}{0pt}{}{i_2^{*}}{i_{12}^{*}}}\right)=c_2\left(s^{*}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{s^{*}}{{\sigma }_2{r}_1\left(s^{*}\right)}}\right)>0.$$

Finally,

$$r_{12}\left(s^{*}\right)={\displaystyle \frac{{\gamma }_1{\sigma }_1{r}_2\left(s^{*}\right)c_1\left(s^{*}\right)+{\gamma }_2{\sigma }_2{r}_1\left(s^{*}\right)c_2\left(s^{*}\right)}{d_{12}}}>0.$$

Thus all reconstructed coordinates are positive.    □

Theorem 5

(Ordered invasion and feasible RUR root imply coexistence for the Gavish model [Ord]). Consider the Gavish model (6), with

$${\beta }_i,{\sigma }_i,{\eta }_i,{\gamma }_i,\mu >0,i=1,2.$$

Assume first that

$$s_1<s_2.$$

If

$$R_2^{E_1}>1$$

and $P_G$ has a root

$$s^{*}\in (s_{\mathrm{\Delta }},s_1),$$

then

$$\exists E_{12}\in {\mathcal{R}}_{>0}^n.$$

Symmetrically, if

$$s_2<s_1,R_1^{E_2}>1,$$

and $P_G$ has a root

$$s^{*}\in (s_{\mathrm{\Delta }},s_2),$$

then

$$\exists E_{12}\in {\mathcal{R}}_{>0}^n.$$

Proof. 

Assume $s_1<s_2$. Then $\overline{s}=s_1$. By hypothesis there exists

$$s^{*}\in (s_{\mathrm{\Delta }},s_1)=(s_{\mathrm{\Delta }},\overline{s})$$

such that

$$P_G\left(s^{*}\right)=0.$$

Lemma 8 gives a positive reconstructed equilibrium, hence

$$E_{12}\in {\mathcal{R}}_{>0}^n.$$

The case $s_2<s_1$ is symmetric.    □

Theorem 5 assumes existence of a feasible root of $P_G$ in $(s_{\mathrm{\Delta }},\overline{s})$; it does not address how many such roots there are. Since $P_G$ has degree 4 (one factor of $\mathrm{\Delta }$, itself quadratic, times the degree-2 polynomial $F_G$; the exact degree depends on which terms of $F_G$ survive the substitution), multiplicity of feasible coexistence equilibria — and hence the possibility of saddle-node/backward-type bifurcations in s — is not excluded by Lemma 7 alone: that lemma pins down only the sign of $P_G$ at the two endpoints $0,s_1,s_2$, which by the intermediate value theorem guarantees an odd number of roots (counted with sign changes) in each sign-change subinterval, not a unique one.

Question 2

(Root count for the coexistence polynomial). Determine explicit semialgebraic parameter conditions under which the Gavish coexistence polynomial $P_G$ has

1.

exactly one root in the feasible interval $(s_{\mathrm{\Delta }},\overline{s})$;

2.

more than one root in this interval;

3.

a multiple root in this interval, characterized by

$$P_G\left(s\right)=0,P_G^{\prime }\left(s\right)=0.$$

A sufficient uniqueness criterion would already be obtained by proving

$$P_G^{\prime }\left(s\right)<0,s\in (s_{\mathrm{\Delta }},\overline{s}).$$

More generally, the number of feasible roots may be determined by a Sturm sequence for $P_G$ relative to the algebraic endpoints $s_{\mathrm{\Delta }},\overline{s}$.    □

Remark 9.

A multiple feasible root corresponds algebraically to a fold of the coexistence-equilibrium branch, provided the reconstruction map is regular and the usual parameter transversality condition holds. This interior saddle-node mechanism is distinct from the boundary transcritical relays of Section 6.

Solving Problem 2 is left open; the tools available are the discriminant $Disc\left(P_G\right)$, the resultant $Res(P_G,P_G^{\prime })$, and a Sturm sequence for $P_G$ on $(s_{\mathrm{\Delta }},\overline{s})$, none of which we have computed here. This is a natural next step: it connects invasion ($R_j^{E_i}>1$, Lemma 7) to coexistence existence (Theorem 5) to coexistence multiplicity (Problem 2) to possible saddle-node/backward bifurcations of the coexistence branch itself, a genuinely different bifurcation phenomenon from the boundary transcritical relays of Section 6.

6. Relay Mechanisms and Coexistence Results for Eco-Epidemic Models

This section develops boundary transcritical relays. Such phenomena are repeatedly observed in explicit epidemic and ecological models: a resident equilibrium loses stability precisely when a new equilibrium appears on a larger invariant face.

The key observation is that this coincidence is already visible in competitive exclusion partitions. In the simplest examples, the same equality of reproduction functions simultaneously determines loss of transversal stability, appearance of a successor equilibrium, and change of the relay inequality.

The relay mechanism naturally separates into two parts. First, one has a threshold theorem identifying invasion thresholds with Perron eigenvalue crossings of the transversal Metzler block. Second, one has a local bifurcation theorem showing that, under suitable nondegeneracy conditions, such a crossing generates a unique equilibrium branch on the larger invariant face. Together, these results explain why invasion thresholds, transversal stability thresholds, equilibrium-creation thresholds, and relay thresholds coincide.

Theorem 6

(reproduction-function threshold principle). Let $E\left(\lambda \right)$ be a smooth resident equilibrium branch contained in an invariant face. Let S be a minimal siphon absent at $E\left(\lambda \right)$, and let $J_S^{\perp }\left(E\left(\lambda \right)\right)$ denote the transversal Jacobian block in the variables of S. Assume that $J_S^{\perp }\left(E\left(\lambda \right)\right)$ admits a regular splitting

$$J_S^{\perp }\left(E\left(\lambda \right)\right)=F_S\left(E\left(\lambda \right)\right)-V_S\left(E\left(\lambda \right)\right),$$

where

$$F_S\left(E\left(\lambda \right)\right)\ge 0,V_S\left(E\left(\lambda \right)\right)$$

is a nonsingular M-matrix. Assume further that $J_S^{\perp }\left(E\left(\lambda \right)\right)$ is irreducible. Define

$$K_S\left(E\left(\lambda \right)\right)=F_S\left(E\left(\lambda \right)\right)V_S{\left(E\left(\lambda \right)\right)}^{-1},$$

and

$$R_S\left(E\left(\lambda \right)\right)=\rho \left(K_S\left(E\left(\lambda \right)\right)\right).$$

Thus $R_S\left(E\left(\lambda \right)\right)$ is the reproduction function of the absent siphon S, evaluated at the resident equilibrium $E\left(\lambda \right)$. Let $\alpha \left(\lambda \right)$ denote the Perron eigenvalue of $J_S^{\perp }\left(E\left(\lambda \right)\right)$. Then

$$R_S\left(E\left(\lambda \right)\right)<1⟺\alpha \left(\lambda \right)<0,$$

$$R_S\left(E\left(\lambda \right)\right)=1⟺\alpha \left(\lambda \right)=0,$$

and

$$R_S\left(E\left(\lambda \right)\right)>1⟺\alpha \left(\lambda \right)>0.$$

Consequently

$$R_S\left(E\left({\lambda }_{*}\right)\right)=1$$

is equivalent to loss of transversal hyperbolicity of the resident branch. Moreover, at ${\lambda }_{*}$, the Perron eigenvalue $\alpha \left({\lambda }_{*}\right)=0$ is simple, $J_S^{\perp }\left(E\left({\lambda }_{*}\right)\right)$ has a one-dimensional kernel, and the corresponding right and left Perron vectors $w\gg 0$, $\pi \gg 0$ are unique up to positive scaling.

Proof. 

The next-generation matrix theorem applied to the regular splitting

$$J_S^{\perp }\left(E\left(\lambda \right)\right)=F_S\left(E\left(\lambda \right)\right)-V_S\left(E\left(\lambda \right)\right)$$

gives

$$s\left(J_S^{\perp }\left(E\left(\lambda \right)\right)\right)<0⟺\rho \left(K_S\left(E\left(\lambda \right)\right)\right)<1,$$

$$s\left(J_S^{\perp }\left(E\left(\lambda \right)\right)\right)=0⟺\rho \left(K_S\left(E\left(\lambda \right)\right)\right)=1,$$

and

$$s\left(J_S^{\perp }\left(E\left(\lambda \right)\right)\right)>0⟺\rho \left(K_S\left(E\left(\lambda \right)\right)\right)>1.$$

Since $J_S^{\perp }\left(E\left(\lambda \right)\right)$ is Metzler [11], its spectral abscissa is its Perron eigenvalue $\alpha \left(\lambda \right)$. Substituting

$$R_S\left(E\left(\lambda \right)\right)=\rho \left(K_S\left(E\left(\lambda \right)\right)\right)$$

gives the three stated equivalences.

It remains to establish simplicity of $\alpha \left({\lambda }_{*}\right)=0$. Since $J_S^{\perp }\left(E\left({\lambda }_{*}\right)\right)$ is irreducible and Metzler, the Perron–Frobenius theorem for irreducible Metzler matrices asserts that the spectral abscissa is a simple eigenvalue with strictly positive right and left eigenvectors $w\gg 0$, $\pi \gg 0$, unique up to positive scaling. In particular $ker{J}_S^{\perp }\left(E\left({\lambda }_{*}\right)\right)$ is one-dimensional, spanned by w.    □

Theorem 7

(local boundary transcritical relay and successor existence [loc]). Let

$$E\left(\lambda \right)=\left(y\right(\lambda ),0)$$

be a smooth resident equilibrium branch, and let S be a minimal siphon absent at $E\left(\lambda \right)$. Let z denote the variables of S, so that the full state is $x=(z,y)$. Write

$$\dot{y}=f_y(y,z,\lambda ),\dot{z}=M(y,z,\lambda )z,$$

where λ is a scalar unfolding parameter.

Assume that $f_y,M$ are of class $C^2$, that for some ${\lambda }_{*}$,

$$R_S\left(E\left({\lambda }_{*}\right)\right)=1,$$

and that $M(y_{*},0,{\lambda }_{*})$, with $y_{*}=y\left({\lambda }_{*}\right)$, has a simple Perron eigenvalue 0 with right and left Perron eigenvectors

$$w\gg 0,\pi \gg 0,\mathrm{normalized}\mathrm{so}\mathrm{that}{\pi }^{\top }w=1,$$

and that every other eigenvalue of $M(y_{*},0,{\lambda }_{*})$ has strictly negative real part. (This normalization fixes the otherwise free positive scaling of $w,\pi$ independently; $a_{\sigma }\left(E\right),b_{\sigma }\left(E\right)$ below are the standard simple-eigenvalue-perturbation coefficients ${\pi }^{\top }(·)w$, which equal the true derivatives of the Perron eigenvalue only under ${\pi }^{\top }w=1$ — for unnormalized eigenvectors the formulas acquire a factor $1/\left({\pi }^{\top }w\right)$. The ratio $a_{\sigma }\left(E\right)/b_{\sigma }\left(E\right)$, the sign condition ${\alpha }_{\mathrm{res}}\left(\lambda \right)b_{\sigma }\left(E\right)<0$, and the branch $x=z\left(\lambda \right)w$ itself are unaffected by this choice: rescaling $w\to cw$, $\pi \to \pi /c$ rescales $a_{\sigma }\to a_{\sigma }$, $b_{\sigma }\to c{b}_{\sigma }$, and $z\left(\lambda \right)\to z\left(\lambda \right)/c$, leaving the physical branch $z\left(\lambda \right)w$ and all sign statements invariant.) Put

$$J_{yy}=D_y{f}_y(y_{*},0,{\lambda }_{*}),J_{yz}=D_z{f}_y(y_{*},0,{\lambda }_{*}).$$

Assume:

(tgS)

Tangential stability:

$$\alpha \left(J_{yy}\right)<0.$$

(trC)

Transversal crossing:

$$a_{\sigma }\left(E\right)={\pi }^{\top }\left(D_{y}M\xi +{\partial }_{\lambda }M\right)w\ne 0,\xi =-J_{yy}^{-1}{\partial }_{\lambda }{f}_y.$$

(nzB)

Nonzero branch coefficient:

$$b_{\sigma }\left(E\right)={\pi }^{\top }\left(D_{y}M\eta +D_{z}Mw\right)w\ne 0,\eta =-J_{yy}^{-1}{J}_{yz}w.$$

Then there exists a unique local nonresident equilibrium branch

$$E_S\left(\lambda \right)=(y_S\left(\lambda \right),z_S\left(\lambda \right))$$

meeting the resident branch at ${\lambda }_{*}$, and it has the expansion

$$z_S\left(\lambda \right)=\tau \left(\lambda \right)w+o\left(\tau \left(\lambda \right)\right),y_S\left(\lambda \right)=y_{*}+\tau \left(\lambda \right)\eta +o\left(\tau \left(\lambda \right)\right),$$

where

$$\tau \left(\lambda \right)=-{\displaystyle \frac{a_{\sigma }\left(E\right)}{b_{\sigma }\left(E\right)}}(\lambda -{\lambda }_{*})+o(\lambda -{\lambda }_{*}).$$

If $\alpha \left(\lambda \right)$ denotes the critical transversal eigenvalue along the resident branch, then

$$\alpha \left(\lambda \right)=a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})+o(\lambda -{\lambda }_{*}).$$

Moreover, the successor equilibrium on the adjacent face exists locally in the positive region precisely on the side where

$$\tau \left(\lambda \right)>0,$$

equivalently

$$\alpha \left(\lambda \right)b_{\sigma }\left(E\right)<0.$$

Thus, under the nondegeneracy assumptions

$$a_{\sigma }\left(E\right)\ne 0,b_{\sigma }\left(E\right)\ne 0,$$

the local event is a generic boundary transcritical bifurcation. It is forward iff

$$b_{\sigma }\left(E\right)<0,$$

in which case the successor equilibrium exists on the side where the resident branch is transversally unstable. It is backward iff

$$b_{\sigma }\left(E\right)>0,$$

in which case the successor equilibrium exists on the side where the resident branch is still transversally stable. The cases $a_{\sigma }\left(E\right)=0$ and $b_{\sigma }\left(E\right)=0$ are degenerate and are treated in §Section 6.1.

Moreover, the two branchesexchange stabilityat ${\lambda }_{*}$. If ${\alpha }_S\left(\lambda \right)$ denotes the critical transversal eigenvalue along the successor branch $E_S\left(\lambda \right)$, then

$${\alpha }_S\left(\lambda \right)=-a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})+o(\lambda -{\lambda }_{*}),$$

so that, provided all noncritical eigenvalues remain in the open left half-plane along both branches, exactly one of $\alpha \left(\lambda \right)$, ${\alpha }_S\left(\lambda \right)$ is negative on each side of ${\lambda }_{*}$: the successor branch is transversally stable precisely where the resident branch is unstable, and conversely.

Proof. 

Since $\alpha \left(J_{yy}\right)<0$, the matrix $J_{yy}$ is invertible. The implicit-function theorem applied to

$$f_y(y,0,\lambda )=0$$

gives a unique smooth resident branch near ${\lambda }_{*}$. Differentiating the resident equilibrium equation at ${\lambda }_{*}$ gives

$$J_{yy}{y}^{\prime }\left({\lambda }_{*}\right)+{\partial }_{\lambda }{f}_y=0,$$

hence

$$y^{\prime }\left({\lambda }_{*}\right)=\xi =-J_{yy}^{-1}{\partial }_{\lambda }{f}_y.$$

Along the resident branch, the critical transversal eigenvalue is the simple Perron eigenvalue of $M\left(y\right(\lambda ),0,\lambda )$. The standard simple-eigenvalue perturbation formula gives

$${\alpha }^{\prime }\left({\lambda }_{*}\right)={\pi }^{\top }\left(D_{y}M{y}^{\prime }\left({\lambda }_{*}\right)+{\partial }_{\lambda }M\right)w={\pi }^{\top }\left(D_{y}M\xi +{\partial }_{\lambda }M\right)w=a_{\sigma }\left(E\right).$$

Since $\alpha \left({\lambda }_{*}\right)=0$, it follows that

$$\alpha \left(\lambda \right)=a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})+o(\lambda -{\lambda }_{*}).$$

We now seek nonresident equilibria in the form

$$z=\tau w+o\left(\tau \right),y=y_{*}+\tau \eta +o\left(\tau \right).$$

The first-order terms in the resident equation give

$$J_{yy}\eta +J_{yz}w=0,$$

and therefore

$$\eta =-J_{yy}^{-1}{J}_{yz}w.$$

Substituting

$$z=\tau w+o\left(\tau \right),y=y_{*}+\tau \eta +o\left(\tau \right),$$

into

$$M(y,z,\lambda )z=0$$

and projecting with ${\pi }^{\top }$, using

$$M(y_{*},0,{\lambda }_{*})w=0,{\pi }^{\top }M(y_{*},0,{\lambda }_{*})=0,$$

gives

$$\tau \left(a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})+b_{\sigma }\left(E\right)\tau +o\left(\right|\lambda -{\lambda }_{*}|+|\tau \left|\right)\right)=0.$$

The factor $\tau =0$ is the resident branch. Since

$$a_{\sigma }\left(E\right)\ne 0,b_{\sigma }\left(E\right)\ne 0,$$

the implicit-function theorem applied to the second factor yields a unique nonresident branch satisfying

$$\tau \left(\lambda \right)=-{\displaystyle \frac{a_{\sigma }\left(E\right)}{b_{\sigma }\left(E\right)}}(\lambda -{\lambda }_{*})+o(\lambda -{\lambda }_{*}).$$

Since $w\gg 0$, the condition $z_S\left(\lambda \right)\gg 0$ is equivalent, for sufficiently small $|\lambda -{\lambda }_{*}|$, to $\tau \left(\lambda \right)>0$. Combining the expansion of $\tau \left(\lambda \right)$ with the expansion of $\alpha \left(\lambda \right)$ gives

$$\tau \left(\lambda \right)>0⟺-{\displaystyle \frac{\alpha \left(\lambda \right)}{b_{\sigma }\left(E\right)}}>0⟺\alpha \left(\lambda \right)b_{\sigma }\left(E\right)<0.$$

This proves the asserted local existence condition for the successor equilibrium and the forward/backward relay alternative.

For the stability-exchange claim, the reduced dynamics on the one-dimensional center manifold parametrized by $\tau$ is, to leading order, the transcritical normal form

$$\dot{\tau }\approx G(\tau ,\lambda )\tau =a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})\tau +b_{\sigma }\left(E\right){\tau }^2,$$

since $\alpha \left(\lambda \right)={\partial }_{\tau }{\left[G(\tau ,\lambda )\tau \right]|}_{\tau =0}=a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})$ recovers the resident eigenvalue expansion already established. The critical eigenvalue along the successor branch is the derivative of the same normal form evaluated at $\tau =\tau \left(\lambda \right)\ne 0$:

$${\alpha }_S\left(\lambda \right)={\partial }_{\tau }\left[a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})\tau +b_{\sigma }\left(E\right){\tau }^2\right]{|}_{\tau =\tau \left(\lambda \right)}=a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})+2b_{\sigma }\left(E\right)\tau \left(\lambda \right)=-a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})+o(\lambda -{\lambda }_{*}).$$

Since the critical eigenvalue 0 is simple and every other eigenvalue of $M(y_{*},0,{\lambda }_{*})$ has negative real part, continuity of the spectrum keeps all noncritical eigenvalues in the open left half-plane for $\lambda$ near ${\lambda }_{*}$ along both branches; hence $\mathrm{sign}\left(\alpha \right)$ and $\mathrm{sign}\left({\alpha }_S\right)=-\mathrm{sign}\left(\alpha \right)$ determine the full transversal stability of each branch, proving the stated exchange of stability.    □

Remark 10

(Scope of the stability-exchange conclusion). The proof above computes only the critical eigenvalue: it shows that the two branches exchange stabilityin the critical (z-)direction. Full local asymptotic stability of a branch requires in addition that the full Jacobian at that branch point be Hurwitz, which needs:

(a)

the resident tangential block $J_{yy}$ is Hurwitz (this is hypothesis (tgS));

(b)

all noncritical eigenvalues of $M(y_{*},0,{\lambda }_{*})$ remain in the open left half-plane for λ near ${\lambda }_{*}$ along both branches (established above by continuity of the spectrum, since 0 is simple);

(c)

the tangential Jacobian $J_{yy}$ evaluated along thesuccessorbranch, at $(y_S\left(\lambda \right),z_S\left(\lambda \right))$, also remains Hurwitz for λ near ${\lambda }_{*}$.

Condition (c) is not automatic from (a) alone, since the full Jacobian at a successor-branch point is not block-triangular in $(y,z)$ once $z_S\left(\lambda \right)\ne 0$: the block ${\partial }_y{\left(Mz\right)|}_{z=z_S\left(\lambda \right)}$ is generally nonzero off the resident branch. It does follow by continuity: at $\lambda ={\lambda }_{*}$ the branches coincide ($z_S\left({\lambda }_{*}\right)=0$), where the full Jacobian is exactly block lower-triangular with blocks $J_{yy}(y_{*},{\lambda }_{*})$ and $M(y_{*},0,{\lambda }_{*})$, whose spectrum lies in the closed left half-plane with a single eigenvalue at 0. Since $z_S\left(\lambda \right)\to 0$ as $\lambda \to {\lambda }_{*}$, the full Jacobian along the successor branch converges to this same block-triangular limit, so by continuity of eigenvalues its strictly negative eigenvalues (in particular those coming from $J_{yy}$) remain strictly negative for λ sufficiently close to ${\lambda }_{*}$, giving (c). Thus, precisely stated: the critical eigenvalue changes sign with opposite orientation on the two branches; hence, under (a)–(c) — the persistence of the stable complementary spectrum on both branches, which holds locally near ${\lambda }_{*}$ by the continuity argument above — the branches exchange full local asymptotic stability, not merely transversal stability.

6.1. Forward, Backward, and Degenerate Boundary Bifurcations

The invasion equality

$$R_S\left(E\left({\lambda }_{*}\right)\right)=1$$

identifies loss of transversal hyperbolicity, but does not by itself determine either the existence or the direction of a successor equilibrium branch. Under the hypotheses of Theorem 7, the reduced equilibrium equation is

$$\tau \left(a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})+b_{\sigma }\left(E\right)\tau +o\left(|\lambda -{\lambda }_{*}|+|\tau |\right)\right)=0.$$

The factor $\tau =0$ represents the resident branch, while the nonresident branch satisfies

$$\tau \left(\lambda \right)=-{\displaystyle \frac{a_{\sigma }\left(E\right)}{b_{\sigma }\left(E\right)}}(\lambda -{\lambda }_{*})+o(\lambda -{\lambda }_{*}).$$

Since

$$\alpha \left(\lambda \right)=a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})+o(\lambda -{\lambda }_{*}),$$

the positivity condition for the successor branch is

$$\tau \left(\lambda \right)>0⟺\alpha \left(\lambda \right)b_{\sigma }\left(E\right)<0.$$

Theorem 8

(Classification at an invasion threshold [Cla]). At a threshold $R_S\left(E\left({\lambda }_{*}\right)\right)=1$, under the standing hypotheses of Theorem 6 (regular splitting, irreducibility), the local behavior is classified by which of the three conditions (tgS), (trC), (nzB) of Theorem 7 hold:

Condition	Local classification
$a_{\sigma }\left(E\right)\ne 0,b_{\sigma }\left(E\right)<0$	Forward boundary transcritical relay: successor branch exists exactly where the resident becomes unstable.
$a_{\sigma }\left(E\right)\ne 0,b_{\sigma }\left(E\right)>0$	Backward boundary transcritical bifurcation: successor branch exists while the resident is still stable (bistability).
$a_{\sigma }\left(E\right)=0$	Tangential threshold contact / higher-order crossing: $\alpha$ is tangent to 0 at ${\lambda }_{*}$; classification needs the first nonzero derivative $a_q$, $q\ge 2$, of $\alpha$.
$a_{\sigma }\left(E\right)\ne 0,b_{\sigma }\left(E\right)=0$	Degenerate boundary transcritical bifurcation: with first nonzero higher-order coefficient $c_q$ ($q\ge 2$), ${\tau }^q\propto -a_{\sigma }\left(E\right)(\lambda -{\lambda }_{*})/c_q$; for $q=2$ this is a one-sided square-root branch (no $\tau \mapsto -\tau$ symmetry survives positivity, so it is not a true pitchfork).
$0\in \sigma \left(J_{yy}\right)$ or complex pair	(tgS) fails: tangential–transversal codimension-two interaction (fold–transcritical if $J_{yy}$ has a zero eigenvalue, or transcritical–Hopf if a purely imaginary pair coincides with the transversal crossing); the scalar reduction of Theorem 7 no longer suffices.
$dimker{J}_S^{\perp }\left(E\left({\lambda }_{*}\right)\right)>1$	Simplicity fails: several invasion directions become critical simultaneously, giving a multiple relay pointwith several successor branches rather than one.

When (tgS), (trC), (nzB) all hold, the relay is generic and forward or backward exactly according to the sign of $b_{\sigma }\left(E\right)$; no single degenerate equality by itself proves backward bifurcation.

Remark 11

(Meaning of the relay threshold principle). Theorem 6 and Theorem 7 play different roles. Theorem 6 is a threshold theorem:

$$R_S\left(E\left({\lambda }_{*}\right)\right)=1⟺\alpha \left({\lambda }_{*}\right)=0.$$

Theorem 7 is the local existence theorem for the successor equilibrium. Under(tgS),(trC)and(nzB), the threshold produces a unique local branch $E_S\left(\lambda \right)$, and the condition

$$\alpha \left(\lambda \right)b_{\sigma }\left(E\right)<0$$

is exactly the local positivity condition for this successor branch.

Thus an invasion threshold is the spectral signature of a relay, while Theorem 7 gives the additional bifurcation mechanism which turns the spectral threshold into an adjacent equilibrium.

The equality

$$R_S\left(E\left({\lambda }_{*}\right)\right)=1$$

therefore detects only the spectral threshold. A genuine generic relay additionally requires

$$a_{\sigma }\left(E\right)\ne 0,b_{\sigma }\left(E\right)\ne 0.$$

The strict sign $b_{\sigma }\left(E\right)>0$ yields a backward boundary transcritical bifurcation, whereas the equalities $a_{\sigma }\left(E\right)=0$ or $b_{\sigma }\left(E\right)=0$ indicate higher-order degeneracies rather than backward bifurcation.

Remark 12

(relation with competitive exclusion partitions). In a two-resident competitive exclusion partition, the relay thresholds are written as invasion equalities

$$R_q\left(E_p\right)=1,R_p\left(E_q\right)=1,$$

where $E_p$ is the resident equilibrium with block p present and $R_q\left(E_p\right)$ is the reproduction function of the absent block q evaluated at $E_p$.

The general threshold chain is

$$R_q\left(E_p\right)=1⟺s\left(J_{q|p}^{\perp }\right)=0⟺\mathrm{loss}\mathrm{of}\mathrm{transversal}\mathrm{hyperbolicity}\mathrm{of}{E}_p.$$

Under the nondegeneracy assumptions of Theorem 7, this is also equivalent locally to the birth or death of the adjacent successor equilibrium.

In the rank-one case one further obtains the explicit reduction

$$R_q\left(E_p\right)=1⟺R_p\left(E_0\right)=R_q\left(E_0\right).$$

Thus the equality of reproduction functions at the DFE is the explicit rank-one form of the more general resident-face invasion equality.

§Section 5.3 also develops a boundary-factorization coexistence criterion for a class of two-strain models with essentially rank-one siphon matrices and a single susceptible variable. Unlike the previous results, which are primarily geometric or dynamical, this criterion requires an additional algebraic ingredient. After elimination, the equilibrium problem reduces to a scalar RUR polynomial satisfying several special properties. In this setting, invasibility of the less fit strain by the fitter resident is sufficient for coexistence (Theorem 4 of §Section 5.3); the converse (that coexistence forces invasibility) remains conjectural, even in the permanent-immunity special case (Conjecture 1).

6.2. Computational Realization by Conic Optimization

The relay conditions of this section are sign conditions on the Perron root of the siphon block: S invades at E when $s\left(M(y_{*},0,\lambda )\right)>0$ and cannot when $s\le 0$, with the threshold $R_S=1$ at $s=0$. Computing a Perron root is an eigenvalue problem, but certifying its sign is not: for a Metzler matrix, $s\left(M\right)\le 0$ is equivalent to the linear feasibility problem

$$\exists w\gg 0:Mw\le 0,$$

by the Collatz–Wielandt characterization. The difference matters because a feasibility problem returns a witness, and one witness can be made to serve a whole parameter region — which is exactly what the chamber structure of a relay graph requires.

Theorem 9

(Uniform noninvasion over a monotone parameter box [UniBox])). [][UniBox] Let $\mathrm{\Lambda }={\prod }_k[{\lambda }_k^{-},{\lambda }_k^{+}]$ be a box, and let $M\left(\lambda \right)$ be Metzler and irreducible for every $\lambda \in \mathrm{\Lambda }$, with every entry nondecreasing in every ${\lambda }_k$. Then the following are equivalent:

(i)

$s\left(M\right(\lambda \left)\right)\le 0$ for every $\lambda \in \mathrm{\Lambda }$;

(ii)

$s\left(M\left({\lambda }^{+}\right)\right)\le 0$ at the single corner ${\lambda }^{+}$;

(iii)

there is one $w\gg 0$ with $M\left(\lambda \right)w\le 0$ for every $\lambda \in \mathrm{\Lambda }$.

Such a w is produced by a single linear feasibility problem posed at ${\lambda }^{+}$, and it certifies noninvasion of S at every resident equilibrium whose parameters lie in Λ. The certificate is exact, not conservative: it exists precisely on the true noninvasion region, not on a proper subset of it.

Proof. (i)⇒(ii) is immediate. For (ii)⇒(iii), take w to be the right Perron vector of $M\left({\lambda }^{+}\right)$, which is strictly positive by irreducibility, so that $M\left({\lambda }^{+}\right)w=s\left(M\left({\lambda }^{+}\right)\right)w\le 0$. For any $\lambda \in \mathrm{\Lambda }$ monotonicity gives $M\left(\lambda \right)\le M\left({\lambda }^{+}\right)$ entrywise, hence $M\left(\lambda \right)w\le M\left({\lambda }^{+}\right)w\le 0$ because $w\gg 0$. For (iii)⇒(i), $M\left(\lambda \right)w\le 0$ with $w\gg 0$ bounds the Collatz–Wielandt ratio: $s\left(M\left(\lambda \right)\right)={min}_{v\gg 0}{max}_i{\left(M\left(\lambda \right)v\right)}_i/v_i\le {max}_i{\left(M\left(\lambda \right)w\right)}_i/w_i\le 0$. The final claim is the equivalence itself: (iii) holds exactly when (i) does.    □

Monotonicity is the biologically ordinary case — transmission, progression and contact coefficients enter the infected block with nonnegative signs, so raising them can only help the siphon invade — and it is what makes the corner test valid. Without it (iii) still implies (i), so a certificate found anywhere is still sound; only the converse, the guarantee that one corner LP finds a certificate whenever one exists, uses the ordering.

A two-stage instance. Let $S=\{E,I\}$ be an exposed/infectious siphon block over a resident susceptible density y,

$$M\left(y\right)=\left(\begin{array}{cc}-(\sigma +\mu )& \beta y\\ \sigma & -(\gamma +\mu )\end{array}\right),\sigma =\frac{4}{5},\mu =\frac{1}{10},\gamma =\frac{1}{2},\beta =\frac{7}{5}.$$

M is Metzler and irreducible, and every entry is nondecreasing in y, so Theorem  applies on $\mathrm{\Lambda }=[0,\overline{y}]$. Here

$$R_S\left(y\right)={\displaystyle \frac{\beta \sigma y}{(\sigma +\mu )(\gamma +\mu )}},{\overline{y}}^{*}={\displaystyle \frac{(\sigma +\mu )(\gamma +\mu )}{\beta \sigma }}={\displaystyle \frac{27}{56}},$$

and the corner LP returns, at $\overline{y}={\overline{y}}^{*}$, the exact certificate

$$w=\frac{1}{7}{(3,4)}^{\top },M\left({\overline{y}}^{*}\right)w=0,$$

the Perron vector, with $w_1/w_2=(\gamma +\mu )/\sigma =3/4$. For $\overline{y}<27/56$ the same w gives $M\left(y\right)w\le 0$ throughout $[0,\overline{y}]$; for $\overline{y}>27/56$ the LP is infeasible, and its infeasibility is the statement that no Volterra-type linear certificate can exist there, consistent with $R_S>1$. One vector of two rational numbers therefore certifies noninvasion over an entire interval of resident densities, in place of an eigenvalue computation at each one.

Code and check. The feasibility problem is

    from scipy.optimize import linprog

    linprog(c=zeros(n), A_ub=M, b_ub=zeros(n),

            A_eq=ones((1,n)), b_eq=[1], bounds=[(1e-6,None)]*n)

with M evaluated at the corner ${\lambda }^{+}$; success returns w, infeasibility certifies $s>0$ somewhere in the box. The equivalence of Theorem  was checked on 4000 random monotone Metzler families of sizes 2 to 4, comparing (i) sampled on a grid, (ii) at the corner and (iii) by the LP, with no disagreement, and in every feasible case the returned w was re-verified to satisfy $M\left(\lambda \right)w\le 0$ across the box.

What the conic formulation contributes here is the passage from a pointwise spectral condition to a single certificate valid on a region. The relay graph is constant on each chamber of parameter space; the LP supplies, per chamber, one strictly positive vector that proves it, and proves it in a form checkable by hand. The same device appears in the chemostat setting in [61], where the resource box plays the role of $\mathrm{\Lambda }$.

6.3. Open Problem: Boundary Relay Geometry and Multistationarity

Beyond its local content, the relay mechanism suggests a global strategy for multistationarity. Boundary equilibria solve the smaller, generically simpler equilibrium equations attached to each siphon face, while interior equilibria require solving the full system; the collection of relay hypersurfaces $R_{\Sigma }\left(E\right)=1$ attached to all resident siphon equilibria partitions parameter space by invasion structure, raising the question of how much this boundary geometry constrains the number and location of strictly positive equilibria. Table 2 summarizes what the relay and rank-one coexistence theorems already establish versus what remains open.

Conjecture 2

(Relay-generation conjecture for positive equilibria). Consider a smooth chemical ODE depending on parameters $\lambda$, with a finite lattice of forward-invariant siphon faces, all equilibria isolated except at ordinary local bifurcation values, and every codimension-one boundary loss of hyperbolicity generated by a simple transversal Metzler eigenvalue. Then every connected component of the set of strictly positive equilibria in state–parameter space has a closure intersecting a boundary equilibrium branch through a finite chain

$$E_{{\Sigma }_0}⟶E_{{\Sigma }_1}⟶\dots ⟶E_{{\Sigma }_q},{\Sigma }_0⊋{\Sigma }_1⊋\dots ⊋{\Sigma }_q=⌀,$$

with $E_{{\Sigma }_j}\in relint{F}_{{\Sigma }_j}$ and each arrow generated at a parameter value with $R_{{\Sigma }_j\setminus {\Sigma }_{j+1}}\left(E_{{\Sigma }_j}\right)=1$ together with the local nondegeneracy conditions of Theorem 7. Equivalently, in the boundary relay graph ${\mathcal{G}}_{\mathrm{relay}}$ (vertices the boundary equilibria $E_{\Sigma }$, edge $E_{\Sigma }\to E_{{\Sigma }^{\prime }}$ whenever ${\Sigma }^{\prime }⊊\Sigma$ and the relay hypotheses hold at $R_{\Sigma \setminus {\Sigma }^{\prime }}\left(E_{\Sigma }\right)=1$), every connected component of positive equilibria has closure intersecting a directed path in ${\mathcal{G}}_{\mathrm{relay}}$ terminating at the DFE vertex ⌀. (The weaker connectivity form of this conjecture only asks that the closure intersect ${\mathcal{G}}_{\mathrm{relay}}$ somewhere, without requiring one specific terminating chain; a counterexample to either form is an isola of strictly positive equilibria disconnected from every vertex and edge of ${\mathcal{G}}_{\mathrm{relay}}$.)    □

Remark 13

(Scope). The conjecture does not assert that every positive equilibrium is created directly from the boundary, nor that boundary relays alone determine its global continuation: after entering the interior, an equilibrium branch may undergo saddle-node, Hopf, or other interior bifurcations. The claim is only that each connected positive branch has a boundary origin through a finite sequence of siphon-face invasions — for example, whether every positive equilibrium connects to the boundary this way, whether disconnected relay chains force disconnected positive branches, whether multiple positive equilibria require several relay hypersurfaces to intersect, and whether backward relays are necessary for some forms of multistationarity, all remain open.

6.4. Relays and LAS–CEP Table for the Feng–Gavish Model

Theorem 10

(LAS-CEP for the Feng–Gavish model [LAS]). Let $R_1:=R_1\left(E_0\right)$, $R_2:=R_2\left(E_0\right)$. Away from equality thresholds, the boundary LAS partition is:

• if $R_1<1$ and $R_2<1$, then DFE is LAS;
• if $R_1>1$ and $R_2<1$, then $E_1$ is LAS;
• if $R_2>1$ and $R_1<1$, then $E_2$ is LAS;
• if $R_1>R_2>1$ and $R_2^{E_1}<1$, then $E_1$ is LAS;
• if $R_2>R_1>1$ and $R_1^{E_2}<1$, then $E_2$ is LAS;
• if $R_2^{E_1}>1$ and $R_1^{E_2}>1$, all boundary equilibria are transversally unstable, the system is uniformly persistent, and $E_{12}\in {\mathbb{R}}_{>0}^n$ exists.

Proof. 

At $E_0$ the Jacobian is block lower-triangular with transversal Metzler blocks $M_1(s_0,0)$, $M_2(s_0,0)$. With the usual regular splitting $M_k=F_k-V_k$, Perron–Frobenius gives

$$sign\alpha \left(M_k\left(E_0\right)\right)=sign(R_k\left(E_0\right)-1).$$

Hence DFE is LAS iff $R_1<1$ and $R_2<1$.

If $R_1>1$ and $R_2<1$, then $E_1$ exists uniquely on $F_1$ and the transversal block at $E_0$ for ${\Sigma }_1$ is unstable, while ${\Sigma }_2$ is stable; $E_2$ does not exist. At $E_1$, the missing block is ${\Sigma }_2$, with transversal matrix $M_2(s_1,r_1^{E_1})$. Again

$$sign\alpha \left(M_2\left(E_1\right)\right)=sign(R_2^{E_1}-1).$$

The tangential block at $E_1$ is Hurwitz whenever $E_1$ exists, by the characteristic-polynomial factorization of [5,47]. Hence $E_1$ is LAS iff $R_1>1$ and $R_2^{E_1}<1$. The same argument, mutatis mutandis with indices $1\leftrightarrow 2$, gives: $E_2$ is LAS iff $R_2>1$ and $R_1^{E_2}<1$. The cases $R_1>R_2>1$ and $R_2>R_1>1$ with the respective invasibility number below 1 fall into these two cases. If both $R_2^{E_1}>1$ and $R_1^{E_2}>1$, both one-strain equilibria are transversally unstable and no boundary equilibrium is LAS; persistence and existence of $E_{12}$ then follow from Corollary 1 and Theorem 5.    □

$$\begin{array}{|c|}\hline \mathrm{LAS}-\mathrm{CEP}\mathrm{is}\mathrm{a}\mathrm{boundary}\mathrm{local}-\mathrm{stability}\mathrm{partition},\mathrm{not}\mathrm{a}\mathrm{global}\mathrm{attractor}\mathrm{partition}.\\ \hline\end{array}$$

Theorem 10 classifies which boundary equilibrium (if any) is locally asymptotically stable, and identifies the regime in which persistence and $E_{12}$ existence follow; it is not a complete competitive-exclusion partition of every attractor, since in the last regime the interior equilibrium $E_{12}$ may itself be unstable, there may be several interior equilibria (cf. the open multiplicity question of Question Section 5.4), or the attractor may be periodic rather than an equilibrium. We reserve the unqualified term “complete exclusion partition” for a future theorem establishing global convergence, and use “LAS-CEP” throughout to mean exactly the boundary local-stability statement proved above.

Remark 14

(Local relay calculation at $E_1$). For the boundary transition $E_1\to E_{12}$, the invading block is $x_{{\Sigma }_2}={(i_2,i_{12})}^{\top }$ and

$${\dot{x}}_{{\Sigma }_2}=M_2(s,r_1)x_{{\Sigma }_2},M_2(s,r_1)={\beta }_2\left({\displaystyle \genfrac{}{}{0pt}{}{s}{{\sigma }_2{r}_1}}\right)\left(1{\eta }_2\right)-v_2{I}_2,$$

so ${\alpha }_2(s,r_1)={\beta }_2(s+{\eta }_2{\sigma }_2{r}_1)-v_2$ and $R_2^{E_1}=({\beta }_2/v_2)(s_1+{\eta }_2{\sigma }_2{r}_1^{E_1})$. At criticality $R_2^{E_1}=1$, put $q_2=s_1+{\eta }_2{\sigma }_2{r}_1^{E_1}=v_2/{\beta }_2$. A normalised Perron pair for $M_2(y_{*},0,{\lambda }_{*})$ is

$$w=\left({\displaystyle \genfrac{}{}{0pt}{}{s_1}{{\sigma }_2{r}_1^{E_1}}}\right)\gg 0,\pi ={\displaystyle \frac{1}{q_2}}\left({\displaystyle \genfrac{}{}{0pt}{}{1}{{\eta }_2}}\right)\gg 0,{\pi }^{\top }w=1.$$

Take $\lambda ={\beta }_2$ as the unfolding parameter. Since ${\partial }_{{\beta }_2}{M}_2=\left({\displaystyle \genfrac{}{}{0pt}{}{s}{{\sigma }_2{r}_1}}\right)\left(1{\eta }_2\right)$ and $E_1$ does not depend on ${\beta }_2$, we have $\xi =0$ and

$$a={\pi }^{\top }\left({\partial }_{{\beta }_2}{M}_2\right)w=(s_1+{\eta }_2{\sigma }_2{r}_1^{E_1}){\pi }^{\top }w=q_2>0.$$

Since $M_2$ depends only on the resident variables $(s,r_1)$, $D_x{M}_2=0$, so $b={\pi }^{\top }\left(D_y{M}_2\eta \right)w$ with $\eta =-J_{yy}^{-1}{J}_{yx}w$. The only resident equations containing $x_{{\Sigma }_2}={(i_2,i_{12})}^{\top }$ are $\dot{s}$ and ${\dot{r}}_1$, giving

$$J_{yx}w=-{\displaystyle \frac{{\beta }_2}{q_2}}\left(\begin{array}{c}{s}_1(s_1+{\eta }_2{\sigma }_2{r}_1^{E_1})\\ {\sigma }_2{r}_1^{E_1}(s_1+{\eta }_2{\sigma }_2{r}_1^{E_1})\\ 0\dots \end{array}\right)=-{\beta }_{2}p,p={(s_1,{\sigma }_2{r}_1^{E_1},0,\dots )}^{\top }.$$

Hence $\eta ={\beta }_2{J}_{yy}^{-1}p$. Since $D_y{M}_{2}v={\beta }_2\left({\displaystyle \genfrac{}{}{0pt}{}{v_s}{{\sigma }_2{v}_{r_1}}}\right)\left(1{\eta }_2\right)$,

$$b={\beta }_2^2({\widehat{p}}_s+{\eta }_2{\sigma }_2{\widehat{p}}_{r_1}),\widehat{p}=J_{yy}^{-1}p.$$

Sign of $\widehat{p}$.The resident block $y=(s,i_1,i_{21},r_1,r_{12})$ at $E_1$ isnotMetzler in this ordering (e.g. ${\partial }_s{\dot{ı}}_1={\beta }_1{i}_1^{E_1}>0$ but ${\partial }_{i_1}\dot{s}=-{\beta }_1{s}_1<0$, an antagonistic off-diagonal pair that no diagonal $\pm 1$ signature can repair, since such a similarity multiplies $J_{ij}$ and $J_{ji}$ by thesamefactor $d_i{d}_j\in \{\pm 1\}$ and so cannot change their relative sign). The Hurwitz property of $J_{yy}$ therefore doesnotby itself give $J_{yy}^{-1}\le 0$, and $\widehat{p}=J_{yy}^{-1}p$ must be computed directly.

Since $r_2^{E_1}=0$ and $i_{21}$ only self-decays at this linearization (${\partial }_s{\dot{ı}}_{21}={\partial }_{i_1}{\dot{ı}}_{21}={\beta }_1{\sigma }_1{r}_2^{E_1}=0$), row $i_{21}$ of $J_{yy}$ is $(0,0,-v_1,0,0)$; since $p_{i_{21}}=0$, this forces ${\widehat{p}}_{i_{21}}=0$, and then row $r_{12}=(0,0,{\gamma }_1,0,-(\mu +{\theta }_{12}))$ with $p_{r_{12}}=0$ forces ${\widehat{p}}_{r_{12}}=0$ as well: both decouple. The remaining $(s,i_1,r_1)$-block of $\widehat{p}=J_{yy}^{-1}p$ solves

$$\left(\begin{array}{ccc}-(\mu +{\beta }_1{i}_1^{E_1})& -{\beta }_1{s}_1& {\theta }_1\\ {\beta }_1{i}_1^{E_1}& {\beta }_1{s}_1-v_1& 0\\ 0& {\gamma }_1& -(\mu +{\theta }_1)\end{array}\right)\left(\begin{array}{c}{\widehat{p}}_s\\ {\widehat{p}}_{i_1}\\ {\widehat{p}}_{r_1}\end{array}\right)=\left(\begin{array}{c}{s}_1\\ 0\\ {\sigma }_2{r}_1^{E_1}\end{array}\right),$$

using ${\beta }_1{s}_1=v_1$ (definition of $s_1$), so the $(i_1,i_1)$ entry is exactly 0. Row $i_1$ then reads ${\beta }_1{i}_1^{E_1}{\widehat{p}}_s=0$; since $i_1^{E_1}>0$ (because $E_1$ exists), this forces

$${\widehat{p}}_s=0$$

exactly. Rows s and $r_1$ then reduce, writing $d_1:=\mu +{\theta }_1$, to the $2\times 2$ system

$$-v_1{\widehat{p}}_{i_1}+{\theta }_1{\widehat{p}}_{r_1}=s_1,{\gamma }_1{\widehat{p}}_{i_1}-d_1{\widehat{p}}_{r_1}={\sigma }_2{r}_1^{E_1},$$

whose solution is

$${\widehat{p}}_{r_1}={\displaystyle \frac{{\sigma }_2{r}_1^{E_1}+{\gamma }_1{s}_1/v_1}{{\gamma }_1{\theta }_1/v_1-d_1}}.$$

The numerator is a sum of strictly positive terms. The denominator satisfies, using $v_1={\gamma }_1+\mu$,

$${\displaystyle \frac{{\gamma }_1{\theta }_1}{v_1}}-d_1={\displaystyle \frac{{\gamma }_1{\theta }_1-(\mu +{\theta }_1)({\gamma }_1+\mu )}{v_1}}={\displaystyle \frac{-\mu ({\gamma }_1+\mu +{\theta }_1)}{v_1}}<0$$

strictly. Hence ${\widehat{p}}_{r_1}<0$, and since ${\widehat{p}}_s=0$,

$$b={\beta }_2^2\left({\widehat{p}}_s+{\eta }_2{\sigma }_2{\widehat{p}}_{r_1}\right)={\beta }_2^2{\eta }_2{\sigma }_2{\widehat{p}}_{r_1}<0.$$

(Equivalently: $\widehat{p}=\eta /{\beta }_2$ is exactly the first-order shift of the resident branch as the invading coordinate grows, so $b={\beta }_{2}d{q}_2{/dz|}_0$ with $q_2=s+{\eta }_2{\sigma }_2{r}_1$ — the direct computation above is a rigorous instance of the informal principle that invasion by strain 2 depletes the effective resource $q_2$ available to it.) Since $a>0$ and $b<0$, the branch satisfies $z>0$ when ${\beta }_2>{\beta }_2^{*}$, i.e. when $R_2^{E_1}>1$: aforward relay, consistent with Table 3. The calculation at $E_2$ is symmetric.

Remark 15

(Features of the Feng–Gavish model).

1.

The stability domain of $E_1$ has the form $R_1>1$, $R_2^{E_1}<1$, with the symmetric statement for $E_2$. This was proved in [5,47] by factoring the full Jacobian at the boundary equilibrium; in the relay formulation it follows from the Metzler transversal block.

2.

In the persistence region of Table 3, all boundary equilibria are unstable. Existence of $E_{12}$ is supplied by Theorem 5, not by the local relay theorem alone.

3.

Metzler transversality alone does not exclude a boundary Hopf bifurcation: it only governs thetransversalblock, and a boundary equilibrium could in principle lose stability through a conjugate eigenvalue pair of itstangentialblock instead. For the Feng–Gavish model specifically, this possibility is excluded because the tangential block at $E_1$ and $E_2$ is Hurwitz throughout their existence regions (item 1 above, citing [5,47]), so the only remaining boundary stability change is the real Perron eigenvalue of the missing-strain Metzler block. Consequently, for this model, any boundary loss of stability occurs through that real root, and a Hopf bifurcation, if present, must occur at the interior equilibrium $E_{12}$. This conclusion uses the model-specific tangential-Hurwitz fact together with transversal Metzler-ness; it is not a general consequence of Metzler transversality by itself.

6.5. What the Relay Table Proves, and What it Does Not: The Three-Strain Case

For a three-strain model, let $E_A$ denote the equilibrium with resident strain set $A\subset \{1,2,3\}$, and write $R_j^{E_A}:=R_j\left(E_A\right)$ for the invasibility number of strain j into the resident equilibrium $E_A$. The local relay theorem applies only to distance-one edges of the siphon lattice,

$$E_A⟶E_{A\cup \left\{j\right\}},j\notin A,$$

provided conditions (tgS), (trC), and (nzB) hold at $E_A$. Each such edge asserts the local existence and direction of a bifurcating branch, not the global existence of $E_{A\cup \left\{j\right\}}$.

The logical gap between levels grows with dimension, for three distinct reasons:

1.

One-strain level.$E_i$ exists iff $R_i^{E_0}>1$. This is a genuine local relay: the bifurcating branch is one-dimensional, its direction is given by the right Perron eigenvector of $J_{{\Sigma }_i}^{\perp }\left(E_0\right)$, and global existence follows from forward invariance of $F_i$ together with the Perron root crossing.

2.

Two-strain level. The relay theorem at $E_i$ gives a branch pointing into the interior of $F_j$ (strain j absent), but does not construct $E_{ij}$. Existence of $E_{ij}$ on the face $F_k=\{x_{{\Sigma }_k}=0\}$ is a two-dimensional fixed-point problem, solved here by reducing to a scalar RUR polynomial $P\left(s\right)$ and applying the intermediate value theorem after establishing the boundary factorization $P\left(s_i\right)\propto -(R_j^{E_i}-1)$.

3.

Three-strain level. Even if every pair $(E_i,E_{ij})$ satisfies the ordered invasion condition, existence of $E_{123}$ in $\mathrm{relint}\mathrm{\Gamma }$ requires solving three coupled fixed-point equations after imposing the three rank-one degeneracy conditions. The resulting system generally does not reduce to a single scalar polynomial, and positivity of the reconstructed root is not automatic. Pairwise invasion inequalities are natural candidate boundary conditions for constructing a three-strain equilibrium. They are not sufficient in general. Their necessity requires additional monotonicity or branch-connectivity hypotheses and should not be asserted without proof: a positive three-strain equilibrium need not, in full generality, imply that each strain can invade every lower-level resident equilibrium. We therefore separate three distinct questions — the local relay condition at each face, global continuation of a branch across a face, and existence of the full interior equilibrium — rather than conflate them.

Table 4. Logical status of existence statements in multi-strain relay tables. Each row requires strictly more than the row above it.

Equilibrium	Residents	Existence condition	Logical status	Proof tool
$E_0$	none	always	definition	conservation / forward invariance
$E_i$	strain i	$R_i^{E_0}>1$	local relay (proved)	Perron root crossing + forward invariance
$E_{ij}$	strains $i,j$	ordered invasion: $R_j^{E_i}>1$, $s_i<s_j$	face coexistence (proved)	scalar RUR polynomial + IVT
$E_{123}$	strains $1,2,3$	pairwise invasions $R_j^{E_i}>1$	candidate; sufficiency and necessity both open	open: three-strain elimination + positivity

Table 4. Logical status of existence statements in multi-strain relay tables. Each row requires strictly more than the row above it.

Equilibrium	Residents	Existence condition	Logical status	Proof tool
$E_0$	none	always	definition	conservation / forward invariance
$E_i$	strain i	$R_i^{E_0}>1$	local relay (proved)	Perron root crossing + forward invariance
$E_{ij}$	strains $i,j$	ordered invasion: $R_j^{E_i}>1$, $s_i<s_j$	face coexistence (proved)	scalar RUR polynomial + IVT
$E_{123}$	strains $1,2,3$	pairwise invasions $R_j^{E_i}>1$	candidate; sufficiency and necessity both open	open: three-strain elimination + positivity

Thus relay tables must be read hierarchically: each level is proved by a different method, and the tools used at level k do not extend automatically to level $k+1$. In particular,

$$\begin{array}{cc}\hfill \mathrm{local}\mathrm{boundary}\mathrm{instability}\mathrm{at}\mathrm{every}{E}_A& \neg \Rightarrow \mathrm{existence}\mathrm{of}{E}_{123}.\hfill \end{array}$$

and filling this gap for three or more strains remains an open problem. A sufficient condition in the spirit of Section 5.4 would be a three-strain analogue of the RUR reduction: a scalar polynomial $P_{123}\left(s\right)$ whose boundary values factor through the pairwise invasion deficits, and whose roots in the feasible region reconstruct a positive equilibrium. Whether such a reduction exists for rank-one three-strain models is an open question.

The scalar RUR/IVT argument of Section 5.4 is, in effect, one-dimensional degree theory: a sign change of $P\left(s\right)$ on an interval is exactly a nonzero mod-2 (or, tracking orientation, integer) degree of the map $s\mapsto P\left(s\right)$ relative to the two endpoints. This suggests a higher-dimensional formulation that does not presuppose a scalar reduction at all.

Proof 

(Relay degree conjecture). Let A be a resident strain set and suppose every missing adjacent strain has positive invasion index at the corresponding boundary equilibrium. Determine conditions under which the Brouwer degree of the interior-equilibrium map on the face generated by A is nonzero. For two strains, the boundary factorization of $P\left(s\right)$ (Lemma 7) computes this degree by a sign change. For three strains, a nonzero multidimensional degree may replace the unavailable scalar RUR polynomial $P_{123}\left(s\right)$ discussed above.    □

This appears a more promising route to the three-strain case than insisting that every three-strain model reduce to one scalar polynomial.

6.6. Numerical Experiment for Hopf Bifurcation

All computations for the two strains Feng-Chung-Gavish model may be found in the GavScan.wl file, together with a numerical scan of the parameter partition (Figure ); no Hopf bifurcations were observed in the range studied. This scan is evidence, not a theorem, and should be read as such: absence of Hopf points in a finite sample does not exclude them elsewhere in parameter space. A structural (symbolic) test that would upgrade this scan to a theorem, or to a rigorous algebraic candidate condition, is recorded as an open computational problem in the conclusion (Question Section 7).

Figure 3. Stability partition scan from GavScan.wl. The numerical classification (coloured points) agrees with the symbolic LAS-CEP derived from the four reproduction and invasibility numbers.

7. Conclusions

We developed a framework for two-strain epidemic models resting on siphon geometry, next-generation matrices, and boundary transcritical relays, yielding algorithmic relay tables and a boundary LAS-CEP partition for the Feng–Gavish model.

Table 5. Synthesis: what is proved in this paper versus what remains open, organized by the logical chain siphon geometry → invasion thresholds → local relays → persistence or coexistence.

Question	Proved here	Still open
Local invasion	$R_{\sigma }\left(E\right)-1$ has the sign of the transversal Perron root (Theorem 6)	Reducible or multiple simultaneously-critical transversal blocks (Theorem 8)
Local successor branch	$a_{\sigma }\left(E\right)b_{\sigma }\left(E\right)\ne 0$ gives a unique forward/backward relay with stability exchange (Theorem 7)	Degenerate classification when $a_{\sigma }\left(E\right)=0$ or $b_{\sigma }\left(E\right)=0$ (§Section 6.1)
Persistence	Acyclic boundary Morse decomposition criterion (Theorem 1, Corollary 1)	Efficient verification of acyclicity for many strains
Coexistence	Boundary-factorization sufficient criterion (Theorem 4, Theorem 5)	Necessity (Conjecture Section 5.3) and root multiplicity (Question Section 5.4)
Global dynamics	Boundary local-stability partition (Theorem 10)	Global asymptotic stability, interior Hopf, cycles, and interior multistationarity (Conjecture Section 6.3)

Table 5. Synthesis: what is proved in this paper versus what remains open, organized by the logical chain siphon geometry → invasion thresholds → local relays → persistence or coexistence.

Question	Proved here	Still open
Local invasion	$R_{\sigma }\left(E\right)-1$ has the sign of the transversal Perron root (Theorem 6)	Reducible or multiple simultaneously-critical transversal blocks (Theorem 8)
Local successor branch	$a_{\sigma }\left(E\right)b_{\sigma }\left(E\right)\ne 0$ gives a unique forward/backward relay with stability exchange (Theorem 7)	Degenerate classification when $a_{\sigma }\left(E\right)=0$ or $b_{\sigma }\left(E\right)=0$ (§Section 6.1)
Persistence	Acyclic boundary Morse decomposition criterion (Theorem 1, Corollary 1)	Efficient verification of acyclicity for many strains
Coexistence	Boundary-factorization sufficient criterion (Theorem 4, Theorem 5)	Necessity (Conjecture Section 5.3) and root multiplicity (Question Section 5.4)
Global dynamics	Boundary local-stability partition (Theorem 10)	Global asymptotic stability, interior Hopf, cycles, and interior multistationarity (Conjecture Section 6.3)

Question 4

(Structural test for Hopf bifurcation in the Feng–Gavish coexistence branch). The numerical scan of the Feng–Gavish LAS-CEP relay table (GavScan.wl, Figure ??) found no Hopf bifurcations in the range studied, but this is evidence, not a theorem. Carry out the following structural (symbolic) test, which was not executed here. At the coexistence equilibrium, form the characteristic polynomial

$${\chi }_{12}\left(\lambda \right)={\lambda }^n+a_1{\lambda }^{n-1}+\dots +a_n.$$

(i) Reduce its coefficients modulo the coexistence ideal of Section 5.4. (ii) Express the coefficients as rational functions of the scalar coordinate s via the RUR reconstruction of Lemma 6. (iii) Compute the Hurwitz determinants ${\mathrm{\Delta }}_1\left(s\right),\dots ,{\mathrm{\Delta }}_{n-1}\left(s\right)$ of ${\chi }_{12}$. (iv) Determine whether any ${\mathrm{\Delta }}_k\left(s\right)$ changes sign on $s\in (s_{\mathrm{\Delta }},\overline{s})$. (v) If some ${\mathrm{\Delta }}_k$ can vanish there, solve the joint system $P_G\left(s\right)=0,{\mathrm{\Delta }}_k\left(s\right)=0$ for candidate Hopf parameter values. (vi) Verify transversality and the sign of the first Lyapunov coefficient numerically only after a symbolic candidate has been isolated this way. Carried through, this programme would turn the present numerical scan into either a theorem excluding Hopf under explicit parameter inequalities, or a rigorous algebraic candidate condition for where a Hopf bifurcation could occur. □

Future work should focus on extending the RUR reduction to three-strain models, resolving Question 4 and Question 2, and developing algorithmic siphon-based tools for automatic stability classification in higher-dimensional epidemiological networks.

7.1. Some Notations, and a Perron-Volterra Lyapunov Function Ansatz

For rank one models, we will obtain below GAS partitions of the parameter space using Lyapunov functions (16) which combine Volterra entropy terms for the resident variables with Perron-weighted linear functionals for the invading variables, whose weights obtained from the rank one decompositions of the new infections matrices $F_{\sigma }$. Defining the family of Perron-Volterra candidate Lyapunov functions for this model uses several concepts:

1.

siphons (a chemical reaction network and Petri nets concept), which yield forward-invariant coordinate faces for chemical ODEs;

2.

the disease-free siphon, which is the union of minimal siphons;

3.

the fact that Jacobians on siphon faces have a triangular block form, cf. [10];

4.

the fact that a transversal Jacobian block on a siphon face is Metzler [11], which puts under spotlight the roles of its Perron eigenvectors;

5.

the Perron-Frobenius decomposition of the transversal Jacobian on the DFE siphon face $J_0^{\perp }$, which defines the “invading strains", with transversal blocks $J_{\sigma }^{\perp }$;

6.

rank one decompositions $J_{\sigma }^{\perp }=F_{\sigma }-V_{\sigma }$, with the “ new infections" matrix $F_{\sigma }$ having rank one, an assumption often verified in ME models.

Some basic notations are summarized in the table below:

Notation for rank one models (dependent on the regular splitting)

Symbol	Meaning
$J_{\sigma }^{\perp }=F_{\sigma }-V_{\sigma }$	transversal Jacobian on the siphon face $\sigma$, written in some regular splitting form
$F_{\sigma },V_{\sigma }$	new-infections matrix (often of rank-one) and transition matrix associated with the regular splitting of $J_{\sigma }^{\perp }$
$w_k,{\pi }_k$	right and left Perron vectors of new-infections strain k; ${\pi }_k$ is a probability vector
${\mathcal{R}}_k={\pi }_k^{\top }{V}_k^{-1}{w}_k$	per-susceptible reproduction coefficient of block k
$R_k\left(s\right)$	invasion function=Perron eigenvalue of block k at susceptible level s
$R_k=R_k\left(s_0\right)$	basic reproduction number of block k at the DFE
${\tilde{\pi }}_k^{\top }={\pi }_k^{\top }{V}_k^{-1}/{\mathcal{R}}_k$	normalized left Perron eigenvector of $K_k=F_k{V}_k^{-1}$; equivalently, normalized left nullvector of $s_k{F}_k-V_k$ at $s_k=1/{\mathcal{R}}_k$, with ${\tilde{\pi }}_k^{\top }{w}_k=1$
${\tilde{w}}_{k,j}={\tilde{\pi }}_{k,j}{w}_{k,j}$	probability weights for Jensen’s inequality (${\tilde{w}}_{k,j}>0,{\sum }_j{\tilde{w}}_{k,j}=1$)
$E_p=(s_p,0,\dots ,0,z_p^{*},0,\dots ,0)$	boundary equilibrium where block p is the unique resident block
$z_k$	state vector of block k
$z_k^{*}$	equilibrium value of block k at $E_k$
${\tilde{i}}_k={\pi }_k^{\top }{z}_k$	Perron aggregate infective level of block k
${\Omega }_p$	positively invariant region where $E_p$ is the resident equilibrium
$s_0$	disease-free susceptible level
$s_p=1/{\mathcal{R}}_p$	susceptible coordinate of the boundary equilibrium $E_p$
$u=s/s_p$	normalized susceptible variable
$y_j=z_{p,j}/z_{p,j}^{*}$	normalized resident coordinates
$G\left(u\right)=u-1-lnu$	Volterra entropy
$\tilde{G}(s,s_p)=s_{p}G\left({\displaystyle \frac{s}{s_p}}\right)$	two variable Volterra entropy
$H_p$	resident Perron–Volterra entropy defined in (16)
$L_p=H_p+{\sum }_{k\ne p}{\tilde{\pi }}_k^{\top }{z}_k$	Perron–Volterra Lyapunov function for the equilibrium $E_p$
${\displaystyle \sum _{k\ne p}\left(\frac{1}{{\mathcal{R}}_p}-\frac{1}{{\mathcal{R}}_k}\right){\tilde{i}}_k}$	canonical invader contribution to the derivative of $L_p$
${\tilde{H}}_p(u,y)$	resident entropy production term for $H_p$, defined in (17), used for proving Lyapunov inequality

$$H_p=s_{p}G\left({\displaystyle \frac{s}{s_p}}\right)+{\displaystyle \sum _{j=1}^{n_p}}{\tilde{\pi }}_{p,j}{z}_{p,j}^{*}G\left({\displaystyle \frac{z_{p,j}}{z_{p,j}^{*}}}\right)=\tilde{G}(s,s_p)+{\tilde{\pi }}_p^{\top }\tilde{G}(z_p,z_p^{*}),$$

(16)

where

$$\tilde{G}(z_p,z_p^{*})={\left(z_{p,j}^{*}G\left({\displaystyle \frac{z_{p,j}}{z_{p,j}^{*}}}\right)\right)}_{j=1}^{n_p}.$$

The resident entropy production term in the derivative of the resident entropy $H_p$ is:

$${\tilde{H}}_p(u,y)=\left(1-{\displaystyle \frac{1}{u}}\right)(1-uY)+{\displaystyle \sum _{j=1}^{n_p}}{\tilde{w}}_{p,j}\left(1-{\displaystyle \frac{1}{y_j}}\right)(uY-y_j),Y={\displaystyle \sum _{j=1}^{n_p}}{\tilde{w}}_{p,j}{y}_j$$

(17)

7.2. The weighted Perron–Volterra Ansatz

The $L_p$ of (16) fixes a single overall scalar for the whole function. The canonical invader directions ${\tilde{\pi }}_k^{\top }={\pi }_k^{\top }{V}_k^{-1}/{\mathcal{R}}_k$ ($k\ne p$) remain unchanged, but each invading block may be given its own independent positive multiplier $a_k$, alongside a multiplier b on the resident entropy $H_p$ itself:

$$L_p^{b,\left(a\right)}=b{H}_p+\sum _{k\ne p}{a}_k{\tilde{\pi }}_k^{\top }{z}_k,b>0,a={\left(a_k\right)}_{k\ne p}\gg 0.$$

(18)

Here a denotes the whole tuple ${\left(a_k\right)}_{k\ne p}$ of invader multipliers, one per block other than p; the index k ranges over $k\ne p$ only inside the defining sum, never in the superscript $\left(a\right)$ itself. The Ansatz $L_p=H_p+{\sum }_{k\ne p}{\tilde{\pi }}_k^{\top }{z}_k$ of (16) is the special case $b=1$, $a=(1,\dots ,1)$. Which tuples $(b,a)$ keep $L_p^{b,\left(a\right)}$ a genuine Lyapunov function – the admissible weight cone at $E_p$ – is investigated, for the papers using this notation file, wherever they discuss the admissible cone explicitly.