# Transmission of multiple pathogens across species

Clotilde Djuikem  Julien Arino

###### Abstract

We analyse a model that describes the propagation of many pathogens within and
between many species. A branching process approximation is used to compute the
probability of disease outbreaks. Special cases of aquatic environments with
two host species and one or two pathogens are considered both analytically and
computationally.

Keywords: multiple species–multiple pathogens, branching process
approximation, introductions

##  1 Introduction

The ranges of species are continuously changing [23, 39]. However, the process
has accelerated in recent years because of climate change [9, 35]. Regardless
of what is driving their evolution, a consequence of the modification of
ranges is more frequent interactions between species that did not use to
interact or interacted quite infrequently.

This has a wide variety of consequences. Competition for resources is
modified, if an invading species is, for instance, using the same resource as
a a resident one. This is thought to be one of the main drivers of species
evolution [36]. Range shifting can also lead to the introduction into
ecosystems of pathogens from which they were absent, when species whose range
now includes these ecosystems become more frequent there [13]. Introductions
of pathogens due to range shifting is also very similar to what happens when
human populations encroach into the ranges of species, which has led to an
increasing number of spillover events.

In both cases, some of the populations involved may be hosts to a wide variety
of pathogens. Understanding a situation with different pathogens and different
species is therefore important.

The specific motivation for the present work comes from the observation that
salmonids are observed increasingly frequently the Mackenzie River, in the
western Canadian arctic. Of interest to collaborators from Fisheries and
Oceans Canada (see _Acknowledgements_) is the fact that these vagrant salmon
species spend most of their lives in distant ecosystems, where they can
acquire pathogens that are, to this point, mostly absent from the Mackenzie
River aquatic ecosystem. When they are collocated in that ecosystem, those
vagrant species can in turn transmit those novel pathogens to resident
species.

While pathogens abound in terrestrial ecosystems, the situation is even more
pronounced in aquatic ecosystems, where numerous pathogens are present [11,
45, 46]. Viruses, for instance, are estimated to be the most abundant
“lifeforms” in the oceans, representing over 90% of the nucleic-acid-
containing particles and about 5% of the biomass there [42]. Many aquatic
pathogens infect fish species, so that the invading species mentioned earlier
may be coming into contact with a wide variety of pathogens prior to their
entering a novel ecosystem.

Our aim is therefore to establish models to help understand the introduction
of pathogens in species from which they were absent up to that point, when
these species come into contact with other species potentially bearing the
pathogen. The model in this paper is a simplified model and serves to set the
general setting in which we operate. We use a simple SLIR model, whose
dynamics in a single location and single population is well understood, but
assume that there are multiple species of hosts as well as multiple pathogen
species. We also assume that there is no coinfection with multiple pathogens.

The article is organised as follows. In Section 2, an ODE multi-species
epidemiological model is introduced, followed by its continuous time Markov
chain (CTMC) equivalent. The section also presents an analysis of both the ODE
and CTMC models, the latter using a branching process approximation to compute
the probability of a disease outbreak. Section 3 focuses on the case of two
species and one pathogen, with example scenarios corresponding to three
different fish viruses investigated numerically, as well as a particular case
focusing on introduction of a pathogen by a species in which it is endemic. In
Section 4, the case of two species and two pathogens is discussed.

##  2 The general model

Consider $P$ populations. These populations could be the same or different
species, the important feature being that they be distinguishable according to
some criterion. In the sequel, we use both terms, _species_ and _population_ ,
interchangeably. Within and between these populations, $V$ pathogens can
propagate. Specifically, each population is described by an SLIR epidemic
model, where susceptible individuals in a given population can become infected
by any of the $V$ viruses if they come into contact with an individual
infected by it, regardless of the population that individual belongs to.
Further, we assume that coinfection does not occur, i.e., once infected by a
given pathogen, an individual cannot acquire infection from another pathogen.

###  2.1 Formulation of the deterministic model

For $p=1,\ldots,P$, denote $S_{p}$ the number of individuals susceptible to
infection in the $p$th population. Individuals of species $p$ may become
infected by any of the $V$ pathogens present when they come across an
individual infectious with that pathogen. The parameter describing the rate at
which contacts between a susceptible from species $p$ and an infectious
individual from species $q=1,\ldots,P$ infected with pathogen $v=1,\ldots,V$,
results in new infections, is $\beta_{pqv}$, i.e., in words, the parameter
$\beta$ has indices

$\beta_{\text{who becomes infected, who infects, with which pathogen}}.$

Upon infection, individuals of species $p=1,\ldots,P$ infected by pathogen
$v=1,\ldots,V$ become latently infected, with numbers denoted $L_{pv}$. We do
not consider coinfection with multiple pathogens; once an individual is
contaminated with any of the viruses, they cannot become infected by any other
pathogen.

After an incubation period of mean duration $1/\varepsilon_{pv}$ time units,
individuals of species $p$ infected with pathogen $v$ become infectious to
others. The number of such infectious individuals is denoted $I_{pv}$.
Finally, after an average $1/\gamma_{pv}$ time units spent infectious with the
pathogen, individuals recover and move to the $R_{p}$ compartment. At this
point, the pathogen they were infected with is ignored as it is not relevant
to the problem under consideration. Regarding species demography, birth into
population $p=1,\ldots,P$ occurs at the fixed rate $b_{p}$, while death occurs
in all compartments at the _per capita_ rate $d_{p}$.

Taking all this into account, we have a group model, with dynamics of the
different states governed for population $p=1,\ldots P$ and pathogen
$v=1,\ldots,V$ by the following ordinary differential equations:

$\displaystyle\dot{S}_{p}=b_{p}-\left(\sum_{q=1}^{P}\sum_{v=1}^{V}\beta_{pqv}I_{qv}+d_{p}\right)S_{p},$
(1a)
$\displaystyle\dot{L}_{pv}=\sum_{q=1}^{P}\beta_{pqv}I_{qv}S_{p}-(\varepsilon_{pv}+d_{p})L_{pv},$
(1b)
$\displaystyle\dot{I}_{pv}=\varepsilon_{pv}L_{pv}-(\gamma_{pv}+d_{p})I_{pv},$
(1c) $\displaystyle\dot{R}_{p}=\sum_{v=1}^{V}\gamma_{pv}I_{pv}-d_{p}R_{p}.$
(1d)

System (1) is considered with nonnegative initial conditions. To avoid a
trivial case, it is assumed that $L_{pv}+I_{pv}>0$ for at least one
$(p,v)\in\\{1,\ldots,P\\}\times\\{1,\ldots,V\\}$. The total size of each
population $p$ for $i=1,\ldots,P$ is given by:

$N_{p}=S_{p}+\sum_{v=1}^{V}L_{pv}+\sum_{v=1}^{V}I_{pv}+R_{p}.$ (2)

Variable | Meaning
---|---
$S_{p}$ | Susceptible individual in population $p$
$L_{pv}$ | Latent individual in population $p$ infected by virus $v$
$I_{pv}$ | Individual in population $p$ infectious with virus $v$
$R_{p}$ | Recovered individual in population $p$
Table 1: State variables and their meaning.

###  2.2 Notation

Equations (1b) and (1c) involve two different indices. Analysis of the system
often requires to list these indices. To simplify presentation, for given
symbols $X$ and $Y$, we use the notation

$\left\llbracket X_{p}\right\rrbracket=X_{1},\ldots,X_{P},$ $\left\llbracket
X_{pv}\right\rrbracket=X_{11},X_{12},\ldots,X_{1V},X_{21},X_{22},\ldots,X_{2V},\ldots,X_{P1},X_{P2},\ldots,X_{PV}$

and

$\left\llbracket
X_{pv}+Y_{p}\right\rrbracket=X_{11}+Y_{1},X_{12}+Y_{1},\ldots,X_{1V}+Y_{1},\\\
X_{21}+Y_{2},X_{22}+Y_{2},\ldots,X_{2V}+Y_{2},\ldots,\\\
X_{P1}+Y_{P},X_{P2}+Y_{P},\ldots,X_{PV}+Y_{P}.$

Thus, when multiple indices are present, we present indices as would the row-
first enumeration of indices of the entries of a $P\times V$ matrix. Note that
the assumption is that indices $p$ and $v$ are reserved, respectively, for
population and virus species indices and therefore run in $1,\ldots,P$ and
$1,\ldots,V$.

###  2.3 Basic analysis of the deterministic model

The disease-free equilibrium (DFE) of system (1) is

$\bm{E}_{0}^{\eqref{sys:ODE-multi-species}}=\left(\left\llbracket
S_{p}^{0}\right\rrbracket,0_{\mathbb{R}^{P(2V+1)}}\right),$ (3)

where $S_{p}^{0}=b_{p}/d_{p}$ for $p=1,\ldots,P$.

To determine the matrices used in the computation of the basic reproduction
number using the next generation matrix method of [44], order infected
variables by type: $\left\llbracket L_{pv}\right\rrbracket,\left\llbracket
I_{pv}\right\rrbracket$. Then the non-negative $2PV\times 2PV$-matrix
$\mathbf{G}$ has block form

$\mathbf{G}=\begin{bmatrix}0&\mathbf{G}_{12}\\\ 0&0\\\ \end{bmatrix},$ (4)

where the $PV\times PV$-matrix $\mathbf{G}_{12}$ is itself a block matrix,
with each $V\times V$ sized block taking the form, for
$p,q\in\\{1,\ldots,P\\}$,

$\mathbf{G}_{12}^{pq}=S_{p}^{0}\mathsf{diag}(\beta_{pq1},\ldots,\beta_{pqV}).$
(5)

The matrix $\mathbf{W}$ is a non-negative $2PV\times 2PV$-matrix and has block
form

$\mathbf{W}=\begin{bmatrix}\mathbf{W}_{11}&\mathbf{0}\\\
-\mathbf{W}_{21}&\mathbf{W}_{22}\\\ \end{bmatrix},$ (6)

with the $PV\times PV$-sized block taking the form

$\mathbf{W}_{11}=\mathsf{diag}(\llbracket\varepsilon_{pv}+d_{p}\rrbracket),\mathbf{W}_{21}=\mathsf{diag}(\llbracket\varepsilon_{pv}\rrbracket)\text{
and }\mathbf{W}_{22}=\mathsf{diag}(\llbracket\gamma_{pv}+d_{p}\rrbracket).$
(7)

The basic reproduction number of (1) is then the spectral radius of
$\mathbf{GW}^{-1}$ and is given by $\mathcal{R}_{0}^{\eqref{sys:ODE-multi-
species}}=\rho(\mathbf{GW^{-1}})$. Since $\mathbf{W}$ is block lower
triangular,

$\mathbf{W^{-1}}=\begin{bmatrix}\mathbf{W}_{11}^{-1}&\mathbf{0}\\\
\mathbf{W}_{22}^{-1}\mathbf{W}_{21}\mathbf{W}_{11}^{-1}&\mathbf{W}_{22}^{-1}\end{bmatrix},$

whence, from the form of $\mathbf{G}$, we obtain

$\mathcal{R}_{0}^{\eqref{sys:ODE-multi-
species}}=\rho(\mathbf{G}_{12}\mathbf{W}_{22}^{-1}\mathbf{W}_{21}\mathbf{W}_{11}^{-1}).$

The matrices (7) are diagonal and the structure of $\mathbf{G}_{12}$ is
relatively simple; it is therefore possible to further simplify the expression
above. We obtain

$\mathcal{R}_{0}^{\eqref{sys:ODE-multi-
species}}=\max_{v=1,\ldots,V}\mathcal{R}_{0}^{v},$ (8)

where

$\mathcal{R}_{0}^{v}=\rho(\mathcal{B}_{v}),$

is the basic reproduction number of virus $v=1,\ldots,V$ and $\mathcal{B}_{v}$
is a $P\times P$-matrix defined as

$\mathcal{B}_{v}=[\mathcal{B}_{pqv}]_{pq},$ (9)

where $\mathcal{B}_{pqv}$ denotes the element in the $p$-th row and $q$-th
column of the matrix $\mathcal{B}_{v}$ and is given by

$\mathcal{B}_{pqv}=\frac{\beta_{pqv}\varepsilon_{qv}S_{p}^{0}}{(\varepsilon_{qv}+d_{q})(\gamma_{qv}+d_{q})},$

for $p,q=1,\ldots,P,\;v=1,\ldots,V$.

From [44, Theorem 2], we deduce the following result concerning the local
asymptotic stability of the disease-free equilibrium
$\bm{E}_{0}^{\eqref{sys:ODE-multi-species}}$.

######  Lemma 1.

The disease-free equilibrium $\bm{E}_{0}^{\eqref{sys:ODE-multi-species}}$ of
(1) is locally asymptotically stable if $\mathcal{R}_{0}^{\eqref{sys:ODE-
multi-species}}<1$ and unstable if $\mathcal{R}_{0}^{\eqref{sys:ODE-multi-
species}}>1$.

Remark that in the absence of interaction between the populations, i.e., when
$\beta_{p\ell v}=0$ if $p\neq\ell$, the basic reproduction number in each
population $p=1,\ldots,P$ is given by

$\mathcal{R}_{0p}=\max_{v=1,\ldots,V}(\mathcal{R}_{0p}^{v}),\text{ where
}\mathcal{R}_{0p}^{v}=\frac{\beta_{ppv}\varepsilon_{pv}S_{p}^{0}}{(\varepsilon_{pv}+d_{p})(\gamma_{pv}+d_{p})}.$
(10)

Various forms of the reproduction numbers appear in the remainder of the text.
To clarify, we list those here.

  * •

$\mathcal{R}_{0p}$ denotes the basic reproduction number of species $p$ in the
presence of multiple pathogens $V$, excluding other species.

  * •

$\mathcal{R}_{0p}^{v}$ denotes the basic reproduction number of species $p$ in
the presence of a single pathogen $v$, excluding other species.

  * •

$\mathcal{R}_{0}^{v}$ denotes the basic reproduction number of a single
pathogen $v$ across $P$ interacting species.

  * •

$\mathcal{R}_{0}$ denotes the basic reproduction number of multiple pathogens
$V$ across multiple interacting species $P$, i.e, for the entire System (1).

Asymptotic stability is in fact global for a given pathogen when
$\mathcal{R}_{0v}<1$, as established in the following theorem.

######  Theorem 2.

If $\mathcal{R}_{0v}^{\eqref{sys:ODE-multi-species}}<1$, then the DFE
$\bm{E}_{0v}^{\eqref{sys:ODE-multi-species}}$ is globally asymptotically
stable (GAS) in $\Omega$, where

$\Omega=\Biggl{\\{}(\llbracket S_{p}\rrbracket,\llbracket
L_{pv}\rrbracket,\llbracket I_{pv}\rrbracket,\llbracket
R_{p}\rrbracket)\in\mathbb{R}^{2P(V+1)}:\\\
N_{p}=S_{P}+\sum_{v=1}^{V}(L_{pv}+I_{pv})+R_{p}\leq\frac{b_{p}}{d_{p}};\quad
p=1,\ldots,P\Biggr{\\}}.$

This result is proved in Appendix B.

###  2.4 The continuous time Markov chain model

By using the continuous-time Markov chain (CTMC) model related to the
deterministic model (1), we gain a deeper understanding of complex disease
dynamics and explore various scenarios that the deterministic models cannot
capture. In particular, the CTMC model allows tracking discrete counts of
individuals, allowing transitions toward a state where the disease is
eradicated. This property does not exist in ordinary differential equation
(ODE) models, where such states are typically approached only as limits.

The CTMC model related to the deterministic model (1),

$\mathbf{X}_{t}=\left(\llbracket S_{p}(t)\rrbracket,\llbracket
L_{pv}(t)\rrbracket,\llbracket I_{pv}(t)\rrbracket,\llbracket
R_{p}(t)\rrbracket\right),\quad t\in\mathbb{R}_{+},$ (11a) is characterised by
transition probabilities from state $\mathbf{k}$ to state $\mathbf{j}$,
$\mathbb{P}(\mathbf{X}(t+\Delta
t)=\mathbf{j}\mid\mathbf{X}(t)=\mathbf{k})=\sigma(\mathbf{k},\mathbf{j}),$
(11b)

with transition rates $\sigma(\mathbf{k},\mathbf{j})$ given in Table 2.

Event | Transition | Transition rate
---|---|---
Birth of $S_{p}$ | $S_{p}\rightarrow S_{p}+1$ | $b_{p}$
Natural death of $S_{p}$ | $S_{p}\rightarrow S_{p}-1$ | $d_{p}S_{p}$
Natural death of $L_{pv}$ | $L_{pv}\rightarrow L_{pv}-1$ | $d_{p}L_{pv}$
Natural death of $I_{pv}$ | $I_{pv}\rightarrow I_{pv}-1$ | $d_{p}I_{pv}$
Natural death of $R_{p}$ | $R_{p}\rightarrow R_{p}-1$ | $d_{p}R_{p}$
Infection of $S_{p}$ by $I_{qv}$ | $S_{p}\rightarrow S_{p}-1,L_{pv}\rightarrow L_{pv}+1$ | $\beta_{pqv}I_{qv}S_{p}$
End of incubation of $L_{pv}$ | $L_{pv}\rightarrow L_{pv}-1,I_{pv}\rightarrow I_{pv}+1$ | $\varepsilon_{pv}L_{pv}$
Recovery of $I_{pv}$ | $I_{pv}\rightarrow I_{pv}+1,R_{p}\rightarrow R_{p}+1$ | $\gamma_{pv}I_{pv}$
Table 2: Events, transitions $\mathbf{k}\to\mathbf{j}$ and transition rates
$\sigma(\mathbf{k},\mathbf{j})$ of the general CTMC model (11).

###  2.5 Branching process approximation of the CTMC

Branching process approximations (BPA) simplifies the analysis of CTMC by
approximating them with simpler branching processes. By considering only the
number of offspring generated by each state transition, BPA reduce the
complexity associated with modelling and analysing CTMC. We use multitype
branching process approximation (MBPA) to approximate the CTMC (11) near the
disease-free equilibrium, where the term _multitype_ is used to indicate that
there are different types of offspring.

The reasoning below follows [2]. Let $Z=(Z^{L},Z^{I})$ be a MBPA of the CTMC
$\mathbf{X}_{t}$ and $\mathbf{u}=(\llbracket
u^{\ell}_{pv}\rrbracket,\llbracket u^{i}_{pv}\rrbracket)\in\mathbb{R}^{2PV}$.
Define the offspring probability-generating function (p.g.f.)

$\mathbf{F}(\mathbf{u})=\left(\llbracket
f_{pv}^{\ell}(\mathbf{u})\rrbracket,\llbracket
f_{pv}^{i}(\mathbf{u})\rrbracket\right),$ (12a) where $\displaystyle
f_{pv}^{\ell}(\mathbf{u})=\frac{\varepsilon_{pv}u_{pv}^{i}+d_{p}}{\varepsilon_{pv}+d_{p}},$
(12b) $\displaystyle
f_{pv}^{i}(\mathbf{u})=\frac{\left(\displaystyle\sum_{q=1}^{P}\beta_{pqv}S_{q}^{0}u_{qv}^{\ell}\right)u_{pv}^{i}+\gamma_{pv}+d_{p}}{\Lambda_{pv}},$
(12c)

with

$\Lambda_{pv}=\sum_{q=1}^{P}\beta_{pqv}S_{q}^{0}+\gamma_{pv}+d_{p}.$

The following result then holds.

######  Theorem 3.

The probability of extinction in the multitype branching process with
generating function (12) is given by

$\displaystyle\mathbb{P}_{\text{ext}}$
$\displaystyle=\prod_{p=1}^{P}\prod_{v=1}^{V}\left(\frac{\varepsilon_{pv}z_{pv}^{i}+d_{p}}{\varepsilon_{pv}+d_{p}}\right)^{\ell_{pv0}}(z_{pv}^{i})^{i_{pv0}},$
(13a) $\displaystyle\mathbb{P}_{\text{outbreak}}$
$\displaystyle=1-\mathbb{P}_{\text{ext}},$ (13b)

where

$\bm{z}:=\left(\llbracket z_{pv}^{\ell}\rrbracket,\llbracket
z_{pv}^{i}\rrbracket\right)$

is a fixed point on $\displaystyle[0,1]^{2PV}$ of the p.g.f. (12) and
$\llbracket L_{pv}(0)\rrbracket=\llbracket\ell_{pv0}\rrbracket$, $\llbracket
I_{pv}(0)\rrbracket=\llbracket i_{pv0}\rrbracket$ is the initial condition.
The following alternative holds:

  * •

if $\mathcal{R}_{0}^{\eqref{sys:ODE-multi-species}}\leq 1$, then
$\bm{z}=\mathbf{1}$, i.e., $\mathbb{P}_{\text{ext}}=1$;

  * •

if $\mathcal{R}_{0}^{\eqref{sys:ODE-multi-species}}>1$, then additionally to
$\bm{z}=\bm{1}$, there is a unique vector $\mathbf{0}<\bm{z}<\mathbf{1}$ such
that $\mathbf{F}(\bm{z})=\bm{z}$.

###### Proof.

In the offspring probability-generating function (12), $f_{pv}^{\ell}$ and
$f_{pv}^{i}$ are, for $p=1,\ldots,P$ and $v=1,\ldots,V$, the offspring p.g.f.
associated to latent ($\ell$) and infected ($i$) individuals, respectively.
Let $r^{\ell}_{pv}$ and $r^{i}_{pv}$ be the number of offspring born from a
parent who is $\ell$ and $i$, respectively. Then p.g.f. (12b) and (12c) take
the form, for $k\in\\{\ell,i\\}$,

$f^{k}_{pv}(\mathbf{u})=\sum_{\llbracket r^{\ell}_{pv}\rrbracket,\llbracket
r^{i}_{pv}\rrbracket=0}^{\infty}\mathbb{P}(\llbracket r_{pv}^{k}\rrbracket)\
(u_{11}^{\ell})^{r^{\ell}_{11}}\cdots\ (u_{PV}^{\ell})^{r^{\ell}_{PV}}\
(u_{11}^{i})^{r^{i}_{11}}\ \cdots\ (u_{PV}^{i})^{r^{i}_{PV}},$ (14)

where $\mathbb{P}(\llbracket r_{pv}\rrbracket)=\mathbb{P}(\llbracket
r^{\ell}_{pv}\rrbracket,\llbracket r^{i}_{pv}\rrbracket)$ is the probability
that an individual of type $pv$ gives birth to $r^{k}_{j}$ individuals of type
$j$, for $j\in(\llbracket\ell_{pv}\rrbracket,\llbracket i_{pv}\rrbracket)$ and
$k\in\\{\ell,i\\}$.

These probabilities are computed using the infinitesimal transition
probabilities of the CTMC (11) (Table 2) and give the expressions (12b) and
(12c).

The probability of disease extinction or disease outbreak relative to the CTMC
$\mathbf{X}_{t}$ can then be approximated by the extinction probability of the
MBPA with generating function (12). For a small number of latent and
infectious individuals and under the assumptions above about the branching
process approximation, the probabilities of disease extinction and disease
outbreak, given $\llbracket
L_{pv}(0)\rrbracket=\llbracket\ell_{pv0}\rrbracket$, $\llbracket
I_{pv}(0)\rrbracket=\llbracket i_{pv0}\rrbracket$ and
$\mathcal{R}_{0}^{\eqref{sys:ODE-multi-species}}>1$, are

$\displaystyle\mathbb{P}_{\text{ext}}$
$\displaystyle=\prod_{p=1}^{P}\prod_{v=1}^{V}(z_{pv}^{\ell})^{\ell_{pv0}}(z_{pv}^{i})^{i_{pv0}}$
(15a) $\displaystyle\mathbb{P}_{\text{outbreak}}$
$\displaystyle=1-\mathbb{P}_{\text{ext}},$ (15b)

where $\bm{z}:=\left(\llbracket z_{pv}^{\ell}\rrbracket,\llbracket
z_{pv}^{i}\rrbracket\right)$ is a fixed point of the p.g.f. (12) on
$\displaystyle[0,1]^{2PV}$ such that $\mathbf{0}\leq\bm{z}<\mathbf{1}$,
denoting $\mathbf{0}=\mathbf{0}_{2PV}$ and $\mathbf{1}$ the $2PV$ vector of
all ones. It is clear that (12b) implies that
$z_{pv}^{\ell}=(\varepsilon_{pv}z_{pv}^{i}+d_{p})/(\varepsilon_{pv}+d_{p})$,
where $z_{pv}^{i}$ is the fixed point of (12c). As a consequence, the
probabilities in (15) defined for $\mathcal{R}_{0}^{\eqref{sys:ODE-multi-
species}}>1$ become those in (13). Proof of the existence part is deferred to
Appendix C. ∎

In MBPA, processes either reach zero or approach infinity. The probability of
extinction is interpreted in our model as the probability of a minor epidemic,
while an outbreak is the establishment of the pathogen. The latter is only
temporary since the MBPA does not address the situation far from the DFE.

##  3 Case of one pathogen and two species

To get better insight into the behaviour of the system in a tractable case, we
consider the case with $P=2$ species and $V=1$ pathogen. We first specialise
the model and results of mathematical analysis to this special case (Section
3.1), then consider numerically four specific transmission scenarios, which we
also summarise by indicating the type of propagation taking place between
populations $P_{1}$ and $P_{2}$.

  1. 1.

Pathogen propagation within and between species (Section 3.2,
$P_{1}\leftrightarrow P_{2}$).

  2. 2.

The pathogen is transmitted to both species, but one species can only infect
members of its species (Section 3.3, $P_{1}\to P_{2},P_{2}\not\to P_{1}$).

  3. 3.

Both species can acquire the pathogen, but one of the two species does not
become a transmitter (Section 3.4, $P_{1}\to P_{2},P_{2}\not\to$).

  4. 4.

The pathogen is established at an endemic level in one species and is absent
from the other species (Section 3.5, $P_{1}^{\star}\to P_{2}$).

###  3.1 The models and their basic analysis

Because $V=1$, the second index of $\varepsilon$ and $\gamma$ and the third
index of $\beta$ always equals 1. To simplify notation, in the remainder of
Section 3, we drop this superfluous index and write $\varepsilon_{p}$,
$\gamma_{p}$ and $\beta_{pq}$, for $p,q=1,\ldots,P$.

####  3.1.1 The ODE model when $P=2$ and $V=1$

Setting $P=2$ and $V=1$ in (1) gives

$\displaystyle\dot{S}_{1}=b_{1}-\beta_{11}S_{1}I_{1}-\beta_{12}S_{1}I_{2}-d_{1}S_{1}$
(16a)
$\displaystyle\dot{L}_{1}=\beta_{11}S_{1}I_{1}+\beta_{12}S_{1}I_{2}-(\varepsilon_{1}+d_{1})L_{1}$
(16b) $\displaystyle\dot{I}_{1}=\varepsilon_{1}L_{1}-(\gamma_{1}+d_{1})I_{1}$
(16c) $\displaystyle\dot{R}_{1}=\gamma_{1}I_{1}-d_{1}R_{1}$ (16d)
$\displaystyle\dot{S}_{2}=b_{2}-\beta_{21}S_{2}I_{1}-\beta_{22}S_{2}I_{2}-d_{2}S_{2}$
(16e)
$\displaystyle\dot{L}_{2}=\beta_{21}S_{2}I_{1}+\beta_{22}S_{2}I_{2}-(\varepsilon_{2}+d_{2})L_{2}$
(16f) $\displaystyle\dot{I}_{2}=\varepsilon_{2}L_{2}-(\gamma_{2}+d_{2})I_{2}$
(16g) $\displaystyle\dot{R}_{2}=\gamma_{2}I_{2}-d_{2}R_{2},$ (16h)

where $\Phi_{1}=\beta_{11}I_{1}+\beta_{12}I_{2}$ and
$\Phi_{2}=\beta_{12}I_{1}+\beta_{22}I_{1}$. Recall that for the analysis, the
system is considered with solutions in the order
$S_{1},S_{2},L_{1},L_{2},I_{1},I_{2},R_{1},R_{2}$.

The analysis in Section 2.3 carries through here, taking into account that
since $V=1$, a few adaptations have to be made to the terms and matrices
defined there. For instance, the disease-free equilibrium (DFE) of (16) is

$\bm{E}_{0}^{\eqref{sys:ode-2p1v}}=\left(\left\llbracket
S_{p}^{0}\right\rrbracket,\bm{0}_{\mathbb{R}^{3P}}\right)=\left(\frac{b_{1}}{d_{1}},\frac{b_{2}}{d_{2}},\bm{0}_{\mathbb{R}^{6}}\right),$
(17)

infected variables are $\left\llbracket L_{p}\right\rrbracket,\left\llbracket
I_{p}\right\rrbracket$, $P\times P$ matrix $\mathbf{G}_{12}$ is not a block
matrix but instead has entry $(p,q)$ equal to $S_{p}^{0}\beta_{pq}$ and, since
$P=2$,

$\mathbf{G}_{12}=\begin{pmatrix}\beta_{11}S_{1}^{0}&\beta_{12}S_{1}^{0}\\\
\beta_{21}S_{2}^{0}&\beta_{22}S_{2}^{0}\end{pmatrix}.$

Blocks in the matrix $\mathbf{W}$ take the form

$\mathbf{W}_{11}=\mathsf{diag}(\llbracket\varepsilon_{p}+d_{p}\rrbracket),\mathbf{W}_{21}=\mathsf{diag}(\llbracket\varepsilon_{p}\rrbracket)\text{
and }\mathbf{W}_{22}=\mathsf{diag}(\llbracket\gamma_{p}+d_{p}\rrbracket),$

so that in the case $P=2$ under consideration,

$\mathbf{W}=\begin{pmatrix}\varepsilon_{1}+d_{1}&0&0&0\\\
0&\varepsilon_{2}+d_{2}&0&0\\\ -\varepsilon_{1}&0&\gamma_{1}+d_{1}&0\\\
0&-\varepsilon_{2}&0&\gamma_{2}+d_{2}\end{pmatrix}.$

It follows that the basic reproduction number of (16) is

$\mathcal{R}_{0}^{\eqref{sys:ode-2p1v}}=\frac{\beta_{11}\kappa_{1}S_{1}^{0}+\beta_{22}\kappa_{2}S_{2}^{0}+\sqrt{\left(\beta_{11}\kappa_{1}S_{1}^{0}-\beta_{22}\kappa_{2}S_{2}^{0}\right)^{2}+4\beta_{12}\beta_{21}\kappa_{1}\kappa_{2}S_{1}^{0}S_{2}^{0}}}{2},$
(18)

where, for $p=1,2$,

$\kappa_{p}=\frac{\varepsilon_{p}}{(\varepsilon_{p}+d_{p})(\gamma_{p}+d_{p})}.$

Results of Lemma 1 carry forward to the local asymptotic stability or
instability of (17) based on the value of
$\mathcal{R}_{0}^{\eqref{sys:ode-2p1v}}$ as defined by (18).

####  3.1.2 The CTMC model when $P=2$ and $V=1$

In this case of one virus and two species, the CTMC takes the form

$\mathbf{X}_{t}=(S_{1}(t),S_{2}(t),L_{1}(t),L_{2}(t),I_{1}(t),I_{2}(t),R_{1}(t),R_{2}(t)),\quad
t\in\mathbb{R}_{+}$ (19)

and is characterised by the transition rates in Table 3.

Event ($p=1,2$) | Transition | Transition rates
---|---|---
Birth of $S_{p}$ | $S_{p}\rightarrow S_{p}+1$ | $b_{p}$
Natural death of $S_{p}$ | $S_{p}\rightarrow S_{p}-1$ | $d_{p}S_{p}$
Natural death of $L_{p}$ | $L_{p}\rightarrow L_{p}-1$ | $d_{p}L_{p}$
Natural death of $I_{p}$ | $I_{p}\rightarrow I_{p}-1$ | $d_{p}I_{p}$
Natural death of $R_{p}$ | $R_{p}\rightarrow R_{p}-1$ | $d_{p}R_{p}$
Infection of $S_{p}$ by $I_{q}$ | $S_{p}\rightarrow S_{p}-1,L_{p}\rightarrow L_{p}+1$ | $\beta_{pq}S_{p}I_{q}$
End of incubation of $L_{p}$ | $L_{p}\rightarrow L_{p}-1,I_{p}\rightarrow I_{p}+1$ | $\varepsilon_{p}L_{p}$
Recovery in $I_{p}$ | $L_{p}\rightarrow L_{p}-1,I_{p}\rightarrow I_{p}+1$ | $\gamma_{p}I_{p}$
Table 3: Reaction rates used to determine transition probabilities for
2-species and 1-pathogen CTMC model.

####  3.1.3 Branching process approximation

Theorem 3 is specialised to the $P=2$, $V=1$ case by letting
$Z=(L_{1},L_{2},I_{1},I_{2})$ be a MBPA of the CTMC defined in (19), with
infected types $\ell_{10}$, $\ell_{20}$, $i_{10}$ and $i_{20}$. The p.g.f.
(12) takes here the form, for
$\mathbf{u}=(u_{1}^{\ell},u_{2}^{\ell},u_{1}^{i},u_{2}^{i})$,

$\mathbf{F}(\mathbf{u})=(f_{1}^{\ell}(\mathbf{u}),f_{2}^{\ell}(\mathbf{u}),f_{1}^{i}(\mathbf{u}),f_{2}^{i}(\mathbf{u})),$
(20)

where

$\displaystyle f_{1}^{\ell}(\mathbf{u})$
$\displaystyle=\frac{\varepsilon_{1}u_{1}^{i}+d_{1}}{\varepsilon_{1}+d_{1}}$
(21a) $\displaystyle f_{2}^{\ell}(\mathbf{u})$
$\displaystyle=\frac{\varepsilon_{2}u_{2}^{i}+d_{2}}{\varepsilon_{2}+d_{2}}$
(21b) $\displaystyle f_{1}^{i}(\mathbf{u})$
$\displaystyle=\frac{(\beta_{11}S_{1}^{0}u_{1}^{\ell}+\beta_{12}S_{2}^{0}u_{2}^{\ell})u_{1}^{i}+\gamma_{1}+d_{1}}{\Lambda_{1}}$
(21c) $\displaystyle f_{2}^{l}(\mathbf{u})$
$\displaystyle=\frac{(\beta_{21}S_{1}^{0}u_{1}^{\ell}+\beta_{22}S_{2}^{0}u_{2}^{\ell})u_{2}^{i}+\gamma_{2}+d_{2}}{\Lambda_{2}},$
(21d)

with

$\Lambda_{p}=\sum_{q=1}^{P}\beta_{pq}S_{q}^{0}+\gamma_{p}+d_{p},\quad p=1,2.$

Solving the equation $\mathbf{F}(\bm{z})=\bm{z}$ in the present case involves
finding $(z_{1},z_{2},z_{3},z_{4})$ such that

$\displaystyle
z_{1}=\frac{\varepsilon_{1}z_{3}+d_{1}}{\varepsilon_{1}+d_{1}},\;z_{2}=\frac{\varepsilon_{2}z_{4}+d_{2}}{\varepsilon_{2}+d_{2}}$
(22a)
$\displaystyle(\beta_{11}S_{1}^{0}z_{1}+\beta_{12}S_{2}^{0}z_{2})z_{3}+\gamma_{1}+d_{1}=(\beta_{11}S_{1}^{0}+\beta_{12}S_{2}^{0}+\gamma_{1}+d_{1})z_{3}$
(22b)
$\displaystyle(\beta_{21}S_{1}^{0}z_{1}+\beta_{22}S_{2}^{0}z_{2})z_{4}+\gamma_{2}+d_{2}=(\beta_{21}S_{1}^{0}+\beta_{22}S_{2}^{0}+\gamma_{2}+d_{2})z_{4}.$
(22c)

This is easily done numerically in applications and is known by Theorem 3 to
have a unique solution in $[\bm{0},\mathbf{1})$ when $\mathcal{R}_{0}>1$.

###  3.2 Infectious Hematopoietic Necrosis ($P_{1}\leftrightarrow P_{2}$)

Initially observed at fish hatcheries in Oregon and Washington in the 1950s
[37], Infectious Hematopoietic Necrosis (IHN) is a viral disease that affects
various species of salmonids, including Sockeye Salmon (_Oncorhynchus nerka_)
and Chum Salmon (_Oncorhynchus keta_). The causative agent of IHN is the
Infectious Hematopoietic Necrosis virus (IHNV), which belongs to the
_Rhabdoviridae_ family. This virus primarily targets the hematopoietic
tissues, leading to severe anaemia and necrosis. The incubation period ranges
from 5 to 45 days [40]. Clinical signs of IHN include lethargy, darkening of
skin colour, haemorrhages in various organs and eventual death. Infected fish
may display reduced swimming ability and impaired feeding behaviour due to
anaemia caused by red blood cell destruction [47]. IHN can have a significant
economic impact on fish farms that raise young rainbow trout or salmon, with
mortality rates reaching 90-95% in highly susceptible fish species [14, 40].
The virus can be transmitted horizontally through direct contact or vertically
from infected parents to their offspring. Waterborne transmission is also
possible, particularly in crowded aquaculture settings.

Assume Chum Salmon is species 1 and Sockeye Salmon is species 2. IHN can
spread within and between these two species. From (22), the fixed point is
solution to

$\displaystyle
z_{1}=\frac{\varepsilon_{1}z_{3}+d_{1}}{\varepsilon_{1}+d_{1}},\,z_{2}=\frac{\varepsilon_{2}z_{4}+d_{2}}{\varepsilon_{2}+d_{2}}$
(23a)
$\displaystyle\mathcal{R}_{01}z_{3}^{2}-\left(\mathcal{R}_{01}+\frac{\beta_{11}S_{2}^{0}}{\gamma_{1}+d_{1}}(1-z_{4})+1\right)z_{3}+1=0$
(23b)
$\displaystyle\mathcal{R}_{02}z_{4}^{2}-\left(\mathcal{R}_{02}+\frac{\beta_{22}S_{1}^{0}}{\gamma_{2}+d_{2}}(1-z_{3})+1\right)z_{4}+1=0.$
(23c)

Computing exact expressions of $z_{3}$ and $z_{4}$ is not easy, but from
Theorem 3, this fixed point exists. The probabilities of IHN extinction and
INH outbreak are:

$\displaystyle\mathbb{P}_{\text{ext}}^{\text{IHN}}$
$\displaystyle=\left\\{\begin{aligned}
&\left(\frac{\varepsilon_{1}z_{3}+d_{1}}{\varepsilon_{1}+d_{1}}\right)^{\ell_{10}}\left(\frac{\varepsilon_{2}z_{4}+d_{2}}{\varepsilon_{2}+d_{2}}\right)^{\ell_{20}}z_{3}^{i_{10}}z_{4}^{i_{20}},&\mathcal{R}_{0}^{\eqref{sys:ode-2p1v}}>1\\\
&1,&\mathcal{R}_{0}^{\eqref{sys:ode-2p1v}}<1,\end{aligned}\right.$
$\displaystyle\mathbb{P}_{\text{outbreak}}^{\text{IHN}}$
$\displaystyle=1-\mathbb{P}_{\text{ext}}^{\text{IHN}}$ (24)

We perform a sensitivity analysis of the probability of an outbreak of the INH
virus, assessing the impact of each parameter. Note that here and throughout
the computational work, we assume reproduction numbers larger than 1 unless
otherwise specified.

Chum salmon has a lifespan of 3 to 6 years, during which females lay between
2,000 and 4,000 eggs [31]. Similarly, sockeye salmon has a lifespan of 4 to 5
years, with females laying between 2,000 and 4,500 eggs [32]. Assuming an 80%
hatch rate, the birth rates for Chum (species 1) and Sockeye (species 2) in
(16) are within the range of [355.55, 711.11] and [355.55, 800] per year,
respectively. In a study [15], the incidence of infection after release ranged
from 0% to 20%. These percentages do not indicate transmission rates. However,
when considering the transmission rates of (16), we compute their values using
the corresponding basic reproduction number. Refer to Table 4 for parameter
values ranges used in the sensitivity analysis. Moreover, with an expression
of the basic reproduction number of species 1 given by $\mathcal{R}_{01}$, the
transmission rates are given by

$\displaystyle\beta_{11}$
$\displaystyle=\frac{\mathcal{R}_{01}(\varepsilon_{1}+d_{1})(\gamma_{1}+d_{1})}{\varepsilon_{1}S_{1}^{0}}$
(25a) and $\displaystyle\beta_{22}$
$\displaystyle=\frac{\mathcal{R}_{02}(\varepsilon_{2}+d_{2})(\gamma_{2}+d_{2})}{\varepsilon_{2}S_{2}^{0}}.$
(25b)

We assume that the infection rate of species 2 by species 1 and species 1 by
species 2 is five times higher than the infection rate within species 1 and
species 2, respectively, i.e., $\beta_{21}=5\beta_{11}$ and
$\beta_{12}=5\beta_{22}$.

| Meaning |  Range/day |
---|---|---|---
$b_{1}$ | Birth rate species 1 | $[1,5]$ | [31]
$b_{2}$ | Birth rate species 2 | $[1,5]$ | [32]
$\varepsilon_{1}$ | Incubation rate of 1 | $[0.02,0.2]$ | [40]
$\varepsilon_{2}$ | Incubation rate of 2 | $[0.02,0.2]$ | [40]
$\gamma_{1}$ | Recovery rate of 1 | $[0.1,0.33]$ | [26]
$\gamma_{2}$ | Recovery rate of 2 | $[0.1,0.33]$ | [26]
$d_{1}$ | Mortality rate of 1 | [$\frac{1}{6}$, $\frac{1}{2}$]$\times$ $\frac{1}{365}$] |
$d_{2}$ | Mortality rate of 2 | [$\frac{1}{5}$, $\frac{1}{2}$]$\times$ $\frac{1}{365}$] |
Table 4: Parameter ranges and values for IHN transmission between Chum Salmon
(species 1) and Sockeye Salmon (species 2).

Figure 1 presents a sensitivity analysis of the probability of disease
outbreak, illustrating the significant impact of incubation rates on the
probability of an IHN outbreak. Also important are demographic parameters
$b_{i}$ and $d_{i}$.

Figure 1: PRCC of the probability (24) of IHN outbreak for different initial
conditions $y_{0}=(l_{10},l_{20},i_{10},i_{20})\in(e_{1},e_{2},e_{3},e_{4})$.
The range of parameter values remaining in Table 4.

###  3.3 Transmission from wild to farmed fish ($P_{1}\to P_{2},P_{2}\not\to
P_{1}$)

Interactions between wild and farmed fish populations are complex and crucial
aspects of ecological and economic landscapes. In the scenario where wild fish
could penetrate the space where fish are farmed, the potential transmission of
diseases between these populations becomes a significant concern [4]. Note
that contamination of wild fish by farm fish is also of concern [21, 25].

The situation we consider here has wild fish (species 1) able to introduce
diseases to farmed populations, while farmed fish (species 2) cannot transmit
diseases to wild fish, because they are raised in captivity. This asymmetry in
transmission dynamics implies that (19) is considered with $\beta_{22}\neq 0$,
while $\beta_{12}=0$.

It follows that the fixed point relations (22) take the form

$\displaystyle
z_{1}=\frac{\varepsilon_{1}z_{3}+d_{1}}{\varepsilon_{1}+d_{1}},\quad
z_{2}=\frac{\varepsilon_{2}z_{4}+d_{2}}{\varepsilon_{2}+d_{2}},\quad
z_{3}=\frac{1}{\mathcal{R}_{01}},$ (26a)
$\displaystyle\mathcal{R}_{02}z_{4}^{2}-\left(\mathcal{R}_{02}+c+1\right)z_{4}+1=0.$
(26b)

where we have denoted
$c=\beta_{21}(\mathcal{R}_{01}-1)(\gamma_{1}+d_{1})/(\beta_{11}(\gamma_{2}+d_{2}))$
and

$\mathcal{R}_{02}=\frac{\beta_{22}\varepsilon_{2}S_{2}^{0}}{(\varepsilon_{2}+d_{1})(\gamma_{2}+d_{1})}$
(27)

the basic reproduction for the pathogen in species 2 _in the absence of
contact_ with species 1.

To solve (26b), we first consider the discriminant and using the fact the
$c>0$,

$D=\left(\mathcal{R}_{02}+c+1\right)^{2}-4\mathcal{R}_{02}>\left(\mathcal{R}_{02}+1\right)^{2}-4\mathcal{R}_{02}=(\mathcal{R}_{02}-1)^{2}>0.$

Therefore, (26b) has two real solutions $z_{4}^{-}$ and $z_{4}^{+}$ given by

$z_{4}^{-}=\frac{\left(\mathcal{R}_{02}+c+1\right)-\sqrt{D}}{2\mathcal{R}_{02}}\quad\text{and}\quad
z_{4}^{+}=\frac{\left(\mathcal{R}_{02}+c+1\right)+\sqrt{D}}{2\mathcal{R}_{02}}.$

Given that $\mathcal{R}_{02}>1$, both $z_{4}^{-}$ and $z_{4}^{+}$ are
positive, $z_{4}^{+}>z_{4}^{-}$, and $z_{4}^{+}z_{4}^{-}=1/\mathcal{R}_{02}$.
Moreover, since $D>(\mathcal{R}_{02}-1)^{2}$, we obtain

$z_{4}^{+}>1+\frac{c}{2\mathcal{R}_{02}}\Longrightarrow
z_{4}^{+}z_{4}^{-}>z_{4}^{-}\Longrightarrow
z_{4}^{-}<\frac{1}{\mathcal{R}_{02}}.$

Thus, $z_{4}^{+}>1$, $z_{4}^{-}<1$ and solution to the fixed point problem has

$z_{4}:=z_{4}^{-}=\frac{\left(\mathcal{R}_{02}+c+1\right)-\sqrt{\left(\mathcal{R}_{02}+c+1\right)^{2}-4\mathcal{R}_{02}}}{2}.$
(28)

Then the probability of extinction of the disease in a wild-to-farmed ($W\to
F$) context

$\mathbb{P}_{\text{ext}}^{W\to F}=\left\\{\begin{aligned}
&z_{1}^{\ell_{10}}z_{2}^{\ell_{20}}z_{3}^{i_{10}}z_{4}^{i_{20}}&\mathcal{R}_{01},\mathcal{R}_{02}>1\\\
&1,&\mathcal{R}_{01},\mathcal{R}_{02}\leq 1,\end{aligned}\right.$

where $z_{1},z_{2}$ and $z_{3}$ are given by (26a) and $z_{4}$ by (28);
furthermore,

$\mathbb{P}_{\text{outbreak}}^{W\to F}=1-\mathbb{P}_{\text{ext}}^{W\to F}.$
(29)

For the sensitivity analysis of (29), we utilized the assumed parameter values
provided in Table 5. The parameter ranges for $\beta_{11}$ and $\beta_{21}$
were computed using the relations given in (25a) and (25b).

Parameter | Range | Value
---|---|---
$\mathcal{R}_{0}$ | $[0.1,5]$ | -
$b_{1}$ | $[2,20]$ | 8
$b_{2}$ | $[2,20]$ | 10
$\varepsilon_{1}$ | $[0.01,0.2]$ | 0.05
$\varepsilon_{2}$ | $[0.01,0.2]$ | 0.05
$\gamma_{1}$ | $[0.01,0.1]$ | 0.02
$\gamma_{2}$ | $[0.01,0.1]$ | 0.02
$d_{1}$ | [$\frac{1}{10}$, $\frac{1}{2}$]$\times$ $\frac{1}{365}$ | -
$d_{2}$ | [$\frac{1}{10}$, $\frac{1}{2}$]$\times$ $\frac{1}{365}$ | -
Table 5: Parameter ranges and values for transmission from wild to farmed
fish.

Sensitivity analysis of disease outbreak probability
$\mathbb{P}_{\text{outbreak}}^{W\to F}$ (29) (Figure 2) shows that birth
($b_{1}$), death ($d_{1}$) and recovery ($\gamma_{1}$) rates of species 1 are
key drivers of the probability of an outbreak, regardless of initial
conditions. Interestingly, for the incubation rates ($\varepsilon_{1}$ and
$\varepsilon_{2}$) of species 1 and 2, respectively, there is a notable
difference in impact based on the initial condition. The PRCC is 0.3 and 0.22
for the initial condition starting with one latent individual in Species 1 and
Species 2, respectively, while it is very small when initially starting with
infected individuals.

In contrast, factors such as infection rates $\beta$ of both species and death
rate of species 2 exert minimal influence on the probability of an outbreak.
Moreover, specific cases involving parameters $\gamma_{2}$ and $d_{1}$ exhibit
interesting patterns depending on where the infection starts. For instance,
these parameters show a positive impact when the infection originates within
species 1 but have a negative impact if it starts in species 2.

Figure 2: PRCC of the probability (29) of disease outbreak
$\mathbb{P}_{\text{outbreak}}^{W\to F}$ for four different initial conditions
$y_{0}=(\ell_{10},\ell_{20},i_{10},i_{20})\in(e_{1},e_{2},e_{3},e_{4})$.
Parameter values ranges in Table 5.

Focusing on sensitivity to mortality rates, we observe that effects can be
both positive and negative. To investigate this further, we consider in Figure
3 the probability $\mathbb{P}_{\text{outbreak}}^{W\to F}$ of an outbreak as a
function of the mortality rates $d_{1}$ and $d_{2}$ of species 1 and 2,
respectively. When the infections start with one infectious individual in
species 1, the probability of an outbreak is between 0.55 and 1. When the
disease starts with one infectious individual of species 2, this range narrows
to between 0.9 and 1. For initial infections in species 1, it is primarily the
mortality rate of species 2 that impacts the value of this probability. In
contrast, when the infection is initiated by species 2, both mortality rates
play a role in determining the outcome. Specifically, as the mortality rate
$d_{1}$ of species 1 increases, there is a decrease in the probability of an
outbreak, while an increase in $d_{2}$ leads to a higher probability value.
Thus, initial infections in different species can lead to varying
probabilities of outbreak.

Figure 3: Probability (29) of disease outbreak as a function of mortality
rates of wild and farmed fish. For two initial conditions
$y_{0}=(\ell_{10},\ell_{20},i_{10},i_{20})\in(e_{3},e_{4})$. Parameter values
in Table 4, with $\beta_{11}=\beta_{22}=10^{-5}$, $\beta_{21}=5\beta_{11}$.

###  3.4 Viral Hemorrhagic Septicaemia ($P_{1}\to P_{2},P_{2}\not\to$)

Viral Hemorrhagic Septicaemia (VHS) is a highly contagious and deadly disease
that affects various species of fish in both cultured and wild populations. It
is prevalent in freshwater and marine environments across several regions of
the Northern Hemisphere [16, 29, 43]. VHS is caused by the Viral Hemorrhagic
Septicemia Virus (VHSV) and was first isolated in Alaska from skin lesion
material of two Pacific cod _Gadus macrocephalus_ [28]. VHSV is known for
causing hemorrhagic septicemia, which leads to severe internal bleeding and
organ damage in infected fish. The symptoms vary depending on the species and
the stage of infection. Some common signs include lethargy, loss of appetite,
abnormal swimming behaviour and external haemorrhaging. Over 100 species of
freshwater and marine fish have been reported to be naturally or
experimentally susceptible to VHSV [10]. These species include the Pacific
herring _Clupea pallasii_ [24, 27].

We consider the latter as species 1 in (16). We suppose that the latent stage
is an enzootic stage, the infectious stage is the combination of disease
amplification and outbreak stage and the recovery stage combines recovery and
refractory as suggested by [30]. Now, while VHSV can infect a wide range of
fish, not all infected individuals show signs of clinical disease and not all
are capable of transmitting the virus to others [33]. So we consider as
species 2 the hybrid fry, Steelhead trout (Salmo gairdneri)$\times$Coho Salmon
(Oncorhynchus kisutch), which can get infected, but whose ability to transmit
VHS has not been demonstrated [33], further simplifying the situation by
making the assumption that this hybrid fry cannot transmit VHS. As a
consequence, $\beta_{121}=0$ and $\beta_{221}=0$. Then, the fixed point
relations (22) become

$\displaystyle
z_{1}=\frac{\varepsilon_{1}z_{3}+d_{1}}{\varepsilon_{1}+d_{1}},\;z_{2}=\frac{\varepsilon_{2}z_{4}+d_{2}}{\varepsilon_{2}+d_{2}}$
(30a) $\displaystyle\mathcal{R}_{01}z_{3}^{2}-(\mathcal{R}_{01}+1)z_{3}+1=0$
(30b)
$\displaystyle\beta_{21}S_{1}^{0}z_{1}z_{4}+\gamma_{2}+d_{2}=(\beta_{21}S_{1}^{0}+\gamma_{2}+d_{2})z_{4}$
(30c)

where we have denoted

$\mathcal{R}_{01}=\frac{\beta_{11}\varepsilon_{11}S_{1}^{0}}{(\varepsilon_{11}+d_{1})(\gamma_{11}+d_{1})}$
(31)

the basic reproduction for species 1 and the virus _in the absence of contact_
with species 2.

Solving (30) yields

$(1,1,1,1)\;\text{and}\;\left(z_{1},z_{2},z_{3},z_{4}\right)$

where

$z_{1}=\frac{\varepsilon_{1}z_{3}+d_{1}}{\varepsilon_{1}+d_{1}},\;z_{2}=\frac{\varepsilon_{2}z_{4}+d_{2}}{\varepsilon_{2}+d_{2}},\;z_{3}=\frac{1}{\mathcal{R}_{01}},\;z_{4}=\frac{1}{c+1}$
(32)

with
$c=\beta_{21}(\mathcal{R}_{01}-1)(\gamma_{1}+d_{1})/(\beta_{11}(\gamma_{2}+d_{2}))$.

Therefore, the probability of extinction (no outbreak) is one if
$\mathcal{R}_{01}\leq 1$, but less than one if $\mathcal{R}_{01}>1$. Given the
initial conditions $L_{1}(0)=\ell_{10}$, $I_{1}(0)=i_{10}$,
$L_{2}(0)=\ell_{20}$, and $I_{2}(0)=i_{20}$, it follows from the independent
branching process approximation that the probabilities of the disease
extinction and the disease outbreak are

$\displaystyle\mathbb{P}_{\text{ext}}^{\text{VHS}}$
$\displaystyle=\left\\{\begin{aligned}
&z_{1}^{\ell_{10}}z_{2}^{\ell_{20}}z_{3}^{i_{10}}z_{4}^{i_{20}},&\mathcal{R}_{01}>1\\\
&1,&\mathcal{R}_{01}\leq 1,\end{aligned}\right.$
$\displaystyle\mathbb{P}_{\text{outbreak}}^{\text{VHS}}$
$\displaystyle=1-\mathbb{P}_{\text{ext}}^{\text{VHS}},$ (33)

where the expressions of $z_{1},z_{2},z_{3}$ and $z_{4}$ are given in (32).

On average, a female Pacific herring lays 20,000 eggs each year [1]. If we
suppose that the viable percentage of the population is around 60%, then
$b_{1}=12000/year\approx 33/day$. In laboratory experiments with infected
herring, it was observed that shed VHSV can be identified in water within 4-5
days after exposure (PE), preceding the onset of host mortality due to the
disease. The peak of viral shedding occurs between 6-10 days PE [17]. We then
consider the incubation period of 5-100 days. The duration of recovery varied
depending on the phase of the epizootic. During the acute phase, which
occurred around day 13 post-exposure, virus loads in tissues were
significantly higher compared to the recovery phase, which spanned days 30 to
42 [20]; for sensitivity analysis of this duration is taken between 20 and 100
days. We suppose hybrid fry averaging 12,000 eggs per year with 80% viability,
giving $b_{2}=9600/year=26.3/day$. For the sensitivity analysis of the
probability $\mathbb{P}_{\text{outbreak}}^{\text{VHS}}$ of VHS outbreak, the
range of parameter values are given in Table 6, (25a) and (25b) is used to
compute infection rates.

| Meaning | Range | Values | Ref
---|---|---|---|---
$\mathcal{R}_{01}$ | Basic reproduction number | $[0.1,5]$ |  |
$b_{1}$ | Birth rate Pacific herring | $[10,50]$ | 33 | [1]
$b_{2}$ | Birth rate hybrid fry | $[10,30]$ | 26.3 |
$\varepsilon_{1}$ | Incubation rate of Pacific herring | $[0.01,0.2]$ | 0.02 | [17]
$\varepsilon_{2}$ | Incubation rate of hybrid fry | $[0.01,0.2]$ | 0.02 | [17]
$\gamma_{1}$ | Recovery rate of Pacific herring | $[0.01,0.05]$ | 0.03 | [20]
$\gamma_{2}$ | Recovery rate of hybrid fry | $[0.01,0.05]$ | 0.03 | [20]
$d_{1}$ | Mortality rate of species 1 | [$\frac{1}{15}$, $\frac{1}{2}$]$\times$ $\frac{1}{365}$ | 3.$10^{-4}$ | [1, 41]
$d_{2}$ | Mortality rate of species 2 | [$\frac{1}{15}$, $\frac{1}{2}$]$\times$ $\frac{1}{365}$ | 3.$10^{-4}$ |
Table 6: Parameter range and values for (16) with two species, Pacific herring
(species 1) and hybrid fry (species 2) and one pathogen, VHS.

In Figure 4, a sensitivity analysis of the probability of disease outbreak
$\mathbb{P}_{\text{outbreak}}^{\text{VHS}}$ is presented for four different
initial conditions. The birth rate $b_{1}$ of species 1 stands out as having a
significant impact on the probability of an outbreak regardless of the initial
condition. When the infection begins with one latent and one infectious
individual in species 2, parameters $b_{1}$, $\gamma_{1}$, $\gamma_{2}$, and
$d_{1}$ all exhibit a similar level of impact on the probability of an
outbreak. The mortality rate $d_{1}$ of species 1 shows varying effects on the
probability of an outbreak depending on where the disease originates. When the
infection starts in species 2, $d_{1}$ has a significantly positive impact on
the likelihood of an outbreak. However, this impact becomes negative when the
disease originates in species 1.

Amongst all parameters analysed, certain factors show no discernible impact on
the probability of an outbreak. This includes parameter $b_{2}$, as well as
transmission rates $\beta_{11}$ and $\beta_{21}$.

These findings highlight which variables have little to no effect on
influencing whether a disease outbreak will occur. Moreover, the influence of
different initial conditions on outcomes is a crucial aspect to consider.

Figure 4: PRCC of the probability of disease outbreak (33) for different
initial conditions
$y_{0}=(\ell_{10},\ell_{20},i_{10},i_{20})\in(e_{1},e_{2},e_{3},e_{4})$.
Parameter values ranges in Table 6.

###  3.5 The pathogen is endemic in one species ($P_{1}^{\star}\to P_{2}$)

The problem motivating this study concerns vagrant species coming in contact
with resident species while bearing a pathogen the resident species has not
been exposed to yet.

To model this situation, we artificially impose that species 1 be at an
endemic equilibrium while species 2 be at a disease-free equilibrium. A
situation with mixed equilibria of this type is not possible in (16), so we
“cheat”: we start by assuming that both species are not in contact and that
species 1 is at the endemic equilibrium
$E_{1}^{\star}:=(S_{1}^{\star},L_{1}^{\star},I_{1}^{\star},R_{1}^{\star})$, in
which

$E_{1}^{\star}=\left(\frac{S_{1}^{0}}{\mathcal{R}_{01}},\;\frac{(\mathcal{R}_{01}-1)(\gamma_{1}+d_{1})d_{1}}{\beta_{11}\varepsilon_{1}},\;\frac{(\mathcal{R}_{01}-1)d_{1}}{\beta_{11}},\;\frac{(\mathcal{R}_{01}-1)\gamma_{1}}{\beta_{11}}\right),$
(34)

where $\mathcal{R}_{01}$ given by (31) is the basic reproduction for the
pathogen in species 1 in the absence of contact with species 2. We suppose
that this determines the dynamics of the pathogen in species 1, which we now
ignore except insofar as the infecting potential of the $I_{1}^{\star}$
individuals from species 1 infectious with the pathogen potentially coming
into contact with susceptible individuals from species 2. We then consider the
second population, focusing on conditions leading to the disease becoming
established there.

####  3.5.1 The ODE introduction model

Given a prevalence of infection $I_{1}^{\star}$ in species 1, the infection
dynamics in species 2 is governed by

$\displaystyle\dot{S}_{2}$
$\displaystyle=b_{2}-\beta_{21}I_{1}^{\star}S_{2}-\beta_{22}S_{2}I_{2}-d_{2}S_{2}$
(35a) $\displaystyle\dot{L}_{2}$
$\displaystyle=\beta_{21}I_{1}^{\star}S_{2}+\beta_{22}S_{2}I_{2}-(\varepsilon_{2}+d_{2})L_{2}$
(35b) $\displaystyle\dot{I}_{2}$
$\displaystyle=\varepsilon_{2}L_{2}-(\gamma_{2}+d_{2})I_{2}$ (35c)
$\displaystyle\dot{R}_{2}$ $\displaystyle=\gamma_{2}I_{2}-d_{2}R_{2},$ (35d)

considered with nonnegative initial conditions.

Model (35) is not a classic immigration model, since the term
$\beta_{21}I_{1}^{\star}$ is factor of $S_{2}$, not a constant. However,
similarly to immigration models, the term $\beta_{21}I_{1}^{\star}S_{2}$ in
(35b) precludes the existence of a disease-free equilibrium for (35). As a
consequence, no basic reproduction number can be computed for (35). The method
of [3] is not applicable here since immigration occurs at the _per capita_
rate $\beta_{21}I_{1}^{\star}$, not at a constant rate. The term
$\beta_{21}I_{1}^{\star}$ cannot be considered as encoding horizontal
transmission either, at least not in the usual sense, since inflow into (35b)
is function of $S_{2}$, not $I_{2}$.

As a consequence, a thorough analysis of properties of (35) is difficult and
beyond the scope of this work. Instead, we focus on simple properties as well
as computational considerations.

System (35) admits a unique equilibrium, an endemic equilibrium taking the
form
$E_{2}^{\star}:=(S_{2}^{\star},L_{2}^{\star},I_{2}^{\star},R_{2}^{\star})$,
where

$S_{2}^{\star}=\frac{\varepsilon_{21}+d_{2}}{\beta_{21}I_{1}^{\star}+\beta_{22}I_{2}^{\star}}L_{2}^{\star},\;L_{2}^{\star}=\frac{\gamma_{21}+d_{2}}{\varepsilon_{21}}I_{2}^{\star},\;R_{2}^{\star}=\frac{\gamma_{2}I_{2}^{\star}}{d_{2}}$

and $I_{2}^{\star}$ is a root of the second order polynomial

$\beta_{22}(I_{2}^{\star})^{2}-((\mathcal{R}_{02}-1)d_{2}-\beta_{21}I_{1}^{\star})I_{2}^{\star}-\frac{\beta_{21}d_{2}I_{1}^{\star}\mathcal{R}_{02}}{\beta_{22}}=0,$
(36)

where $\mathcal{R}_{02}$ given in (27) is the basic reproduction for species 2
and the virus _in the absence of contact_ with species 1.

Since the coefficients $\beta_{22}>0$ and
$\beta_{21}\mathcal{R}_{02}I_{1}^{\star}/\beta_{22}>0$, Descartes rule of
signs implies that the polynomial (36) has a unique positive solution. The
expression of this solution is given by

$I_{2}^{\star}=\frac{(\mathcal{R}_{02}-1)d_{2}-\beta_{21}I_{1}^{\star}+\sqrt{((\mathcal{R}_{02}-1)d_{2}-\beta_{21}I_{1}^{\star})^{2}+4\mathcal{R}_{02}d_{2}\beta_{21}I_{1}^{\star}}}{2\beta_{22}}.$
(37)

Figure 5: Equilibrium prevalence of infection $I_{2}^{\star}$ in population 2
as a function of the reproduction number $\mathcal{R}_{02}$ for the pathogen
in species 2. The different curves correspond to different values of
prevalence $I_{1}^{\star}$ in population 1.

In Figure 5, we show how (equilibrium) prevalence of infection in species 2
depends on prevalence of infection in species 1. For low values of the
reproduction number in species 2, the situation is quite dependent on the
prevalence of infection in the “introducing species”, but as the reproduction
number increases, this dependence diminishes to the point of the situation
becoming indistinguishable for large values of $\mathcal{R}_{02}$.

####  3.5.2 The CTMC introduction model

Note that it is not possible here to use a multitype branching process
approximation, for roughly the same reasons that a basic reproduction number
does not exist for the deterministic model (35). Indeed, while branching
processes incorporating _immigration_ exist [19, 34], they assume that
immigration is a process that is independent from the population immigrants
are joining. In our model, this is not true, since “immigration” depends on
the susceptible population $S_{2}$. So we focus here on the computational
analysis of the CTMC. In all results presented here, 10,000 simulations
(realisations) of the CTMC associated to (35) are used for each data point.

Figure 6: Percentage of 10K realisations in which spread by species 1
(introductions) and species 2 occurs. Also shown is the percentage of
realisations where introductions by species 1 are followed by extinctions.

Figure 6 shows the percentage of 10K simulations in which at least one
transmission occurs from species 1 to species 2 (red curve) and from species 2
to species 2 (blue curve). Simulations assume that $\mathcal{R}_{02}=1.5$ and
are run for 90 days. Note that the latter type of transmission always requires
that introduction by species 1 has taken place. Additionally, we show the
percentage of realisations in which spread by species 1 occurs followed by
extinctions, where we characterise an extinction of the infection in species 2
as a moment when the total number of infected in species 2,
$L_{2}(t)+I_{2}(t)$, is zero after having been positive.

This illustrates a very important part of the introduction process. In keeping
with the terminology in [6, 5], let us assess success of an introduction from
the perspective of the pathogen. In view of Figure 6, what drives successful
introductions is the size of the introduction, i.e., the so-called _inoculum
size_. However, before it becomes established in species 2 (the consequence of
which is shown in a deterministic context by Figure 5), the infection in
species 2 must “survive” the _stochastic phase_ of the epidemic; see, e.g.,
the one location case in [8]. This illustrates that spillover events are often
unsuccessful, as observed in the zoonotic case with bats [38].

Figure 7: Distribution of the times at which the first infection event occurs
in species 2 following contact with an infected individual from species 1
(red) and species 2 (blue), for different values of $I_{1}^{\star}$. Boxplots
within the violins show the median, interquartile range and whiskers extents.

To better understand this issue, consider Figure 7, where we show violin plots
of the distribution of times at which the first infection in species 2 arises
stemming from contact with, respectively, species 1 and species 2. These
values are from the same simulations as used in Figure 6 and thus represent
the percentages shown there of 10K simulations. First, consider infections
with source species 1, i.e., introductions into species 2. We observe that as
the prevalence $I_{1}^{\star}$ in species 1 increases, the time to first
introduction progressively diminishes, with the interquartile range covering
smaller and smaller values. Now consider the timing of infections originating
from species 2, i.e., taking place after the infection has become somewhat
established in species 2. There, we observe that as the prevalence
$I_{1}^{\star}$ increases, times to the first transmission decrease, but not
as acutely as those to the first introduction. This confirms our earlier
observation that introduction and establishment, although evidently
correlated, are not as tightly tied as can be expected.

##  4 Case of two pathogens and two species

Here the dynamics with $P=2$ species and $V=2$ pathogens are considered. While
the cases in Section 3 are more tractable mathematically, the situation here
is more realistic. In practice, the collaboration motivating this work is
interested in over a dozen pathogens potentially infecting four fish species.
We do not consider such a general situation here, but illustrate the
computational complexities that arise even when $P=V=2$. One particularly
interesting feature is the existence of mixed equilibria, i.e., equilibria in
which one of the pathogens is present and the other absent.

###  4.1 Deterministic model and basic analysis

$\displaystyle\dot{S}_{1}=b_{1}-\left(\sum_{q=1}^{2}\sum_{v=1}^{2}\beta_{1qv}I_{qv}+d_{1}\right)S_{1},$
$\displaystyle\dot{S}_{2}$
$\displaystyle=b_{2}-\left(\sum_{q=1}^{2}\sum_{v=1}^{2}\beta_{2qv}I_{qv}+d_{2}\right)S_{2},$
(38a)
$\displaystyle\dot{L}_{11}=\sum_{q=1}^{2}\beta_{1q1}I_{q1}S_{1}-(\varepsilon_{11}+d_{1})L_{11},$
$\displaystyle\dot{L}_{21}$
$\displaystyle=\sum_{q=1}^{2}\beta_{2q1}I_{q1}S_{2}-(\varepsilon_{21}+d_{2})L_{21},$
(38b)
$\displaystyle\dot{L}_{12}=\sum_{q=1}^{2}\beta_{1q2}I_{q2}S_{1}-(\varepsilon_{12}+d_{1})L_{12},$
$\displaystyle\dot{L}_{22}$
$\displaystyle=\sum_{q=1}^{2}\beta_{2q2}I_{q2}S_{2}-(\varepsilon_{22}+d_{2})L_{22},$
(38c)
$\displaystyle\dot{I}_{11}=\varepsilon_{11}L_{11}-(\gamma_{11}+d_{1})I_{11},$
$\displaystyle\dot{I}_{21}$
$\displaystyle=\varepsilon_{21}L_{21}-(\gamma_{21}+d_{2})I_{21},$ (38d)
$\displaystyle\dot{I}_{12}=\varepsilon_{12}L_{12}-(\gamma_{12}+d_{1})I_{12},$
$\displaystyle\dot{I}_{22}$
$\displaystyle=\varepsilon_{22}L_{22}-(\gamma_{22}+d_{2})I_{22},$ (38e)
$\displaystyle\dot{R}_{1}=\sum_{v=1}^{2}\gamma_{1v}I_{1v}-d_{1}R_{1},$
$\displaystyle\dot{R}_{2}$
$\displaystyle=\sum_{v=1}^{2}\gamma_{2v}I_{2v}-d_{2}R_{2}.$ (38f)

The disease-free equilibrium of (38) is

$\bm{E}_{0}^{\eqref{sys:ode-2p2v}}=(S_{1}^{0},S_{2}^{0},0_{\mathbb{R}^{10}}),\;\text{
with }S_{1}^{0}=\frac{b_{1}}{d_{1}}\text{ and }S_{2}^{0}=\frac{b_{2}}{d_{2}}.$
(39)

To compute the basic reproduction number $\mathcal{R}_{0}$, observe that the
matrices $\mathbf{G}_{12}$ and $\mathbf{W}$ derived in Section 2.3 take here
the form

$\mathbf{G}_{12}=\begin{bmatrix}\mathbf{G}_{12}^{11}&\mathbf{G}_{12}^{12}\\\
\mathbf{G}_{12}^{21}&\mathbf{G}_{12}^{22}\\\ \end{bmatrix},$

where $\mathbf{G}_{12}^{11}=S_{1}^{0}\mathsf{diag}(\beta_{111},\beta_{112})$,
$\mathbf{G}_{12}^{12}=S_{1}^{0}\mathsf{diag}(\beta_{121},\beta_{122})$,
$\mathbf{G}_{12}^{21}=S_{2}^{0}\mathsf{diag}(\beta_{211},\beta_{212})$ and
$\mathbf{G}_{12}^{22}=S_{2}^{0}\mathsf{diag}(\beta_{221},\beta_{222})$ are
diagonal matrices, and $S_{1}^{0}$ and $S_{2}^{0}$ are scalars. The matrix
$\mathbf{W}$ is block lower triangular, with blocks

$\mathbf{W_{11}}=\mathsf{diag}(d_{1}+\varepsilon_{11},d_{1}+\varepsilon_{12},d_{2}+\varepsilon_{21},d_{2}+\varepsilon_{22}),$
$\mathbf{W_{21}}=\mathsf{diag}(\varepsilon_{11},\varepsilon_{12},\varepsilon_{21},\varepsilon_{22})$

and

$\mathbf{W_{22}}=\mathsf{diag}(d_{1}+\gamma_{11},d_{1}+\gamma_{12},d_{2}+\gamma_{21},d_{2}+\gamma_{22}).$

The basic reproduction number of system (38) following the formula in equation
(8) is given by:

$\mathcal{R}_{0}^{\eqref{sys:ode-2p2v}}=\max(\mathcal{R}_{0}^{1},\mathcal{R}_{0}^{2})$

where $\mathcal{R}_{0}^{1}=\rho(\mathcal{B}_{1})$ and
$\mathcal{R}_{0}^{2}=\rho(\mathcal{B}_{2})$, with the matrices
$\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ obtained by using a form of a matrix
in equation (9), and then

$\mathcal{B}_{1}=\displaystyle\begin{pmatrix}\frac{\beta_{111}\varepsilon_{11}S_{1}^{0}}{(\varepsilon_{11}+d_{1})(\gamma_{11}+d_{1})}&\frac{\beta_{121}\varepsilon_{21}S_{1}^{0}}{(\varepsilon_{21}+d_{2})(\gamma_{21}+d_{2})}\\\
\frac{\beta_{211}\varepsilon_{11}S_{2}^{0}}{(\varepsilon_{11}+d_{1})(\gamma_{11}+d_{1})}&\frac{\beta_{221}\varepsilon_{21}S_{2}^{0}}{(\varepsilon_{21}+d_{2})(\gamma_{21}+d_{2})}\end{pmatrix},\mathcal{B}_{2}=\displaystyle\begin{pmatrix}\frac{\beta_{112}\varepsilon_{12}S_{1}^{0}}{(\varepsilon_{12}+d_{1})(\gamma_{12}+d_{1})}&\frac{\beta_{122}\varepsilon_{22}S_{1}^{0}}{(\varepsilon_{22}+d_{2})(\gamma_{22}+d_{2})}\\\
\frac{\beta_{212}\varepsilon_{12}S_{2}^{0}}{(\varepsilon_{12}+d_{1})(\gamma_{12}+d_{1})}&\frac{\beta_{222}\varepsilon_{22}S_{2}^{0}}{(\varepsilon_{22}+d_{2})(\gamma_{22}+d_{2})}\end{pmatrix}$
(40)

Note that the result provided by using
$\mathcal{R}_{0}^{\eqref{sys:ode-2p2v}}$ does not show the whole picture.
Indeed, one interesting characteristic of (1) is that the viruses function in
a disconnected way. This can be inferred from the reducibility of the system
discussed in Appendix A.

######  Theorem 4.

Consider (38) with $\mathcal{R}_{01}>1$ and $\mathcal{R}_{02}\leq 1$. Then the
DFE of (38) is unstable and consists of a mixed equilibrium wherein pathogen 1
is present at an endemic level and pathogen 2 is absent. Stability of the
pathogen-2-free equilibrium is global and asymptotic with respect to the
pathogen-2 subsystem.

The existence part of the proof of this result is shown in Appendix D. Global
asymptotic stability of the pathogen-2-free equilibrium follows directly from
Theorem 2 and reducibility of the system.

###  4.2 Branching process approximation

Let $Z=(L_{11},L_{12},L_{21},L_{22},I_{11},I_{12},I_{21},I_{22})$ be the
multitype branching process approximation of CTMC $X^{2,2}(t)$ with infected
types $\ell_{11}$, $\ell_{12}$, $\ell_{21}$, $\ell_{22}$, $i_{11}$, $i_{12}$,
$i_{21}$ and $i_{22}$. The p.g.f. for
$\mathbf{u}=(u_{11}^{\ell},u_{12}^{\ell},u_{21}^{\ell},u_{22}^{\ell},u_{11}^{i},u_{21}^{i},u_{21}^{i},u_{22}^{i})$
is defined as

$\mathbf{F}(\mathbf{u})=(f_{11}^{\ell}(\mathbf{u}),f_{12}^{\ell}(\mathbf{u}),f_{21}^{\ell}(\mathbf{u}),f_{22}^{\ell}(\mathbf{u}),f_{11}^{i}(\mathbf{u}),f_{21}^{i}(\mathbf{u}),f_{21}^{i}(\mathbf{u}),f_{22}^{i}(\mathbf{u})),$
(41)

where, for $p,q=1,2$,

$\displaystyle f_{pq}^{\ell}(\mathbf{u})$
$\displaystyle=\frac{\varepsilon_{pq}u_{pq}^{i}+d_{p}}{\varepsilon_{pq}+d_{p}},$
(42a) $\displaystyle f_{pq}^{i}(\mathbf{u})$
$\displaystyle=\frac{\left(\sum_{k=1}^{2}\beta_{pkq}S_{k}^{0}u_{kq}^{\ell}\right)u_{pq}^{i}+\gamma_{pq}+d_{p}}{\Lambda_{pq}}.$
(42b)

The Jacobian matrix then takes the form

$DF(u)=\begin{bmatrix}\begin{array}[]{c:c}0&\mathbb{M}_{12}\\\
\hdashline\mathbb{M}_{21}&\mathbb{M}_{22}\\\ \end{array}\end{bmatrix},$

where

$\mathbb{M}_{12}=\mathsf{diag}\left(\frac{\varepsilon_{11}}{\varepsilon_{11}+d_{1}},\frac{\varepsilon_{12}}{\varepsilon_{12}+d_{1}},\frac{\varepsilon_{21}}{\varepsilon_{21}+d_{2}},\frac{\varepsilon_{22}}{\varepsilon_{22}+d_{2}}\right),$
$\mathbb{M}_{21}=\begin{pmatrix}\frac{\beta_{111}S_{1}^{0}u_{11}^{i}}{\Lambda_{11}}&0&\frac{\beta_{121}S_{2}^{0}u_{11}^{i}}{\Lambda_{11}}&0\\\
0&\frac{\beta_{112}S_{1}^{0}u_{12}^{i}}{\Lambda_{12}}&0&\frac{\beta_{122}S_{2}^{0}u_{12}^{i}}{\Lambda_{12}}\\\
\frac{\beta_{211}S_{1}^{0}u_{21}^{i}}{\Lambda_{21}}&0&\frac{S_{2}^{0}\beta_{221}u_{21}^{i}}{\Lambda_{21}}&0\\\
0&\frac{\beta_{212}S_{1}^{0}u_{22}^{i}}{\Lambda_{22}}&0&\frac{\beta_{222}S_{2}^{0}u_{22}^{i}}{\Lambda_{22}}\end{pmatrix}$

and

$\mathbb{M}_{22}=\mathsf{diag}\Biggl{(}\frac{\beta_{111}S_{1}^{0}u_{11}^{\ell}+\beta_{121}S_{2}^{0}u_{21}^{\ell}}{\Lambda_{11}},\frac{\beta_{112}S_{1}^{0}u_{12}^{\ell}+\beta_{122}S_{2}^{0}u_{22}^{\ell}}{\Lambda_{12}},\\\
\frac{\beta_{211}S_{1}^{0}u_{11}^{\ell}+\beta_{221}S_{2}^{0}u_{21}^{\ell}}{\Lambda_{21}},\frac{\beta_{212}S_{1}^{0}u_{12}^{\ell}+\beta_{222}S_{2}^{0}u_{22}^{\ell}}{\Lambda_{22}}\Biggr{)}.$

Theorem 3 applies here. Given $L_{11}(0)=\ell_{110}$, $L_{12}(0)=\ell_{120}$,
$L_{21}(0)=\ell_{210}$, $L_{22}(0)=\ell_{220}$, $I_{11}(0)=i_{110}$,
$I_{12}(0)=i_{120}$, $I_{21}(0)=i_{210}$ and $I_{22}(0)=i_{220}$, it follows
from the independent branching process approximation that the probabilities of
extinction and disease outbreak are:

$\mathbb{P}_{\text{ext}}^{\eqref{eq:pdf-2p2v}}=\left\\{\begin{aligned}
&(z_{11}^{\ell})^{\ell_{110}}(z_{12}^{\ell})^{\ell_{120}}(z_{21}^{\ell})^{\ell_{210}}(z_{22}^{\ell})^{\ell_{220}}(z_{11}^{i})^{i_{110}}(z_{12}^{i})^{i_{120}}(z_{21}^{i})^{i_{210}}(z_{22}^{i})^{i_{220}},&\mathcal{R}_{0}^{\eqref{sys:ode-2p2v}}>1\\\
&1,&\mathcal{R}_{0}^{\eqref{sys:ode-2p2v}}<1,\end{aligned}\right.$
$\mathbb{P}_{\text{outbreak}}^{\eqref{eq:pdf-2p2v}}=1-\mathbb{P}_{\text{ext}}^{\eqref{eq:pdf-2p2v}}$
(43)

##  5 Discussion

We formulated an SLIR model for the spread of $V$ pathogens between and within
$P$ species. The model was first formulated as a set of $2P(1+V)$ ordinary
differential equations, of which we considered elementary properties: basic
reproduction number $\mathcal{R}_{0}$ and local asymptotic stability as a
function of the value of the reproduction number. Global asymptotic stability
when $\mathcal{R}_{0}<1$ was also established. We then considered the
corresponding continuous-time Markov chain (CTMC) model, which provides a
better tool to study the behaviour of the system close to the disease-free
equilibrium (DFE), which is our main interest here. To do so, we employed a
multitype branching process approximation of the CTMC near the DFE, obtaining
an expression for the probability of an outbreak when $\mathcal{R}_{0}>1$.
This probability was interpreted, in the context of our model, as the
probability that the pathogen becomes established in the population, at least
temporarily. (The result is local and does not address the proper
establishment at an endemic level.)

The case of a single pathogen spreading between two species was then
investigated computationally. A metapopulation (spatial) version of this
situation with a slightly different underlying SLIR model was investigated
both mathematically and computationally in [7], with, using the notation here,
$P$ species present. However, focusing on just two species as we did here
allows to get a better understanding of the processes. To this effect, in
particular, we investigated the sensitivity of the system to its parameters in
the case of three viruses affecting fish leading to very different
transmission scenarios. These highlighted in particular the important role
played by demographic parameters. A fourth, more abstract case concerned
introduction of a pathogen in a population by another in which the pathogen is
present at an endemic level.

One interesting feature of the model is that despite its complication,
pathogens function more or less independently from one another. This was shown
in the case $P=V=2$, which we considered next. We conjecture that the
following natural extension of Theorem 4 holds.

######  Conjecture 5.

Consider (1) in which pathogens $a=1,\ldots,k$ have $\mathcal{R}_{0a}<1$ and
pathogens $e=k+1,\ldots,V$ have $\mathcal{R}_{0e}>1$, for some
$k\in\\{1,\ldots,n\\}$. Then (1) has an unstable mixed equilibrium in which
pathogens $1,\ldots,k$ are absent and pathogens $k+1,\ldots,V$ are present at
an endemic level.

By Theorem 2, those pathogens that are absent would naturally be globally
asymptotically stably so. This highlights a limitation of the model: because
we assume that species dynamics is independent of the pathogens and that
coinfections cannot occur, a situation as described by Theorem 4 or Conjecture
5 is possible. While not necessarily unrealistic, taking into account more
advanced interactions between pathogens, or effects of pathogens onto their
host species, could be interesting and lead to wholly different results.
Competition effects between species could also be incorporated and would also
likely lead to different results.

Another interesting variation could involve considering an _epidemic_ model.
In its current form, the model is an _endemic_ model, with the basic
reproduction numbers distinguishing between a situation where the disease is
absent and one where the disease is present at an endemic level. In this case,
of course, the pathogen always becomes extinct; it is not clear, though, if a
situation similar to that in Theorem 4 and Conjecture 5, in this case with
some pathogens undergoing an epidemic and others not, would hold.

### Acknowledgements

The authors thank Sharon Clouthier, Karen Dunmall and Collin Ghallagher, from
the Freshwater Institute (Department of Fisheries and Oceans Canada, Winnipeg)
for initiating the project and helpful discussions. JA acknowledges support
from NSERC through the Discovery Grants program.

##  Appendix A Normal forms of matrices

When considering the global asymptotic stability of the DFE in Appendix B or
existence of threshold behaviour in Appendix C, we observe that the matrices
involved are reducible. While the presentation used in the body of the paper
is the most natural, the proofs in these appendices require to use the normal
form of the matrices involved. Further understanding of the reproduction
number is also gained by using the normal form. Let us illustrate this using
matrix $A_{22}$ in Appendix B. The permutation matrix obtained in the process
is the same for all four matrices we consider here.

It is clear that matrix $A_{22}$ as given by (47) is reducible. Indeed,
consider the weighted loop-directed graph having $A_{22}$ as its adjacency
matrix. Diagonal entries of the two main diagonal blocks correspond to loops.
The blocks $\bm{\Psi}(\mathbf{x}_{S},\mathbf{x}_{I})$ and
$\mathsf{diag}(\llbracket\varepsilon_{pv}\rrbracket)$, on the other hand, show
that vertices form $PV$ strongly connected components, with each component
comprising $2$ vertices. Consider for instance the $(1,1)$ entry in
$\bm{\Psi}$ and the $(1,1)$ entry in
$\mathsf{diag}(\llbracket\varepsilon_{pv}\rrbracket)$. They establish that
vertex 1 is connected to vertex $PV+1$ and that vertex $PV+1$ is connected to
vertex 1. All entries in these matrices define similar pairs of vertices,
giving the $PV$ strong components. Ordering vertices (and corresponding matrix
entries) so that vertices in a strong connected component are listed
consecutively, it is then easy to find the permutation matrix $\bm{\Pi}$ such
that

$\bm{\Pi}\
A_{22}(\mathbf{x}_{S},\mathbf{x}_{I})\bm{\Pi}^{T}=\bigoplus_{\llbracket
pv\rrbracket}\begin{pmatrix}-(\varepsilon_{pv}+d_{p})&\sum_{q=1}^{P}\beta_{qpv}S_{q}\\\
\varepsilon_{pv}&-(\gamma_{pv}+d_{p})\end{pmatrix}.$ (44)

Applying $\bm{\Pi}$ to $D\mathbf{F}$ given by (49) in Appendix C, we find

$\bm{\Pi}\ D\mathbf{F}(\mathbf{u}^{\ell},\mathbf{u}^{i})\
\bm{\Pi}^{T}=\bigoplus_{\llbracket
pv\rrbracket}\begin{pmatrix}0&\dfrac{\varepsilon_{pv}}{\varepsilon_{pv}+d_{p}}\\\
S_{p}^{0}\dfrac{\beta_{pqv}u_{pv}^{i}}{\Lambda_{pv}}&\dfrac{\sum_{q=1}^{P}S_{q}^{0}u_{qv}^{\ell}}{\Lambda_{pv}}\end{pmatrix}.$
(45)

##  Appendix B Proof of Theorem 2

To prove global asymptotic stability of the disease-free equilibrium of (1)
when $\mathcal{R}_{0}<1$, we use a result of Kamgang and Sallet [22, Theorems
4.3 and 4.5]. For the convenience of the reader, we recall this result but
with notation adapted to the problem under consideration here.

Let $\mathbf{x}_{S}=(\llbracket S_{p}\rrbracket,\llbracket
R_{p}\rrbracket)^{T}\in\mathbb{R}^{2P}$ and $\mathbf{x}_{I}=(\llbracket
L_{pv}\rrbracket,\llbracket I_{pv}\rrbracket)^{T}\in\mathbb{R}^{2PV}$ be,
respectively, the vectors of non-infected compartments and infected
compartments. Denote $\mathbf{x}_{S}^{0}=(\llbracket
S_{p}^{0}\rrbracket,\llbracket 0\rrbracket)^{T}\in\mathbb{R}^{2P}$ the part of
the disease-free equilibrium of (1) corresponding to $\mathbf{x}_{S}$, i.e.,
the DFE is $(\mathbf{x}_{S}^{0},\bm{0}_{2PV})$. Rewrite (1) in the following
compact form,

$\displaystyle\dot{\mathbf{x}}_{S}$
$\displaystyle=A_{11}(\mathbf{x}_{S},\mathbf{x}_{I})(\mathbf{x}_{S}-\mathbf{x}_{S}^{0})+A_{12}(\mathbf{x}_{S},\mathbf{x}_{I})\mathbf{x}_{I},$
(46) $\displaystyle\dot{\mathbf{x}}_{I}$
$\displaystyle=A_{22}(\mathbf{x}_{S},\mathbf{x}_{I})\mathbf{x}_{I},$

with

$\displaystyle A_{11}(\mathbf{x}_{S},\mathbf{x}_{I})$
$\displaystyle=\begin{pmatrix}-\mathsf{diag}\left(\llbracket
d_{p}\rrbracket\right)&\bm{0}\\\ \bm{0}&-\mathsf{diag}\left(\llbracket
d_{p}\rrbracket\right)\\\ \end{pmatrix}\in\mathbb{R}^{2P\times 2P},$
$\displaystyle A_{12}(\mathbf{x}_{S},\mathbf{x}_{I})$
$\displaystyle=\begin{pmatrix}\bm{0}&-\bm{\Phi}(\mathbf{x}_{S},\mathbf{x}_{I})\\\
\bm{0}&\bm{\Gamma}(\mathbf{x}_{S},\mathbf{x}_{I})\\\
\end{pmatrix}\in\mathbb{R}^{2P\times 2PV},$ $\displaystyle
A_{22}(\mathbf{x}_{S},\mathbf{x}_{I})$
$\displaystyle=\begin{pmatrix}-\mathsf{diag}\left(\llbracket\varepsilon_{pv}+d_{p}\rrbracket\right)&\bm{\Psi}(\mathbf{x}_{S},\mathbf{x}_{I})\\\
\mathsf{diag}\left(\llbracket\varepsilon_{pv}\rrbracket\right)&-\mathsf{diag}\left(\llbracket\gamma_{pv}+d_{p}\rrbracket\right)\end{pmatrix}\in\mathbb{R}^{2PV\times
2PV},$ (47)

where $\Phi(\mathbf{x}_{S},\mathbf{x}_{I})$ is a $P\times PV$-matrix whose
$j$th row is a $PV$-vector with entries $\llbracket
S_{j}\beta_{jpv}\rrbracket$, $\bm{\Gamma}(\mathbf{x}_{S},\mathbf{x}_{I})$ is a
$P\times PV$-matrix with $j$th row having $V$ nonzero entries
$\gamma_{j1},\ldots,\gamma_{jV}$ in columns $(j-1)V+1$ to $jV$ and
$\bm{\Psi}(\mathbf{x}_{S},\mathbf{x}_{I})$ is a diagonal $PV\times PV$-matrix,

$\bm{\Psi}(\mathbf{x}_{S},\mathbf{x}_{I})=\mathsf{diag}\left(\left\llbracket\sum_{q=1}^{P}\beta_{qpv}S_{q}\right\rrbracket\right),$

with the enumerator running over indices $p=1,\ldots,P$ and $v=1,\ldots,V$.

As discussed in Appendix A, it is clear that matrix $A_{22}$ as given by (47)
is reducible. So we apply the method described in that Appendix and instead
work with the similar normal form matrix (44). We can apply [22, Theorem 4.3]
to each of the $PV$ blocks in (44) or apply [22, Theorem 4.5], which considers
the reducible case. We use a combination of the two results.

Let $\Omega\subset\mathbb{R}^{2P}_{+}\times\mathbb{R}^{2PV}_{+}$ be the set
defined in the statement of Theorem 2. If the following five conditions hold
true, then [22, Theorem 4.5] establishes that the disease-free equilibrium
$(\mathbf{x}_{S}^{0},\bm{0}_{2PV})$ is globally asymptotically stable when
$\mathcal{R}_{0}<1$.

C1

System (46) is defined on a positively invariant set $\Omega$ of the
nonnegative orthant and dissipative on $\Omega$.

C2

Subsystem
$\dot{\mathbf{x}}_{S}=A_{11}(\mathbf{x}_{S},\bm{0})(\mathbf{x}_{S}-\mathbf{x}_{S}^{0})$
is globally asymptotically stable at the disease-free equilibrium
$\mathbf{x}_{S}^{0}$ on the canonical projection of $\Omega$ on
$\mathbb{R}^{2P}_{+}$.

C3

Matrix $\tilde{A}_{22}(\mathbf{x}_{S},\mathbf{x}_{I})$ given by (44) is block
upper triangular, with each diagonal block $\hat{A}_{22}^{pv}$ being Metzler
and irreducible for any $\mathbf{x}=(\mathbf{x}_{S},\mathbf{x}_{I})\in\Omega$.

C4

For each $p,v$, there exists an upper-bound matrix $\hat{A}_{22}^{pv}$ for
$M=\\{A_{22}^{pv}(\mathbf{x})\in\mathbb{R}^{2\times 2};\mathbf{x}\in\Omega\\}$
with the property that either $\hat{A}_{22}^{pv}\not\in M$ or if
$\hat{A}_{22}^{pv}\in M$, (i.e., $\hat{A}_{22}=\max_{\Omega}M$), then for any
$\hat{\mathbf{x}}\in\Omega$ such that $\hat{A}_{22}=A_{22}(\hat{\mathbf{x}})$,
$\hat{\mathbf{x}}\in\mathbb{R}^{2P}_{+}\times\\{0\\}^{2PV}$.

C5

The spectral abscissa of the matrix $\hat{A}_{22}^{pv}$ verifies
$\sigma(\hat{A}_{22}^{pv})\leq 0$ when $\mathcal{R}_{0}^{\eqref{sys:ODE-multi-
species}}\leq 1$.

The right-hand side of (46) is of class $C^{1}$ on the open set
$\mathbb{R}^{2P}_{+}\times\mathbb{R}^{2PV}_{+}$, so solutions are defined.
Solutions in $\Omega$ remain in $\Omega$. Furthermore, extending $\Omega$ as

$\Omega_{\delta}=\Biggl{\\{}(\llbracket S_{p}\rrbracket,\llbracket
L_{pv}\rrbracket,\llbracket I_{pv}\rrbracket,\llbracket
R_{p}\rrbracket)\in\mathbb{R}^{2P(V+1)}:\\\
N_{p}=S_{P}+\sum_{v=1}^{V}(L_{pv}+I_{pv})+R_{p}\leq\frac{b_{p}}{d_{p}}+\delta;\quad
p=1,\ldots,P\Biggr{\\}},$

we have that solutions with initial conditions in $\mathbb{R}^{2P(V+1)}$
eventually enter $\Omega_{\delta}$, for any $\delta>0$. As a consequence, (46)
is dissipative and C1 is satisfied. Then note that in the original form, the
model without disease is, for $p=1,\ldots,P$,

$\displaystyle\dot{S}_{p}$ $\displaystyle=b_{p}-d_{p}S_{p},$ (48a)
$\displaystyle\dot{R}_{p}$ $\displaystyle=-d_{p}R_{pv}.$ (48b)

Thus it is clear that for each $p=1,\ldots,P$, $S_{p}(t)\to b_{d}/d_{p}$ and
$R_{p}(t)\to 0$ regardless of initial conditions, meaning that condition C2
holds. Matrix (44) is in block-diagonal form and as a consequence, is block
upper triangular. Recall that a matrix is Metzler if its offdiagonal entries
are nonnegative, so C3 holds.

To verify that C4 holds, remark that for solutions in $\Omega$, the maximal
value for a given $S_{p}$ is attained when $S_{p}=N_{p}$, i.e., at the
disease-free equilibrium $\mathbf{x}_{S}^{0}$. In other words, for
$\tilde{\mathbf{x}}\in\mathbb{R}^{2P}_{+}\times\\{0\\}^{2PV}$,
$\tilde{\mathbf{x}}_{S}\leq\mathbf{x}_{S}^{0}$. Thus, the upper-bound matrix
is
$\hat{A}_{22}=\hat{A}_{22}(\mathbf{x}_{S}^{0},\bm{0})=\mathbf{G}-\mathbf{W}$.
Note that the result can also be formulated using blocks $\hat{A}_{22}^{pv}$
as in (44) and selecting the relevant components in
$\bm{\Pi}(\mathbf{G}-\mathbf{W})\bm{\Pi}^{T}$, but we do need the form using
the unreduced matrix $\hat{A}_{22}$ to show that condition C5 holds.

Indeed, to show C5, return to the unreduced form (47) and note that we have
$\sigma(\hat{A}_{22})=\sigma(\mathbf{G}-\mathbf{W})$, with the matrices as
defined in Section 2.3 and satisfying the conditions of [44, Theorem 2]. The
proof of that theorem establishes that $\mathcal{R}_{0}<1$, i.e., the spectral
radius $\rho(\mathbf{G}\mathbf{W}^{-1})<1$, is equivalent to the spectral
abscissa $\sigma(\mathbf{G}-\mathbf{W})<0$. As a consequence, when
$\mathcal{R}_{0}<1$, then $\sigma(\hat{A}_{22})<0$ and the same is true for
each diagonal block $\hat{A}_{22}^{pv}$ in the matrix in normal form (44),
since it is similar to (47).

Since conditions C1–C5 are satisfied, the proof is done.

##  Appendix C Proof of threshold behaviour in Theorem 3

To continue the proof of Theorem 3, we now need to consider the existence of
fixed points of the p.g.f. $\mathbf{F}$. The result will follow from
application of the Threshold Theorem [2] together with [18, Theorem 7.1],
provided we can show that the branching process is positive and regular.

The Jacobian of (12) is the $2PV\times 2PV$ block matrix

$D\mathbf{F}\left(\mathbf{u}^{\ell},\mathbf{u}^{i}\right)=\begin{pmatrix}\begin{array}[]{cc}\bm{0}&\mathbb{M}_{12}\\\
\mathbb{M}_{21}&\mathbb{M}_{22}\\\ \end{array}\end{pmatrix},$ (49)

where each block has size $PV\times PV$. First,

$\mathbb{M}_{12}=\mathsf{diag}\left(\left\llbracket\frac{\varepsilon_{pv}}{\varepsilon_{pv}+d_{p}}\right\rrbracket\right),\;\mathbb{M}_{22}=\mathsf{diag}\left(\left\llbracket\frac{\sum_{q=1}^{P}S_{q}^{0}u_{qv}^{\ell}}{\Lambda_{pv}}\right\rrbracket\right).$

Then, the $PV\times PV$-matrix $\mathbb{M}_{21}$ is itself a block matrix,
with each $V\times V$ sized block taking the form, for
$p,q\in\\{1,\ldots,P\\}$,

$\mathcal{K}_{pq}=\displaystyle
S_{q}^{0}\mathsf{diag}\left(\left\llbracket\frac{\beta_{pqv}u_{pv}^{i}}{\Lambda_{pv}}\right\rrbracket\right).$

As established in Appendix A, (49) is reducible. Using the similarity
transformation in Appendix A, we derive the normal form of $D\mathbf{F}$,
(45). As in the proof in Appendix B, it is useful to use both the original
matrix (49) and its normal form (45).

  1. (i)

It is clear that (49) is such that if $\mathbf{x},\mathbf{y}\in[0,1)^{2PV}$
are such that $\mathbf{x}\leq\mathbf{y}$, one has $D\mathbf{F}(\mathbf{x})\leq
D\mathbf{F}(\mathbf{y})$. Indeed, terms $\mathbf{u}^{\ell}$ and
$\mathbf{u}^{i}$ appear as sums in the numerators of the expressions involving
them. Furthermore, $\mathbf{F}(\bm{0})>\bm{0}$, i.e., it is a nonnegative
matrix with some positive entries. This implies that the multitype branching
processes are not singular [12, Theorem 2.3].

  2. (ii)

The matrix of first moments is $\mathbb{M}=D\mathbf{F}(\mathds{1}_{2PV})$,
where $\mathds{1}_{2PV}$ is the unit column vector of size $2PV$. In the
transformed matrix (45), diagonal blocks of $\mathbb{M}$ take the form

$\begin{pmatrix}0&\dfrac{\varepsilon_{pv}}{\varepsilon_{pv}+d_{p}}\\\
S_{p}^{0}\dfrac{\beta_{pqv}}{\Lambda_{pv}}&\dfrac{\sum_{q=1}^{P}S_{q}^{0}}{\Lambda_{pv}}\end{pmatrix}$

and are therefore irreducible (and even primitive). Consequently, the matrix
of first moments $\mathbb{M}$ is block-primitive.

From $(i)$ and $(ii)$, we conclude that the branching process is positive and
regular. As a consequence, applying the Threshold Theorem [2] together with
[18, Theorem 7.1] to each of the diagonal blocks in the matrix in normal form
(45), gives the threshold behaviour, with existence of a fixed point
$(0,0)<(z_{pv}^{\ell},z_{pv}^{i})<(1,1)$ additionally to
$(z_{pv}^{\ell},z_{pv}^{i})=(1,1)$ when the process is supercritical. Putting
things together, under the conditions of Theorem 3, there exists an additional
fixed point $\bm{0}<\bm{z}<\bm{1}$ when $\mathcal{R}_{0}>1$.

##  Appendix D Existence of a mixed equilibrium for (38)

Suppose that the conditions of Theorem 4 are satisfied: $P=V=2$, virus 2 is at
the disease-free equilibrium (DFE), i.e., $L_{12}=L_{22}=I_{12}=I_{22}=0$ and
$\mathcal{R}_{02}<1$. We seek equilibria of (38) with positive values for
$L_{11}^{\star}$, $I_{21}^{\star}$, $I_{11}^{\star}$ and $I_{21}^{\star}$,
under the assumption that $\mathcal{R}_{01}>1$.

Substituting the DFE of species 2 into (38), we obtain from (38c) that

$S_{1}^{\star}=\frac{\varepsilon_{11}+d_{1}}{\beta_{111}I_{11}^{\star}+\beta_{121}I_{21}^{\star}}L_{11}^{\star},\quad
S_{2}^{\star}=\frac{\varepsilon_{21}+d_{2}}{\beta_{211}I_{11}^{\star}+\beta_{221}I_{21}^{\star}}L_{21}^{\star},$

while (38e) gives

$L_{11}^{\star}=\frac{\gamma_{11}+d_{1}}{\varepsilon_{11}}I_{11}^{\star},\quad
L_{21}^{\star}=\frac{\gamma_{21}+d_{2}}{\varepsilon_{21}}I_{21}^{\star}$

and, finally, from (38f),

$R_{1}^{\star}=\frac{\gamma_{11}}{d_{1}}I_{11}^{\star},\quad
R_{2}^{\star}=\frac{\gamma_{21}}{d_{2}}I_{21}^{\star}.$

Since the total population of each species is governed, for $i=1,2$, by
$\dot{N}_{i}=b_{i}-d_{i}N_{i}$, at an equilibrium,
$b_{i}-d_{i}N_{1}^{\star}=0$ and thus, when (38) is at an equilibrium with
virus 2 at the DFE, one has

$b_{1}-d_{1}\left(S_{1}^{\star}+L_{11}^{\star}+I_{11}^{\star}+R_{1}^{\star}\right)=0\quad\text{and}\quad
b_{2}-d_{2}\left(S_{2}^{\star}+L_{21}^{\star}+I_{21}^{\star}+R_{2}^{\star}\right)=0.$

Expressing all terms as functions of $I_{i1}$, $i=1,2$, and using the
expressions of $S_{1}$, $L_{11}$, $R_{1}$, $S_{2}$, $L_{21}$ and $R_{2}$ gives

$\displaystyle\beta_{111}(I_{11}^{\star})^{2}-\left[d_{1}\left(\mathcal{R}_{01}-1\right)-\beta_{121}I_{21}^{\star}\right]I_{11}^{\star}-\frac{\beta_{121}}{\beta_{111}}d_{1}\mathcal{R}_{01}I_{21}^{\star}=0,$
(50a)
$\displaystyle\beta_{221}(I_{21}^{\star})^{2}-\left[d_{2}\left(\mathcal{R}_{02}-1\right)-\beta_{211}I_{11}^{\star}\right]I_{21}^{\star}-\frac{\beta_{211}}{\beta_{221}}d_{2}\mathcal{R}_{02}I_{11}^{\star}=0.$
(50b)

To simplify computations, let us denote $x=I_{11}^{\star}$ and
$y=I_{21}^{\star}$. Then (50) can be written as

$\displaystyle\Gamma_{1}(x,y)$ $\displaystyle:=\beta_{111}x^{2}+\beta_{121}xy-
d_{1}\left(\mathcal{R}_{01}-1\right)x-\frac{\beta_{121}}{\beta_{111}}d_{1}\mathcal{R}_{01}y=0,$
(51a) $\displaystyle\Gamma_{2}(x,y)$
$\displaystyle:=\beta_{211}xy+\beta_{221}y^{2}-\frac{\beta_{211}}{\beta_{221}}d_{2}\mathcal{R}_{02}x-d_{2}\left(\mathcal{R}_{02}-1\right)y=0.$
(51b)

Both $\Gamma_{1}$ and $\Gamma_{2}$ are conic sections. Their discriminants
$-\beta_{121}^{2}$ and $-\beta_{211}^{2}$ are both negative, hence they are
both hyperbolas. As (51a) has no second degree $y$ monomial, one of its
asymptotes is vertical. Likewise, since (51b) has no second degree $x$
monomial, one of its asymptotes is horizontal. If both of these asymptotes
intersect the positive (or first) quadrant
$Q_{1}=\mathbb{R}_{+}\times\mathbb{R}_{+}$, then $\Gamma_{1}$ and $\Gamma_{2}$
intersect a single time there. The situation is shown in Figures 8 and 9.

Figure 8: Situation leading to the existence of a mixed equilibrium. Red
curve: $\Gamma_{1}$. Blue curve: $\Gamma_{2}$. The centres of the hyperbolas
are shown as triangles of corresponding colours, as are the vertical and
horizontal asymptotes relevant to the problem. The shaded box shows the
possible range of values of the endemic component of the mixed equilibrium.

Let us show that this is indeed the case. There are two ingredients:

  1. 1.

The curves $\Gamma_{1}$ and $\Gamma_{2}$ intersect $Q_{1}$.

  2. 2.

The vertical asymptote of $\Gamma_{1}$ and the horizontal asymptote of
$\Gamma_{2}$ intersect $Q_{1}$.

First, note that the point of intersection $(x,y)=(0,0)$ is obvious since
neither (51a) nor (51b) have terms of degree 0. Now consider the $x$\- and
$y$-intercepts of $\Gamma_{1}$ and $\Gamma_{2}$. For $\Gamma_{1}$,
$x$-intercepts satisfy

$\Gamma_{1}(x,0)=\left(\beta_{111}x-d_{1}(\mathcal{R}_{01}-1)\right)x,$

i.e., $x=0$ and $x=d_{1}(\mathcal{R}_{01}-1)/\beta_{111}>0$ by the assumption
$\mathcal{R}_{01}>1$. For $y$-intercepts,

$\Gamma_{1}(0,y)=-\frac{\beta_{121}}{\beta_{111}}d_{1}\mathcal{R}_{01}y,$

i.e., $\Gamma_{1}$ intercepts the $y$-axis only at the origin. For
$\Gamma_{2}$, $x$-intercepts are given by

$\Gamma_{2}(x,0)=-\frac{\beta_{211}}{\beta_{221}}d_{2}\mathcal{R}_{02}x,$

giving only the origin as $x$-intercept, while $y$-intercepts are given by

$\Gamma_{2}(0,y)=\left(\beta_{221}y-d_{2}(\mathcal{R}_{02}-1)\right)y,$

i.e., $y$-intercepts are $y=0$ and $y=d_{2}(\mathcal{R}_{02}-1)/\beta_{221}<0$
by assumption.

We can therefore establish that $\Gamma_{1}$ and $\Gamma_{2}$ intersect
$Q_{1}$. Indeed, from the gradients $\nabla\Gamma_{1}(x,y)$ and
$\nabla\Gamma_{2}(x,y)$, we deduce that vectors tangent to $\Gamma_{1}$ and
$\Gamma_{2}$ are

$\displaystyle
T_{1}(x,y)=\left(\beta_{121}x-\frac{\beta_{121}}{\beta_{111}}d_{1}\mathcal{R}_{01},d_{1}(\mathcal{R}_{01}-1)-2\beta_{111}x-\beta_{121}y\right),$
$\displaystyle
T_{2}(x,y)\left(2\beta_{221}y+\beta_{211}x-d_{2}(\mathcal{R}_{02}-1),\frac{\beta_{211}}{\beta_{221}}d_{2}\mathcal{R}_{02}-\beta_{211}y\right).$

At the origin,
$T_{1}(0,0)=(-\beta_{121}d_{1}\mathcal{R}_{01}/\beta_{111},d_{1}(\mathcal{R}_{01}-1))$
has signs $(-,+)$. This means that $\Gamma_{1}$ “moves” through origin from
the second quadrant $Q_{2}=\mathbb{R}_{-}\times\mathbb{R}_{+}$ when $x<0$ to
the fourth quadrant $Q_{4}=\mathbb{R}_{+}\times\mathbb{R}_{-}$ when $x>0$. On
the other hand,
$T_{2}(0,0)=(-d_{2}(\mathcal{R}_{02}-1),\beta_{211}d_{2}\mathcal{R}_{02}/\beta_{221})$
has signs $(+,+)$, implying that left of the $y$-axis, $\Gamma_{2}$ is in the
third quadrant $Q_{3}=\mathbb{R}_{-}\times\mathbb{R}_{-}$, while it is in
$Q_{1}$ when $x>0$.

It remains to show that the vertical and horizontal asymptotes of $\Gamma_{1}$
and $\Gamma_{2}$, respectively, intersect the positive quadrant $Q_{1}$.

The centre of $\Gamma_{1}$ is
$(d_{1}\mathcal{R}_{01}/\beta_{111},-d_{1}(\mathcal{R}_{01}+1)/\beta_{121})\in
Q_{4}$ (red triangle in Figure 8). Since $\Gamma_{1}$ intersects $Q_{2}$, the
asymptote to $\Gamma_{1}\cap Q_{2}$ has negative slope. It follows that the
vertical asymptote to $\Gamma_{1}$ is the one to $\Gamma_{1}\cap Q_{1}$ and
therefore, intersects $Q_{1}$. Reasoning similarly, observe that since the
centre
$(-d_{2}(\mathcal{R}_{02}+1)/\beta_{211},d_{2}\mathcal{R}_{02}/\beta_{221})$
of $\Gamma_{2}$ lies in $Q_{2}$ (blue triangle in Figure 8) and $\Gamma_{2}$
intersects $Q_{4}$, the asymptote to $\Gamma_{2}\cap Q_{4}$ has negative slope
and that to the part of $\Gamma_{2}\cap Q_{1}$ is horizontal.

As a consequence, there is a point in the interior of the positive quadrant
$Q_{1}$ where $\Gamma_{1}$ intersects $\Gamma_{2}$. More precisely, remark
that since the centres of the hyperbola lie on the vertical and horizontal
asymptotes, respectively, the $x$-coordinate of the point of intersection
cannot be larger that the $x$-component of the centre of $\Gamma_{1}$ and its
$y$-coordinate cannot exceed the $y$-component of the centre of $\Gamma_{2}$.
This means that the endemic equilibrium belongs to the box

$\left(0,\frac{d_{1}\mathcal{R}_{01}}{\beta_{111}}\right]\times\left(0,\frac{d_{2}\mathcal{R}_{02}}{\beta_{221}}\right]$
(52)

shown shaded in Figures 8 and 9.

Figure 9: Situation leading to the existence of a mixed equilibrium. Zoom on
Figure 8 focusing on the positive quadrant. $\Gamma_{1}$ and its related
features is shown in red, while $\Gamma_{2}$ is in blue.

## References

  * [1] Alaska Department of Fish and Game.  Pacific herring (Clupea pallasii), n.d.  Accessed on February 2, 2024.
  * [2] L.J.S. Allen and P. van den Driessche.  Relations between deterministic and stochastic thresholds for disease extinction in continuous- and discrete-time infectious disease models.  Mathematical Biosciences, 243(1):99–108, May 2013.
  * [3] R.M. Almarashi and C.C. McCluskey.  The effect of immigration of infectives on disease-free equilibria.  Journal of Mathematical Biology, 79:1015–1028, 2019.
  * [4] P. Arechavala-Lopez, P. Sanchez-Jerez, J.T. Bayle-Sempere, I. Uglem, and I. Mladineo.  Reared fish, farmed escapees and wild fish stocks—a triangle of pathogen transmission of concern to mediterranean aquaculture management.  Aquaculture Environment Interactions, 3(2):153–161, 2013.
  * [5] J. Arino, N. Bajeux, S. Portet, and J. Watmough.  Quarantine and the risk of COVID-19 importation.  Epidemiology & Infection, 148:e298, 2020.
  * [6] J. Arino, P.-Y. Boëlle, E.M. Milliken, and S. Portet.  Risk of covid-19 variant importation - how useful are travel control measures?  Infectious Disease Modelling, 6:875–897, 2021.
  * [7] J. Arino, J.R. Davis, D. Hartley, R. Jordan, J. Miller, and P. van den Driessche.  A multi-species epidemic model with spatial migration.  Mathematical Medicine and Biology, 22(2):129–142, 2005.
  * [8] J. Arino and E.M. Milliken.  Effect of mobility on the stochastic phase of an epidemic.  Bulleting of Mathematical Biology, 84, 2022.
  * [9] K.E. Atkins and J.M.J. Travis.  Local adaptation and the evolution of species’ ranges under climate change.  Journal of Theoretical Biology, 266(3):449–457, 2010.
  * [10] W.N. Batts, J. Lovy, R.G. Getchell, M. Faisal, I. Standish, J.V. Warg, N.B.D. Phelps, G. Glenney, and J.R. Winton.  Viral Hemorrhagic Septicemia.  USGS Western Fisheries Research Center, 2020.
  * [11] Ø. Bergh, K.Y. Børsheim, G. Bratbak, and M. Heldal.  High abundance of viruses found in aquatic environments.  Nature, 340(6233):467–468, 1989.
  * [12] A. Berman and R.J. Plemmons.  Nonnegative matrices in the mathematical sciences, January 1979.
  * [13] C.J. Carlson, G.F. Albery, C. Merow, C.H. Trisos, C.M. Zipfel, E.A. Eskew, K.J. Olival, N. Ross, and S. Bansal.  Climate change increases cross-species viral transmission risk.  Nature, 607(7919):555–562, 2022.
  * [14] P. Dixon, R. Paley, R. Alegria-Moran, and B. Oidtmann.  Epidemiological characteristics of infectious hematopoietic necrosis virus (IHNV): a review.  Veterinary Research, 47(1):1–26, 2016.
  * [15] J.S. Foott, D. Free, T. McDowell, K.D. Arkush, and R.P. Hedrick.  Infectious Hematopoietic Necrosis Virus transmission and disease among juvenile Chinook salmon exposed in culture compared to environmentally relevant conditions.  San Francisco Estuary and Watershed Science, 4(1):2, 2006.
  * [16] N. Gagné, A.-M. MacKinnon, L. Boston, B. Souter, M. Cook-Versloot, S. Griffiths, and G. Olivier.  Isolation of Viral Haemorrhagic Septicaemia Virus from Mummichog, Stickleback, Striped bass and Brown trout in eastern Canada.  Journal of Fish Diseases, 30:213–223, 2007.
  * [17] K.A. Garver and L.M. Hawley.  Characterization of viral haemorrhagic septicaemia virus (VHSV) to inform pathogen transfer risk assessments in British Columbia.  Technical Report DFO Can. Sci. Advis. Sec. Res. Doc., Fisheries and Oceans Canada, Canadian Science Advisory Secretariat, 2021.
  * [18] T.E. Harris.  The Theory of Branching Processes, volume 119 of Grundlehren der mathematischen Wissenschaften.  Springer Berlin, Heidelberg, 1963.
  * [19] C.R. Heathcote.  A branching process allowing immigration.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 27(1):138–143, 1965.
  * [20] P.K. Hershberger, J.L. Gregg, C.A. Grady, R.M. Collins, and J.R. Winton.  Susceptibility of three stocks of Pacific herring to Viral Hemorrhagic Septicemia.  Journal of Aquatic Animal Health, 22(1):1–7, 2010.
  * [21] L.-H. Johansen, I. Jensen, H. Mikkelsen, P.-A. Bjørn, P.A. Jansen, and Ø. Bergh.  Disease interaction and pathogens exchange between wild and farmed fish populations with special reference to norway.  Aquaculture, 315(3-4):167–186, 2011.
  * [22] J.-C. Kamgang and G. Sallet.  Computation of threshold conditions for epidemiological models and global stability of the disease-free equilibrium (DFE).  Mathematical Biosciences, 213(1):1–12, 2008.
  * [23] M. Kirkpatrick and N.H. Barton.  Evolution of a species’ range.  The American Naturalist, 150(1):1–23, 1997.
  * [24] R.M. Kocan, P.K. Hershberger, N.E. Elder, and J.R. Winton.  Epidemiology of viral hemorrhagic septicemia among juvenile Pacific herring and Pacific sand lances in Puget Sound, Washington.  Journal of Aquatic Animal Health, 13:77–85, 2001.
  * [25] M. Krkosek, J.S. Ford, A. Morton, S. Lele, R.A. Myers, and M.A. Lewis.  Declining wild salmon populations in relation to parasites from farm salmon.  Science, 318(5857):1772–1775, 2007.
  * [26] S.E. LaPatra, L. Barone, G.R. Jones, and L.I. Zon.  Effects of infectious Hematopoietic Necrosis Virus and infectious Pancreatic Necrosis Virus infection on hematopoietic precursors of the zebrafish.  Blood Cells, Molecules, and Diseases, 26(5):445–452, 2000.
  * [27] T.R. Meyers, S. Short, K. Lipson, W.N. Batts, J.R. Winton, J. Wilcock, and E. Brown.  Association of viral hemorrhagic septicemia virus with epizootic hemorrhages of the skin in Pacific herring Clupea harengus pallasi from Prince William Sound and Kodiak Island, Alaska, USA.  Diseases of Aquatic Organisms, 19:27–37, 1994.
  * [28] T.R. Meyers, J. Sullivan, E. Emmenegger, J. Follett, S. Short, W.N. Batts, and J.R. Winton.  Identification of viral hemorrhagic septicemia virus isolated from Pacific cod Gadus macrocephalus in Prince William Sound, Alaska, USA.  Diseases of Aquatic Organisms, 12:167–175, 1992.
  * [29] T.R. Meyers and J.R. Winton.  Viral Hemorrhagic Septicemia Virus in North America.  Annual Review of Fish Diseases, 5:3–24, 1995.
  * [30] M.G. Mitro and A.L. White.  Viral Hemorrhagic Septicemia and Freshwater Fisheries: The State of the Science.  Technical Bulletin No. 196 Department of Natural Resources, 2008.
  * [31] National Oceanic and Atmospheric Administration (NOAA).  Chum salmon.  NOAA Fisheries website, 2023.  Last updated on 12/22/2023.
  * [32] National Oceanic and Atmospheric Administration (NOAA).  Sockeye salmon.  NOAA Fisheries website, 2023.  Last updated on 12/22/2023.
  * [33] W.M. Ord, M. Le Berre, and P. de Kinkelin.  Viral Hemorrhagic Septicemia: Comparative susceptibility of Rainbow trout (salmo gairdneri) and Hybrids (s. gairdneri $\times$ oncorhynchus kisutch) to experimental infection.  Journal of the Fisheries Research Board of Canada, 33(5):1205–1208, 1976.
  * [34] A.G. Pakes.  Branching processes with immigration.  Journal of Applied Probability, 8(1):32–42, 1971.
  * [35] C. Parmesan.  Ecological and evolutionary responses to recent climate change.  Annu. Rev. Ecol. Evol. Syst., 37:637–669, 2006.
  * [36] B.L. Phillips, G.P. Brown, and R. Shine.  Life-history evolution in range-shifting populations.  Ecology, 91(6):1617–1627, 2010.
  * [37] R.R. Rucker, W.J. Whipple, J.R. Parvin, and C.A. Evans.  A contagious disease of salmon possibly of virus origin.  1953\.
  * [38] C.A. Sánchez, H. Li, K.L. Phelps, C. Zambrana-Torrelio, L.-F. Wang, P. Zhou, Z.-L. Shi, K.J. Olival, and P. Daszak.  A strategy to assess spillover risk of bat SARS-related coronaviruses in Southeast Asia.  Nature Communications, 13(1):4380, 2022.
  * [39] J.P. Sexton, P.J. McIntyre, A.L. Angert, and K.J. Rice.  Evolution and ecology of species range limits.  Annu. Rev. Ecol. Evol. Syst., 40:415–436, 2009.
  * [40] A.R. Spickler.  Infectious hematopoietic necrosis, 2007.  Retrieved on January 19, 2024.
  * [41] K.D.E. Stokesbury, J. Kirsch, E.V. Patrick, and B.L/ Norcross.  Natural mortality estimates of juvenile Pacific herring (Clupea pallasi) in Prince William Sound, Alaska.  Canadian Journal of Fisheries and Aquatic Sciences, 59(3):424–433, 2002.
  * [42] C.A. Suttle.  Marine viruses—major players in the global ecosystem.  Nature Reviews Microbiology, 5(10):801–812, 2007.
  * [43] R. Takano, K. Mori, T. Nishizawa, M. Arimoto, and K. Muroga.  Isolation of viruses from wild Japanese flounder, Paralichthys olivaceus.  Fish Pathology, 36:153–160, 2001.
  * [44] P. van den Driessche and J. Watmough.  Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission.  Mathematical Biosciences, 180(1):29–48, 2002.
  * [45] K.E. Wommack and R.R. Colwell.  Virioplankton: viruses in aquatic ecosystems.  Microbiology and Molecular Biology Reviews, 64(1):69–114, 2000\.
  * [46] A.P. Wyn-Jones and J. Sellwood.  Enteric viruses in the aquatic environment.  Journal of Applied Microbiology, 91(6):945–962, 2001.
  * [47] C.Y. Yong, H.K. Ong, H.C. Tang, S.K. Yeap, A.R. Omar, K.L. Ho, and W.S. Tan.  Infectious Hematopoietic Necrosis Virus: Advances in diagnosis and vaccine development.  PeerJ, 7:e7151, 2019.