# Algorithmic approach for an unique definition of the next generation matrix

Florin Avram1, Rim Adenane2, Lasko Basnarkov3, Matthew Johnston4

###### Abstract

The basic reproduction number $R_{0}$ is a concept which originated in
population dynamics, mathematical epidemiology, and ecology, and is closely
related to the mean number of children in branching processes (reflecting the
fact that the phenomena of interest are well approximated by branching
processes, at their inception). Despite the very extensive literature around
$R_{0}$ for deterministic epidemic models, we believe there are still aspects
which are not fully understood. Foremost is the fact that $R_{0}$ is not a
function of the original ODE model, unless we include in it also a certain
$(F,V)$ gradient decomposition, which is not unique. This is related to the
specification of the "infected compartments", which is also not unique. A
second interesting question is whether the extinction probabilities of the
natural continuous time Markovian chain approximation of an ODE model around
boundary points (disease-free equilibrium and invasion points) are also
related to the $(F,V)$ gradient decomposition, as suggested below in Section
2.4.

We offer below three new contributions to the literature:

1) We offer a universal algorithmic definition of a $(F,V)$ gradient
decomposition (and hence of the resulting $R_{0}$), which requires a minimal
input from the user, namely the specification of an admissible set of
disease/infection variables. We also present examples where other choices may
be more reasonable, with more terms in $F$, or more terms in $V$. This trade-
off is explained briefly in section 2.2, Remark 5, and further examined in
examples in Sections 5.1, 8.2.3.

2) We glean out from the works of Bacaer a fixed point equation (8) for the
extinction probabilities of a stochastic model associated to a deterministic
ODE model, which may be expressed in terms of the $(F,V)$ decomposition. The
fact that both $R_{0}$ and the extinction probabilities are functions of
$(F,V)$ underlines the centrality of this pair, which may be viewed as more
fundamental than the famous next generation matrix $FV^{-1}$. The equation (8)
may be rarely solved explicitly; however, even when this is not the case,
useful quasi-explicit solutions may be provided via "rational univariate
representations" (for which an algorithm which works often is also provided).

3) We suggest introducing a new concept of sufficient/minimal
disease/infection set (sufficient for determining $R_{0}$). More precisely,
our universal recipe of choosing "new infections" once the "infections" are
specified suggests focusing on the choice of the latter, which is also not
unique. The maximal choice of choosing all compartments which become $0$ at
the given boundary point seems to always work, but is the least useful for
analytic computations, therefore we propose to investigate the minimal one. As
a bonus, this idea seems useful for understanding the Jacobian factorization
approach for computing $R_{0}$ (see below). We view this as a method for
obtaining an approximation, which we show to always yield upper or lower
bounds of the true $R_{0}$, depending on whether $R_{0}\leq 1$ or not. This
raises interesting questions of determining conditions on the epidemic model
which ensure tightness of our bound, of getting better bounds when tightness
doesn’t hold, and seems to deserve further work, only the surface of which is
skimmed below.

Last but not least, we offer Mathematica scripts and implement them for a
large variety of examples, which illustrate that our recipe offers always
reasonable results, but that sometimes other reasonable $(F,V)$ decompositions
are available as well.

1 Laboratoire de Mathématiques Appliquées, Université de Pau, 64000, Pau,
France<EMAIL_ADDRESS>
2 Laboratoire des Equations aux dérivées partielles, Algébre et Géométrie
spectrales, département des Mathématiques, Université Ibn-Tofail, 14000,
Kenitra, Maroc<EMAIL_ADDRESS>
3 Faculty of Computer Science and Engineering, Ss. Cyril and Methodius
University in Skopje, 1000 Skopje, North Macedonia;
<EMAIL_ADDRESS>
4 Department of Mathematics + Computer Science Lawrence Technological
University, 21000 W 10 Mile Rd, Southfield, MI 48075, USA<EMAIL_ADDRESS>
Keywords: deterministic epidemic model, disease-free equilibrium, stability
threshold, basic reproduction number, $(F,V)$ gradient decomposition, next
generation matrix, Jacobian approach, CTMC stochastic model associated to a
deterministic epidemic model, probability of extinction, rational univariate
representation

###### Contents

  1. 1 Introduction
  2. 2 A bird’s eye view of mathematical epidemiology: the disease-free equilibrium, the next generation matrix, and an algorithmic definition of a stability threshold associated to the basic reproduction number $R_{0}$
    1. 2.1 The disease free equilibrium (DFE)
    2. 2.2 $(F,V)$ gradient decompositions, the next generation matrix, $R_{0}$, and a simple recipe for computing them
    3. 2.3 An algorithmic $F-V$ decomposition
    4. 2.4 A multi-dimensional birth and death CTMC process associated to a $(F,V)$ decomposition, its branching process approximation, and the Bacaer equation for the probability of extinction
    5. 2.5 The Jacobian factorization bound
    6. 2.6 The “rational univariate representation" (RUR) and the reduced order quasi -stationary approximation
  3. 3 $R_{0}$ and extinction probabilities for the SEIR epidemic model
  4. 4 Rank one host-only models with pathogen, and $R_{0}$ readable from the flow-chart
    1. 4.1 The SEIARW model with “catalyzing pathogen" of [32] has rank one next generation matrix and $R_{0}=R_{J}$
  5. 5 Target-infection-virus models
    1. 5.1 Two admisible $(F,V)$ decompositions and $R_{0}$’s for the three dimensional model of [48]
    2. 5.2 Two distinct approximate extinction probabilities, one for each admisible $(F,V)$ decomposition for the model of [48]
  6. 6 Multi-strain host only models
    1. 6.1 The two-strain SIS tuberculosis model of [56, Sec4.4]
    2. 6.2 The minimal disease set of the multi-strain host only dengue model with antibody-dependent enhancement (ADE) [4]
    3. 6.3 Effects of single-strain vaccination on the dynamics of a multi-strain host only dengue model with ADE
      1. 6.3.1 The Jacobian $R_{J}(X)$ is the max of two polynomials, obtained using a minimal disease set
      2. 6.3.2 The endemic boundary equilibrium $E_{i}$ exist iff $R_{i}(E_{0})>1$
      3. 6.3.3 The recipe next generation matrix $R_{0}$ and the Jacobian factorization one coincide
      4. 6.3.4 The invasion reproduction number of $E_{i}$ is given by $R_{j}(E_{i})$
  7. 7 Vector-host models
    1. 7.1 The Jacobian $R_{0}$ is the square of the recipe NGM $R_{0}$ for the dengue vector-host model without demography of [16]
    2. 7.2 The two groups model in [42, (5.8)] does not obey a square relation
  8. 8 Multi-strain vector-host models
    1. 8.1 A two-strain vector-host model of Feng and Velasco-Hernández [28], where the square relation holds for the basic reproduction number
      1. 8.1.1 Invasion numbers of [28]
    2. 8.2 The Dengue- Zika model with coinfection and ADE [47]
      1. 8.2.1 The dengue only resident fixed point $E_{d}$
      2. 8.2.2 The Zika only resident fixed point $E_{z}$
      3. 8.2.3 The dengue invasion reproduction number (IRN) and two possible $(F,V)$ decompositions
  9. 9 Conclusions
  10. 10 Appendix: the implementation of the Jacobian factorization approach

##  1 Introduction

Motivation. Mathematical epidemiology has started by proposing simple models
for specific epidemics, and computing explicitly certain important
characteristics like the basic reproduction number and the final size; for
example the SIR model was introduced, among other concepts, in the celebrated
"A contribution to the mathematical theory of epidemics" [39]. The most
fundamental and actually only general result of the field, due to Diekmann,
Heesterbeek, Van den Driesche and Watmough, is expressing the disease-free
equilibrium stability domain in terms of $R_{0}$, which is defined as the
Perron Frobenius eigenvalue of a certain $(F,V)$ gradient decomposition (this
is presented in detail in Section 2.2). But, since the $(F,V)$ decomposition
is not unique, it seems to us that the question of what is $R_{0}$ still
deserves further discussion.

On the other hand, one may note that nowadays, mathematical epidemiologists
typically either restrict themselves to low-dimensional models, resolved
symbolically, even by hand, or consider very complicated models which are
resolved only numerically, for particular values gleaned from the medical
literature. Missing from here are moderately complex models, which may be
solved partly symbolically for any values of the parameters, but where the use
of computer algebra systems (CAS) is either indispensable or greatly
facilitating. Even in the case of papers belonging to this level – see for
example [47], the role of the CAS is deemphasized. Our paper is also an
attempt to cast the CAS as one of the main heroes of our story.

Our main result. We provide below, for the first time, a universal recipe for
choosing a natural $(F,V)$ gradient decomposition, which only requires
specifying the disease compartments (a subset of those which are zero for the
boundary point under consideration) (informally, these are not far
conceptually from the so called fast components of singular perturbation
theory). This decomposition is useful both for determining $R_{0}$ and for
computing the extinction probabilities of an associated stochastic model. We
identify also examples in which the $(F,V)$ decomposition is not unique, and
in which choosing another decomposition with $F$ of lower rank may be
beneficial for simplifying the $R_{0}$ formula.

First restriction (among others to follow). In this paper, we will restrict to
mathematical epidemiology models for which there exist at least two possible
special fixed states. The first, the disease-free equilibrium (DFE),
corresponds to the elimination of all possible compartments involving
sickness, and will be assumed to be unique. Typically, this point is locally
stable only for certain values of the parameters. Outside this domain it is
typically replaced by another fixed point, which will be called “endemic" if
all its components are positive, and "resident boundary point" otherwise.

Importantly, the stability of the DFE may be related to the historically
famous basic reproduction number and net reproduction rate–see (1). These
pillar concepts in population dynamics, mathematical epidemiology, virology,
ecology, etc, were already introduced by the father of mathematical demography
Lotka – see [41, 26], and also the introduction of the book [13].

A bit of history of the net reproduction rate $\mathcal{R}$, and its evolution
into the mathematical concept of basic reproduction number/stability threshold
$R_{0}$. Loosely speaking, in the case of only one infectious class, the net
reproduction rate $\mathcal{R}$ describes the expected number of secondary
cases which one infected case would produce in a homogeneous, completely
susceptible population, during the lifetime of the infection. This description
is especially relevant at the start of an epidemics, when the dynamics is well
approximated by that of a branching process (a fact which goes back to
Bartlett and Kendall – see for example [24, 38]). The main characteristic of a
branching process is the “fertility", i.e. the expected number of descendants
one individual produces in the next generation. As a consequence, the
branching result insuring extinction when the fertility is less than one
translates in epidemiology into local stability results of the disease free
equilibrium involving $\mathcal{R}$.

The reproduction number $\mathcal{R}$ intervened already, in a particular
case, in the foundational paper “A contribution to the mathematical theory of
epidemics" [39], which showed that:

  1. 1.

The condition

$R_{0}<1,\quad\text{where }R_{0}=\;\mathsf{s}_{dfe}\mathcal{R},$ (1)

implies local stability of the DFE. Here $\mathcal{R}$ is the net reproduction
rate (number of secondary infections produced by one infectious individual),
and $\;\mathsf{s}_{dfe}$ is the fraction of susceptibles at the DFE.

  2. 2.

The condition

$R_{0}>1$

implies instability of the DFE.

With more infectious classes, one deals at the inception of an epidemics with
approximate multi-class branching processes, whose stability is determined by
a “next generation matrix " (NGM) –see section 2.2.

The “Jacobian approach" for computing $R_{0}$. For big size problems, this
approach is doomed to fail symbolically, since it is equivalent to the Routh-
Hurwicz conditions (RH), which rarely succeed symbolically beyond dimension
$4$ (also, RH is irrelevant numerically, since the eigenvalues themselves are
just as easy to compute). Therefore, we studied below a variant, the “Jacobian
factorization approach", which focuses on an approximation, which we show to
yield always upper or lower bounds of the NGM $R_{0}$, depending on whether
$R_{0}\leq 1$ or not – see Theorem 1. Several questions around this bound are
scattered below in Sections 6.3, 8.1, 8.2.

Note, as mentioned in [34], that an example where the Jacobian method does not
yield $R_{0}$ is offered in Diekmann & Heesterbeek (2000, Exe 5.43), and that
Roberts & Heesterbeek (2003), suggest that when threshold parameters
determined from the Jacobian do not have the biological interpretation of
dominant eigenvalue of the next generation matrix, then they should not be
called basic reproductive ratios, nor denoted $R_{0}$ (we follow their
suggestion and use the notation $R_{J}$ in this case).

The dilemma of the several different methods for computing $R_{0}$ has been
discussed in many papers, see for example [34, 40]. But this is a direct
consequence of the non-uniqueness of the $(F,V)$ decomposition.

Deterministic or stochastic models? Most of the mathematical epidemiology
papers belong exclusively to one of these two paradigms. However, any
deterministic model may also be viewed as a stochastic continuous time Markov
field (CTMC), evolving on the integers. One interesting CTMC, which seems not
to have been discussed before, is presented in section 2.4.

Contents. Our paper is structured as follows. Section 2 recalls the definition
of the DFE and provides our algorithmic definition of the (F,V) decomposition,
in the form of a Mathematica script, as well as a discussion – see Remark 5–
of why other decompositions might turn out useful. This section also provides
a new equation (8) for computing extinction probabilities for associated
continuous time Markov chain models in terms of the (F,V) decomposition, and
shows that the Jacobian factorization approach yields upper bounds and lower
bounds for NGM $R_{0}$’s, in section 10.

We turn then to a series of examples, chosen to help investigating what may be
the major open problem in the field nowadays, which, in our opinion, is
relating on one hand $R_{0}$, and on the other hand the extinction
probabilities – see below– and the duration of minor epidemics [8, 9, 58, 53,
44], which is not further touched here.

Some hint that such a connection might exist may be gleaned from the SIR case
studied by Whittle, and the SEIR case, recalled in section 3.

Section 4.1 presents a host only model, with a single susceptible class and an
$F$ matrix of rank one, where the formula of $R_{0}$ may be “guessed by
inspection" of the flow chart. This kind of examples have kept alive the hope
of “interpretable $R_{0}$ formulas", as illustrated in other recent papers –
see for example [32, 51]. These authors start by presenting simple cases, and
propose then algorithms for more complex cases based on the graph structure of
the flow chart, which in our opinion are not sufficiently detailed or
documented. While it may well be that tools like Petri nets, as proposed in
the second paper, will one day succeed for resolving flow chart with certain
structures, this does not seem to have happened yet. Also, for models with
next generation matrix of high-rank, the lack of simple formulas for $R_{0}$
and of "simple biological interpretations" is natural to be expected; simple
formulas for the spectral radius can only be a consequence of a simple graph
structure which has not been pinpointed yet.

Sections 5.1, 7.1, offer two examples in which several $R_{0}$ formulas were
offered in the literature, but we are at a loss of how to choose among them.
In the first case (a virus-tumor model), the recipe $R_{0}$ is simpler than
its competitor, but in the second case (a vector-host model), it is more
complicated.

Section 6.2 shows that the boundary equilibria, and the (invasion)
reproduction numbers may be easily computed with our scripts; to illustrate
this, we use a two-strain host only model from [42, Ch.8], where our recipe
NGM yields the same answer as that given by the Jacobian factorization.

Section 6.3 offers another two-strain host only example, this time including
also vaccination, in which our recipe NGM yields again the same answer as that
given by the Jacobian factorization.

Section 7.2 offers an example from the textbook of [42] in which the square
relation stops holding.

Sections 8.1, 8.2 offer yet more examples, this time in the two-strain vector-
host context, in which our recipe NGM yields an $R_{0}$ formula which is
precisely the square root of that given by the Jacobian approach. Note that
here the first of the three elegant relations concerning the invasion numbers
from Section 6.3 – see Remark 20–holds, but the other two seem to break down.

The last subsection provides, for the invasion numbers, a second example where
another choice of $R_{0}$ may be more reasonable, on the grounds of leading to
a simpler answer (but the admissibility requirement forces then extra
assumptions on the parameters).

##  2 A bird’s eye view of mathematical epidemiology: the disease-free
equilibrium, the next generation matrix, and an algorithmic definition of a
stability threshold associated to the basic reproduction number $R_{0}$

###  2.1 The disease free equilibrium (DFE)

The DFE may be defined as a “maximal boundary state", and may be found by
identifying a maximal sub-system of the ODE epidemic model which factors

$\bm{{i}}^{\prime}=M\bm{{i}},$ (2)

where prime denotes derivative with respect to time, and $M$ is a matrix may
depend on $\bm{{i}},$ but may not explode in the domain of interest, which we
will take for simplicity to be $\mathbb{R}_{+}^{n}$.

###### Remark 1

One fixed point of this system is $\bm{{i}}=0$. This motivates us to call the
components $\bm{{i}}$ disease or infectious states. The set of all its indices
will be denoted by $\mathcal{I}$. Note that specifying ‘$\mathcal{I}$ induces
a partition of both the coordinates and the equations of our original system
into infection (eliminable) components, and “non-infection" (the others).

The eventual other fixed points may be found by solving $M=0$ together with
the other non-infection equations under the condition $\bm{{i}}=0$.

In this paper we will assume uniqueness of the DFE, at least after excluding
biologically irrelevant fixed points, like an unreachable origin.

We end this section with the very elementary script that implements this. Note
that any ODE model “mod" (like SIR, etc…), is a pair mod= (dyn,X) consisting
of a vector field “dyn" and a list of variables “X", and that to find any
boundary fixed point it suffices to know the set of indices “inf" where it is
$0$, so that we solve the system “dyn==0" under the condition “X[[inf]]->0".
But, since sometimes only numeric solutions are possible, our DFE Mathematica
script below has also an optional numerical condition parameter “cn", which is
taken by default as the empty set.

     DFE[mod_,inf_,cn_:{}]:=Module[{dyn,X},
      dyn=mod[[1]]/.cn;X=mod[[2]];
      Solve[Thread[dyn==0]/.Thread[X[[inf]]->0],X]];

For the non-Mathematica users, only the Solve command is relevant, the others
being just Mathematica implementation details.

###  2.2 $(F,V)$ gradient decompositions, the next generation matrix, $R_{0}$,
and a simple recipe for computing them

From now on, the infection equations (2) will be rewritten as

$\bm{{i}}^{\prime}=\mathcal{F}-{\mathcal{V}}=(F-V)\bm{{i}}.$ (3)

Of course, such a decomposition is not unique, but we will also ask, following
[24, 56, 54], that $F$, the gradient of $\mathcal{F}$, is a matrix with
nonnegative elements, and $-V$, the gradient of $-{\mathcal{V}}$ is a
Markovian generating matrix (i.e. a matrix with non-negative off-diagonal
elements and non-positive row-sums). Conceptually, $F$ models inputs to the
disease compartments from outside ("new infections"), and $-V$ models
transfers between the disease compartments. Still, a priori, the decomposition
(3) is not unique.

###### Example 1

Let us illustrate this via a SIR example with superinfection, in which the
classes S and R play symmetric roles, inspired by the works of [7, 36, 1]

$\displaystyle\begin{cases}\;\mathsf{s}^{\prime}(t)&=b_{s}-\;\mathsf{s}(t)\left[\beta_{s}\;\mathsf{i}(t)(1+\xi\;\mathsf{i}(t))+d_{s}\right]+i_{s}\;\mathsf{i}(t)+\gamma_{r}\;\mathsf{r}(t)\\\
\;\mathsf{i}^{\prime}(t)&=\;\mathsf{i}(t)\Big{[}\left[\beta_{s}\;\mathsf{s}(t)+\beta_{r}\;\mathsf{r}(t)\right](1+\xi\;\mathsf{i}(t))\Big{]}-d_{i}\;\mathsf{i}(t)\\\
\;\mathsf{r}^{\prime}(t)&=b_{r}-\;\mathsf{r}(t)\left[\beta_{r}\;\mathsf{i}(t)(1+\xi\;\mathsf{i}(t))+d_{r}\right]+i_{r}\;\mathsf{i}(t)+\gamma_{s}\;\mathsf{s}(t)\end{cases}.$

Here, the only infection equation, the second, is already written in a
decomposed form
$\mathcal{F}-{\mathcal{V}},{\mathcal{V}}=d_{i}\;\mathsf{i}(t)$, and
$F=\Big{[}\left[\beta_{s}\;\mathsf{s}(t)+\beta_{r}\;\mathsf{r}(t)\right](1+\xi\;\mathsf{i}(t))\Big{]}+\xi\;\mathsf{i}(t)\left[\beta_{s}\;\mathsf{s}(t)+\beta_{r}\;\mathsf{r}(t)\right]$.

Note that for the application of the next generation matrix method we must
plug finally $\;\mathsf{i}=0$; therefore, the second term in $F$, due to
"superinfection", is irrelevant for this purpose.

###### Remark 2

The possible non-uniqueness of the decomposition brings us to a delicate point
in mathematical epidemiology. Anticipating a bit, since $R_{0}$ is the Perron
Frobenius eigenvalue of $FV^{-1},$ strictly speaking, $R_{0}$ is not
determined just by an ODE epidemical model, but also by the $(F,V)$ gradient
decomposition. If we want that an ODE epidemical model determines uniquely an
$R_{0}$, we must include in the definition of an ODE epidemical model also the
$(F,V)$ gradient decomposition we adopt.

###### Definition 1

An ODE epidemic model is an ODE dynamical model in which a certain subset of
equations, usually called "disease/infection" equations, and to be referred
from now on also as zeroable set, admits at least one admissible decomposition
(3), with $(\mathcal{F},{\mathcal{V}})$ satisfying the conditions (A1-A5) of
[54].

###### Remark 3

Note that (3) is the most common model used in population dynamics. This makes
natural to define informally ODE epidemic models as population dynamics models
(3), with extra equations modeling interactions with the non-disease
compartments, which admits at least one admissible decomposition.

###### Remark 4

The definition of ODE epidemic models above is imprecise, since it does not
list all the requirements we must put on an ODE model. Some reasonable
restrictions are

  1. 1.

Essentially nonnegative processes having a non-empty set of disease classes,
so we deal with an epidemic (note however that we define disease classes in
the sense of classes which satisfy (2), which excludes for example importation
models).

  2. 2.

Processes with a unique DFE, at least after excluding biologically irrelevant
fixed points, like an unreachable origin.

  3. 3.

The local stability domain of the DFE is non-empty, and not the full set.

  4. 4.

The dynamical system has polynomial coefficients, to be able to take advantage
of the remarkable symbolic computation tools available for this class.

We make these assumptions because they are satisfied by most mathematical
models which have already been used for modelling real life biological
phenomena. However, these assumptions might not be enough, and further ones
might be necessary for obtaining the currently missing precise definition of
“real life ODE mathematical epidemiology models".

###### Remark 5

Admissible decompositions need not be unique. A priori, one may "move terms
from F to V", to lower its rank and simplify the formula for $R_{0}$, and also
"move off-diagonal terms from V to F", which enlarges the domain of parameters
which ensure that $V^{-1}$ has positive terms. There is a tradeoff between
these two possible moves, since simplicity of $R_{0}$ comes at the cost of
extra assumptions on the parameters. Our universal decomposition seems to
strike a balance between the two directions.

###### Remark 6

It was emphasized from the outstart – see for example [22, 40, 18, 55], that
an ODE mathematical epidemiology model might have several “admissible
decompositions", which might yield distinct next generation matrices and
distinct $R_{0}$’s.

This was viewed as a richness rather than a default. A similar point of view
is advocated in the recent paper [21], which argues that the class of all
admissible NGM’s associated to an ODE epidemiology model must be studied as an
ensemble, rather than focusing on individual representatives.

For any admissible decomposition, Diekmann, Heesterbeek, Van den Driesche and
Watmough established the following celebrated DFE stability theorem:

###### Proposition 1

For any admissible decomposition $(F,V)$, let

$R_{0}=\rho(FV^{-1})$

denote the Perron Frobenius eigenvalue of the next generation matrix. Then,
the DFE is unstable on $R_{0}>1,$ and locally stable on $R_{0}<1$ [24, 23, 56,
54].

For a recent historical overview of $R_{0}$, next generating matrices, and
their calculation in many examples, we refer the reader to the delightful
paper [16].

Unfortunately, the standard definition of a next generation matrix (and hence
of $R_{0}$) involves concepts like “new infections", which were defined in the
original papers based on epidemiological considerations, and require therefore
the intervention of a human expert. This had created the impression that this
method can not be encapsulated into a computer program. However, we offer and
implement below a simple algorithmic definition, based only on the structure
of the system and of the “infectious/disease equations".

Our proposal is to use a special F-V decomposition with F constructed as the
positive part of all the interactions in the disease equations which involve
both disease compartments and input/susceptible ones. The latter are defined
as the complement of the disease compartments, after the possible removal of
output compartments, which may be specified as deterministic functions of the
other compartments (i.e. may be computed, once the other compartments are
known). Note now that the concept of "positive part of the interactions" may
be hard to pinpoint mathematically, but useful enough to have been implemented
in CAS’s (Mathematica, Maple, Sage, etc); this made us adopt the following
definition:

###### Definition 2

For a given zeroable set, an admissible $(F,V)$ gradient decomposition (3) is
one where $F$, the gradient of $\mathcal{F}$, does not contain in its expanded
form syntactic minuses in its CAS representation, and where $V$, the gradient
of ${\mathcal{V}}$ is such that $-V$ is a “sub-generating matrix", under the
assumption of nonnegativity of all the model parameters. Note that this last
condition implies, but is not equivalent to the existence of $V^{-1}$ and to
the nonnegativity of all its elements.

The problem of whether the $R_{0}$ of the decomposition provided satisfies
under certain conditions the stability theorem of Van den Driessche and
Watmough is still open; therefore it should be viewed for now just as a recipe
that works well in simple cases.

After lots of experimenting, we have found only few cases – see for example
Section 8.2 where the recipe NGM has a serious competitor; it is for computing
the invasion reproduction number for a two-strain vector-host model with
altered infectivity for co-infected vectors, and with ADE (antibody-dependent
enhancement).

###  2.3 An algorithmic $F-V$ decomposition

We complement now the famous $\mathcal{F}-{\mathcal{V}}$ “equations
decomposition" and next generation matrix method of [24, 56, 54] by an
algorithmic $F-V$ decomposition.

  1. 1.

The user supplies the model “mod" (a pair containing the RHS of the dynamical
system, and its variables), and the indices “inf" of the disease (or
infectious) variables; the indices of the other compartments are denoted by
“infc".

  2. 2.

Subsequently, the Jacobian of the infectious equations $M$ with respect to the
corresponding variables is computed.

  3. 3.

Define the interaction terms as terms in $M$ which contain variables
$s\in\text{infc}$, and which, if positive, must end up in $F$. Their
complement, denoted by $V1$, will form part of $V$.

  4. 4.

A first guess for $F$, $F1$ is constructed as the complement of $V1$. It
contains all the interaction terms (which involve both disease and susceptible
compartments).

  5. 5.

$F$ is obtained by retaining only the positive part of the matrix $F1$, i.e.
the terms which do not contain syntactic minuses. 444we use the simplest
algebraic representation of the equations, and do not study the effect which
algebraic manipulations introducing minuses might have. Finally, $V1$ is
increased to $V$, which is the complement of $F$.

  6. 6.

The script outputs {M,V1,F1,F,V,K}.

    NGM[mod_,inf_]:=Module[{dyn,X,infc,M,V,F,F1,V1,K},
       dyn=mod[[1]];X=mod[[2]];
       infc=Complement[Range[Length[X]],inf];
       M=Grad[dyn[[inf]],X[[inf]]]
       (*The jacobian of the infectious equations*);
       V1=-M/.Thread[X[[infc]]->0]
       (*V1 is a first guess for V, retains all gradient terms which
       disappear when the non infectious components are null*);
       F1=M+V1
       (*F1 is a first guess for F, containing all other gradient terms*);
       F=Replace[F1, _. _?Negative -> 0, {2}];
       (*all terms in F1 containing minuses are set to 0*);
       V=F-M;
       K=(F . Inverse[V])/.Thread[X[[inf]]->0]//FullSimplify;
    {M,V1,F1,F,V,K}]

Note that our NGM script requires a minimal input from the user, just the
specification of the disease compartments; there is no need to specify “new
infections".

The results of this decomposition seem to yield correct results in all the
examples from the literature we checked. We would like to add that for
dynamical systems satisfying the four conditions in the remark 4, this
decomposition yields “admissible gradient decompositions", in the sense that
$V^{-1}$ will contain only non-negative terms, and that it is furthermore
obtainable from an equations decomposition which is admissible in the sense of
[54] (see Definition 1), and yields therefore the correct stability domain.

###### Remark 7

Note that the “Replace" command in the script uses the powerful Mathematica
capability of applying a “rule" to parts of an “expression", specified by
“levelspec", and that it was furnished to us by the user Michael E2 in

    https://mathematica.stackexchange.com/questions/286500/
    how-to-set-to-0-all-terms-in-a-matrix-which-contain-a-minus
    /287406?noredirect=1#comment715559_287406

Finally, let us discuss an alternative possible implementation. We could just
provide NGM with the right-hand side of the differential equations, compute
the steady states, specify one of them, and then define the infected classes
as the components with zeros.

However, this would be unpractical, since for the majority of the models with
explicit DFE the other fixed points are either not explicit, or require very
long execution times. It is therefore much simpler to have the user help the
AI by providing it with $\mathcal{I}$, which leads immediately to the matrix
$M$. Essentially, we jump directly to the factorization (2) of the infected
equations, postponing solving the non-infection variables later.

###  2.4 A multi-dimensional birth and death CTMC process associated to a
$(F,V)$ decomposition, its branching process approximation, and the Bacaer
equation for the probability of extinction

The works of Kendall and Bartlett suggest that ODE epidemic models may be
associated to corresponding birth-and-death CTMC processes, and then
approximated further by branching process.

Citing [31] :"It has been noted by Bartlett (1955), p. 129, that for an
epidemic in a large population, the number of susceptibles may, at least in
the early stages of an outbreak, be regarded as approximately constant at its
initial value and that this approximation will continue to hold throughout the
course of an epidemic, provided that the final epidemic size is small relative
to the total susceptible population. Thus the general epidemic process may be
approximated by a simple birth-and-death process."

To make this more precise, a $(F,V)$ decomposition (3) determines a natural
associated multi-dimensional birth and death CTMC process, by fixing the
values of the non-disease variables, so that the matrices $(F,V)$ depend only
on $\bm{{i}}$, and interpreting the transition rates between compartments as
rates of BD transitions.If the CTMC has rates which are linear in the disease
variables, one may associate to it a branching process, and take advantage of
the well-known equation for extinction probabilities. This procedure has been
detailed in previous works like [8, 9, 58, 53, 44], and used to approximate
extinction and invasion probabilities, as well as the duration of minor
epidemics. If the CTMC has rates which are super-linear in the disease
variables, a further approximation of ignoring the higher power terms in
$\bm{{i}}$ is necessary. At the end, this results in assuming that the
matrices $(F,V)$ are constant (do not depend on $\bm{{i}}$).

Let us illustrate this philosophy on the famous SIR example. However, in line
with our interest in this paper and getting a bit ahead of ourselves, we will
only look at an "disease process" of the infected, with the other components
fixed. The state space of the process will thus be $\mathbb{N}$. We note this
is similar in spirit with the "slow-fast/singular perturbation" technique of
considering only variables whose lifetime is short, and fixing the other
variables whose lifetime is longer, and is in fact the idea behind the famous
next generation matrix approach.

###### Example 2

The "SIR" disease process (i.e. defined on the disease compartments) is
$i^{\prime}=(\beta s-\gamma)i$. The natural SIR/linear CT birth and death
disease stochastic process (DSP) is a Markov process $X_{t}\in\mathbb{N}$ with
generating operator on the set of functions $f:\mathbb{N}->\mathbb{N}$ defined
by

$\mathcal{G}f(i)=\beta si(f(i+1)-f(i))+\gamma i(f(i-1)-f(i))=Af(i),$ (4)

and corresponding semi-infinite generator matrix

$A=\begin{pmatrix}-\beta&\beta&0&\cdots&0\\\
\gamma&-\beta-\gamma&\beta&\cdots&0\\\ 0&\gamma&-\beta-\gamma&\beta&\ddots\\\
\vdots&\ddots&\ddots&\ddots&\beta\\\
0&0&\cdots&\gamma&-\beta-\gamma\end{pmatrix}.$ (5)

###### Remark 8

We recall for the benefit of readers who have not been exposed to the
(immense) literature on Markov processes that the behavior of expectations of
this class of stochastic processes always involves one deterministic operator
$A$, the generator of the Markovian semigroup, which acts on a space of
"appropriate functions" on the state space (4), and where "appropriate" may be
skipped in simple cases like ours (5). The essential thing to notice here is
that our Markov generator operator $A$ is completely defined by the rates,
just like its "mean-field" deterministic ODE. Thus, from the practical point
of view of estimating rates, we have added nothing to the parameters of the
ODE model (as would be the case with other stochastic processes involving
Brownian motion, etc). We have only modified the state space and the operator;
however, this way, phenomena which are invisible in the continuous mean-field
limit, become relevant.

Finally, for readers puzzled by the question of where is the randomness hidden
in the deterministic operator (4), we mention that this arrives via two
Poisson processes describing the times when the process jumps up and down,
respectively, and we refer to the literature for more details.

This process converges to $\infty$ (i.e. is non-recurrent) or to a stationary
distribution iff $R_{0}:=\frac{\beta s}{\gamma}$ is strictly bigger than $1$,
or strictly smaller than $1$, respectively. The probability of
"extinction/absorbtion into $0$", when starting the process with $j$ infected
are

$p(j)=q^{j},\quad q=\begin{cases}1&R_{0}<1\\\ \frac{\gamma}{\beta
s}=\frac{1}{R_{0}}&R_{0}\geq 1\end{cases}.$ (6)

This result may be found for example in the textbook [19] (it is, up to
technical difficulties caused by the non-compact state space, the simplest
illustration of the fact that solutions of "Dirichlet problems" of the form
$p(j)=E_{X_{0}=j}[g(X_{\tau})],$ where $\tau$ is the exit time from a domain
$S$, must solve $\mathcal{G}p=0$ and $p=g$ on the boundary of $S$).

The expected time to extinction when starting the process with is $j$ infected
and when $R_{0}<1$ may be found using the fact that solutions of "Poisson
problems" of the form $T(j)=E_{X_{0}=j}[\int_{0}^{\tau}h(X(s)ds],$ must solve

$\begin{cases}\mathcal{G}T+h=0\\\ T=0\text{ on the boundary of }S\end{cases}.$

Another interesting quantity is the expected time to extinction when
$R_{0}>1,$ in the case that extinction occurs. This "Dirichlet-Poisson
problem" may be written as

$T(j)=E_{X_{0}=j}[g(X_{\tau})\int_{0}^{\tau}h(X(s)ds],$

where $h=1$, and $g$ is the indicator of extinction occurring. Such
expectations must solve

$\begin{cases}\mathcal{G}T+hp=0\\\ T=0\text{ on the boundary of
}S\end{cases},$

where $p$ is the solution of the Dirichlet problem with boundary value $g$.

For SIR, we must solve respectively

$\begin{cases}\beta sx(T(x+1)-T(x))+\gamma
x(T(x-1)-T(x))+1=0,T(0)=0,T(K)=0,K->\infty,&\text{when }R_{0}<1\\\ \beta
sx(T(x+1)-T(x))+\gamma
x(T(x-1)-T(x))+q^{x}=0,T(0)=0,T(K)=0,K->\infty,&\text{when }R_{0}\geq
1\end{cases}.$

(7)

These two equations may be solved explicitly. The limits are quite challenging
even with Mathematica, as shown in the Appendix 10. We are able to recover and
generalize the results of [57] (see also [53, (10)]) when $j\geq 1$ for the
first problem, but not for the second one.

The Bacaer equation. One missing aspect in the previous works however is
characterizing the extinction probabilities via one final equation, without
going through the discretization procedure employed in [8, 9, 58, 53, 44], and
solving each example individually. Interestingly, such an equation in terms of
$(F,V)$ decompositions was provided by Griffiths in [31], except that this
paper considers only BD’s with no transfers.

We review now the work of [14] (who were motivated by analyzing the case of
periodic steady solutions), but on the way spelled out also the simple
equation (8) below. To each fixed values for the disease variables, one may
associate to a $(F,V)$ decomposition $\bm{{i}}^{\prime}=(F-V)\bm{{i}}$ a
"multi-dimensional birth and death process" (BD), with birth rates given by
$F$, and with transfer and death rates given by $-V$.
444$\bm{{i}}^{\prime}=(F-V)\bm{{i}}$ are precisely the mean-field equation for
the multi-dimensional birth and death process; this is precisely [31, (6)],
under the extra condition that we assume, that the immigration vector into the
disease compartments is $0$. In fact, the $-V$ matrix by itself generates an
absorbing CTMC (and the $F$ matrix models roughly inputs to be fed into this
absorbing CTMC). This observation explains that an ODE mathematical
epidemiology model has associated to it both a birth and death process, as
well as a "death and transfer only" absorbtion CTMC –see Remark 22 for an
example. If furthermore $(F,V)$ are independent of $\bm{{i}},$ we are dealing
with a branching process (approximation).

A useful fact to recall is that the probabilities of extinction of a multi-
variate discrete time branching process are of the form

$\mathbb{P}_{0}=q_{1}^{j_{1}}q_{2}^{j_{2}}...,$

where $q=(q_{j},j=1,...,J)$, and $J$ is the number of disease compartments,
and where $q_{j}$ satisfy the "Bacaer equation"

$(q^{t}\circ F)*q+(1-q)\circ V-q*f=0\Leftrightarrow
q_{j}=\frac{\sum_{k=1}^{J}(1-q_{k})V_{kj}}{f_{j}-\sum_{k=1}^{J}q_{k}F_{kj}},$
(8)

where $*$ denotes coordinate-wise product, dot product is denoted by $\circ$,
and $f_{j}=\sum_{k}F_{k,j}$. This equation is new, but it may be inferred from
[31, (9)], [14, (11)] and [45, (5.3)] (after some changes of variables).

For the SIR process for example, (8) becomes

$(q-1)\beta sq+(1-q)\gamma=(q-1)(\beta sq-\gamma)=0,$

with the two roots $q=1$ and $q=\frac{\gamma}{\beta s}=\frac{1}{R_{0}}$,
recovering Whittle’s result (the two roots yield the correct result when
$R_{0}$ is strictly smaller than $1$ and strictly bigger than $1$,
respectively).

We will check below that (8) also recovers other explicit particular cases
offered in the literature, like SEIR [8], [9, sec.4], SIV [43], etc

###  2.5 The Jacobian factorization bound

Note first the following elementary fact:

###### Lemma 1

A sufficient (but not necessary) condition for a polynomial with real
coefficients and positive leading term to admit a positive root is that
$c_{0}<0$, where $c_{0}$ is the constant term of the polynomial.

For polynomials of degree $1$, this condition is also necessary. This converse
result may be strenghtened to “Descartes type polynomials".

###### Definition 3

We will say that a parametric polynomial with real coefficients, whose
constant coefficient may change sign, but whose all other coefficients are
“sign definite", and of the same sign (which may be supposed w.l.o.g. to be
$+$), is of Descartes type.

As an immediate consequence of Descartes’s rule of signs, it follows that

###### Lemma 2

A sufficient and necessary condition for a Descartes polynomial with a
positive leading term to admit a positive root is that $c_{0}<0$, where
$c_{0}$ is the constant term of the polynomial.

###### Remark 9

Note the immense simplification with respect to the Routh-Hurwitz conditions,
when we need to establish the existence of a positive root, for a Descartes
type polynomial .

We believe that “the mystery of the success of the Jacobian factorization
approach" comes from the fact that “simple epidemic models" often feature
Descartes type polynomials. However, this leaves us with many further
questions, like when does this happen and what to do when it doesn’t.

The Jacobian factorization approach consists in:

  1. 1.

Putting all the rational factors of the characteristic polynomial of the
Jacobian, in a form normalized to have positive leading term, assuming they
are sign definite (if this is not the case, this approach does not work, but
may be generalized).

  2. 2.

Removing all linear factors with eigenvalues which are negative.

  3. 3.

For all remaining factors $F_{i}$ for which $c_{0}^{(i)}<0$ may hold for
certain parameter values, rewrite this inequality into the form

$c_{0}^{(i)}=c_{+}-c_{-}=c_{+}(1-R_{J}^{(i)})<0\Leftrightarrow
R_{J}^{(i)}:=\frac{c_{-}}{c_{+}}>1,$

where $c_{+},c_{-}$ are the positive and negative parts of the expanded form
of $c_{0}^{(i)}.$

  4. 4.

Define the "Jacobian factorization $R_{0}$"

$R_{J}=\max_{i}[R_{J}^{(i)}].$ (9)

###### Theorem 1

A) In the instability domain, $R_{J}$ is a lower bound for $\inf_{F\text{
admissible}}R_{F}$.

B) In the stability domain, $R_{J}$ is an upper bound for $\sup_{F\text{
admissible}}R_{F}$.

Proof: A) Fix any admissible $F$ and let $R_{F}$ be its associated NGM
$R_{0}$. Then

$R_{J}>1\Leftrightarrow\exists i:R_{J}^{(i)}>1\Leftrightarrow\exists
i:c_{0}^{(i)}<0\Longrightarrow\text{ DFE instability }\Leftrightarrow
R_{F}>1.$

Thus

$R_{J}>1\Longrightarrow R_{F}>1\Leftrightarrow R_{J}\leq R_{F},$ (10)

and the result follows.

B) Similar proof.

Conjecture: We conjecture that if all the factors $F_{i}$ are Descartes
polynomials, then $R_{J}=R_{F}$ for any admissible decomposition $(F,V)$, and
will denote the resulting object by $R_{0}$.

Open question 1: Under what conditions do our NGM $R_{0}$ and our Jacobian
$R_{J}$ coincide?

The implementation of the Jacobian factorization approach is provided in
Section 10.

###  2.6 The “rational univariate representation" (RUR) and the reduced order
quasi -stationary approximation

Hundreds of mathematical epidemiology papers have already employed the idea of
reducing the fixed point system to one scalar equation in one of the disease
variables, via rational substitutions for the other variables. We note that
this is a particular case of the so-called "rational univariate
representation" (RUR), but for Mathematica users this is irrelevant, since RUR
is not implemented currently, and we had to write our own script, included
below, in which the user chooses in a system the variable he wants to restrict
to.

The current code for this reduction to one equation algorithm is

    RUR[mod_, ind_, cn_ : {}] := Module[{dyn, X, par, eq, elim},
          dyn = mod[[1]]; X = mod[[2]]; par = mod[[3]];
          elim = Complement[Range[Length[X]],ind];
          eq = Thread[dyn == 0];
          ratsub = Solve[Drop[eq, ind], X[[elim]]][[1]];
          pol =
           Collect[GroebnerBasis[Numerator[Together[dyn /. cn]],
             Join[par, X[[ind]]], X[[elim]]], X[[ind]]];
              {ratsub, pol}
        ]

###### Remark 10

The command which does the essential work is ”GroebnerBasis". When ”ind" is a
set with just one component, this reduces the system to a polynomial in this
variable. Alternatively, this could be achieved by plugging the results of
”ratsub" into the system.

The script above works directly for models with demographics, but must be
modified for “conservation systems", where the fixed points are only
determined by adding the total mass conservation equation to the fixed point
equations.

This script may also be used for order reduction, in the spirit of the (
quasi-steady-state assumption) QSSA method in biochemistry, and of the recent
epidemiology paper [37]. We illustrate this for the simplest SIR example.

###### Example 3

For the SIR process $(S(t),I(t),R(t),t\geq 0)$ with linear birth rates
$b_{s},b_{r}$ for susceptible and recovered, the system for the fractions
$\;\mathsf{s}(t)=\frac{S(t)}{N},\;\mathsf{i}(t)=\frac{I(t)}{N},\;\mathsf{r}(t)=\frac{R(t)}{N},N=S+I+R$
is:

$\begin{cases}\;\mathsf{s}^{\prime}(t)=b_{s}-\beta\;\mathsf{s}(t)\;\mathsf{i}(t)+\gamma_{r}\;\mathsf{r}(t)-d_{s}\;\mathsf{s}(t),&d_{s}=\gamma_{s}+\mu\\\
\;\mathsf{i}^{\prime}(t)=\beta\;\mathsf{s}(t)\;\mathsf{i}(t)-d_{i}\;\mathsf{i}(t),&d_{i}=\gamma_{i}+\mu+\delta\\\
\;\mathsf{r}^{\prime}(t)=b_{r}+\gamma_{i}\;\mathsf{i}(t)+\gamma_{s}\;\mathsf{s}(t)-d_{r}\;\mathsf{r}(t),&d_{r}=\gamma_{r}+\mu\end{cases}$
(11)

The DFE is:
$(\frac{b_{r}\gamma_{r}+b_{s}\left(\mu+\gamma_{r}\right)}{\mu\left(\mu+\gamma_{r}+\gamma_{s}\right)},0,\frac{b_{r}\left(\mu+\gamma_{s}\right)+b_{s}\gamma_{s}}{\mu\left(\mu+\gamma_{r}+\gamma_{s}\right)})$.

The rational substitution with respect to $i$ obtained via RUR is:

$\left(s\to\frac{\gamma_{r}\left(b_{r}+b_{s}+i\gamma_{i}\right)+\mu
b_{s}}{\beta
i\left(\mu+\gamma_{r}\right)+\mu\left(\mu+\gamma_{r}+\gamma_{s}\right)},r\to\frac{b_{r}\left(\beta
i+\mu+\gamma_{s}\right)+b_{s}\gamma_{s}+i\gamma_{i}\left(\beta
i+\mu+\gamma_{s}\right)}{\beta
i\left(\mu+\gamma_{r}\right)+\mu\left(\mu+\gamma_{r}+\gamma_{s}\right)}.\right)$

Note this reduces to the DFE when $i=0$.

The reduced approximate model obtained via RUR is:

$\displaystyle
i^{\prime}=i\left[a_{0}-a_{1}i\right],a_{1}=\beta\left(\mu\gamma_{i}+(\delta+\mu)\left(\mu+\gamma_{r}\right)\right),$
$\displaystyle
a_{0}=\beta\left(b_{r}\gamma_{r}+b_{s}\left(\mu+\gamma_{r}\right)\right)-\mu\left(\mu+\gamma_{r}+\gamma_{s}\right)\left(\delta+\gamma_{i}+\mu\right)$
$\displaystyle=\mu\left(\mu+\gamma_{r}+\gamma_{s}\right)\left(\delta+\gamma_{i}+\mu\right)\left(\;\mathsf{s}_{dfe}\mathcal{R}-1\right)=\mu\left(\mu+\gamma_{r}+\gamma_{s}\right)\left(\delta+\gamma_{i}+\mu\right)\left(R_{0}-1\right).$

This has an explicit (rather formidable) analytic solution, provided in the
Mathematica file.

Figure 1: Illustration of the asymptotic convergence of $i(t)$ towards the
endemic value $\frac{a_{0}}{a_{1}}=.092$, both for the full SIR model and its
approximation.

One may notice that for the chosen numerical illustration, the plots of $i$
and its approximation converge towards the same value, but differ sharply for
the chosen numeric values as far as shape is concerned.

We mention finally the possibility to develop yet another possible algorithm
for computing a “bifurcation $R_{0}$", suggested by the example above, which
is based on the known fact that this parameter is expected to produce
bifurcations at $R_{0}=1$.

The steps are:

  1. 1.

Factor out the variable in the scalar polynomial of the reduced model (always
possible iff this is a disease variable).

  2. 2.

Write the free coefficient of the divided polynomial as
$F(R_{0})=G(R_{0})(R_{0}-1)$, where $F(R)$ is rational (always possible due to
the known bifurcation at $R_{0}=1$).

  3. 3.

Identify a factor which is linear in susceptible variables like
$\;\mathsf{s}_{dfe}$, etc, and write it as a difference of positive and
negative terms. Upon normalizing one of them to one, the other will be
$R_{0}$, or $1/R_{0}$.

##  3 $R_{0}$ and extinction probabilities for the SEIR epidemic model

The SEIR process $(S(t),E(t),I(t),R(t),t\geq 0)$ adds to the SIR model the
class $E$(exposed). The model for the fractions
$\;\mathsf{s}(t)=\frac{S(t)}{N},\;\mathsf{e}(t)=\frac{E(t)}{N},\;\mathsf{i}(t)=\frac{I(t)}{N},\;\mathsf{r}(t)=\frac{R(t)}{N},N=S+E+I+R$
is:

$\begin{cases}\;\mathsf{s}^{\prime}(t)=b_{s}-\beta\;\mathsf{s}(t)\;\mathsf{i}(t)+\gamma_{r}\;\mathsf{r}(t)-d_{s}\;\mathsf{s}(t),&d_{s}=\gamma_{s}+\mu\\\
\;\mathsf{e}^{\prime}(t)=\beta\;\mathsf{s}(t)\;\mathsf{i}(t)-\gamma_{e}\;\mathsf{e}(t)-d_{e}\;\mathsf{e}(t),&d_{e}=\gamma_{e}+\mu\\\
\;\mathsf{i}^{\prime}(t)=\gamma_{e}\;\mathsf{e}(t)-d_{i}\;\mathsf{i}(t),&d_{i}=\gamma_{i}+\mu+\delta\\\
\;\mathsf{r}^{\prime}(t)=b_{r}+\gamma_{i}\;\mathsf{i}(t)+\gamma_{s}\;\mathsf{s}(t)-d_{r}\;\mathsf{r}(t),&d_{r}=\gamma_{r}+\mu\end{cases}$

This is both a textbook model, and one for which answers to many open
questions (concerning for example the emergence of chaos under stochastic and
periodic perturbations) are still awaited – see for example [25, 50, 30].

The DFE of (3) is:
$(\frac{b_{r}\gamma_{r}+b_{s}\left(\mu+\gamma_{r}\right)}{\mu\left(\mu+\gamma_{r}+\gamma_{s}\right)},0,0,\frac{b_{r}\left(\mu+\gamma_{s}\right)+b_{s}\gamma_{s}}{\mu\left(\mu+\gamma_{r}+\gamma_{s}\right)})$.
The decomposition matrices and basic reproduction number are:

$F=\left(\begin{array}[]{cc}0&\beta s\\\ 0&0\\\
\end{array}\right),V=\left(\begin{array}[]{cc}\gamma_{e}+\mu&0\\\
-\gamma_{e}&\delta+\gamma_{i}+\mu\\\ \end{array}\right),R_{0}=\frac{\beta
s\gamma_{e}}{\left(\gamma_{e}+\mu\right)\left(\delta+\gamma_{i}+\mu\right)}.$

The associated disease stochastic process ${X=(e,i)}\in\mathbb{N}^{2}$ has
generating operator

$\scalebox{0.9}{$\mathcal{G}=\beta
si(f(x+e_{1})-f(x))+\gamma_{e}e(f(x+tr)-f(x))+d_{e}e(f(x+e_{3})-f(x))+d_{i}i(f(x+e_{2})-f(x))$},$

where $x=(e,i),e_{1}=(1,0)=-e_{3},tr=(-1,1),e_{2}=(0,-1)$.

The extinction probabilities obtained by solving (8) are

$\begin{cases}q_{i}=1,q_{e}=1&\text{when }R_{0}<1\\\
q_{i}=\frac{1}{R_{0}},q_{e}=\frac{\mu}{d_{e}}+\frac{\gamma_{e}}{d_{e}}\frac{1}{R_{0}}&\text{when
}R_{0}\geq 1.\end{cases}$

This checks with the particular case in [8], where $\;\mathsf{s}_{dfe}=1$.

###### Remark 11

It is not clear intuitively why separating the transition rates into those of
$F$ (which increase the norm of $x$) and those of $V$ (which do not increase
the norm of $x$) should matter for determining the extinction probabilities,
as happens in (8). This seems an interesting question.

##  4 Rank one host-only models with pathogen, and $R_{0}$ readable from the
flow-chart

###  4.1 The SEIARW model with “catalyzing pathogen" of [32] has rank one next
generation matrix and $R_{0}=R_{J}$

[32] attempted to offer a “ definition-based method" for “computing $R_{0}$ of
dynamic models of single host species, which is mutually coherent with the
next-generation method (NGM)" (and somewhat unclear for "computing $R_{0}$ for
a population with multi-group models"). Unfortunately, these authors do not
seem aware of the fact that all the single host species they examined have
next generation matrix of rank one, and that in this case there exists a
simple general formula [10, 11, 12], which is also related to the definition-
based method of [20].

We review now the SEIAR model (susceptibles, exposed, infected, asymptomatic
and recovered), to which [32] add also a pathogen compartment $W$, resulting
in the SEIARW model.

$\begin{cases}\begin{cases}e^{\prime}=s\left(a\beta_{a}+i\beta_{i}+w\beta_{w}\right)-ed_{e},&d_{e}=e_{i}+e_{a}+\mu\\\
i^{\prime}=ee_{i}-id_{i},&d_{i}=\gamma_{i}+\mu+\delta,\\\
a^{\prime}=ee_{a}-ad_{a},&d_{a}=\gamma_{a}+\mu,\\\
w^{\prime}=a\varepsilon_{a}+i\varepsilon_{i}-wd_{w}\\\
r^{\prime}=a\gamma_{a}+\gamma_{i}i-\mu r\end{cases}\\\
s^{\prime}=\Lambda-s\left(a\beta_{a}+i\beta_{i}+w\beta_{w}+\mu\right).\end{cases}$
(12)

In matrix form, the disease equations are:

$\displaystyle\ \begin{pmatrix}e^{\prime}\\\ i^{\prime}\\\ a^{\prime}\\\
w^{\prime}\\\
r^{\prime}\end{pmatrix}=\begin{pmatrix}-d_{e}&s\beta_{i}&s\beta_{a}&s\beta_{w}&0\\\
e_{i}&-d_{i}&0&0&0\\\ e_{a}&0&-d_{a}&0&0\\\
0&\varepsilon_{i}&\varepsilon_{a}&-d_{w}&0\\\
0&\gamma_{i}&\gamma_{a}&0&-\mu\end{pmatrix}\begin{pmatrix}e\\\ i\\\ a\\\ w\\\
r\end{pmatrix}$

In the absence of pathogen, SEIAR is a rank one "generalized stage-structured
infectious disease model" as revealed by its
${F=}\left(\begin{array}[]{ccc}0&s\beta_{i}&s\beta_{a}\\\ 0&0&0\\\ 0&0&0\\\
\end{array}\right)$ and by its $V=\left(\begin{array}[]{ccc}d_{e}&0&0\\\
-e_{i}&d_{i}&0\\\ -e_{a}&0&d_{a}\\\ \end{array}\right)$ matrix, which is
triangular (compare to [9, Sec. 3]).

The $R_{0}$ has a very intuitive, and easily explainable form:

$R_{0}=\frac{\;\mathsf{s}_{dfe}}{d_{e}}\left[\beta_{i}\frac{e_{i}}{d_{i}}+{\beta_{a}}\frac{d_{a}}{e_{a}}\right]$
(13)

(compare to [9, Sec. 3], to see the general pattern for more stages).

###### Remark 12

Note that this result may be obtained with $\mathcal{I}=\\{e,i,a,r\\}$, and
also with $\mathcal{I}=\\{e,i,a\\},$ which raises the question of defining the
concept of minimal or "sufficient disease" set, in such a way that it allows
deriving both $R_{0}$ and the extinction probabilities.

After the addition of the catalyzing pathogen, the SEIARW is still a rank one
"generalized stage-structured infectious disease model", but the $R_{0}$ is
less intuitive

$R_{0}=\frac{\;\mathsf{s}_{dfe}}{d_{e}}\left[\beta_{i}\frac{e_{i}}{d_{i}}+{\beta_{a}}\frac{d_{a}}{e_{a}}+\frac{\beta_{w}}{d_{w}}\left(\varepsilon_{i}\frac{e_{i}}{d_{i}}+\varepsilon_{a}\frac{e_{a}}{d_{a}}\right)\right];$
(14)

still, it may be read out of the flow chart “almost by inspection" (see also
[32] for an algorithm computing this).

Figure 2: Flow chart corresponding to the SEIARW model (12).

Despite the fact that the characteristic polynomial is not of Descartes type,
all our three $R_{0}$ recipes yield the above result. We provide now details
for the NGM method. After removing the compartment $r$ (since it does not
appear in the other equations), the calls “inf=Range[4];DFE[SEIARW,inf];
NGM[SEIARW,inf]" of our scripts yield that the DFE is

$\left\\{s\to\frac{\Lambda}{\mu},e\to 0,i\to 0,a\to 0,w\to 0\right\\}$

and

$F=s\left(\begin{array}[]{cccc}0&\beta_{i}&\beta_{a}&\beta_{w}\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ \end{array}\right)=s\begin{pmatrix}1\\\ 0\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}0&\beta_{i}&\beta_{a}&\beta_{w}\end{pmatrix},V=\left(\begin{array}[]{cccc}d_{e}&0&0&0\\\
-e_{i}&d_{i}&0&0\\\ -e_{a}&0&d_{a}&0\\\
0&-\varepsilon_{i}&-\varepsilon_{a}&d_{w}\\\ \end{array}\right).$ (15)

Here the dominant eigenvalue of $K=FV^{-1}$, that of the transpose

$K^{t}=s\left(\begin{array}[]{cccc}\frac{\beta_{i}d_{w}e_{i}d_{a}+\beta_{a}e_{a}d_{w}d_{i}+\beta_{w}\left(e_{i}\varepsilon_{i}d_{a}+e_{a}\varepsilon_{a}d_{i}\right)}{d_{e}d_{i}d_{a}d_{w}}&\frac{\beta_{i}d_{w}+\varepsilon_{i}\beta_{w}}{d_{w}d_{i}}&\frac{\beta_{a}d_{w}+\varepsilon_{a}\beta_{w}}{d_{w}d_{a}}&\frac{\beta_{w}}{d_{w}}\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ \end{array}\right),$

the rank 1 formula of [10, 11, 12]

$(1,0,0,0).(V^{t})^{-1}.(0,\beta_{i},\beta_{a},\beta_{w})^{t},$

as well as the Jacobian factorization confirm all the result (14) of [32].

##  5 Target-infection-virus models

###  5.1 Two admisible $(F,V)$ decompositions and $R_{0}$’s for the three
dimensional model of [48]

The three dimensional model of [48, 1] is:

$\begin{cases}x^{\prime}=\Lambda(x)-\beta
xv-\beta_{xy}xy,\quad\Lambda(x)=\mu_{x}(x_{dfe}-x)\\\
\begin{pmatrix}y^{\prime}\\\ v^{\prime}\end{pmatrix}=\begin{pmatrix}\beta
xv+\beta_{xy}xy\\\ 0\end{pmatrix}-\begin{pmatrix}\mu_{y}y\\\
\mu_{v}v+\beta_{xv}xv+\beta_{yv}yv-b\mu_{y}y\end{pmatrix},\end{cases}$ (16)

where we represented already the infectious equations as a difference of "new
(positive) infection" terms and "transfers". The DFE is $x=x_{dfe},y=0,v=0$.

This reduces to the case with 0-delays in [60, 5.1] when the rate of viruses
moving into a healthy cell $\beta_{xv}$ and the rate of viruses moving into an
infected cell $\beta_{yv}$ are both $0$, and to the case in [2] when
$\beta_{yv}=0=\beta_{xy}$ (the latter is the cell to cell infection rate), and
$\beta_{xv}=\beta$.

The gradient of the infectious equations is

$M=\left(\begin{array}[]{cc}x\beta_{xy}-\mu_{y}&\beta x\\\
b\mu_{y}-v\beta_{yv}&-\mu_{v}-x\beta_{xv}-y\beta_{yv}\end{array}\right)$ (17)

Calling our NGM script with “inf={2,3}" yield’s [48]’s result: namely,

$F=x_{dfe}\left(\begin{array}[]{cc}\beta_{xy}&\beta\\\ 0&0\\\
\end{array}\right),-V=\left(\begin{array}[]{cc}-\mu_{y}&0\\\
b\mu_{y}&-\mu_{v}-x_{dfe}\beta_{xv}\end{array}\right),$

the next generation matrix (NGM) of the infectious coordinates at the DFE

$K=\left(\begin{array}[]{cc}\frac{\beta
bx}{\mu_{v}+x\beta_{xv}}+\frac{x\beta_{xy}}{\mu_{y}}&\frac{\beta
x}{\mu_{v}+x\beta_{xv}}\\\ 0&0\\\ \end{array}\right),$

and that the DFE is Lyapunov-Malkin stable when $R_{0}$ defined in

$R_{0}=\frac{x_{dfe}}{x_{c}}+\beta
b\frac{x_{dfe}}{\mu_{v}+x_{{dfe}}\beta_{xv}},\;\;x_{c}:=\frac{\mu_{y}}{\beta_{xy}},$
(18)

is smaller than $1$, and unstable when $R_{0}>1$.

The Jacobian factorization provides the same formula, despite the fact that
the characteristic polynomial is of Descartes type only conditionally, when
$\beta_{xv}\geq\beta_{xy}$.

###### Remark 13

Interestingly, another admissible decomposition
$\mathcal{F}=\begin{pmatrix}\beta xv+\beta_{xy}xy\\\
beta\mu_{y}y\end{pmatrix}$, appears in an earlier version of [48] at
https://people.clas.ufl.edu/pilyugin/files/cosner60-dcdsB.pdf :

$F=\left(\begin{array}[]{cc}x\beta_{xy}&\beta x\\\ b\mu_{y}&0\\\
\end{array}\right),\quad V=\left(\begin{array}[]{cc}\mu_{y}&0\\\
0&\mu_{v}+\beta_{xv}x_{dfe}\\\
\end{array}\right),K=\left(\begin{array}[]{cc}\frac{x_{dfe}}{x_{c}}&\frac{\beta
x_{dfe}}{\mu_{v}+\beta_{xv}x_{dfe}}\\\ b&0\\\ \end{array}\right)$ (19)

This second decomposition yields a different $R_{0}$:

$R_{0}^{\prime}=\frac{x_{dfe}}{2x_{c}}\left(1+\sqrt{1+\frac{4\beta
bx_{c}^{2}}{x_{dfe}(x_{dfe}\beta_{xv}+\mu_{v})}}\right).$ (20)

Furthermore, this early version also shows that the two decompositions have
the same stability domain for the DFE, which may be reexpressed as

$R_{0}=K_{1,1}+K_{1,2}K_{2,1}=\frac{x_{dfe}}{x_{c}}+b\beta\frac{x_{dfe}}{\mu_{v}+\beta_{xv}x_{dfe}}<1.$
(21)

We note that this equivalence follows also by applying the first criterion in
[59] (when the characteristic polynomial, given here by
$\lambda^{n}-a_{1}\lambda^{n-1}-a_{2}\lambda^{n-1}-...$ has all coefficients
non-negative, than $\sum_{i}a_{i}$ may be used as threshold parameter instead
of $R_{0}$), with $n=3,a_{1}=K_{1,1},a_{2}=K_{1,2}K_{2,1},a_{3}=0$.

###### Remark 14

Note the second decomposition has one more non-zero term in $F$, which does
not appear in ours, since we view it as a transfer and not as an interaction.
We see here an excellent example of non-uniqueness, where one must chose
between an answer with $F$ of lower rank and a simpler $R_{0}$ formula, but
which is valid only under certain conditions (that the non-diagonal term
$b\mu_{y}$ in $V$ is small enough), and an answer with simpler $V$, which
requires less assumptions on the parameters, but yields a more complicated
$R_{0}$.

###### Remark 15

The domain of stability, in terms of the parameters. As an aside, it is easy
to show that (21) is equivalent to

$x_{dfe}<x_{c},b<b_{0}:=\frac{\mu_{v}+\beta_{xv}x_{dfe}}{\beta
x_{dfe}}\left(1-\frac{x_{dfe}}{x_{c}}\right),$ (22)

where $b_{0}$ is the solution of the equation $R_{0}(b)=1$. Thus, stability of
the DFE is equivalent to both $x_{dfe}$ and the “burst parameter" $b$ being
small enough.

We offer now a third gradient decomposition, which turns out to be
inadmissible sometimes, but yields again our recipe’s $R_{0}$. Taking
$F=\left(\begin{array}[]{cc}x_{dfe}\beta_{xy}&0\\\ b\mu_{y}&0\\\
\end{array}\right)$ yields

$V=F-M=\left(\begin{array}[]{cc}\mu_{y}&-\beta x_{dfe}\\\
0&\mu_{v}+x_{dfe}\beta_{xv}\\\
\end{array}\right),V^{-1}=\left(\begin{array}[]{cc}\frac{1}{\mu_{y}}&\frac{\beta
x_{dfe}}{\mu_{y}(\mu_{v}+x_{dfe}\beta_{\text{xv}})}\\\
0&\frac{1}{\mu_{v}+x_{dfe}\beta_{\text{xv}}}\\\ \end{array}\right).$

Note that $V$ is a subgenerating matrix only if
$x_{dfe}\leq\frac{\mu_{y}}{\beta}$.

However
$K=\left(\begin{array}[]{cc}\frac{x_{dfe}\beta_{xy}}{\mu_{y}}&\frac{\beta
x_{dfe}^{2}\beta_{xy}}{\mu_{v}\mu_{y}+x_{dfe}\beta_{xv}\mu_{y}}\\\
b&\frac{\beta bx_{dfe}}{\mu_{v}+x_{dfe}\beta_{xv}}\\\ \end{array}\right)$
yields the correct

$R_{0}=\max\left(0,\frac{\beta
bx_{dfe}}{\mu_{v}+x_{dfe}\beta_{xv}}+\frac{x_{dfe}\beta_{xy}}{\mu_{y}}\right).$

In the current example, the RUR algorithm works as well. The difference of the
two positive terms is

$\beta
bx_{dfe}\mu_{y}+x_{dfe}\beta_{{xy}}(\mu_{v}+x_{dfe}\beta_{{xv}})-\mu_{y}(\mu_{v}+x_{dfe}\beta_{{xv}})=\mu_{y}(\mu_{v}+x_{dfe}\beta_{{xv}})(R_{0}-1)$

for both choices $y$ and $v$ as scalar variable, and the appropriate cosmetics
recovers the recipe NGM $R_{0}$.

This example illustrates the fact that sometimes several admissible and even
conditionally non-admissible decompositions, as well as other approaches, may
lead to the same $R_{0}$.

###  5.2 Two distinct approximate extinction probabilities, one for each
admisible $(F,V)$ decomposition for the model of [48]

The extinction probabilities of the stochastic model are of course unique. We
may use the result of Bacaer’s formula as approximations. In this interesting
example, we find out that both $(F,V)$ decompositions yield reasonable
results. This suggests that we have not one, but two deterministic
epidemiologic approximations for a single stochastic model. This strengthens
our point of view that a deterministic epidemiologic model must include a
specification of the $(F,V)$ decomposition.

The respective results we got are:

  1. 1.

For the first decomposition, the extinction probabilities obtained by solving
(8) are

  2. 2.

For the second decomposition, the extinction probabilities obtained by solving
(8) are:

In a numeric instance, we found the two results reasonably close to each
other.

##  6 Multi-strain host only models

Multi-strain diseases are diseases that consist of several strains, or
serotypes. One interesting thing about multi-strain models is that besides the
DFE we have new boundary points which are relevant epidemiologically, in which
one subset of strains $A$ is present (“resident"). We have then a natural
coexistence of several “$\mathcal{R}$ thresholds":

  1. 1.

$R_{A}$ is the bifurcation threshold at which the DFE stops being stable, when
the only compartments present are those of $A$.

  2. 2.

${\mathcal{R}_{A}}$ is the bifurcation threshold at which the boundary point
$E_{A}$ starts existing (in the presence of the $A^{c}$ compartments).

  3. 3.

$R_{A^{c},A}$ is the bifurcation threshold at which the boundary point $E_{A}$
stops being stable, i.e. when the $A^{c}$ compartments invade the $A$
compartments.

Note that for two strains already, we have at least two new thresholds,
$R_{21},R_{12},$ which together with $R_{0}$ and the thresholds $R_{1}$,
$R_{2}$ of the individual strains divide the line into $6$ regions with
different stability properties. Studying the relations between the various
thresholds in parameter space is quite a challenging topic. However, their
calculation is a priori of the same level of difficulty as for the DFE.

###  6.1 The two-strain SIS tuberculosis model of [56, Sec4.4]

The model presented here is a limiting case of that presented in the next
section, obtained when the transition rates $\gamma_{1},\gamma_{2}$ converge
to $\infty$. It also generalizes the two-strain SIS tuberculosis model of [56,
Sec4.4] by allowing cross-infections in both directions

$\displaystyle\begin{cases}i_{1}^{\prime}=i_{1}\left(i_{2}\left(\nu_{2}-\nu_{1}\right)+\beta_{1}s-\sigma_{1}-b\right)=i_{1}\left(i_{2}\left(\nu_{2}-\nu_{1}\right)+\beta_{1}s-d_{1}\right),\\\
i_{2}^{\prime}=i_{2}\left(i_{1}\left(\nu_{1}-\nu_{2}\right)+\beta_{2}s-\sigma_{2}-b\right)=i_{2}\left(i_{1}\left(\nu_{1}-\nu_{2}\right)+\beta_{2}s-d_{2}\right),\\\
s^{\prime}={b}-s\left(\beta_{1}i_{1}+\beta_{2}i_{2}+b\right)+i_{1}\sigma_{1}+i_{2}\sigma_{2},\end{cases}$

where we put $d_{1}=\sigma_{1}-b,d_{2}=\sigma_{2}-b$ in the first two
equations, to simplify their notation (the last equation may be removed, since
$s=1-i_{1}-i_{2}$).

Noting that the first two equations factor yields the following three boundary
steady states, where $\mathbf{x}=(i_{1},i_{2},s)$:

$\displaystyle\mathbf{x}_{0}$ $\displaystyle=$
$\displaystyle\left(0,0,1\right)$ $\displaystyle\mathbf{x}_{1}$
$\displaystyle=$
$\displaystyle\displaystyle{\left(1-\mathcal{R}_{1}^{-1},0,\mathcal{R}_{1}^{-1}\right)}$
$\displaystyle\mathbf{x}_{2}$ $\displaystyle=$
$\displaystyle\displaystyle{\left(0,1-\mathcal{R}_{2}^{-1},\mathcal{R}_{2}^{-1}\right)},$

where we put

$\mathcal{R}_{1}=\displaystyle{\frac{\beta_{1}}{{b+\sigma_{1}}}},\mathcal{R}_{2}=\displaystyle{\frac{\beta_{2}}{{b+\sigma_{2}}}}.$

The _disease free steady state_ $\mathbf{x}_{0}$ exists for all parameter
values, while the _original strain only steady state_ $\mathbf{x}_{1}$ is
physically relevant if and only if $\mathcal{R}_{1}>1$ and the _emerging
strain only steady state_ $\mathbf{x}_{2}$ is physically relevant if and only
if $\mathcal{R}_{2}>1$.

There may also be a fourth non-negative coexistence equilibrium (COE), given
by

$\begin{cases}i_{1}=\frac{\beta_{1}d_{2}-\beta_{2}d_{1}-\left(\nu_{1}-\nu_{2}\right)\left(\beta_{2}-d_{2}\right)}{\left(\nu_{1}-\nu_{2}\right)\left(\beta_{1}-\beta_{2}+\nu_{1}-\nu_{2}\right)}\\\
i_{2}=\frac{d_{1}\left(\beta_{2}-\nu_{1}+\nu_{2}\right)-\beta_{1}\left(d_{2}-\nu_{1}+\nu_{2}\right)}{\left(\nu_{1}-\nu_{2}\right)\left(\beta_{1}-\beta_{2}+\nu_{1}-\nu_{2}\right)}\\\
s=1-i_{1}-i_{2}.\end{cases}$ (23)

Note this depends only on $\nu_{1}-\nu_{2},$ which shows that the case
$\nu_{1}=0$ considered in [56, Sec4.4] is not that restrictive 444However, the
appearance of $\nu_{1}-\nu_{2}$ in the denominator suggests limiting diffusion
phenomena, which may be worth studying in their own right. . In this case, the
COE point simplifies to:

$\begin{cases}i_{1}=\frac{d_{2}\left(d_{1}\left(\mathcal{R}_{1}-\mathcal{R}_{2}\right)+\nu\left(\mathcal{R}_{2}-1\right)\right)}{\nu\left(-d_{1}\mathcal{R}_{1}+d_{2}\mathcal{R}_{2}+\nu\right)}\\\
i_{2}=\frac{d_{1}\left(d_{2}\left(\mathcal{R}_{2}-\mathcal{R}_{1}\right)+\nu(1-\mathcal{R}_{1})\right)}{\nu\left(-d_{1}\mathcal{R}_{1}+d_{2}\mathcal{R}_{2}+\nu\right)}\\\
s=1-i_{1}-i_{2}\end{cases},$ (24)

which is positive iff $\mathcal{R}_{2}>1$ and the following conditions hold

$\begin{cases}\mathcal{R}_{1}>\frac{\nu+\mathcal{R}_{2}d_{2}}{\nu+d_{2}},\;0<\nu<\frac{d_{1}(\mathcal{R}_{2}-\mathcal{R}_{1})}{\mathcal{R}_{2}-1},\;\mbox{or},\\\
\mathcal{R}_{1}<\frac{\nu+\mathcal{R}_{2}d_{2}}{\nu+d_{2}},\left(0<d_{1}<\nu(1-\frac{1}{\mathcal{R}_{2}})\;\mbox{or}\;d_{1}>\nu(1-\frac{1}{\mathcal{R}_{2}}),\nu<\frac{d_{1}(\mathcal{R}_{2}-\mathcal{R}_{1})}{1-\mathcal{R}_{2}}\right).\end{cases}$
(25)

We give now some details of the NGM implementation for the three boundary
points. Recall that the idea is to project the ODE at each boundary point on
the $0$ coordinates (or some subset), while fixing the other coordinates. We
must compute therefore new $(F,V)$ pairs at each boundary point, since the
respective zero coordinates are different.

  1. 1.

At the DFE, the zero coordinates are $\left\\{i_{1},i_{2}\right\\}$, and so
$\mathcal{I}=\\{1,2\\}$.

Our script yields the expected result

$R_{0}=Max\left[\frac{\beta_{2}s_{dfe}}{\sigma_{2}+b},\frac{\beta_{1}s_{dfe}}{\sigma_{1}+b}\right]=Max\left[R_{1},R_{2}\right],R_{i}=\;\mathsf{s}_{dfe}\mathcal{R}_{i}=\mathcal{R}_{i},i=1,2.$

  2. 2.

At $\mathbf{x}_{2}$, $\mathcal{I}=\\{1\\}$, and

$R_{12}=\frac{\mathcal{R}_{1}}{\mathcal{R}_{2}}+\frac{\left(\nu_{2}-\nu_{1}\right)\left(1-\mathcal{R}_{2}^{-1}\right)}{b+\sigma_{1}}.$

When $\nu_{1}=0,\nu_{2}=\nu,$ we recover the result [56, (18)]
$R_{12}=\frac{\mathcal{R}_{1}}{\mathcal{R}_{2}}+\frac{\nu}{b+\sigma_{1}}\left(1-\mathcal{R}_{2}^{-1}\right).$

This implies that stability holds iff $\mathcal{R}_{2}>1$ and
$\mathcal{R}_{1}$ is not too big, more precisely:

$R_{12}<1\Leftrightarrow\mathcal{R}_{1}<\mathcal{R}_{2}+\frac{\nu\left(1-\mathcal{R}_{2}\right)}{b+\sigma_{1}}.$
(26)

For a sanity check, we will derive the stability condition of the point
$\mathbf{x}_{2}$ also by the direct Jacobian approach. The Jacobian at
$\mathbf{x}_{2}$ is

$\left(\begin{array}[]{ccc}\frac{-\beta_{2}\left(b+\nu_{1}-\nu_{2}+\sigma_{1}\right)+\beta_{1}\left(b+\sigma_{2}\right)+\left(\nu_{1}-\nu_{2}\right)\left(b+\sigma_{2}\right)}{\beta_{2}}&0&0\\\
-\frac{\left(\nu_{1}-\nu_{2}\right)\left(b-\beta_{2}+\sigma_{2}\right)}{\beta_{2}}&0&-b+\beta_{2}-\sigma_{2}\\\
\sigma_{1}-\frac{\beta_{1}\left(b+\sigma_{2}\right)}{\beta_{2}}&-b&\sigma_{2}-\beta_{2}\\\
\end{array}\right).$

In the case of [56], the eigenvalues are

$\left\\{-b,-\left(\left(\mathcal{R}_{2}-1\right)\left(b+\sigma_{2}\right)\right),\frac{(b+\sigma_{1})\left(\mathcal{R}_{1}-\mathcal{R}_{2}\right)+\nu\left(\mathcal{R}_{2}-1\right)}{\mathcal{R}_{1}}\right\\}.$

The second eigenvalue is negative iff $\mathcal{R}_{2}>1$, and the third
eigenvalue is negative when

$\left(\mathcal{R}_{1}-\mathcal{R}_{2}\right)+\frac{\nu}{b+\sigma_{1}}\left(\mathcal{R}_{2}-1\right)<0\Leftrightarrow
R_{12}<1\text{ see }\eqref{Vs}.$

  3. 3.

An analog result holds by symmetry at $\mathbf{x}_{1}$, where
$\mathcal{I}=\\{2\\}$, and

$\mathcal{R}_{21}=\frac{\left(\nu_{1}-\nu_{2}\right)\left(\mathcal{R}_{1}-1\right)}{d_{2}\mathcal{R}_{1}}+\frac{\mathcal{R}_{2}}{\mathcal{R}_{1}}.$

We illustrate now via a $i_{1}$ bifurcation diagram that, as natural, when
$\beta_{1}$ is small enough, the $x_{2}$ fixed point is stable, to be replaced
as attractor first by the COE, and finally by the $x_{1}$ fixed point, when
$\beta_{1}$ increases.

Figure 3: $i_{1}$ bifurcation diagram when $\beta_{1}$ varies and
$\nu_{1}=0,\nu_{2}=\nu=3=\beta_{2}=3,b=\sigma_{1}=\sigma_{2}=1,R_{1}=\frac{\beta_{1}}{2},R_{2}=\frac{3}{2},$
so that $x_{2}$ is always positive. Since $R_{0}\geq R_{2}>1$, the DFE is
never stable. Observe the following three regimes: a) until $\beta_{12}=1.5$
defined by equality in
${R_{12}:=\frac{\nu\left(\mathcal{R}_{2}-1\right)}{d_{1}\mathcal{R}_{2}}+\frac{\mathcal{R}_{1}}{\mathcal{R}_{2}}}\leq
1\Leftrightarrow\beta_{12}\leq\frac{\beta_{2}\left(b-\nu+\sigma_{1}\right)}{b+\sigma_{2}}+\nu,$
the only stable solution is $x_{2}$. b) At $\beta_{12}=1.5$, $x_{2}$ becomes
unstable and the coexistence solution becomes nonnegative and stable, until
$\beta_{21}$ defined by
$\mathcal{R}_{21}=\frac{\mathcal{R}_{2}}{\mathcal{R}_{1}}-\frac{\nu\left(\mathcal{R}_{1}-1\right)}{d_{2}\mathcal{R}_{1}}=1\Leftrightarrow\beta_{21}=\frac{\left(b+\sigma_{1}\right)(\beta_{2}+\nu)}{b+\nu+\sigma_{2}}=2.4$.
This is also the first intersection point of the COE and $x_{1}$. For a
numerical check, at $\beta_{1c}=2$, defined by
$\mathcal{R}_{1}=1\Leftrightarrow\beta_{1c}=b+\sigma_{1},$ where the $x_{1}$
solution emerges and is initially unstable, the eigenvalues for the COE are
$\left(-1,-0.333333\pm 0.235702Im\right)$. c) After
$\beta_{1}=\beta_{21}\Leftrightarrow\mathcal{R}_{21}<1$, the $x_{1}$ solution
becomes stable and the COE loses its stability (the latter was checked
numerically). Note that at
$\beta_{2c}=3\Leftrightarrow\mathcal{R}_{1}=R_{12}\Leftrightarrow\beta_{1}=\nu$
there is no stability change: the COE and $x_{1}$ continue to be unstable and
stable, respectively.

(a) $(i_{1},i_{2},s)$-time plot at the point $\beta_{1c}=2$ reveals
convergence towards the COE

(b) $(i_{1},i_{2})-$stream plot

Figure 4: Time and phase-plot at the point $\beta_{1c}=2$ illustrating
convergence towards COE=$\left(i_{1}\to 0.0833333,i_{2}\to 0.166667,s\to
0.75\right)$ .

###  6.2 The minimal disease set of the multi-strain host only dengue model
with antibody-dependent enhancement (ADE) [4]

The ADE (antibody-dependent enhancement) effect, believed to occur for Dengue
and Zika, means that infection with a single serotype is asymptomatic, but
infection with a second serotype may lead to serious illness, accompanied by
greater infectivity. It was first studied mathematically by [29, 49], who
showed that for sufficiently small ADE, the numbers of infectives of each
serotype synchronize, with outbreaks occurring in phase, but when the ADE
increases past a threshold, the system becomes chaotic, and infectives of each
serotype desynchronize (however, certain groupings of the primary and
secondary infectives remain synchronized even in the chaotic regime).
Subsequently, [15] examined the effects of single-strain vaccine campaigns on
the dynamics of an epidemic multi-strain Dengue model. We cite now the
eloquent Dengue description given by these authors:

"What makes modeling the dengue virus so interesting is that it has developed
a sophisticated spreading process. Dengue is known to exhibit as many as four
coexisting serotypes (strains) in a region. Once a person is infected and
recovered from one serotype, they confer life-long immunity from that
serotype. However, the antibodies that the body develops for the first
serotype will not counteract a second infection by a different serotype. In
fact, due to the nature of the disease, the antibodies developed from the
first infection form complexes with the second serotype so that the virus can
enter more cells, increasing viral production. The increased transmission rate
in subsequent infections is known as antibody-dependent enhancement (ADE). ADE
is an alarming evolutionary development in multistrain viruses with respect to
vaccines. An optimal vaccination would need to cover all strains of the
disease at once, or the vaccinations could increase transmission of the
strains not covered. This is particularly dangerous for people who have dengue
because the infections are more severe in individuals who already have dengue
antibodies."

A multi-strain model which adds further compartments allowing for temporary
cross-immunity has been developed in the works of Aguiar, Stollenwerk and Kooi
[5, 4, 6, 52, 3].

In this section we consider a ten compartments asymmetric version of the model
of [4], whose variables, denoted by capital letters, represent:

  1. 1.

$S$ are individuals susceptible to both strains;

  2. 2.

$I_{i}$, for $i,j=1,2$ are individuals infected with strain $i$ and with
temporary cross-immunity to strain $j\neq i$;

  3. 3.

$R_{i}$ are individuals who have recovered from strain $i$, but are not yet
susceptible to the other strain $j$;

  4. 4.

$S_{i}$ are individuals who have recovered from strain $i$, and have become
susceptible to the other strain $j$;

  5. 5.

$Y_{j}=I_{ij}$ are individuals previously infected with strain $i$ and now
immune to it, but reinfected with strain $j,i,j=1,2,i\neq j$;

  6. 6.

$R$, omitted in (27) since they do not feed back to the other components, are
the recovered individuals immune to all the strains.

After denoting by small letters the corresponding proportions, we arrive at:

$\begin{cases}s^{\prime}=\mu-s\left(\beta_{1}i_{1}+\beta_{2}i_{2}+\mu+\beta_{1}y_{1}\phi_{1}+\beta_{2}y_{2}\phi_{2}\right),\\\
i_{1}^{\prime}=\beta_{1}s\left(i_{1}+y_{1}\phi_{1}\right)-i_{1}\left(\gamma_{1}+\mu\right),\\\
r_{1}^{\prime}=\gamma_{1}i_{1}-r_{1}\left(\theta_{1}+\mu\right),\\\
s_{1}^{\prime}=\theta_{1}r_{1}-s_{1}\left(\beta_{2}\alpha_{2}\left(i_{2}+y_{2}\phi_{2}\right)+\mu\right),\\\
{y_{2}^{\prime}}=\beta_{2}\alpha_{2}s_{1}\left(i_{2}+y_{2}\phi_{2}\right)-y_{2}\left(\gamma_{2}+\mu\right),\\\
{i_{2}^{\prime}}=\beta_{2}s\left(i_{2}+y_{2}\phi_{2}\right)-i_{2}\left(\gamma_{2}+\mu\right),\\\
{r_{2}^{\prime}}=\gamma_{2}i_{2}-r_{2}\left(\theta_{2}+\mu\right),\\\
s_{2}^{\prime}=\theta_{2}r_{2}-s_{2}\left(\beta_{1}\alpha_{1}\left(i_{1}+y_{1}\phi_{1}\right)+\mu\right),\\\
y_{1}^{\prime}=\beta_{1}\alpha_{1}s_{2}\left(i_{1}+y_{1}\phi_{1}\right)-y_{1}\left(\gamma_{1}+\mu\right).\end{cases}$
(27)

Besides the DFE where $s=1$ and all the other compartments are $0$, this
system has also two other boundary points. With
$\mathcal{R}_{i}=\frac{\beta_{i}}{\gamma_{i}+\mu}$, these are:

  1. 1.

one with $i_{2}=r_{2}=s_{2}=y_{1}=y_{2}=0$, given by

$\displaystyle
E_{1}=\left(\frac{\mu}{\beta_{1}}\left(\mathcal{R}_{1}-1\right),\frac{\mu\gamma_{1}}{\beta_{1}\left(\alpha_{1}+\mu\right)}(\mathcal{R}_{1}-1),\frac{\alpha_{1}\gamma_{1}}{\beta_{1}\left(\alpha_{1}+\mu\right)}(\mathcal{R}_{1}-1),0,0,0,0,0,\frac{1}{\mathcal{R}_{1}}\right)$

  2. 2.

and one with $i_{1}=r_{1}=s_{1}=y_{1}=y_{2}=0$, given by

$\displaystyle
E_{2}=\left(0,0,0,0,\frac{\mu}{\beta_{2}}\left(\mathcal{R}_{2}-1\right),\frac{\mu\gamma_{2}}{\beta_{2}\left(\alpha_{2}+\mu\right)}(\mathcal{R}_{2}-1),\frac{\alpha_{2}\gamma_{2}}{\beta_{2}\left(\alpha_{2}+\mu\right)}(\mathcal{R}_{2}-1),0,\frac{1}{\mathcal{R}_{2}}\right)$

Thus, $\mathcal{R}_{i},i=1,2$ are the bifurcation values at which these two
boundary points appear.

The maximal disease set contains $I_{i},R_{i},S_{i},Y_{i},i=1,2$. The DFE may
be determined already using the disease set $I_{i},Y_{i},i=1,2$, which has the
advantage of possessing a simple characteristic polynomial with two factors
$R_{1}(X),R_{2}(X)$, which yields:

$R_{J}=\max[R_{1}(X),R_{2}(X)],R_{1}(X)=\frac{\beta_{2}\left(\alpha_{2}s_{1}\phi_{2}+s\right)}{\gamma_{2}+\mu},R_{2}(X)=\frac{\beta_{1}\left(\alpha_{1}s_{2}\phi_{1}+s\right)}{\gamma_{1}+\mu}.$

Also, our scripts find that

$R_{ji}=\;\mathsf{s}_{dfe}\;\mathcal{R}_{j},j\neq i,i=1,2.$ (28)

Finally, applying the NGM script to $E_{i},i=1,2$ yields the elegant relation

$R_{0}=\;\mathsf{s}_{dfe}\max\left[\mathcal{R}_{1},\mathcal{R}_{2}\right]=\max\left[R_{21},R_{12}\right].$
(29)

###### Remark 16

Note the notations $R_{1}(X),R_{2}(X)$, suggesting that we want to view these
as polynomials in the variables of the model, rather than values evaluated at
one of the fixed points.

We end this section by drawing the attention to the object which allowed
computing the key polynomials $R_{1}(X),R_{2}(X)$.

###### Definition 4

A) A minimal disease set $\mathcal{I}$ is a minimal set which still allows the
computation of the DFE, after being set to $0$.

B) The model factors are the factors which may admit positive roots in the
characteristic polynomial of the Jacobian with all variables in $\mathcal{I}$
set to $0$.

###### Remark 17

Assume w.l.o.g. $\mathcal{R}_{1}<\mathcal{R}_{2}$. Two situations may arise:

$\displaystyle\begin{cases}\;\mathsf{s}_{dfe}\mathcal{R}_{1}<\mathcal{R}_{1}<\;\mathsf{s}_{dfe}\mathcal{R}_{2}<\mathcal{R}_{2}\\\
\;\mathsf{s}_{dfe}\mathcal{R}_{1}<\;\mathsf{s}_{dfe}\mathcal{R}_{2}<\mathcal{R}_{1}<\mathcal{R}_{2},\end{cases}$

and in each of them $1$ may lie in any of the partition intervals. This gives
raise to $6$ disjoint cases:

$\begin{cases}\mathcal{R}_{1}<\mathcal{R}_{2}\leq 1&\text{the DFE is the only
boundary equilibrium}\\\
\;\mathsf{s}_{dfe}\mathcal{R}_{1}<\;\mathsf{s}_{dfe}\mathcal{R}_{2}<1<\mathcal{R}_{1}<\mathcal{R}_{2}&\text{both
$E_{1},E_{2}$ exist and are unstable}\\\
\;\mathsf{s}_{dfe}\mathcal{R}_{1}<1<\min[\mathcal{R}_{1},\;\mathsf{s}_{dfe}\mathcal{R}_{2}]<\mathcal{R}_{2}&\text{$E_{1}$
unstable, $E_{2}$ stable}\\\
\;\mathsf{s}_{dfe}\mathcal{R}_{1}<\mathcal{R}_{1}<1<\;\mathsf{s}_{dfe}\mathcal{R}_{2}<\mathcal{R}_{2}&\text{only
$E_{2}$ exists and is stable}\\\
\;\mathsf{s}_{dfe}\max[\mathcal{R}_{1},\;\mathsf{s}_{dfe}\mathcal{R}_{2}]<1<\mathcal{R}_{2}&\text{only
$E_{2}$ exists and is unstable}\\\
1<\;\mathsf{s}_{dfe}\mathcal{R}_{1}<\;\mathsf{s}_{dfe}\mathcal{R}_{2}&\text{competition
between the two stable dominant}\\\ &\text{strains }E_{1},E_{2}.\end{cases}$
(30)

All these cases have been investigated in detail, for a more general model, in
[17], reviewed in the next section; it turns out that the results are fully
determined by the model factors.

Before proceeding, let us give a name to the very interesting structure we
have started to investigate.

###### Definition 5

A Descartes multi strain model of order $M$ is an epidemic model for which the
characteristic polynomial of the Jacobian factors completely over the
rationals as a product of terms, precisely $M$ of which are “Descartes
polynomials". For such models, the Jacobian factorization threshold is defined
as

$R_{J}(X):=\max_{1\leq m\leq M}R_{m}(X).$

One may check that

###### Lemma 3

For Descartes multi strain models of order $K$, the local stability set is a
subset of

$R_{J}(X)\leq 1.$

###### Remark 18

The example of this section is a Descartes two -strain model (since the
characteristic polynomial has only linear factors, precisely $2$ of which have
constant coefficient which may change sign).

###  6.3 Effects of single-strain vaccination on the dynamics of a multi-
strain host only dengue model with ADE

In this section we will show that the mysterious formula (29) continues to
hold under the considerably more complicated two strains model of [17], with
vaccination applied to one strain only. The model studied in [17] is depicted
in Figure 5.

Figure 5: Schematic representation of the infection status due to the
concomitant transmission of viruses $1$ and $2$, considering that the
population is vaccinated against the virus $1$.

This model involves twelve compartments, two of which capture the vaccination
against strain $1$.

  1. 1.

$S=S_{0}$ are individuals susceptible to both strains;

  2. 2.

$I_{i}$, for $i,j=1,2$ are individuals infected with strain $i$ and with
temporary cross-immunity to strain $j\neq i$;

  3. 3.

$C_{i}$ ($R_{i}$ in the original model of [5]) are individuals recovered from
strain $i$ and hence permanently immune to it, and with temporary cross-
immunity to strain $j\neq i$;

  4. 4.

$R_{i}$ ($S_{i}$ in the original model of [5]) are unvaccinated individuals
who have recovered from strain $i$, but have become susceptible to the other
strain $j$;

  5. 5.

$Y_{j}$ ($I_{ij}$ in the original model of [5]) are individuals previously
infected with strain $i$ and immune to it, but reinfected with strain
$j,i,j=1,2,i\neq j$;

  6. 6.

$R=R_{12}$ are individuals immune to all the strains;

  7. 7.

Finally, there are individuals $V$ who are vaccinated against strain $1$ and
still susceptible to strain $2$, and individuals $R_{v1}=S_{v}=Z$ who have
been vaccinated against strain $1$ and subsequently became infected by strain
$2$.

Denote by
$N(t)=S(t)+V(t)+I_{1}(t)+I_{2}(t)+C_{1}(t)+C_{2}(t)+R_{1}(t)+R_{2}(t)+Y_{1}(t)+Y_{2}(t)+S_{v}(t)+R_{12}(t)$
the total population, put $J_{i}=I_{i}+Y_{i},\quad i=1,2$, and assume that the
two forces of infection acting on $S$ are:

$F_{i}=\beta_{i}\frac{J_{i}}{N},$

and that the forces of infection acting on $Y_{i}=S_{i},i=1,2$ are:

$\alpha_{1}\beta_{1}\frac{J_{1}}{N},\alpha_{2}\beta_{2}\frac{J_{2}}{N},\alpha_{v}\beta_{2}\frac{J_{2}}{N},$

where $\alpha_{1},\alpha_{2},\alpha_{v}$ denote decreases or increases of the
susceptibility to secondary infections ($\alpha_{i}>1$ implying an ADE
effect).

The following equations, with appropriate initial conditions, represent the
disease dynamics model:

$\displaystyle\dfrac{dS}{dt}$
$\displaystyle=(1-\xi)\mu-\beta_{1}J_{1}\frac{S}{N}-\beta_{2}J_{2}\frac{S}{N}-\mu
S$ $\displaystyle\dfrac{dI_{1}}{dt}$
$\displaystyle=\beta_{1}J_{1}\frac{S}{N}-(\gamma_{1}+\mu)I_{1}$
$\displaystyle\dfrac{dC_{1}}{dt}$
$\displaystyle=\gamma_{1}I_{1}-(\theta_{1}+\mu)C_{1}$
$\displaystyle\dfrac{dR_{1}}{dt}$
$\displaystyle=\theta_{2}C_{1}-\alpha_{2}\beta_{2}J_{2}\frac{R_{1}}{N}-\mu
R_{1}$ $\displaystyle\dfrac{dY_{2}}{dt}$
$\displaystyle=\alpha_{2}\beta_{2}J_{2}\frac{R_{1}}{N}+\alpha_{v}\beta_{2}J_{2}\frac{S_{v}}{N}-(\gamma_{2}+\mu)Y_{2}$
$\displaystyle\dfrac{dI_{2}}{dt}$
$\displaystyle=\beta_{2}J_{2}\frac{S}{N}-(\gamma_{2}+\mu)I_{2}$
$\displaystyle\dfrac{dC_{2}}{dt}$
$\displaystyle=\gamma_{2}I_{2}-(\theta_{2}+\mu)C_{2}$
$\displaystyle\dfrac{dR_{2}}{dt}$
$\displaystyle=\theta_{1}C_{2}-\alpha_{1}\beta_{1}J_{1}\frac{R_{2}}{N}-\mu
R_{2}$ $\displaystyle\dfrac{dY_{1}}{dt}$
$\displaystyle=\alpha_{1}\beta_{1}J_{1}\frac{R_{2}}{N}-(\gamma_{1}+\mu)Y_{1}$
$\displaystyle\dfrac{dV}{dt}$ $\displaystyle=\xi\mu-(\theta_{v}+\mu)V$
$\displaystyle\dfrac{dS_{v}}{dt}$
$\displaystyle=\theta_{v}V-\alpha_{v}\beta_{2}J_{2}\frac{S_{v}}{N}-\mu S_{v}$
$\displaystyle\dfrac{dR_{12}}{dt}$
$\displaystyle=\gamma_{1}Y_{1}+\gamma_{2}Y_{2}-\mu R_{12}$ (31)

Table 1 summarizes the parameters and compartments of the model.

Table 1: Parameters and compartments of the model.

Parameter | Description (for $i,j=1,2$)
---|---
$\mu$ | Birth rate
$\mu$ | Per capita death rate
$\beta_{i}$ | Transmission rate of virus $i$
$\gamma_{i}$ | Per capita recovery rate of infected people with virus $i$
$\theta_{i}$ | Per capita loss rate of cross-immunity to virus $i$ after previous infection with virus $j$
$\theta_{v}$ | Per capita loss rate of cross-immunity to virus $2$ obtained by vaccination
$\alpha_{i}$ | ADE factor that can alter the susceptibility of unvaccinated individuals to the virus $i$
$\alpha_{v}$ | ADE factor that can alter the susceptibility of vaccinated individuals to virus $2$
$\xi$ | Per capita vaccination rate
Compartments | Description
$S$ | Susceptible individuals to both virus
$V$ | Vaccinated individuals against the virus $1$
$I_{i}$ | Individuals with primary infection by the virus $i$
$C_{i}$ | Individuals recovered from infection with virus $i$ and have cross-immunity to virus $j$
$R_{i}$ | Unvaccinated individuals immune to virus $i$ and susceptible to virus $j$
$Z=S_{v}$ | Individuals vaccinated for virus $1$ and susceptible to virus $2$
$Y_{1}$ | Individuals infected by virus $1$ and recovered and hence immune to virus $2$
$Y_{2}$ | Individuals infected by virus $2$ and immune to virus $1$ either due to recovery or vaccination
$R_{12}$ | Individuals immune to both virus

This system does not have negative cross effects; therefore, it leaves the
non-negative quadrant invariant [33]. It follows from the equations that

$\frac{dN(t)}{dt}=\mu(1-N(t)).$

Therefore,

$\lim_{t\rightarrow+\infty}N(t)=1.$

Assuming $N(0)=1$ implies that $N(t)=1$, for $t\geq 0$. Using this, we may
assume w.l.o.g. that $N=1$, and work with the proportions, to be denoted by
the corresponding lowercase letters.

The only non-zero compartments in the DFE, to be denoted by $E_{0}$, are
easily found to be

$s_{0}=1-\xi,z_{0}=\xi\frac{\theta_{v}}{\mu+\theta_{v}},v_{0}=\xi\frac{\mu}{\mu+\theta_{v}};$

in fact, the last value holds at any fixed point. As known from [17], there
are also two endemic points on the boundary, whose rather complicated formulas
will be given later.

###### Remark 19

From a modeling point of view, this system has crucial parameters like
$\alpha_{v}$ (note that $\alpha_{v}=0$ means perfect vaccination, and
$\alpha_{v}=1$ means that infection by second strain is equally likely for
vaccinated people).

Due to conservation, the system evolves in a compact domain, and so we may
eliminate one compartment, for example $V$, from the analysis. Finally, the
last compartment does not send input to the others, and therefore may also be
disregarded in the analysis.

####  6.3.1 The Jacobian $R_{J}(X)$ is the max of two polynomials, obtained
using a minimal disease set

We may tackle this example via the Jacobian factorization approach, choosing
the minimal disease set $\mathcal{I}=(i_{1},i_{2},y_{1},y_{2})$, just like in
the previous section. Again, the characteristic polynomial of the Jacobian
with the variables in $\mathcal{I}$ set to $0$ factors completely as a product
of linear terms

$(\mu+u)^{5}\left(\gamma_{1}+\mu+u\right)\left(\gamma_{2}+\mu+u\right)\left(\theta_{1}+\mu+u\right)\left(\theta_{2}+\mu+u\right)\left(\mu+u+\theta_{v}\right)\times$
$\displaystyle\left(\gamma_{1}+\mu-\alpha_{1}\beta_{1}r_{2}-\beta_{1}s+u\right)\left(\gamma_{2}+\mu-\alpha_{2}\beta_{2}r_{1}-\beta_{2}s-\beta_{2}z\alpha_{v}+u\right),$

only two of which (the 7’th and 8’th factors) may yield positive eigenvalues.
Both are of Descartes type, and instability may occur iff

$R_{J}(X):=\max[R_{1}(X),R_{2}(X)]=\max[\frac{\beta_{1}\left(\alpha_{1}r_{2}+s\right)}{\gamma_{1}+\mu},\frac{\beta_{2}\left(\alpha_{2}r_{1}+z\alpha_{v}+s\right)}{\gamma_{2}+\mu}]>1.$
(32)

At the DFE, $r_{1}=r_{2}=0,$ and this yields

$R_{J}:=R_{J}(E_{0})=R_{N}=\max[s_{0}\frac{\beta_{1}}{d_{1}},s_{0}\frac{\beta_{2}}{d_{2}}+z_{0}\frac{\beta_{2}\alpha_{v}}{d_{2}}],\quad
d_{1}=\gamma_{1}+\mu,d_{2}=\gamma_{2}+\mu.$ (33)

This expression reveals a pattern similar to (14), with the difference that
the existence of two strains are reflected in the max, and that the second
strain is alimented by two classes of susceptibles, one of which is the people
vaccinated against the first strain.

In addition to the disease-free equilibrium, there might exist two more
equilibriums on the boundary: the endemic equilibrium where there are only
infections by the strain $1$, $E_{1}$, and the endemic equilibrium where there
are only infections by the strain $2$, $E_{2}$, reviewed in the next section.

####  6.3.2 The endemic boundary equilibrium $E_{i}$ exist iff
$R_{i}(E_{0})>1$

At the equilibrium $E_{1}$, the values of $I_{2},C_{2},R_{2},Y_{1},Y_{2}$ and
$R_{12}$ are zero. The coordinates are easily found by the “Solve" command.
Those of $V,Z$ are the same as at the DFE, and the others are:

$\displaystyle
s_{1}=\frac{\gamma_{1}+\mu}{\beta_{1}},\;i_{1}=\frac{\mu}{\beta_{1}}\left[\frac{1-\xi}{s_{1}}-1\right]:=\frac{\mu}{\beta_{1}}\left({\mathcal{R}_{1}}-1\right),\;c_{1}=\frac{\gamma_{1}}{\theta_{1}+\mu}i_{1},\;r_{1}=\frac{\theta_{1}}{\mu}c_{1}$
(34)

where

$\mathcal{R}_{1}=(1-\xi)\frac{\beta_{1}}{\gamma_{1}+\mu}=R_{1}(E_{0})$ (35)

(the endemic equilibrium $E_{1}$ exists if and only if $\mathcal{R}_{1}>1$).

At the equilibrium $E_{2}$, the values of $I_{1},C_{1},R_{1}$ and $Y_{1}$ are
zero, and that of $V$ is the same as at the DFE.

The solutions of the $E_{2}$ system involve all complicated square roots. In
such a case, it is more convenient to replace the “Solve" command by our RUR
algorithm, which requires the user to input a variable to reduce $2$. The
normal choice is $i_{2}$ (which transitions to positive at the bifurcation
value), but here we will use $s$, to check the results of [17], who find,
using as reduction scalar $x={\beta_{2}j_{2}},$ that

$\displaystyle s_{2}=\frac{(1-\xi)\mu}{x+\mu},\quad
i_{2}=\frac{(1-\xi)x\mu}{(x+\mu)(\gamma_{2}+\mu)},\quad
c_{2}=\frac{(1-\xi)x\gamma_{2}\mu}{(x+\mu)(\gamma_{2}+\mu)(\theta_{1}+\mu)},$
$\displaystyle
r_{2}=\frac{(1-\xi)x\gamma_{2}\theta_{1}}{(x+\mu)(\gamma_{2}+\mu)(\theta_{1}+\mu)},\quad
z_{2}=\frac{v\theta_{v}\mu}{(\theta_{v}+\mu)(\alpha_{v}x+\mu)},$ (36)
$\displaystyle
y_{2}=\frac{v\alpha_{v}x\theta_{v}\mu}{(\alpha_{v}x+\mu)(\theta_{v}+\mu)(\gamma_{2}+\mu)},$

and that $x$ is solution of the quadratic equation

$ax^{2}+bx+c=0,\quad\begin{cases}a&=\alpha_{v}\\\
b&=\mu\alpha_{v}\left[1-\frac{\beta_{2}(1-\xi)}{\gamma_{2}+\mu}\right]+\mu\left[1-\frac{\beta_{2}\alpha_{v}\theta_{v}v}{(\gamma_{2}+\mu)(\theta_{v}+\mu)}\right]\\\
c&=\mu^{2}\left(1-\mathcal{R}_{2}\right).\end{cases}$

The equilibrium $E_{2}$ exists iff $\mathcal{R}_{2}>1$, where

$\mathcal{R}_{2}=\frac{\beta_{2}}{\gamma_{2}+\mu}\left[1-\xi+\xi\frac{\alpha_{v}\theta_{v}}{\theta_{v}+\mu}\right]=\frac{\beta_{2}}{\gamma_{2}+\mu}\left[s_{0}+\alpha_{v}z_{0}\right]=R_{2}(E_{0})$
(37)

If $\mathcal{R}_{2}\leq 1$, the fractions in the expression of $b$ must be
smaller than one or equal to one, and it is not possible for both to be one.
Therefore, $b>0$. We also have $c\geq 0$. Since that $a>0$, the equation
(6.3.2) does not have roots with positive real parts. This implies that there
is no endemic equilibrium like $E_{2}$. Thus, in this case, $c<0$. Since the
coefficient $a$ is positive, the equation (6.3.2) has two real roots and only
one of them is positive. In resume, if $\mathcal{R}_{2}>1$, there is a unique
endemic equilibrium where there are infections only by the strain $2$.

####  6.3.3 The recipe next generation matrix $R_{0}$ and the Jacobian
factorization one coincide

This section shows that the polynomials $R_{1}(X),R_{2}(X)$ in this example
may also be obtained via the next generation matrix approach, as eigenvalues
of the $K$ matrix, by a judicious choice of infectious classes.

One may choose as infectious subset the nine compartments that are $0$ in the
limit, but a luckier choice here is the smaller subset
$\mathcal{I}=\\{I_{1},I_{2},Y_{2},Y_{1}\\}$, which has as eigenvalues
precisely the expressions $R_{1}(X),R_{2}(X)$ in (32).

The decomposition matrices are

$V=\begin{pmatrix}\gamma_{1}+\mu&0&0&0\\\ 0&\gamma_{2}+\mu&0&0\\\
0&0&\gamma_{1}+\mu&0\\\
0&0&0&\gamma_{2}+\mu\end{pmatrix},F=\left(\begin{array}[]{cccc}\beta_{1}s&0&0&\beta_{1}s\\\
0&\beta_{2}s&\beta_{2}s&0\\\
0&\beta_{2}z\alpha_{v}&\beta_{2}\;\mathsf{z}\alpha_{v}&0\\\ 0&0&0&0\\\
\end{array}\right)=sB_{0}+zB_{v},$

where $B_{0},B_{v}$ are:

$\displaystyle B_{0}=\left(\begin{array}[]{cccc}\beta_{1}&0&0&\beta_{1}\\\
0&\beta_{2}&\beta_{2}&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{array}\right),\;B_{v}=\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&0\\\
0&\beta_{2}\alpha_{v}&\beta_{2}\alpha_{v}&0\\\ 0&0&0&0\\\ \end{array}\right).$

The explicit non-zero eigenvalues of the next generation matrix
$(\;\mathsf{s}B_{0}+\\\ zB_{v})V^{-1}$ are

$\left(\frac{\beta_{1}\;\mathsf{s}}{\gamma_{1}+\mu},\frac{\beta_{2}\left(\\\
z\alpha_{v}+\;\mathsf{s}\right)}{\gamma_{2}+\mu}\right),$ (39)

confirming the result of the Jacobian method.

Let us note finally that (32), and the result of this section imply the
relation

$R_{0}=\max\left[\mathcal{R}_{1},\mathcal{R}_{2}\right],$ (40)

where $\mathcal{R}_{i},i=1,2$ denote the bifurcation parameters at which the
boundary points $E_{i}$ start to exist.

###### Remark 20

Interestingly, $R_{0}=\max[\mathcal{R}_{1},\mathcal{R}_{2}]$ is the max of two
quantities which satisfy that ${\mathcal{R}_{i}}>1,i=1,2$ are precisely the
domains where endemic points $E_{i}$ containing exactly one of the strains
appear – see (40). This formula, natural in cases where the next generation
matrix has block structure, seems to be a general feature of multi-strain
models, even when the block structure is not apparent.

In the case of this section, there seems to be more specific structure: the
Jacobian factorization approach allows introducing two “Descartes type" (see
Definition 3) factors $R_{i}(X),i=1,2$ of the characteristic polynomial, which
are such that

  1. 1.

The existence conditions for $E_{i}$ may be expressed as
$\mathcal{R}_{i}:=R_{i}(DFE)>1$ – see (35), (37), (42).

  2. 2.

The invasion reproduction numbers may be obtained simply by substituting the
coordinates of the dominance boundary equilibria into the corresponding
factor. More precisely, the invasion number of the fixed point $E_{i}$ for
strain $i$ is given by $R_{ji}=R_{j}(E_{i})$.

Open question 2: Does the relation $R_{0}=\max_{1}^{K}\mathcal{R}_{k}$ hold
for all Descartes multi strain models of order $K$? (recall Definition 5 and
Lemma 3).

####  6.3.4 The invasion reproduction number of $E_{i}$ is given by
$R_{j}(E_{i})$

The invasion reproduction numbers (see for example [27]) may, just as the
basic reproduction number, be calculated using the next generation matrix.

Our script yield quickly that

${\mathcal{R}_{ji}}=R_{j}(E_{i}),i=1,2,j\neq i.$ (41)

Open question 3: Do the formulas connecting (40) and (41) to the Jacobian
factorization

$\begin{cases}\mathcal{R}_{i}=R_{i}(E_{0}),R_{0}=\max\left[\mathcal{R}_{1},\mathcal{R}_{2}\right],\\\
\mathcal{R}_{ji}=R_{j}(E_{i}),\text{ where $R_{i}$ denote polynomials obtained
via}\\\ \text{ the Jacobian factorization approach},\end{cases}$ (42)

hold, for some general class of epidemic models?

C) For “two strain epidemic models", what conditions must be satisfied to
ensure the inequalities $\mathcal{R}_{ji}<\mathcal{R}_{j},i=1,2,j\neq i$?

To resolve these questions, it might be useful to study the three and four
strains generalizations of this problem, and to investigate “non-simple"
multi-strain models (in which the characteristic polynomial contains non
Descartes type polynomials).

##  7 Vector-host models

###  7.1 The Jacobian $R_{0}$ is the square of the recipe NGM $R_{0}$ for the
dengue vector-host model without demography of [16]

[16, (28)] considers a “no demography/conservation" model with 6 compartments,
three of which represent hosts, while the rest represent the vector. Note that
such models with no demography do not have a finite set of fixed points. The
DFE is not unique, it coincides with the initial conditions. However, our
algorithm works just fine. The model, after removing two “R" classes which do
not affect the rest, is:

$\begin{cases}S_{1}^{\prime}=-\frac{\beta_{21}I_{2}S_{1}}{N_{1}}\\\
S_{2}^{\prime}=-\frac{\beta_{12}I_{1}S_{2}}{N_{2}}\\\
\begin{pmatrix}I_{1}^{\prime}\\\
I_{2}^{\prime}\end{pmatrix}=\begin{pmatrix}-\gamma_{1}&\frac{\beta_{21}S_{1}}{N_{1}}\\\
\frac{\beta_{12}S_{2}}{N_{2}}&-\gamma_{2}\end{pmatrix}\begin{pmatrix}I_{1}\\\
I_{2}\end{pmatrix}\end{cases}$ (43)

The call “${\text{inf}}=\\{1,2\\};NGM[Brouwer22,\text{inf}]$" of our script
yields that the decomposition matrices are

$F=\left(\begin{array}[]{cc}0&\frac{\beta_{21}S_{1}}{N_{1}}\\\
\frac{\beta_{12}S_{2}}{N_{2}}&0\\\
\end{array}\right),\;V=\left(\begin{array}[]{cc}\gamma_{1}&0\\\
0&\gamma_{2}\\\ \end{array}\right),$
$K=\left(\begin{array}[]{cc}0&\frac{\beta_{21}S_{1}}{\gamma_{2}N_{1}}\\\
\frac{\beta_{12}S_{2}}{\gamma_{1}N_{2}}&0\\\ \end{array}\right),$

and

$R_{F}=\sqrt{\frac{S_{1}S_{2}\beta_{12}\beta_{21}}{N_{1}N_{2}\gamma_{1}\gamma_{2}}}.$
(44)

After using the fact that the DFE is determined by the initial conditions
$S_{1}=N_{1}$, and $S_{2}=N_{2}$ we obtain the basic reproduction number

$R_{F}=\sqrt{\frac{\beta_{12}\beta_{21}}{\gamma_{1}\gamma_{2}}}$ (45)

of [16, 40].

Here the characteristic polynomial is of Descartes type and the Jacobian
method, and the RUR method, yield both the square of the (modified) formula
(44) $R_{J}={\frac{\beta_{12}\beta_{21}}{\gamma_{1}\gamma_{2}}}$.

###### Remark 21

Note that [16, 35] offers yet another admissible decomposition, based on a
different biological interpretation, with $R_{F}=R_{J}$, and raises the
question of which of the answers is more relevant for a given epidemics.
Deciding this from the ODE model only seems impossible.

###  7.2 The two groups model in [42, (5.8)] does not obey a square relation

The two groups model in [42, (5.8)] defined by

$\displaystyle\begin{cases}S_{1}^{\prime}=-\frac{\beta_{11}I_{1}S_{1}}{N_{1}}-\frac{\beta_{21}I_{2}S_{1}}{N_{1}}+\lambda_{1}-\mu_{1}S_{1}\\\
S_{2}^{\prime}=-\frac{\beta_{12}I_{1}S_{2}}{N_{2}}-\frac{\beta_{22}I_{2}S_{2}}{N_{2}}+\lambda_{2}-\mu_{2}S_{2}\\\
\begin{pmatrix}I_{1}^{\prime}\\\
I_{2}^{\prime}\end{pmatrix}=\begin{pmatrix}-\gamma_{1}-\mu_{1}-\delta_{1}+\frac{\beta_{11}S_{1}}{N_{1}}&\frac{\beta_{21}S_{1}}{N_{1}}\\\
\frac{\beta_{12}S_{2}}{N_{2}}&\frac{\beta_{22}S_{2}}{N_{2}}-\gamma_{2}-\mu_{2}-\delta_{2}\end{pmatrix}\begin{pmatrix}I_{1}\\\
I_{2}\end{pmatrix}\end{cases}$

is not anymore a vector-host model, due to the addition of the "intra-group
contact infection rates" $\beta_{11},\beta_{22}$.

The DFE is
$\left\\{0,0,\frac{\lambda_{1}}{\mu_{1}},\frac{\lambda_{2}}{\mu_{2}}\right\\}$,
and the $R_{N}$ is quite complicated:

$\displaystyle\scalebox{1.0}{$+\frac{\beta_{22}\gamma_{1}N_{1}S_{2}+\beta_{11}\gamma_{2}N_{2}S_{1}+\beta_{22}\delta_{1}N_{1}S_{2}+\beta_{11}\delta_{2}N_{2}S_{1}+\beta_{22}\mu_{1}N_{1}S_{2}+\beta_{11}\mu_{2}N_{2}S_{1}}{2N_{1}N_{2}\left(\gamma_{1}+\delta_{1}+\mu_{1}\right)\left(\gamma_{2}+\delta_{2}+\mu_{2}\right)}$}.$

The Jacobian factorization method yields a different answer, for a
characteristic polynomial which is not of Descartes type, precisely because of
the addition of $\beta_{11},\beta_{22}$.

$R_{J}=\frac{\beta_{22}N_{1}S_{2}\left(\gamma_{1}+\delta_{1}+\mu_{1}\right)+\beta_{11}\gamma_{2}N_{2}S_{1}+\beta_{11}\delta_{2}N_{2}S_{1}+\beta_{11}\mu_{2}N_{2}S_{1}+\beta_{12}\beta_{21}S_{1}S_{2}}{N_{1}N_{2}\left(\gamma_{1}+\delta_{1}+\mu_{1}\right)\left(\gamma_{2}+\delta_{2}+\mu_{2}\right)+\beta_{11}\beta_{22}S_{1}S_{2}}.$

##  8 Multi-strain vector-host models

###  8.1 A two-strain vector-host model of Feng and Velasco-Hernández [28],
where the square relation holds for the basic reproduction number

[28] consider a human population settled in a region where a mosquito
population of the genus Aedes is present and carrier of two strains of the
dengue virus. Let $V_{i},I_{i},Y_{i},i=1,2$ denote the infected mosquitoes,
individuals infected by one strain, and individuals having suffered a
secondary infection, let $N=S+R+\sum_{i=1}^{2}I_{i}+Y_{i}$ denote the total
human population, and let
$B_{1}=\frac{\beta_{1}V_{1}(t)}{c+w_{h}N},B_{2}=\frac{\beta_{2}V_{2}(t)}{c+w_{h}N}$
denote the rates of infections in human hosts produced by the two strains. The
model is defined as follows:

$\displaystyle\begin{cases}S^{\prime}(t)=h-S(t)\left(B_{1}+B_{2}\right)-\mu
S(t),\\\ I_{1}^{\prime}(t)=B_{1}S(t)-\sigma_{2}B_{2}I_{1}(t)-\mu I_{1}(t),\\\
I_{2}^{\prime}(t)=B_{2}S(t)-\sigma_{1}B_{1}I_{2}(t)-\mu I_{2}(t),\\\
Y_{1}^{\prime}(t)=\sigma_{1}B_{1}I_{2}(t)-(e_{1}+\mu+r)Y_{1}(t),\\\
Y_{2}^{\prime}(t)=\sigma_{2}B_{2}I_{1}(t)-(e_{2}+\mu+r)Y_{2}(t),\\\
R^{\prime}(t)=r(Y_{1}(t)+Y_{2}(t))-\mu R(t),\\\ \\\
V_{1}^{\prime}(t)=\alpha_{1}\frac{I_{1}(t)+Y_{1}(t)}{c+w_{v}N}M(t)-\delta
V_{1}(t),\\\
V_{2}^{\prime}(t)=\alpha_{2}\frac{I_{2}(t)+Y_{2}(t)}{c+w_{v}N}M(t)-\delta
V_{2}(t)\\\
M^{\prime}(t)=q-M(t)\left(\alpha_{1}\frac{I_{1}(t)+Y_{1}(t)}{c+w_{v}N}+\alpha_{2}\frac{I_{2}(t)+Y_{2}(t)}{c+w_{v}N}\right)-\delta
M(t).\end{cases}$

The DFE is given by $E_{0}=\left(h/\mu,0,0,0,0,0,0,0,q/\delta\right)$. For the
infectious set $I_{1},I_{2},Y_{1},Y_{2},V_{1},V_{2}$, the $F$ and $V$ matrices
used in the next generation approach are given by

$\displaystyle F=\left(\begin{array}[]{cccccc}0&0&0&0&\beta_{1}s_{dfe}&0\\\
0&0&0&0&0&\beta_{2}s_{dfe}\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
\alpha_{1}M_{dfe}&0&\alpha_{1}M_{dfe}&0&0&0\\\
0&\alpha_{2}M_{dfe}&0&\alpha_{2}M_{dfe}&0&0\\\ \end{array}\right),$
$\displaystyle V=\left(\begin{array}[]{cccccc}\mu&0&0&0&0&0\\\
0&\mu&0&0&0&0\\\ 0&0&e_{1}+r+\mu&0&0&0\\\ 0&0&0&e_{2}+r+\mu&0&0\\\
0&0&0&0&\delta&0\\\ 0&0&0&0&0&\delta\\\ \end{array}\right)$

with $M_{dfe}=q/\delta$. Then,

$FV^{-1}=\left(\begin{array}[]{cccccc}0&0&0&0&\frac{\beta_{1}s_{dfe}}{{\delta}}&0\\\
0&0&0&0&0&\frac{\beta_{2}s_{dfe}}{\delta}\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
\frac{\alpha_{1}M_{dfe}}{{\mu}}&0&\frac{\alpha_{1}M_{dfe}}{{e_{1}+\mu+\xi}}&0&0&0\\\
0&\frac{\alpha_{2}M_{dfe}}{\mu}&0&\frac{\alpha_{2}M_{dfe}}{{e_{2}+\mu+\xi}}&0&0\\\
\end{array}\right)$

We obtain a basic reproduction number which is a max

$\mathcal{R}=max\left(\sqrt{\mathcal{R}_{1}},\sqrt{\mathcal{R}_{2}}\right),\mathcal{R}_{i}:={s_{0}m_{0}}\frac{{{\alpha_{i}\beta_{i}}}}{{\delta\mu}},$
(48)

just like (40), but contains also the extra square roots typical of vector-
host models.

Furthermore, it may be checked that this is precisely the square root of the
answer given by the Jacobian factorization method, which decomposes the
characteristic polynomial of the Jacobian as the product of five linear
factors with negative roots, and two quadratic Descartes type polynomials.

There also two boundary (dominance) equilibria, where only one strain
survives. The non-zero coordinates at the first one, $E_{1}$, are given by

$\alpha_{1}i_{1}=\delta\frac{\mathcal{R}_{1}-1}{m_{0}\beta_{1}/(\mu)+1},\beta_{1}v_{1}=\mu\frac{\mathcal{R}_{1}-1}{s_{0}\alpha_{1}(\delta)+1},s=\mu\frac{\alpha_{1}s_{0}+\delta}{\alpha_{1}\beta_{1}m_{0}+\alpha_{1}\mu},$

with similar formulas holding for the other boundary point $E_{2}$, by
symmetry. Thus, these points become positive precisely when the corresponding
factor of the DFE becomes bigger than $1$, causing instability.

Since we had trouble with computing the invasion reproduction numbers, we
switched to the “simplified model" of [28], in which $M$ is eliminated by
noting that the equation for the total vector population $T=M+V_{1}+V_{2}$ is
$T^{\prime}=q-\delta T$, and that, assuming
$T_{0}=\lim_{t\to\infty}T(t)=q/\delta$, $M$ can be removed from the system by
substituting

$M=q/\delta-V_{1}-V_{2}.$ (49)

As a first consequence of using (49), the $R_{N}$ becomes equal to $R_{J}$.

However, the recipe $R_{0}$ at $E_{1}$ for the natural choice of “inf" is very
complicated, and [28] provide here a laborious local stability analysis, with
complicated result, via the third order Routh-Hurwitz conditions.

We note finally that the characteristic polynomial for $jac(E_{1})$ has two
factors of degree $3$, one of which is Descartes typeand one which is not. The
Descartes type factor yields a polynomial $R_{1}(X)$. Putting this together
with its symmetric $R_{2}(X)$ allows finally defining

$R_{J}(X)=max_{j}[R_{1}(X),R_{2}(X)]=max_{j}[\frac{\alpha_{j}\beta_{j}qs/\delta}{(\beta_{j}v_{j}+\mu)\left(\alpha_{j}i_{j}+\delta\right)+\beta_{j}v_{j}\alpha_{j}s}].$

####  8.1.1 Invasion numbers of [28]

The two-strain vector-host model in [28] admits two boundary equilibria beside
the DFE in which

$S_{1}^{*},S_{2}^{*},I_{1}^{*},I_{2}^{*},V_{1}^{*},V_{2}^{*}$ are the invasion
infection classes. In this case, we consider the subset
$\mbox{in}_{1}=\left(I_{2},Y_{1},Y_{2},V_{2}\right)$ corresponding to the
invasion infection class of $E_{1}$, then

$\displaystyle F=\left(\begin{array}[]{cccc}0&0&0&b_{2}S\\\
b_{1}\sigma_{1}v_{1}&0&0&0\\\ 0&0&0&b_{2}i_{1}\sigma_{2}\\\
a_{2}\left(\frac{q}{\delta}-v_{1}\right)&0&a_{2}\left(\frac{q}{\delta}-V_{1}\right)&0\\\
\end{array}\right),$
$\displaystyle\scalebox{0.85}{$V=\left(\begin{array}[]{cccc}b_{1}\sigma_{1}V_{1}+\mu&0&0&0\\\
0&e_{1}+\mu+\xi&0&0\\\ 0&0&e_{2}+\mu+\xi&0\\\
a_{2}\left(\frac{q}{\delta}-V_{1}\right)-a_{2}\left(\frac{q}{\delta}-V_{1}-V_{2}\right)&0&a_{2}\left(\frac{q}{\delta}-V_{1}\right)-a_{2}\left(\frac{q}{\delta}-V_{1}-V_{2}\right)&a_{2}\left(I_{2}+Y_{2}\right)+\delta\\\
\end{array}\right)$},$ $\displaystyle
K=\left(\begin{array}[]{cccc}0&0&0&\frac{b_{2}S}{\delta}\\\
\frac{b_{1}\sigma_{1}V_{1}}{b_{1}\sigma_{1}V_{1}+\mu}&0&0&0\\\
0&0&0&\frac{b_{2}I_{1}\sigma_{2}}{\delta}\\\
\frac{a_{2}\left(\frac{q}{\delta}-V_{1}\right)}{b_{1}\sigma_{1}V_{1}+\mu}&0&\frac{a_{2}\left(\frac{q}{\delta}-V_{1}\right)}{e_{2}+\mu+\xi}&0\\\
\end{array}\right)$

then the IRN of strain 1 at $E_{1}$ is

$R_{1}=\frac{\sqrt{a_{2}}\sqrt{b_{2}}\sqrt{q}\sqrt{S\left(e_{2}+\mu+\xi\right)}}{\delta\sqrt{\mu}\sqrt{e_{2}+\mu+\xi}}.$

Similarly, we chose the other subset
$\mbox{in}_{1}=\left(I_{1},Y_{1},Y_{2},V_{1}\right)$ corresponding to the
invasion infection class at $E_{2}$, we obtain

$\displaystyle F=\left(\begin{array}[]{cccc}0&0&0&b_{1}S\\\
0&0&0&b_{1}I_{2}\sigma_{1}\\\ b_{2}\sigma_{2}V_{2}&0&0&0\\\
a_{1}\left(\frac{q}{\delta}-V_{2}\right)&a_{1}\left(\frac{q}{\delta}-V_{2}\right)&0&0\\\
\end{array}\right),$
$\displaystyle\scalebox{0.85}{$V=\left(\begin{array}[]{cccc}b_{2}\sigma_{2}V_{2}+\mu&0&0&0\\\
0&e_{1}+\mu+\xi&0&0\\\ 0&0&e_{2}+\mu+\xi&0\\\
a_{1}\left(\frac{q}{\delta}-V_{2}\right)-a_{1}\left(\frac{q}{\delta}-V_{1}-V_{2}\right)&a_{1}\left(\frac{q}{\delta}-V_{2}\right)-a_{1}\left(\frac{q}{\delta}-V_{1}-V_{2}\right)&0&a_{1}\left(I_{1}+Y_{1}\right)+\delta\\\
\end{array}\right)$},$ $\displaystyle
K=\left(\begin{array}[]{cccc}0&0&0&\frac{b_{1}S}{\delta}\\\
0&0&0&\frac{b_{1}I_{2}\sigma_{1}}{\delta}\\\
\frac{b_{2}\sigma_{2}V_{2}}{b_{2}\sigma_{2}V_{2}+\mu}&0&0&0\\\
\frac{a_{1}\left(\frac{q}{\delta}-V_{2}\right)}{b_{2}\sigma_{2}V_{2}+\mu}&\frac{a_{1}\left(\frac{q}{\delta}-V_{2}\right)}{e_{1}+\mu+\xi}&0&0\\\
\end{array}\right)$

then the maximum eigenvalue of K yields the IRN at $E_{2}$ which is

$R_{2}=\frac{\sqrt{a_{1}}\sqrt{b_{1}}\sqrt{q-\delta
v_{2}}\sqrt{I_{2}\sigma_{1}\left(b_{2}\sigma_{2}V_{2}+\mu\right)+e_{1}S+S(\mu+\xi)}}{\delta\sqrt{e_{1}+\mu+\xi}\sqrt{b_{2}\sigma_{2}V_{2}+\mu}}.$

###  8.2 The Dengue- Zika model with coinfection and ADE [47]

The model studied in this paper continues previous papers like Isea and
Lonngren 2016 [35], and Okuneye et al. 2017 [46], most notably by taking into
account the possibility of coinfection and of direct transmission of Zika via
sex (which entails two forces of infection for Zika transmissions in their
flow-chart, and hence an asymmetry in the results).

Introduce the following forces of infection:

$\begin{cases}F_{vd}=\beta_{hd}T_{vd},T_{vd}=I_{vd}+I_{{vc}}\nu_{d},&\text{
dengue vector force}\\\
F_{vz}=\beta_{hz}T_{vz},T_{vz}=I_{vz}+I_{{vc}}\nu_{z},&\text{ zika vector
force}\\\ F_{hz}=\beta_{vz}T_{hz},T_{hz}=I_{z}+I_{c}+J_{z}k_{z},&\text{ zika
human force}\\\ F_{hd}=\beta_{vd}T_{hd},T_{hd}=I_{d}+I_{c}+J_{d}k_{d},&\text{
dengue human force}\\\ F_{s}=\beta_{s}T_{hz}&\text{ zika human to human
force}.\end{cases}$ (56)

Note that $\nu_{d},\nu_{z}$ and $k_{d},k_{z}$ are respectively the parameters
of altered infectivity for co-infected vectors and of ADE, and note that even
when $\nu_{d}=\nu_{z}=1$, the co-infection model is more accurate than
previous works like [28], since it takes into account the existence of doubly
infected vectors $I_{{vc}}$ which influence both chains of infection.

We will consider the model :

$\begin{cases}S_{h}^{\prime}=(N_{h}-S_{h})\mu-
S_{h}\left(F_{vd}+F_{vz}+F_{s}\right),\\\ I_{d}^{\prime}=S_{h}F_{vd}-\rho
I_{d}\left(F_{vz}+F_{s}\right)-I_{d}\left(\gamma_{d}+\mu\right),\\\
I_{z}^{\prime}=S_{h}(F_{vz}+F_{s})-\rho
I_{z}F_{vd}-I_{z}\left(\gamma_{d}+\mu\right),\\\
I_{c}^{\prime}=\rho\left[I_{d}\left(F_{vz}+F_{s}\right)+I_{z}F_{vd}-I_{c}\left(\gamma_{d}+\gamma_{c}\right)\right]-\mu
I_{c},\\\
R_{d}^{\prime}=I_{d}\gamma_{d}-R_{d}\left(F_{vz}+F_{s}+\mu\right),\\\
R_{z}^{\prime}=I_{z}\gamma_{d}-R_{z}\left(F_{vd}+\mu\right),\\\
J_{d}^{\prime}=\rho\gamma_{z}I_{c}+R_{z}\left(F_{vd}-\gamma_{d}-\mu\right),\\\
J_{z}^{\prime}=\rho\gamma_{d}I_{c}+R_{d}\left(F_{vz}+F_{s}-\gamma_{z}-\mu\right),\\\
R^{\prime}=J_{d}\gamma_{d}+J_{z}\gamma_{z}-\mu R,\\\
S_{v}^{\prime}=(N_{v}-S_{v})\mu_{v}-S_{v}(F_{hd}+T_{hz}),\\\
I_{vd}^{\prime}={F_{hd}S_{v}-\rho F_{hz}I_{vd}-I_{vd}\mu_{v}},\\\
I_{vz}^{\prime}={F_{hz}S_{v}-\rho F_{hd}I_{vz}-I_{vz}\mu_{v}},\\\
I_{vc}^{\prime}=\rho\left(F_{hz}I_{vd}+F_{hd}I_{vz}\right)-I_{{vc}}\mu_{v},\end{cases}$
(57)

which generalizes a bit [47] by introducing the parameter $\rho$, whose
purpose is to allow simplifying the model to remove the $I_{c},I_{vc}$
classes, by setting $\rho=0$.

Note that humans are born fully susceptible to dengue and Zika at a rate of
$\mu{N_{h}}$, where $\mu$ is the natural birth/death rate for humans and
$N_{h}$ is the total human population. Susceptible individuals can become
infected with dengue from either a dengue-infected ($I_{vd}$) or coinfected
female mosquito ($I_{vc}$). The mosquito-to-human dengue infection rate is
given by $\beta_{hd}$. This rate is modified by a factor of $\nu_{d}$ to
indicate the altered infectivity of coinfected mosquitoes. Once infected with
dengue, humans can recover or become co-infected with Zika (by a Zika-infected
($I_{vz}$) or coinfected female mosquito ($I_{vc}$), or by sexual transmission
from a Zika-infected ($I_{z}$) or coinfected ($I_{c}$) human) and transition
into the Rd or Ic class, respectively. In a similar manner, fully susceptible
humans become infected with Zika from a mosquito in the $I_{vz}$ or $I_{vc}$
compartment.

The DFE has only non-zero components $S_{v}=N_{v},S_{h}=N_{h}$. Choosing as
infectious set all the compartments except $S_{v},S_{h}$ yields

$R_{0}=\max[\sqrt{\frac{\beta_{{hd}}N_{v}\beta_{{vd}}}{N_{h}\mu_{v}(\gamma_{d}+\mu)}},\frac{\beta_{s}+\sqrt{\beta_{s}^{2}+\frac{4N_{v}\beta_{{hz}}\beta_{{vz}}\left(\mu+\gamma_{z}\right)}{N_{h}\mu_{v}}}}{2\left(\mu+\gamma_{z}\right)}]:=\max[\mathcal{R}_{d},\mathcal{R}_{z}],$
(58)

confirming [47, sec. 4], and also the multi-strain structure we already met in
(40), (48). Furthermore, one may show that
$\mathcal{R}_{d}>1,\mathcal{R}_{z}>1$ are necessary and sufficient conditions
for the existence of the dengue only and Zika only fixed points – see
subsequent sections.

We end this section by reporting on the Jacobian factorizations at $E_{0}$,
when choosing as infectious set

$\mathcal{I}=\left\\{I_{d},I_{z},I_{c},J_{d},J_{z},I_{{vd}},I_{{vz}},I_{{vc}}\right\\}.$

Now the characteristic polynomial has two second order factors:

  1. 1.

One of Descartes type which yields the polynomial
$R_{1}(X)=\frac{\beta_{{hd}}S_{v}\beta_{{vd}}\left(k_{d}R_{z}+S_{h}\right)}{N_{h}^{2}\mu_{v}\left(\gamma_{d}+\mu\right)},$
which generalizes $\mathcal{R}_{d}$, in the sense that
$R_{1}(E_{0})=\mathcal{R}_{d}^{2}$; this raises the question of whether this
is related to the Zika IRN.

  2. 2.

One not of Descartes type, which raises the question of how to exploit non
Descartes type second order factors.

####  8.2.1 The dengue only resident fixed point $E_{d}$

Even though the coordinates of the dengue only resident fixed point $E_{d}$
are pretty simple, obtaining them isn’t. We have an a priori choice of
zeroable set $in_{1^{\prime}}=\\{I_{z},R_{z},J_{z},I_{zv}\\}$ which turns out
to lead to about 2.5 hrs for "Solve" (due to the existence of 4 extra fixed
points, which are non-positive for the numeric values of [47]. After
performing the computation, it turns out that the full zeroable set is
$in_{1}=\\{I_{z},I_{c},R_{z},J_{d},J_{z},R,I_{vz},I_{vc}\\}$. The remaining
set of equations:

$\left(\begin{array}[]{cc}-\gamma_{d}I_{d}-\mu
I_{d}+\frac{I_{{vd}}S_{h}\beta_{{hd}}}{N_{h}}&=0\\\ \gamma_{d}I_{d}-\mu
R_{d}&=0\\\ \frac{I_{d}S_{v}\beta_{{vd}}}{N_{h}}-I_{{vd}}\mu_{v}&=0\\\
\mu\left(N_{h}-S_{h}\right)-\frac{I_{{vd}}S_{h}\beta_{{hd}}}{N_{h}}&=0\\\
-\frac{I_{d}S_{v}\beta_{{vd}}}{N_{h}}+N_{v}\mu_{v}-S_{v}\mu_{v}&=0\\\
\end{array}\right)$

may be easily solved. Besides the DFE, it has one extra fixed point:

$\displaystyle
R_{d}=\frac{\gamma_{d}I_{d}}{\mu},S_{v}=\frac{\mu_{v}N_{h}N_{v}}{\mu_{v}N_{h}+\beta_{vd}I_{d}}=\frac{\mu_{v}\left(\gamma_{d}+\mu\right)\left(\mu
N_{h}+\beta_{{hd}}N_{v}\right)}{\beta_{{hd}}\left[\mu_{v}\left(\gamma_{d}+\mu\right)+\mu\beta_{vd}\right]},$
$\displaystyle
S_{h}=\frac{N_{h}^{2}\left(\mu_{v}\left(\gamma_{d}+\mu\right)+\mu\beta_{{vd}}\right)}{\beta_{{vd}}\left(\mu
N_{h}+\beta_{{hd}}N_{v}\right)},I_{d}=\frac{\mu
N_{h}^{2}\mu_{v}}{\beta_{{vd}}\left(\mu
N_{h}+\beta_{{hd}}N_{v}\right)}\left(\frac{N_{v}\beta_{{hd}}\beta_{{vd}}}{N_{h}\mu_{v}\left(\gamma_{d}+\mu\right)}-1\right),$
$\displaystyle
I_{{dv}}=\beta_{vd}I_{d}\frac{S_{v}}{\mu_{v}N_{h}}=\frac{I_{d}N_{v}\beta_{{vd}}}{I_{d}\beta_{{vd}}+N_{h}\mu_{v}},I_{{vz}}=0,I_{vc}=0.$

The bifurcation value for $E_{d}$ is thus

$\frac{N_{v}\beta_{{hd}}\beta_{{vd}}}{N_{h}\mu_{v}\left(\gamma_{d}+\mu\right)}:=\mathcal{R}_{d}^{2},$

confirming [47, Lem. 1].

The Jacobian factorizations when choosing as infectious set the complement of
$I_{d},R_{d},I_{{vd}},S_{v},S_{h}$ has characteristic polynomial with one non-
Descartes type, third order factor.

####  8.2.2 The Zika only resident fixed point $E_{z}$

Using the full zeroable set given in [47]
$in_{2}=\\{I_{d},I_{c},R_{d},J_{d},J_{z},R,I_{dv},I_{vc}\\}$, yields the set
of equations:

$\left(\begin{array}[]{cc}S_{h}\left(\frac{\beta_{{hz}}I_{{vz}}}{N_{h}}+\frac{I_{z}\beta_{s}}{N_{h}}\right)-I_{z}\gamma_{z}-\mu
I_{z}&=0\\\ I_{z}\gamma_{z}-\mu R_{z}&=0\\\
\frac{I_{z}S_{v}\beta_{{vz}}}{N_{h}}-I_{{vz}}\mu_{v}&=0\\\
\mu\left(N_{h}-S_{h}\right)-S_{h}\left(\frac{\beta_{{hz}}I_{{vz}}}{N_{h}}+\frac{I_{z}\beta_{s}}{N_{h}}\right)&=0\\\
-\frac{I_{z}S_{v}\beta_{{vz}}}{N_{h}}+N_{v}\mu_{v}-S_{v}\mu_{v}&=0\\\
\end{array}\right)$

The Zika only resident fixed point $E_{z}$ satisfies

$\displaystyle
R_{z}=\frac{\gamma_{z}I_{z}}{\mu},S_{v}=\frac{\mu_{v}N_{h}N_{v}}{\mu_{v}N_{h}+\beta_{vz}I_{z}},I_{{dv}}=\beta_{vz}I_{z}\frac{S_{v}}{\mu_{v}N_{h}}=\frac{I_{z}\beta_{{vz}}N_{v}}{\mu_{v}N_{h}+\beta_{vz}I_{z}},I_{{vz}}=0,I_{vc}=0,$
$\displaystyle S_{h}=\frac{\mu
N_{h}^{2}\left(N_{h}\mu_{v}+I_{z}\beta_{{vz}}\right)}{I_{z}\beta_{{vz}}\left(\mu
N_{h}+\beta_{{hz}}N_{v}+I_{z}\beta_{s}\right)+N_{h}\mu_{v}\left(\mu
N_{h}+I_{z}\beta_{s}\right)}=\frac{N_{h}^{2}\left(\mu_{v}\left(\gamma_{d}+\mu\right)+\mu\beta_{{vd}}\right)}{\beta_{{vd}}\left(\mu
N_{h}+\beta_{{hd}}N_{v}\right)},$

where $I_{z}$ is a positive root of the quadratic equation
$aI_{z}^{2}+bI_{z}+c=0$, with coefficients:

$\begin{cases}c=\mu
N_{h}\left(N_{h}\mu_{v}\left(\mu-\beta_{s}+\gamma_{z}\right)-\beta_{{hz}}N_{v}\beta_{{vz}}\right),\\\
b=N_{h}\beta_{s}\mu_{v}\left(\mu+\gamma_{z}\right)+\mu
N_{h}\beta_{{vz}}\left(\mu-\beta_{s}+\gamma_{z}\right)+\beta_{{hz}}N_{v}\beta_{{vz}}\left(\mu+\gamma_{z}\right),\\\
a=\beta_{s}\beta_{{vz}}\left(\mu+\gamma_{z}\right)\end{cases}$

Assume first that $\beta_{s}$ is small enough so that $b>0$; then, this
equation has a unique positive root iff $c<0$, which may be written also as

$\frac{N_{h}\beta_{s}\mu_{v}+\beta_{{hz}}N_{v}\beta_{{vz}}}{N_{h}\mu_{v}(\mu+\gamma_{z})}>1.$
(59)

It is shown in [47, Thm 1] that this is equivalent to $R_{z}>1$ (both
conditions determine the correct stability domain, and both reduce when
$\beta_{s}=0$ to the same answer
$\frac{\beta_{{hz}}N_{v}\beta_{{vz}}}{N_{h}\mu_{v}(\mu+\gamma_{z})}$).

The model of [47] contains several interesting particular cases, to which we
turn next.

####  8.2.3 The dengue invasion reproduction number (IRN) and two possible
$(F,V)$ decompositions

The dengue fixed point has non-zero values $S_{h},S_{v},I_{d},R_{d},I_{dv}$.
Computing the IRN’s requires specifying the "invasion infection classes". [47]
work with a subset of

$in_{2^{\prime}}=\\{I_{d},I_{c},R_{d},J_{d},J_{z},I_{dv},I_{cv},R_{c}\\},$

given by $in_{2}=\\{I_{d},I_{c},J_{d},I_{dv},I_{cv}\\}$.

The resulting recipe $V$ matrix is diagonal, and the recipe $F$ matrix, after
denoting proportions by minuscule letters, is:

$F=\left(\begin{array}[]{ccccc}0&0&0&{s_{h}\beta_{{hd}}}&{\nu_{d}s_{h}\beta_{{hd}}}\\\
{\rho\left(\beta_{{hz}}i_{{zv}}+i_{z}\beta_{s}\right)}&0&0&{\rho\beta_{{hd}}i_{z}}&{\rho\nu_{d}\beta_{{hd}}i_{z}}\\\
0&0&0&{\beta_{{hd}}r_{z}}&{\nu_{d}\beta_{{hd}}r_{z}}\\\
{s_{v}\beta_{{vd}}}&{s_{v}\beta_{{vd}}}&{k_{d}s_{v}\beta_{{vd}}}&0&0\\\
{i_{{zv}}\beta_{{vd}}}&{i_{{zv}}\beta_{{vd}}}&{k_{d}i_{{zv}}\beta_{{vd}}}&{i_{z}\beta_{{vz}}}&0\\\
\end{array}\right)$ (60)

and the spectral radius of the resulting recipe $K$ matrix satisfies a
polynomial equation of degree $4$.

Now [47, Sec 5.1] move two of the $F$ terms in the $V$ matrix, yielding

$F=\left(\begin{array}[]{ccccc}0&0&0&{s_{h}\beta_{{hd}}}&{\nu_{d}s_{h}\beta_{{hd}}}\\\
0&0&0&{\rho\beta_{{hd}}i_{z}}&{\rho\nu_{d}\beta_{{hd}}i_{z}}\\\
0&0&0&{\beta_{{hd}}r_{z}}&{\nu_{d}\beta_{{hd}}r_{z}}\\\
{s_{v}\beta_{{vd}}}&{s_{v}\beta_{{vd}}}&{k_{d}s_{v}\beta_{{vd}}}&0&0\\\
{i_{{zv}}\beta_{{vd}}}&{i_{{zv}}\beta_{{vd}}}&{k_{d}i_{{zv}}\beta_{{vd}}}&0&0\\\
\end{array}\right),$ (61)

with the $-V$ matrix being:

$\left(\begin{array}[]{ccccc}-\gamma_{d}-\mu-{\rho\left(\beta_{{hz}}i_{{zv}}+i_{z}\beta_{s}\right)}&{\rho\left(\beta_{{hz}}i_{{zv}}+i_{z}\beta_{s}\right)}&0&0&0\\\
0&-\rho\left(\gamma_{d}+\gamma_{z}\right)-\mu&\rho\gamma_{z}&0&0\\\
0&0&-\gamma_{d}-\mu&0&0\\\
0&0&0&-{i_{z}\beta_{{vz}}}-\mu_{v}&{i_{z}\beta_{{vz}}}\\\ 0&0&0&0&-\mu_{v}\\\
\end{array}\right)$ (62)

They reduce thus the rank of $K$ to $2$ and getting a simpler $R_{0}$. On the
other hand, their decomposition is admissible only under extra conditions of
the parameters which ensure the non-positivity of the row-sums of $-V$, which
they omit to mention.

###### Remark 22

The associated CTMC is the union of two disjoint generalized Erlangs, on the
host and vector, respectively. These are employed in the
probabilistic/epidemic interpretations in [47].

The probabilistic/epidemic significance of $F$ is better understood after
decomposing this matrix as a sum of matrices of rank $1$ as follows:

$F=\begin{pmatrix}{\beta_{{hd}}}s_{h}\\\ {\rho\beta_{{hd}}i_{z}}\\\
{\beta_{{hd}}k_{d}r_{z}}\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}0&0&0&1&{\nu_{d}}\end{pmatrix}+\begin{pmatrix}0\\\
0\\\ 0\\\ {\beta_{{vd}}}s_{v}\\\
{\beta_{{vd}}}i_{{zv}}\end{pmatrix}\begin{pmatrix}1&1&1&0&0\end{pmatrix}.$
(63)

The column vector are total infectivity rates for the resident compartments,
the row vectors are distribution vectors, and this decomposition yields
immediately both the Diekmann kernel and $R_{0}$ – see [11, 12].

##  9 Conclusions

The possible non-uniqueness of the NGM matrix has not been sufficiently
studied in the literature. Sometimes, like in the example of the last section,
one simplifying choice is justified a posteriori on the grounds of some
interpretability of the results, ignoring the fact that other choices might
lead to even simpler answers, and the fact that a priori there is no reason to
expect simple answers.

We answer to this classic dilemma by showing via numerous examples that the
first “recipe NGM" to come to mind leads quickly to most of the results found
in the literature. The question of whether our recipe may always be associated
to admissible equation decompositions remains open.

We have also examined a variant of the Jacobian approach, a "factorization
Jacobian approach", which draws the attention to certain polynomials with
interesting properties (42), and raises interesting questions – see especially
Section 4.2, Open Question 3. Notably, the relation (40) holds in all the
three "multi-strain" examples we examined, and raises the additional question
of how to define multi-strain models in terms of the dynamical system, to
ensure that this always holds for this class.

##  10 Appendix: the implementation of the Jacobian factorization approach

First, we use a utility which, for a given model, infectious set, and dummy
variable (taken always as $u$, to avoid confusions), outputs the Jacobian at
the DFE, the trace and determinant (for other purposes), the characteristic
polynomial in $u$, the NGM matrix and $R_{F}$.

    JR0[mod_,inf_,u_,cn_:{}]:=
      Module[{dyn,X,par,cinf,cp,cX,jac,tr,det,chp,ngm,K,R0},
        dyn=mod[[1]];X=mod[[2]];par=mod[[3]];
        Print[" dyn=",dyn//FullSimplify//MatrixForm,X,par];
        cinf=Thread[X[[inf]]->0];
        cp=Thread[par>0];cX=Thread[X>0];
        cdfe=Join[DFE[mod,inf],cinf];
        jac=Grad[dyn,X]/.cinf/.cn;
        tr=Tr[jac];
        det=Det[jac];
        chp=CharacteristicPolynomial[jac,u];
        ngm=NGM[mod,inf];
        K=ngm[[6]];
       Print["K=",K//MatrixForm];
       R0=Assuming[Join[cp,cX],Max[Eigenvalues[K]]];
      {chp,R0,K,jac,tr,det}];

Most of the work is done after calling this utility, by another one, JR02.This
splitting of JR0 in two parts is necessary since the detection of the non-sign
definite factors which must be analyzed is easier to perform by eye, than to
program. The JR02 script is:

    JR02[pol_,u_]:=Module[{co,co1,cop,con,R_J},co=CoefficientList[pol,u];
      Print["the  factor ",pol," has degree ",Length[co]-1];
      co1=Expand[co[[1]]* co[[Length[co]]]];
      Print["its leading * constant coefficient product is ",co1];
      cop=Replace[co1, _. _?Negative -> 0, {1}](*level 1 here ?*);
      con=cop-co1;
      Print["R_J is"];
      R_J=con/cop//FullSimplify;
    {R_J,co}
    ]

For a specific “mod", both $R_{0}$’s may be obtained by typing:

    jr = JR0[mod, inf, u];
    chp = jr[[1]] // Factor
    Print["factor is ", pol = chp[[5]]]
    pc = JR02[pol,
       u];(*the script JR02 determines R_J, using the index,
       for example 5, determined by \eye inspection in the previous command*)
    Print["R_J is ", R_J = pc[[1]] // FullSimplify]
    Print["R_N is ", R_N = jr[[2]] // FullSimplify]

### Proof of [57]’s result via Mathematica

  1. 1.

The solution of the first recurrence equation in (7) for the expected time to
extinction of a linear birth and death process with arrival rate $A$ and death
rate $qA$ (relevant when $R_{0}<1$), via Mathematica is:

$\displaystyle\scalebox{0.85}{$\frac{q\left(H_{K}\left(1-q^{j}\right)+H_{j}\left(q^{K}-1\right)+\log\left(\frac{q-1}{q}\right)\left(q^{K}-q^{j}\right)\right)-\left(\left(q^{j}-1\right)\Phi\left(\frac{1}{q},1,K+1\right)\right)+\left(q^{K}-1\right)\Phi\left(\frac{1}{q},1,j+1\right)}{A(q-1)q\left(q^{K}-1\right)}$},$

where $H$ denotes the Harmonic function.

Since Mathematica cannot compute the limit when $K$ converges to infinity
directly, we break the limit into its three parts, and end up with the
following generalization:

Making now $j=1$ yields [57]’s result which is

$\frac{\log(q)-\log(q-1)}{A}.$

  2. 2.

When $R_{0}>1$ we cannot obtain the limit for general $j$. When $j=1$,
similarly with the previous case, the limit is devided into four parts:

Here Mathematica yields that $a_{1}=0$, $a_{2}=-\frac{\log(1-q)}{A}$, the
second being precisely Whittle’s result, but we were unable to confirm with
Mathematica that $a_{3}=a_{4}=0$.

Acknowledgements. We thank Andrew Brouwer, Corey Shanbrom, Matija Vidmar, Dan
Goreac, Andrei Halanay and James Watmough for useful exchanges.

## References

  * [1] Symmetries and normalization in 3-compartment epidemic models. i: The replacement number dynamics.  arXiv:2301.00159, 2022.
  * [2] Rim Adenane, Eric Avila-Vales, Florin Avram, Andrei Halanay, and Angel GC Pérez.  On a three-dimensional and two four-dimensional oncolytic viro-therapy models.  arXiv preprint arXiv:2210.00401, 2022.
  * [3] Ma ́\mathsf{i}ra Aguiar, Vizda Anam, Konstantin B Blyuss, Carlo Delfin S Estadilla, Bruno V Guerrero, Damián Knopoff, Bob W Kooi, Akhil Kumar Srivastav, Vanessa Steindorf, and Nico Stollenwerk.  Mathematical models for dengue fever epidemiology: A 10-year systematic review.  Physics of Life Reviews, 40:65–92, 2022.