A stochastic multi-host model for West Nile virus transmission

Abstract

When initially introduced into a susceptible population, a disease may die out or result in a major outbreak. We present a Continuous-Time Markov Chain model for enzootic WNV transmission between two avian host species and a single vector, and use multitype branching process theory to determine the probability of disease extinction based upon the type of infected individual initially introducing the disease into the population – an exposed vector, infectious vector, or infectious host of either species. We explore how the likelihood of disease extinction depends on the ability of each host species to transmit WNV, vector biting rates on host species, and the relative abundance of host species, as well as vector abundance. Theoretical predictions are compared to the outcome of stochastic simulations. We find the community composition of hosts and vectors, as well as the means of disease introduction, can greatly affect the probability of disease extinction.

KEYWORDS:

• West Nile virus
• continuous-time Markov chain
• multitype branching process
• probability of disease extinction
• host heterogeneity

AMS SUBJECT CLASSIFICATIONS:

• 60J80
• 34A34
• 92D30

1. Introduction

West Nile virus (WNV) is an arbovirus primarily spread between birds and mosquitoes (enzootic transmission). Humans are unable to spread WNV, but remain at risk of infection and complications from the disease. Most infected individuals (around 80%) do not develop any symptoms, but about 1 in 150 infected individuals will develop West Nile encephalitis, a serious illness that causes inflammation of the brain and spinal cord and surrounding tissues [1].

While the large degree of temporal and spatial heterogeneity observed in WNV transmission makes it difficult to predict when and where outbreaks will occur, many factors are known to affect WNV risk [2]. The presence of a WNV-competent mosquito species that feeds on avian hosts is necessary for sustained transmission. The community composition of hosts (the species present and their abundance) is also extremely influential for disease transmission [3], as avian species differ greatly in their ability to contract and spread the virus [4, 5]. Species range from being non-competent or dead-end hosts, like humans, to being highly competent, and increasing host diversity can result in either diluted or amplified transmission, depending on the competence and abundance of the species present [6]. Vector feeding preferences also play an important role in determining the WNV risk of an area [7, 8], as some avian species are fed upon by mosquitoes more or less frequently than would be expected based on their abundance [9].

Many deterministic ordinary differential equation (ODE) models for WNV transmission have been developed and analyzed [8, 10–14]. Mathematical modeling has shown that the risk of an outbreak generally increases with the ratio of vectors to hosts in the population [8, 10, 14, 15]; transmission intensifies when the number of hosts is reduced in the population [14], or with an increase in either the vector carrying capacity or the vector biting rate [13]. Multi-host models have shown increasing host diversity can result in either diluted or amplified transmission, depending on the competence and abundance of the species present along with vector feeding preferences [6, 11, 16].

WNV is likely introduced to new areas by small numbers of infectious hosts or vectors. When disease prevalence is low in either the vector or host population, demographic stochasticity can have a large impact on whether the disease goes extinct before resulting in a major outbreak [17, 18]. Analysis of deterministic models often includes a computation of the basic reproduction number, $R_0$. If $R_0$ is greater than 1, introducing an infected individual into a completely susceptible population results in an outbreak, and if $R_0$ is less than 1, the disease will die out. With stochastic models, however, there is a possibility that the disease will die out before resulting in a major outbreak even when $R_0$ is greater than 1.

A useful stochastic framework for modeling disease dynamics when there are multiple types of infected individuals, such as hosts and vectors, when disease prevalence is low is continuous-time Markov chain (CTMC) modeling and multitype branching process approximation [18–29]. Allen and van den Driessche [30] and Maliyoni [23] both applied this methodology to models for WNV with a single species of host and vector, and found the highest probability of extinction to occur when the disease is introduced by an exposed mosquito and the lowest probability of extinction to occur when the disease is introduced by an infectious mosquito, with the probability of extinction for introduction by an infectious host falling in between. Lloyd et al. [18] determined extinction probabilities in a model for vector-borne disease with m host species and n vector species, with no latent class for the vector. For the case of two host species and a single vector, with hosts differing only in abundance and attractiveness to vectors, they found the probability of disease extinction when introduced by an infectious vector to be maximized when vectors bite hosts in proportion to their abundance. The effect of vector feeding preferences on the probability of disease extinction was dependent on whether the disease was introduced by an infectious vector or a randomly chosen infectious host, as well as the vector abundance.

In this paper we develop a stochastic CTMC model for WNV incorporating two host and one vector species, allowing the two types of hosts to differ in their ability to spread the disease, exposure to vectors, and relative abundance. In Section 2, we develop a deterministic ODE model and compute the basic reproduction number, $R_0$. In Section 3, we derive the CTMC model from the deterministic model along with expressions for the theoretical probability of extinction, ${\mathbb{P}}_0$, when the disease is introduced by an infectious host of species 1, infectious host of species 2, exposed vector, or infectious vector. We compute ${\mathbb{P}}_0$ for each infected type for the baseline parameter set and compare to the approximate probability of extinction, ${\mathbb{P}}_A$, computed from simulations. In Section 4, we explore how the probability of extinction (both theoretical and approximate) depends on factors known to influence the WNV risk of an area, specifically host community composition, vector feeding preferences, and the relative abundance of hosts and vectors. In Section 5 we conclude with a discussion.

2. Deterministic ODE model

In this section we develop a deterministic model for enzootic WNV transmission between two avian host species and a single species of mosquito using ordinary differential equations (ODEs). Hosts may differ in their competence (ability to spread WNV) as well as their exposure to vectors (also referred to throughout as vector feeding preference) and abundance (determined by birth and death rates). A schematic of the model is shown in Figure 1, and the model equations are: (1) $\begin{array}{rl}\frac{d{B}_S}{dt}& ={\mathrm{\Lambda }}_B-a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-{\mu }_B{B}_S\\ \frac{d{B}_I}{dt}& =a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-({\gamma }_B+{\nu }_B+{\mu }_B)B_I\\ \frac{d{B}_R}{dt}& ={\gamma }_B{B}_I-{\mu }_B{B}_R\\ \frac{d{\hat{B}}_S}{dt}& ={\mathrm{\Lambda }}_{\hat{B}}-a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-{\mu }_{\hat{B}}{\hat{B}}_S\\ \frac{d{\hat{B}}_I}{dt}& =a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_I\\ \frac{d{\hat{B}}_R}{dt}& ={\gamma }_{\hat{B}}{\hat{B}}_I-{\mu }_{\hat{B}}{\hat{B}}_R\\ \frac{d{M}_S}{dt}& ={\mathrm{\Lambda }}_M-a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-{\mu }_M{M}_S\\ \frac{d{M}_L}{dt}& =a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-k{M}_L-{\mu }_M{M}_L\\ \frac{d{M}_I}{dt}& =k{M}_L-{\mu }_M{M}_I\end{array}$(1) The number of susceptible, infectious, and recovered birds at time t are denoted by $B_S(t)$, $B_I(t)$, and $B_R(t)$, respectively, for species 1, and ${\hat{B}}_S(t)$, ${\hat{B}}_I(t)$, and ${\hat{B}}_R(t)$, respectively, for species 2. The number of susceptible, latent (or exposed), and infectious mosquitoes at time t are given by $M_S(t),M_L(t)$, and $M_I(t)$, respectively. We assume infectious mosquitoes do not live long enough to recover.

Mosquitoes bite birds at rate a, with total mosquito bites distributed among the two species. As in [11], we denote the fraction of bites going to host species 1 by $\alpha (t)$ and host species 2 by $\hat{\alpha }(t)$, modeled as follows: $\begin{array}{rl}\alpha (t)& =\frac{ϵB_T(t)}{ϵB_T(t)+\widehat{ϵ}{\hat{B}}_T(t)}\\ \hat{\alpha }(t)& =\frac{\widehat{ϵ}{\hat{B}}_T(t)}{ϵB_T(t)+\widehat{ϵ}{\hat{B}}_T(t)}\end{array}$where the subscript T denotes total. The total number of birds of host species 1 is $B_T(t)=B_S(t)+B_I(t)+B_R(t)$, and the total number of birds of host species 2 is ${\hat{B}}_T(t)={\hat{B}}_S(t)+{\hat{B}}_I(t)+{\hat{B}}_R(t)$. The parameters ϵ and $\widehat{ϵ}$ are exposure coefficients for host species 1 and 2, respectively, and assumed to be constant. If $ϵ=\widehat{ϵ}$ then both species are bitten in proportion to their abundance. A higher value of either ϵ or $\widehat{ϵ}$ would denote increased exposure to mosquito bites for one species relative to the other species, or a vector feeding preference for that species.

Figure 1. West Nile virus two-host model schematic. Transmission can occur when a susceptible vector bites an infectious host of species 1 or 2, or when an infectious vector bites a susceptible host of species 1 or 2.

Susceptible birds of host species 1 and 2 are recruited at rate ${\mathrm{\Lambda }}_B$ and ${\mathrm{\Lambda }}_{\hat{B}}$, respectively. WNV transmission may occur when a susceptible bird of either species is bitten by an infectious mosquito. Birds of host species 1 are bitten by infectious vectors at rate $a\alpha M_I$. A fraction of these bites, $\frac{B_S}{B_T}$, are on susceptible hosts, and each bite results in transmission with probability ${\beta }_B$. Similarly, birds of host species 2 are bitten by infectious vectors at rate $a\hat{\alpha }{M}_I$. A fraction of these bites, $\frac{{\hat{B}}_S}{{\hat{B}}_T}$, are on susceptible hosts, and each bite results in transmission with probability ${\beta }_{\hat{B}}$. Infectious birds of species 1 recover at rate ${\gamma }_B$, and infectious birds of species 2 recover at rate ${\gamma }_{\hat{B}}$. Birds exit all classes due to natural mortality at rate ${\mu }_B$ for species 1 and ${\mu }_{\hat{B}}$ for species 2, and infectious birds have increased death rate due to disease ${\nu }_B$ for species 1 and ${\nu }_{\hat{B}}$ for species 2.

Susceptible mosquitoes are recruited at rate ${\mathrm{\Lambda }}_M$. WNV transmission may occur when a susceptible mosquito bites an infectious bird of either species and subsequently moves into the latent class. A susceptible mosquito takes $a\alpha$ and $a\hat{\alpha }$ bites per day on hosts of species 1 and 2, respectively. A fraction ${\displaystyle \frac{B_I}{B_T}}$ and ${\displaystyle \frac{{\hat{B}}_I}{{\hat{B}}_T}}$ of those bites are on infectious hosts for species 1 and 2, respectively, and each bite results in transmission with probability ${\delta }_B$ and ${\delta }_{\hat{B}}$ for species 1 and 2, respectively. Latent mosquitoes move into the infectious class at rate k, where $\frac{1}{k}$ is the extrinsic incubation period of the virus, and mosquitoes in all classes have natural mortality rate ${\mu }_M$.

2.1. The disease-free equilibrium and $R_0$

The disease-free equilibrium (DFE) occurs in the absence of the disease, or when all infected populations are 0. The DFE for model (1(1) $\begin{array}{rl}\frac{d{B}_S}{dt}& ={\mathrm{\Lambda }}_B-a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-{\mu }_B{B}_S\\ \frac{d{B}_I}{dt}& =a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-({\gamma }_B+{\nu }_B+{\mu }_B)B_I\\ \frac{d{B}_R}{dt}& ={\gamma }_B{B}_I-{\mu }_B{B}_R\\ \frac{d{\hat{B}}_S}{dt}& ={\mathrm{\Lambda }}_{\hat{B}}-a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-{\mu }_{\hat{B}}{\hat{B}}_S\\ \frac{d{\hat{B}}_I}{dt}& =a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_I\\ \frac{d{\hat{B}}_R}{dt}& ={\gamma }_{\hat{B}}{\hat{B}}_I-{\mu }_{\hat{B}}{\hat{B}}_R\\ \frac{d{M}_S}{dt}& ={\mathrm{\Lambda }}_M-a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-{\mu }_M{M}_S\\ \frac{d{M}_L}{dt}& =a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-k{M}_L-{\mu }_M{M}_L\\ \frac{d{M}_I}{dt}& =k{M}_L-{\mu }_M{M}_I\end{array}$(1) ) is given by: (2) $(B_S^{\ast },B_I^{\ast },B_R^{\ast },{\hat{B}}_S^{\ast },{\hat{B}}_I^{\ast },{\hat{B}}_R^{\ast },M_S^{\ast },M_L^{\ast },M_I^{\ast })=(\frac{{\mathrm{\Lambda }}_B}{{\mu }_B},0,0,\frac{{\mathrm{\Lambda }}_{\hat{B}}}{{\mu }_{\hat{B}}},0,0,\frac{{\mathrm{\Lambda }}_M}{{\mu }_M},0,0)$(2) With the parameters in Table 1, there are 100 susceptible hosts of each species and 1041.7 susceptible mosquitoes at the DFE. We note that at the DFE, both α and $\hat{\alpha }$ are constant.

Table 1. Parameters for two-species West Nile virus model.

Parameter	Description	Value	Reference
${\mathrm{\Lambda }}_B$ (${\mathrm{\Lambda }}_{\hat{B}}$)	Recruitment rate of avian hosts of species 1 (species 2)	1 birds/day	Assumed
${\mathrm{\Lambda }}_M$	Recruitment rate of mosquitoes	100 mosquitoes/day	Assumed
a	Biting rate on avian hosts	0.33 bites/mosquito/day	[14]
ϵ ($\widehat{ϵ}$)	Host species 1 (species 2) exposure coefficient	varies	[13]
${\mu }_M$	Vector mortality rate	0.096 days${}^{-1}$	[31]
k	Virus extrinsic incubation period	0.106 days${}^{-1}$	[14]
${\gamma }_B$ (${\gamma }_{\hat{B}}$)	Host species 1 (species 2) recovery rate	1/3 days${}^{-1}$	[5]
${\mu }_B$ (${\mu }_{\hat{B}}$)	Host species 1 (species 2) natural mortality rate	0.01 days${}^{-1}$	[8, 32]
${\beta }_B$ (${\beta }_{\hat{B}}$)	Host species 1 (species 2) susceptibility	1	[5]
${\delta }_B$ (${\delta }_{\hat{B}}$)	Host species 1 (species 2) infectivity	0.36	[5]
${\nu }_B$ (${\nu }_{\hat{B}}$)	Host species 1 (species 2) virulence	0.1 days${}^{-1}$	[15]

For model (1(1) $\begin{array}{rl}\frac{d{B}_S}{dt}& ={\mathrm{\Lambda }}_B-a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-{\mu }_B{B}_S\\ \frac{d{B}_I}{dt}& =a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-({\gamma }_B+{\nu }_B+{\mu }_B)B_I\\ \frac{d{B}_R}{dt}& ={\gamma }_B{B}_I-{\mu }_B{B}_R\\ \frac{d{\hat{B}}_S}{dt}& ={\mathrm{\Lambda }}_{\hat{B}}-a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-{\mu }_{\hat{B}}{\hat{B}}_S\\ \frac{d{\hat{B}}_I}{dt}& =a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_I\\ \frac{d{\hat{B}}_R}{dt}& ={\gamma }_{\hat{B}}{\hat{B}}_I-{\mu }_{\hat{B}}{\hat{B}}_R\\ \frac{d{M}_S}{dt}& ={\mathrm{\Lambda }}_M-a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-{\mu }_M{M}_S\\ \frac{d{M}_L}{dt}& =a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-k{M}_L-{\mu }_M{M}_L\\ \frac{d{M}_I}{dt}& =k{M}_L-{\mu }_M{M}_I\end{array}$(1) ), we can compute the basic reproduction number, $R_0$, which gives the number of secondary infections produced by a typical infected individual introduced into a completely susceptible population [33]. If $R_0>1$, there will be an outbreak, while if $R_0<1$ the disease will die out and the population will return to the DFE. We compute the basic reproduction number using the next-generation matrix approach [34], defining the 4-by-4 matrices F and V: (3) $F=\left[\begin{array}{cccc}0& 0& 0& a\alpha {\beta }_B\\ 0& 0& 0& a\hat{\alpha }{\beta }_{\hat{B}}\\ {\displaystyle a\alpha {\delta }_B\frac{M_S^{\ast }}{B_S^{\ast }}}& {\displaystyle a\hat{\alpha }{\delta }_{\hat{B}}\frac{M_S^{\ast }}{{\hat{B}}_S^{\ast }}}& 0& 0\\ 0& 0& 0& 0\end{array}\right]$(3) and (4) $V=\left[\begin{array}{cccc}({\gamma }_B+{\nu }_B+{\mu }_B)& 0& 0& 0\\ 0& ({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}})& 0& 0\\ 0& 0& k+{\mu }_M& 0\\ 0& 0& -k& {\mu }_M\end{array}\right].$(4) The spectral radius of $F{V}^{-1}$, $\rho (F{V}^{-1})$, is the basic reproduction number, so (5) $R_0=\sqrt{\frac{a\alpha {\delta }_B{M}_S^{\ast }}{({\gamma }_B+{\nu }_B+{\mu }_B)B_S^{\ast }}\frac{a\alpha {\beta }_{B}k}{(k+{\mu }_M){\mu }_M}+\frac{a\hat{\alpha }{\delta }_{\hat{B}}{M}_S^{\ast }}{({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_S^{\ast }}\frac{a\hat{\alpha }{\beta }_{\hat{B}}k}{(k+{\mu }_M){\mu }_M}}.$(5) The first term under the square root in (5(5) $R_0=\sqrt{\frac{a\alpha {\delta }_B{M}_S^{\ast }}{({\gamma }_B+{\nu }_B+{\mu }_B)B_S^{\ast }}\frac{a\alpha {\beta }_{B}k}{(k+{\mu }_M){\mu }_M}+\frac{a\hat{\alpha }{\delta }_{\hat{B}}{M}_S^{\ast }}{({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_S^{\ast }}\frac{a\hat{\alpha }{\beta }_{\hat{B}}k}{(k+{\mu }_M){\mu }_M}}.$(5) ) represents transmission between vectors and host species 1, and the second term represents transmission between vectors and host species 2. $R_0$ can be written as $R_0=\sqrt{R_{BV}{R}_{VB}+R_{\hat{B}V}{R}_{V\hat{B}}}$where $R_{BV}=\frac{a\alpha {\delta }_B{M}_S^{\ast }}{({\gamma }_B+{\nu }_B+{\mu }_B)B_S^{\ast }}$is the number of new infections in vectors caused by a single infectious bird of species 1 and $R_{VB}=\frac{k}{k+{\mu }_M}\frac{a\alpha {\beta }_B}{{\mu }_M}$is the number of new infections in host species 1 caused by a single latent vector. There are similar terms for host species 2, where $R_{\hat{B}V}=\frac{a\hat{\alpha }{\delta }_{\hat{B}}{M}_S^{\ast }}{({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_S^{\ast }}$is the number of new infections in vectors caused by a single infectious bird of species 2 and $R_{V\hat{B}}=\frac{k}{k+{\mu }_M}\frac{a\hat{\alpha }{\beta }_{\hat{B}}}{{\mu }_M}$is the number of new infections in host species 2 caused by a single latent vector.

3. Stochastic model

3.1. CTMC model

Here we introduce a stochastic Continuous-Time Markov Chain (CTMC) model for WNV transmission between two host species and a single vector species; we assume this model is time-homogeneous and satisfies the Markov property, which requires that the future state of the process depends only on the current state [23, 35]. State variables and parameters are the same as in the deterministic model (1(1) $\begin{array}{rl}\frac{d{B}_S}{dt}& ={\mathrm{\Lambda }}_B-a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-{\mu }_B{B}_S\\ \frac{d{B}_I}{dt}& =a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-({\gamma }_B+{\nu }_B+{\mu }_B)B_I\\ \frac{d{B}_R}{dt}& ={\gamma }_B{B}_I-{\mu }_B{B}_R\\ \frac{d{\hat{B}}_S}{dt}& ={\mathrm{\Lambda }}_{\hat{B}}-a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-{\mu }_{\hat{B}}{\hat{B}}_S\\ \frac{d{\hat{B}}_I}{dt}& =a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_I\\ \frac{d{\hat{B}}_R}{dt}& ={\gamma }_{\hat{B}}{\hat{B}}_I-{\mu }_{\hat{B}}{\hat{B}}_R\\ \frac{d{M}_S}{dt}& ={\mathrm{\Lambda }}_M-a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-{\mu }_M{M}_S\\ \frac{d{M}_L}{dt}& =a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-k{M}_L-{\mu }_M{M}_L\\ \frac{d{M}_I}{dt}& =k{M}_L-{\mu }_M{M}_I\end{array}$(1) ), where $B_S(t),B_I(t),B_R(t),{\hat{B}}_S(t),{\hat{B}}_I(t),{\hat{B}}_R(t)$, $M_S(t),M_L(t),M_I(t)$ are now integer-valued random variables representing the previously discussed classes at time t. From the deterministic model (1(1) $\begin{array}{rl}\frac{d{B}_S}{dt}& ={\mathrm{\Lambda }}_B-a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-{\mu }_B{B}_S\\ \frac{d{B}_I}{dt}& =a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-({\gamma }_B+{\nu }_B+{\mu }_B)B_I\\ \frac{d{B}_R}{dt}& ={\gamma }_B{B}_I-{\mu }_B{B}_R\\ \frac{d{\hat{B}}_S}{dt}& ={\mathrm{\Lambda }}_{\hat{B}}-a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-{\mu }_{\hat{B}}{\hat{B}}_S\\ \frac{d{\hat{B}}_I}{dt}& =a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_I\\ \frac{d{\hat{B}}_R}{dt}& ={\gamma }_{\hat{B}}{\hat{B}}_I-{\mu }_{\hat{B}}{\hat{B}}_R\\ \frac{d{M}_S}{dt}& ={\mathrm{\Lambda }}_M-a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-{\mu }_M{M}_S\\ \frac{d{M}_L}{dt}& =a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-k{M}_L-{\mu }_M{M}_L\\ \frac{d{M}_I}{dt}& =k{M}_L-{\mu }_M{M}_I\end{array}$(1) ), we create a list of all state transitions; these include hosts and vectors moving through infected classes, such as a host changing from susceptible to infectious or infectious to recovered, or a vector changing from susceptible to latent or latent to infectious, as well as births and deaths of hosts and vectors. The rates at which these transitions occur are computed from model (1(1) $\begin{array}{rl}\frac{d{B}_S}{dt}& ={\mathrm{\Lambda }}_B-a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-{\mu }_B{B}_S\\ \frac{d{B}_I}{dt}& =a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-({\gamma }_B+{\nu }_B+{\mu }_B)B_I\\ \frac{d{B}_R}{dt}& ={\gamma }_B{B}_I-{\mu }_B{B}_R\\ \frac{d{\hat{B}}_S}{dt}& ={\mathrm{\Lambda }}_{\hat{B}}-a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-{\mu }_{\hat{B}}{\hat{B}}_S\\ \frac{d{\hat{B}}_I}{dt}& =a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_I\\ \frac{d{\hat{B}}_R}{dt}& ={\gamma }_{\hat{B}}{\hat{B}}_I-{\mu }_{\hat{B}}{\hat{B}}_R\\ \frac{d{M}_S}{dt}& ={\mathrm{\Lambda }}_M-a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-{\mu }_M{M}_S\\ \frac{d{M}_L}{dt}& =a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-k{M}_L-{\mu }_M{M}_L\\ \frac{d{M}_I}{dt}& =k{M}_L-{\mu }_M{M}_I\end{array}$(1) ) and given in Table 2.

Table 2. Transitions and Rates for CTMC Model.

Event	Change in State	Rate
Recruitment of $B_S$	(1,0,0,0,0,0,0,0,0)	${\mathrm{\Lambda }}_B$
Infection of $B_S$	(−1,1,0,0,0,0,0,0,0)	$a\alpha {\beta }_B{M}_I\frac{B_S}{B_T}$
Recovery of $B_I$	(0,−1,1,0,0,0,0,0,0)	${\gamma }_B{B}_I$
Death of $B_S$	(−1,0,0,0,0,0,0,0,0)	${\mu }_B{B}_S$
Death of $B_I$	(0,−1,0,0,0,0,0,0,0)	$({\mu }_B+{\nu }_B)B_I$
Death of $B_R$	(0,0,−1,0,0,0,0,0,0)	${\mu }_B{B}_R$
Recruitment of ${\hat{B}}_S$	(0,0,0,1,0,0,0,0,0)	${\mathrm{\Lambda }}_{\hat{B}}$
Infection of ${\hat{B}}_S$	(0,0,0,−1,1,0,0,0,0)	$a\hat{\alpha }{\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}$
Recovery of ${\hat{B}}_I$	(0,0,0,0,−1,1,0,0,0)	${\gamma }_{\hat{B}}{\hat{B}}_I$
Death of ${\hat{B}}_S$	(0,0,0,−1,0,0,0,0,0)	${\mu }_{\hat{B}}{\hat{B}}_S$
Death of ${\hat{B}}_I$	(0,0,0,0,−1,0,0,0,0)	$({\mu }_{\hat{B}}+{\nu }_{\hat{B}}){\hat{B}}_I$
Death of ${\hat{B}}_R$	(0,0,0,0,0,−1,0,0,0)	${\mu }_{\hat{B}}{\hat{B}}_R$
Recruitment of $M_S$	(0,0,0,0,0,0,1,0,0)	${\mathrm{\Lambda }}_M$
Infection of $M_S$	(0,0,0,0,0,0,−1,1,0)	$a(\alpha {\delta }_B\frac{B_I}{B_T}+\hat{\alpha }{\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S$
Progression of $M_L$ to $M_I$	(0,0,0,0,0,0,0,−1,1)	$k{M}_L$
Death of $M_S$	(0,0,0,0,0,0,−1,0,0)	${\mu }_M{M}_S$
Death of $M_L$	(0,0,0,0,0,0,0,−1,0)	${\mu }_M{M}_L$
Death of $M_I$	(0,0,0,0,0,0,0,0,−1)	${\mu }_M{M}_I$

Note: For each of ($B_S$, $B_I$, $B_R$, ${\hat{B}}_S$, ${\hat{B}}_I$, ${\hat{B}}_R$, $M_S$, $M_L$, $M_I$), a 1 represents an increase by 1, −1 represents a decrease by 1, and 0 represents no change in the variable's state from time t to time ($t+\mathrm{\Delta }t$).

3.2. Theoretical probability of disease extinction

Galton-Watson branching process theory provides a linear approximation of model dynamics near the disease-free equilibrium provided the susceptible populations are sufficiently large [20]. When the initial infected population is small (here we consider one infected individual in total), the branching process for the total number of infected individuals either reaches zero, resulting in disease extinction, after infecting only a small number of individuals, or grows exponentially, resulting in a major outbreak [17].

We can determine probability generating functions for the distribution of secondary infections for each infected class, called ‘offspring PGFs,’ that are used in determining the probability of disease extinction versus a major outbreak [36]. These offspring PGFs have the general form: (6) $f_i(u_1,\dots ,u_n)=\sum _{k_n=0}^{\mathrm{\infty }}\cdots \sum _{k_1=0}^{\mathrm{\infty }}{P}_i(k_1,\dots ,k_n)u_1^{k_1}\cdots u_n^{k_n}$(6) where $P_i(k_1,\dots ,k_n)$ is the probability that an infected individual of type i creates $k_j$ offspring of type j [17, 20]. We define offspring probability generating functions for the four types of infected individuals below, where $f_1$ is the PGF for infectious bird species 1 ($B_I$), $f_2$ is the PGF for infectious bird species 2 (${\hat{B}}_I$), $f_3$ is the PGF for exposed mosquitoes ($M_L$), and $f_4$ is the PGF for infectious mosquitoes ($M_I$). (7) $\begin{array}{rl}{f}_1(u_1,u_2,u_3,u_4)& =\frac{a\alpha \frac{M_S^{\ast }}{B_T^{\ast }}{\delta }_B{u}_1{u}_3+{\gamma }_B+{\mu }_B+{\nu }_B}{a\alpha {\delta }_B\frac{M_S^{\ast }}{B_T^{\ast }}+{\gamma }_B+{\mu }_B+{\nu }_B}\end{array}$(7) (8) $\begin{array}{rl}{f}_2(u_1,u_2,u_3,u_4)& =\frac{a\hat{\alpha }{\delta }_{\hat{B}}\frac{M_S^{\ast }}{{\hat{B}}_T^{\ast }}{u}_2{u}_3+{\gamma }_{\hat{B}}+{\mu }_{\hat{B}}+{\nu }_{\hat{B}}}{a\hat{\alpha }{\delta }_{\hat{B}}\frac{M_S^{\ast }}{{\hat{B}}_T^{\ast }}+{\gamma }_{\hat{B}}+{\mu }_{\hat{B}}+{\nu }_{\hat{B}}}\end{array}$(8) (9) $\begin{array}{rl}{f}_3(u_1,u_2,u_3,u_4)& =\frac{k{u}_4+{\mu }_M}{k+{\mu }_M}\end{array}$(9) (10) $\begin{array}{rl}{f}_4(u_1,u_2,u_3,u_4)& =\frac{a\alpha {\beta }_B{u}_1{u}_4+a\hat{\alpha }{\beta }_{\hat{B}}{u}_2{u}_4+{\mu }_M}{a\alpha {\beta }_B+a\hat{\alpha }{\beta }_{\hat{B}}+{\mu }_M}\end{array}$(10) An infectious bird of species 1 or 2 can either infect a susceptible mosquito, thereby creating two offspring (itself and the exposed mosquito), or it can die or recover before transmitting the disease, creating zero offspring. An exposed mosquito can either survive the latent period and become infectious, resulting in one offspring, or it can die before becoming infectious, resulting in zero infectious offspring. An infectious mosquito can transmit the disease to a bird of either species 1 or species 2, thereby producing two offspring (itself and the infectious host) or die before transmitting the disease, producing zero offspring.

From the PGFs (7(7) $\begin{array}{rl}{f}_1(u_1,u_2,u_3,u_4)& =\frac{a\alpha \frac{M_S^{\ast }}{B_T^{\ast }}{\delta }_B{u}_1{u}_3+{\gamma }_B+{\mu }_B+{\nu }_B}{a\alpha {\delta }_B\frac{M_S^{\ast }}{B_T^{\ast }}+{\gamma }_B+{\mu }_B+{\nu }_B}\end{array}$(7) )–(10(10) $\begin{array}{rl}{f}_4(u_1,u_2,u_3,u_4)& =\frac{a\alpha {\beta }_B{u}_1{u}_4+a\hat{\alpha }{\beta }_{\hat{B}}{u}_2{u}_4+{\mu }_M}{a\alpha {\beta }_B+a\hat{\alpha }{\beta }_{\hat{B}}+{\mu }_M}\end{array}$(10) ), we can define a 4-by-4 expectation matrix, $\mathbb{M}$, where the (j, i) entry $\frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{u}_j}$ evaluated at the fixed point $(u_1,u_2,u_3,u_4)=(1,1,1,1)$ gives the expected number of new infected individuals of type j created by an infected individual of type i [20]: (11) $\begin{array}{rl}\mathbb{M}& =[\begin{array}{cc}\frac{a\alpha {\delta }_B\frac{M_S^{\ast }}{B_T^{\ast }}}{a\alpha {\delta }_B\frac{M_S^{\ast }}{B_T^{\ast }}+{\gamma }_B+{\mu }_B+{\nu }_B}& 0\\ 0& {\displaystyle \frac{a\hat{\alpha }{\delta }_{\hat{B}}\frac{M_S^{\ast }}{{\hat{B}}_T^{\ast }}}{a\hat{\alpha }{\delta }_{\hat{B}}\frac{M_S^{\ast }}{{\hat{B}}_T^{\ast }}+{\gamma }_{\hat{B}}+{\mu }_{\hat{B}}+{\nu }_{\hat{B}}}}\\ {\displaystyle \frac{a\alpha {\delta }_B\frac{M_S^{\ast }}{B_T^{\ast }}}{a\alpha {\delta }_B\frac{M_S^{\ast }}{B_T^{\ast }}+{\gamma }_B+{\mu }_B+{\nu }_B}}& {\displaystyle \frac{a\hat{\alpha }{\delta }_{\hat{B}}\frac{M_S^{\ast }}{{\hat{B}}_T^{\ast }}}{a\hat{\alpha }{\delta }_{\hat{B}}\frac{M_S^{\ast }}{{\hat{B}}_T^{\ast }}+{\gamma }_{\hat{B}}+{\mu }_{\hat{B}}+{\nu }_{\hat{B}}}}\\ 0& 0\end{array}\\ & \begin{array}{cc}0& {\displaystyle \frac{a\alpha {\beta }_B}{a\alpha {\beta }_B+a\hat{\alpha }{\beta }_{\hat{B}}+{\mu }_m}}\\ 0& {\displaystyle \frac{a\hat{\alpha }{\beta }_{\hat{B}}}{a\alpha {\beta }_B+a\hat{\alpha }{\beta }_{\hat{B}}+{\mu }_m}}\\ 0& 0\\ {\displaystyle \frac{k}{k+{\mu }_m}}& {\displaystyle \frac{a\alpha {\beta }_B+a\hat{\alpha }{\beta }_{\hat{B}}}{a\alpha {\beta }_B+a\hat{\alpha }{\beta }_{\hat{B}}+{\mu }_m}}\end{array}].\end{array}$(11) $\mathbb{M}$ is non-negative and irreducible, and satisfies the relationship $(\mathbb{M}-\mathbb{I})W=F-V$, where $\mathbb{I}$ is the $4\times 4$ identity matrix, F and V are as in (3(3) $F=\left[\begin{array}{cccc}0& 0& 0& a\alpha {\beta }_B\\ 0& 0& 0& a\hat{\alpha }{\beta }_{\hat{B}}\\ {\displaystyle a\alpha {\delta }_B\frac{M_S^{\ast }}{B_S^{\ast }}}& {\displaystyle a\hat{\alpha }{\delta }_{\hat{B}}\frac{M_S^{\ast }}{{\hat{B}}_S^{\ast }}}& 0& 0\\ 0& 0& 0& 0\end{array}\right]$(3) ) and (4(4) $V=\left[\begin{array}{cccc}({\gamma }_B+{\nu }_B+{\mu }_B)& 0& 0& 0\\ 0& ({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}})& 0& 0\\ 0& 0& k+{\mu }_M& 0\\ 0& 0& -k& {\mu }_M\end{array}\right].$(4) ), and $W=\mathrm{diag}(a\alpha {\delta }_B\frac{M_S^{\ast }}{B_T^{\ast }}+{\gamma }_B+{\mu }_B+{\nu }_B,a\hat{\alpha }{\delta }_{\hat{B}}\frac{M_S^{\ast }}{{\hat{B}}_T^{\ast }}+{\gamma }_{\hat{B}}+{\mu }_{\hat{B}}+{\nu }_{\hat{B}},k+{\mu }_M,a\alpha {\beta }_B+a\hat{\alpha }{\beta }_{\hat{B}}+{\mu }_M$). Therefore the Threshold Theorem from [30] holds, guaranteeing $R_0<1(=1,>1)$ if and only if the spectral radius $\rho (\mathbb{M})<1(=1,>1)$.

If $\rho (\mathbb{M})\le 1$, $(1,1,1,1)$ is the only fixed point of the PGFs and the theoretical probability of disease extinction ${\mathbb{P}}_0$ is one. If $\rho (\mathbb{M})>1$, we have a supercritical branching process, and in addition to the fixed point $(1,1,1,1)$ there exists a fixed point $(q_1,q_2,q_3,q_4)$ of the PGFs such that $f_i(q_1,q_2,q_3,q_4)=q_i$ with $q_i<1$ [20, 30]. The probability of disease extinction in the population is then given by (12) ${\mathbb{P}}_0=q_1^{b_i}{q}_2^{{\hat{b}}_i}{q}_3^{m_l}{q}_4^{m_i}$(12) where $b_i=B_I(0)$, ${\hat{b}}_i={\hat{B}}_I(0)$, $m_l=M_L(0)$, $m_i=M_I(0)$, and the probability of a major outbreak is $1-{\mathbb{P}}_0$ [20]. While we are able to compute ${\mathbb{P}}_0$ if the disease is introduced by multiple infected individuals (of the same type or different types), we focus on the case where disease introduction occurs by a single infected individual. The fixed point $(q_1,q_2,q_3,q_4)$ of the PGFs (7(7) $\begin{array}{rl}{f}_1(u_1,u_2,u_3,u_4)& =\frac{a\alpha \frac{M_S^{\ast }}{B_T^{\ast }}{\delta }_B{u}_1{u}_3+{\gamma }_B+{\mu }_B+{\nu }_B}{a\alpha {\delta }_B\frac{M_S^{\ast }}{B_T^{\ast }}+{\gamma }_B+{\mu }_B+{\nu }_B}\end{array}$(7) )–(10(10) $\begin{array}{rl}{f}_4(u_1,u_2,u_3,u_4)& =\frac{a\alpha {\beta }_B{u}_1{u}_4+a\hat{\alpha }{\beta }_{\hat{B}}{u}_2{u}_4+{\mu }_M}{a\alpha {\beta }_B+a\hat{\alpha }{\beta }_{\hat{B}}+{\mu }_M}\end{array}$(10) ) with $q_i<1$ was computed numerically in MATLAB (using fsolve) for the baseline parameters in Table 1, chosen based on the vector Culex pipiens and avian host American robin. The corresponding probabilities of disease extinction, ${\mathbb{P}}_0$, for introduction by a single infectious individual of each type are given in Table 3.

Table 3. Theoretical (${\mathbb{P}}_0$) and approximate (${\mathbb{P}}_A$) probabilities of disease extinction for disease introduction by a single infectious host of species 1 ($b_i$), infectious host of species 2 (${\hat{b}}_i$), exposed vector ($m_l$), and infectious vector ($m_i$) under the baseline parameters in Table 1.

$b_i$	${\hat{b}}_i$	$m_l$	$m_i$	${\mathbb{P}}_0$	${\mathbb{P}}_A$
1	0	0	0	0.745	0.753
0	1	0	0	0.745	0.753
0	0	1	0	0.755	0.769
0	0	0	1	0.533	0.533

3.3. Approximate probability of disease extinction

A numerical approximation of the probability of extinction, denoted ${\mathbb{P}}_A$, can also be computed from simulations of the stochastic CTMC model. Simulations were run in MATLAB and began with a single infected individual introduced into an otherwise susceptible population at the disease free equilibrium. For each simulation, either the number of infected individuals reached zero after only a small total number of infections, or a major outbreak occurred. We classified any simulation with an outbreak size less than $10\mathrm{\%}$ of the total population as disease extinction. For the parameter values in Table 1, this threshold is 124.2 individuals. To find the approximate probability of extinction (${\mathbb{P}}_A$), we computed the fraction of 10, 000 simulations for which the number of infected individuals reaches zero before exceeding an outbreak size of 124.

The theoretical and approximate probabilities of WNV extinction for the CTMC model with the baseline parameter values in Table 1 are given in Table 3 for disease introduction by a single infectious host of species 1 ($b_i$), a single infectious host of species 2 (${\hat{b}}_i$), a single latent (or exposed) mosquito ($m_l$), and a single infectious mosquito ($m_i$). The basic reproduction number $R_0$ for the deterministic model is greater than one ($R_0=1.587$) and does not depend on how the disease is introduced. Simulations of the deterministic model (1(1) $\begin{array}{rl}\frac{d{B}_S}{dt}& ={\mathrm{\Lambda }}_B-a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-{\mu }_B{B}_S\\ \frac{d{B}_I}{dt}& =a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-({\gamma }_B+{\nu }_B+{\mu }_B)B_I\\ \frac{d{B}_R}{dt}& ={\gamma }_B{B}_I-{\mu }_B{B}_R\\ \frac{d{\hat{B}}_S}{dt}& ={\mathrm{\Lambda }}_{\hat{B}}-a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-{\mu }_{\hat{B}}{\hat{B}}_S\\ \frac{d{\hat{B}}_I}{dt}& =a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_I\\ \frac{d{\hat{B}}_R}{dt}& ={\gamma }_{\hat{B}}{\hat{B}}_I-{\mu }_{\hat{B}}{\hat{B}}_R\\ \frac{d{M}_S}{dt}& ={\mathrm{\Lambda }}_M-a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-{\mu }_M{M}_S\\ \frac{d{M}_L}{dt}& =a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-k{M}_L-{\mu }_M{M}_L\\ \frac{d{M}_I}{dt}& =k{M}_L-{\mu }_M{M}_I\end{array}$(1) ) under these baseline parameter values will always result in disease outbreak while simulations of the CTMC model may result in either disease extinction or a major outbreak. For this parameter set, we find the probability of disease extinction is lowest when the disease is introduced by an infectious mosquito and highest when introduced by an exposed mosquito. There is no difference in the probability of extinction when the disease is introduced by a host of species 1 versus species 2, as expected since there is no vector preference and host species are equal in abundance and competence. We find that the approximate probability of extinction, ${\mathbb{P}}_A$, closely matches the theoretical probability of extinction, ${\mathbb{P}}_0$, for introduction by all types of infected individuals.

4. Impact of host and vector community composition on probability of disease extinction

In addition to the type of infected individual introducing the disease, the probability of disease extinction may depend on properties of the host and vector community into which the disease is entering. We next explore how the probability of extinction of WNV, when introduced into a community of hosts and vectors by each type of infected individuals, is affected by differences in competence, exposure to vectors, and relative abundance of the two host species, as well as overall vector abundance. For all results in this section, ${\mathbb{P}}_0$ and ${\mathbb{P}}_A$ are computed as described in Section 3.

4.1. Host competence and exposure

4.1.1. Differential exposure

To explore the effect of differential exposure of host species to vectors, or vector feeding preferences, we set the exposure coefficient of host species 2, $\widehat{ϵ}$, to 1 and allow the exposure coefficient of host species 1, ϵ, to vary between 1 and 10. When the host species are equally abundant, the percentage of bites on species 1 ranges from 50% when $ϵ=1$ and there is no host preference to over 90% when $ϵ=10$ and species 1 is strongly preferred. Figure 2 shows the basic reproduction number, $R_0$, for the deterministic model (1(1) $\begin{array}{rl}\frac{d{B}_S}{dt}& ={\mathrm{\Lambda }}_B-a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-{\mu }_B{B}_S\\ \frac{d{B}_I}{dt}& =a\alpha (t){\beta }_B{M}_I\frac{B_S}{B_T}-({\gamma }_B+{\nu }_B+{\mu }_B)B_I\\ \frac{d{B}_R}{dt}& ={\gamma }_B{B}_I-{\mu }_B{B}_R\\ \frac{d{\hat{B}}_S}{dt}& ={\mathrm{\Lambda }}_{\hat{B}}-a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-{\mu }_{\hat{B}}{\hat{B}}_S\\ \frac{d{\hat{B}}_I}{dt}& =a\hat{\alpha }(t){\beta }_{\hat{B}}{M}_I\frac{{\hat{B}}_S}{{\hat{B}}_T}-({\gamma }_{\hat{B}}+{\nu }_{\hat{B}}+{\mu }_{\hat{B}}){\hat{B}}_I\\ \frac{d{\hat{B}}_R}{dt}& ={\gamma }_{\hat{B}}{\hat{B}}_I-{\mu }_{\hat{B}}{\hat{B}}_R\\ \frac{d{M}_S}{dt}& ={\mathrm{\Lambda }}_M-a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-{\mu }_M{M}_S\\ \frac{d{M}_L}{dt}& =a(\alpha (t){\delta }_B\frac{B_I}{B_T}+\hat{\alpha }(t){\delta }_{\hat{B}}\frac{{\hat{B}}_I}{{\hat{B}}_T})M_S-k{M}_L-{\mu }_M{M}_L\\ \frac{d{M}_I}{dt}& =k{M}_L-{\mu }_M{M}_I\end{array}$(1) ) and both the theoretical probability of disease extinction, ${\mathbb{P}}_0$, and approximate probability of disease extinction, ${\mathbb{P}}_A$, for introduction by the four types of infected individuals (host species 1, host species 2, exposed vector, and infectious vector) as a function of host species 1 exposure, ϵ. We consider three cases: the host species are equally competent (Figure 2(b)), host species 1 is more competent (Figure 2(c)), and host species 2 is more competent (Figure 2(d)). The average competence of the host population remains the same for all cases.

Figure 2. The basic reproduction number, $R_0$, and probabilities of disease extinction are graphed as a function of host species 1 exposure coefficient (ϵ). In (b) ${\delta }_B={\delta }_{\hat{B}}=0.36$, (c) ${\delta }_B=0.56$, ${\delta }_{\hat{B}}=0.16$, and (d) ${\delta }_B=0.16$, ${\delta }_{\hat{B}}=0.56$. Solid curves denote the theoretical probability of disease extinction, ${\mathbb{P}}_0$, and circles denote the approximate probability of disease extinction, ${\mathbb{P}}_A$. All other parameters are as in Table 1.

When host species are equally exposed to vectors ($ϵ=1$), $R_0$ is equal for all cases, but ${\mathbb{P}}_0$ varies. If host species are equally competent, ${\mathbb{P}}_0$ for introduction by either host species is 0.745, in between ${\mathbb{P}}_0$ for introduction by an exposed vector (0.755) and infectious vector (0.533). When species differ in competence, ${\mathbb{P}}_0$ is highest for introduction by the less competent host species (0.874), and ${\mathbb{P}}_0$ for introduction by the more competent host species is 0.665, falling in between ${\mathbb{P}}_0$ for introduction by an exposed vector (0.768) and infectious vector (0.558), both higher than when hosts are equally competent.

As ϵ increases, ${\mathbb{P}}_0$ decreases if the disease is introduced by host species 1 (the more exposed or preferred species) and increases if the disease is introduced by host species 2 (the less exposed or non-preferred species). This trend holds regardless of host competence.

When hosts are equally competent (Figure 2(b)), increasing ϵ results in increased values of $R_0$ along with decreased values of ${\mathbb{P}}_0$ for disease introduction by an exposed vector, infectious vector, or host of species 1, and increased values of ${\mathbb{P}}_0$ for disease introduction by hosts of species 2.

The case where the preferred host is also more competent (Figure 2(c)) results in the highest overall $R_0$ values along with the most extreme values of ${\mathbb{P}}_0$ for introduction by infectious hosts (lowest for introduction by species 1 and highest for species 2) and the lowest probabilities of extinction for introduction by an exposed or infectious vector.

When the preferred host is less competent (Figure 2(d)), $R_0$ first decreases with ϵ, as more bites are taken on the less competent species, before increasing with ϵ as bites become highly concentrated on a single host species. This non-monotonic nature of $R_0$ is reflected in ${\mathbb{P}}_0$ for introduction by an exposed or infectious vector; ${\mathbb{P}}_0$ first increases before eventually decreasing with ϵ. We still see a monotonic decrease in ${\mathbb{P}}_0$ with ϵ for introduction by species 1, the preferred species, and a monotonic increase in ${\mathbb{P}}_0$ with ϵ for the introduction by species 2, the non-preferred species, as in the other two cases.

When hosts are equally competent or species 1 is more competent, ${\mathbb{P}}_0$ is always greater for introduction by host species 2 than host species 1. However, when species 2 is more competent there is a range of ϵ values ($1<ϵ<3.5$) for which ${\mathbb{P}}_0$ is greater for introduction by species 1, the preferred but less competent species, than for introduction by species 2, the more competent but non-preferred species. For all cases, ${\mathbb{P}}_0$ is lowest when the disease is introduced by infectious vectors.

4.1.2. Differential competence

We next explore how $R_0$ and ${\mathbb{P}}_0$ are affected by the competence of the host community. We vary the probability of disease transmission from infectious hosts of species 1 to susceptible vectors, ${\delta }_B$, from 0 to 1 while keeping the infectivity of species 2, ${\hat{\delta }}_B$, at its baseline value of 0.36 (Figure 3). We consider three cases: vectors have no host preference, host species 1 is preferred, and host species 2 is preferred (assuming the preferred species has a 10-fold increase in its exposure coefficient relative to the non-preferred species). Hosts are assumed to be equally abundant, so the average infectivity of the host community increases from 0.18 when ${\delta }_B=0$ to 0.68 when ${\delta }_B=1$.

Figure 3. The basic reproduction number, $R_0$, and the probabilities of disease extinction are graphed as a function of host species 1 infectivity (${\delta }_B$). In (b) $ϵ=\widehat{ϵ}=1$, (c) $ϵ=10$, $\widehat{ϵ}=1$, and (d) $ϵ=1$, $\widehat{ϵ}=10$. Solid curves denote the theoretical probability of disease extinction, ${\mathbb{P}}_0$, and circles denote the approximate probability of disease extinction, ${\mathbb{P}}_A$. All other parameters are as in Table 1.

$R_0$ increases with ${\delta }_B$ for all cases. The lowest overall $R_0$ value (0.20) occurs when species 1 is preferred and non-competent (${\delta }_B=0$), and the greatest overall $R_0$ value (3.41) occurs when species 1 is preferred with bites on an infectious host guaranteed to result in disease transmission (${\delta }_B=1$). When species 2 is preferred, $R_0$ changes minimally with ${\delta }_B$, only increasing from 2.04 for ${\delta }_B=0$ to 2.07 for ${\delta }_B=1$.

When species 1 is preferred and ${\delta }_B$ is small enough, $R_0$ is below 1 and introduction of the disease by any infected individual results in extinction with probability 1 (${\mathbb{P}}_0=1$; Figure 3(c)). Once $R_0$ exceeds 1, ${\mathbb{P}}_0$ decreases with ${\delta }_B$ for introduction by any infected individual regardless of differences in host exposure. Though introduction of the disease by an infectious vector again (as in Figure 2) results in the lowest ${\mathbb{P}}_0$ for almost all combinations of host exposure and infectivity, when species 1 is both preferred and highly competent (${\delta }_B>0.85$) the lowest ${\mathbb{P}}_0$ occurs for introduction by species 1. When species 2 is preferred, ${\delta }_B$ has negligible impact on ${\mathbb{P}}_0$ unless the disease is introduced by species 1.

When hosts are equally exposed, ${\mathbb{P}}_0$ is higher when the disease is introduced by the less competent species (which is species 1 if ${\delta }_B<0.36$ and species 2 if ${\delta }_B>0.36$; Figure 3(b)) than by the more competent species. When there is a vector feeding preference, the highest values of ${\mathbb{P}}_0$ occur when the disease is introduced by the non-preferred species (Figure 3(c,d)).

4.1.3. Differential exposure and competence

We next vary both species 1 infectivity (${\delta }_B$) and exposure (ϵ), again plotting $R_0$ and ${\mathbb{P}}_0$ for introduction by the four types of infected individuals (Figure 4). Species 2 infectivity and exposure remain at baseline values (${\delta }_{\hat{B}}=0.36$, $\widehat{ϵ}=1$). Horizontal cross-sections of Figure 4 at $ϵ=1$ and $ϵ=10$ recover Figure 3(b,c), respectively, and a vertical cross-section at ${\delta }_B=0.36$ recovers Figure 2(b).

Figure 4. The basic reproduction number, $R_0$, and theoretical probability of disease extinction, ${\mathbb{P}}_0$, are graphed as a function of host species 1 infectivity (${\delta }_B$) and exposure coefficient (ϵ). In all plots the black contour line indicates where $R_0=1$. All other parameters are as in Table 1.

The average competence of the host community again increases with ${\delta }_B$. $R_0$ is below 1 for ${\delta }_B$ small enough and ϵ high enough (black contour lines denote the curve where $R_0=1$), so ${\mathbb{P}}_0=1$ for introduction by all types of infected individuals. $R_0$ again increases with ${\delta }_B$ for all values of ϵ, and ${\mathbb{P}}_0$ decreases with ${\delta }_B$ for introduction by all types of infected individuals (once $R_0$ exceeds 1). $R_0$ is maximized when species 1 exposure and competence are maximized ($ϵ=10$, ${\delta }_B=1$), and this is also where ${\mathbb{P}}_0$ is minimized when the disease is introduced by exposed vectors, infectious vectors, or hosts of species 1. When introduced by species 2, however, ${\mathbb{P}}_0$ is minimized when ${\delta }_B=1$ and $ϵ=1$; these parameters maximize the overall community competence as well as the vector biting rate on species 2.

4.2. Relative abundance of hosts and vectors

In this section we vary the percentage of hosts in species 1 versus species 2 while keeping the total host population size constant at 200 birds (${\mathrm{\Lambda }}_B+{\mathrm{\Lambda }}_{\hat{B}}=2$). We compute $R_0$ and ${\mathbb{P}}_0$ for introduction by each of the 4 types of infected individuals as the percentage of hosts in species 1 increases from 1 to 99 (${\mathrm{\Lambda }}_B$ increases from 0.02 to 1.98, with ${\mathrm{\Lambda }}_{\hat{B}}=2-{\mathrm{\Lambda }}_B$), and the exposure coefficient of species 1 (ϵ) varies from 1 to 10 (with $\widehat{ϵ}=1$). Results are shown for three cases: host species equally competent (Figure 5, left column), species 1 more competent (Figure 5, middle column), and species 2 more competent (Figure 5, right column).

Figure 5. The basic reproduction number, $R_0$, and the theoretical probability of disease extinction, ${\mathbb{P}}_0$, are graphed as a function of host species 1 exposure coefficient (ϵ) and the percentage of hosts in species 1. In the left column: ${\delta }_B=0.36,{\delta }_{\hat{B}}=0.36$. In the middle column: ${\delta }_B=0.56,{\delta }_{\hat{B}}=0.16$. In the right column: ${\delta }_B=0.16,{\delta }_{\hat{B}}=0.56$. All other parameters are as in Table 1.

4.2.1. Hosts equally competent

When hosts are equal in competence (Figure 5, left column) and exposure to vectors ($ϵ=1$), changing the relative abundance of hosts of each species does not affect $R_0$ or ${\mathbb{P}}_0$. $R_0$ increases with ϵ, and is maximized when the exposure coefficient of species 1 is maximized ($ϵ=10$) and there is a small percentage of hosts in species 1 receiving a large percentage of vector bites.

${\mathbb{P}}_0$ for introduction by species 1 decreases with ϵ and increases with the percentage of hosts in species 1, for $ϵ>1$. While the percentage of vector bites on species 1 increases with the relative abundance of species 1, the per capita biting rate on individuals of species 1 decreases. As the number of hosts in species 1 increases, it becomes less likely that a mosquito will bite the host that is infectious rather than another host of species 1, thus increasing the probability of disease extinction. When the disease is introduced by species 2, the non-preferred species, ${\mathbb{P}}_0$ increases with ϵ. ${\mathbb{P}}_0$ for introduction by species 2 also increases as the percentage of hosts in species 1 increases (and the percentage of hosts in species 2 decreases), for $ϵ>1$.

When disease introduction occurs via an exposed or infectious vector, ${\mathbb{P}}_0$ decreases with ϵ. When species 1 is preferred, ${\mathbb{P}}_0$ has a non-monotonic response to increasing the percentage of hosts in species 1 as seen for $R_0$, first decreasing and then increasing.

4.2.2. Hosts species 1 more competent

When species 1 is more competent than species 2 (Figure 5, middle column), the average host competence increases with the percentage of hosts in species 1. $R_0$ is again maximized when there are a small number of highly exposed hosts of species 1. The maximum value of $R_0$ is greater compared to the case where hosts are equally competent, even though the average host competence is lower, since the hosts receiving a high number of bites are those with higher competence.

When there is no host preference ($ϵ=1$), ${\mathbb{P}}_0$ decreases as the percentage of hosts in species 1 increases (along with the average competence of the host community) for introduction by all types of infectious individuals. ${\mathbb{P}}_0$ again decreases with ϵ when the disease is introduced by species 1 and increases with ϵ when the disease is introduced by species 2, as long as the percentage of hosts in species 1 is high enough. If the percentage of hosts in species 1 is low enough, ${\mathbb{P}}_0$ decreases with ϵ when introduced by species 2, the non-preferred and less competent species. For introduction by both host species 1 and 2, ${\mathbb{P}}_0$ has a non-monotonic response to increasing the percentage of hosts in species 1 when there is a strong enough preference for species 1, first decreasing and then increasing.

4.2.3. Hosts species 2 more competent

When species 1 is less competent than species 2 (Figure 5, right column), the average host competence now decreases with the percentage of hosts in species 1. The maximum value of $R_0$ is lower than for the case where hosts are equally competent even though the average host competence is higher, as the hosts being bitten at a higher rate are those with lower competence. As in Figure 2(d), the effect of increasing ϵ is non-monotonic for $R_0$ (as well as ${\mathbb{P}}_0$ for disease introduction by exposed or infectious vectors), since increasing the biting rate on less competent hosts first dilutes transmission, and eventually concentrating bites on a single species amplifies transmission.

4.2.4. Vector abundance

Finally, we explore how increasing vector abundance, and the vector-host ratio, affects the probability of disease extinction (Figure 6). The vector birthrate ${\mathrm{\Lambda }}_M$ is varied between its baseline value of 100 and 2500 mosquitoes per day, resulting in the equilibrium number of vectors in the population varying from $1,042-26,042$, or from about $5-130$ times the number of hosts in the population. $R_0$ and ${\mathbb{P}}_0$ for disease introduction by the 4 types of infected individuals are shown in Figure 6 for three cases: no host preference, host species 1 preferred, and host species 2 preferred. For all cases, hosts are equally abundant with species 1 more competent than species 2.

Figure 6. The basic reproduction number, $R_0$, and the probabilities of disease extinction are graphed as a function of vector birthrate, ${\mathrm{\Lambda }}_M$. In (b) $ϵ=\widehat{ϵ}=1$, (c) $ϵ=10$, $\widehat{ϵ}=1$, and (d) $ϵ=1$, $\widehat{ϵ}=10$. Solid curves denote the theoretical probability of disease extinction, ${\mathbb{P}}_0$, and circles denote the approximate probability of disease extinction, ${\mathbb{P}}_A$. For all plots, ${\delta }_B=0.56$ and ${\hat{\delta }}_B=0.16$. All other parameters are as in Table 1.

For all cases, $R_0$ increases with ${\mathrm{\Lambda }}_M$ (and vector abundance), reaching its highest values when the more competent species is preferred and lowest values when the less competent species is preferred. The probability of disease extinction decreases with ${\mathrm{\Lambda }}_M$ for introduction by all types of infected individuals. As vector abundance increases, ${\mathbb{P}}_0$ saturates for introduction by an exposed or infectious vector, approaching the same limiting values regardless of vector feeding preferences. At higher levels of vector abundance, ${\mathbb{P}}_0$ is more sensitive to changes in ${\mathrm{\Lambda }}_M$ when the disease is introduced by an infectious host compared to an exposed or infectious vector.

At the baseline value ${\mathrm{\Lambda }}_M=100$, ${\mathbb{P}}_0$ is lowest when the disease is introduced by an infectious vector, for all assumptions regarding vector feeding preferences. However, once ${\mathrm{\Lambda }}_M$ is large enough, the lowest ${\mathbb{P}}_0$ occurs for introduction by an infectious host. When there is no host preference, ${\mathbb{P}}_0$ for introduction by species 1 (the more competent species) and species 2 (the less competent species) drop below introduction by an infectious vector once ${\mathrm{\Lambda }}_M$ reaches 228 and 1177, respectively. When species 1, the more competent species, is also preferred, ${\mathbb{P}}_0$ for introduction by host species 1 drops below ${\mathbb{P}}_0$ for introduction by an infectious vector for a lower value of ${\mathrm{\Lambda }}_M$ (152), but ${\mathbb{P}}_0$ for introduction by species 2 remains above introduction by an infectious vector for the entire range of ${\mathrm{\Lambda }}_M$ considered. When species 2, the less competent species, is preferred, ${\mathbb{P}}_0$ is lower for introduction by species 2 than for the more competent species 1. ${\mathbb{P}}_0$ is lowest for introduction by species 2 once ${\mathrm{\Lambda }}_M$ exceeds 562 and ${\mathbb{P}}_0$ for introduction by species 1 also drops below ${\mathbb{P}}_0$ for infectious vectors once ${\mathrm{\Lambda }}_M$ is 1949.

5. Discussion

The relationship between $R_0$ and the probability of disease extinction is complex for diseases with multiple types of infected individuals [18, 20, 30]. For the standard SIR model, with a single class of infectious individual, the probability the disease will go extinct quickly after introduction, prior to causing a major outbreak, is $\frac{1}{R_0}$ [20, 30] and decreases monotonically with $R_0$. For the multi-host WNV model considered in this paper, there are four types of infected individuals – infectious birds of species 1 and 2, exposed mosquitoes, and infectious mosquitoes – and therefore four possible ways the disease could be introduced into a susceptible population. When vector preference for one host species increases, $R_0$ also increases as long as the preferred species is also the more competent species, or if both species are equally competent. When the preferred species is less competent, we see a non-monotonic effect of increasing host preference on $R_0$; $R_0$ first decreases, as more bites are taken on the less competent species, then increases as bites become more concentrated on the preferred hosts. When the disease is introduced by an exposed or infectious vector, we largely see an inverse relationship between the probability of disease extinction in the stochastic model and $R_0$ for the deterministic model. When the disease is introduced by a host, however, the probability of disease extinction can increase or decrease with $R_0$ depending on whether the species of bird introducing the disease is more competent, more abundant, and/or is bitten by mosquitoes at a higher rate than other avian species in the community. As preference for one species becomes stronger, the probability of disease extinction typically decreases if the disease is introduced by the preferred species, but increases if the disease is introduced by the non-preferred species.

Allen and van den Driessche [30] and Maliyoni [23] previously determined the probability of extinction for stochastic CTMC models of WNV transmission between a single host and vector, finding disease introduction by exposed vectors to have the highest probability of extinction and disease introduction by infectious vectors to have the lowest, for the parameter values considered. Our results agree in that the probability of extinction is always higher when the disease is introduced by an exposed vector compared to when the disease is introduced by an infectious vector, since the exposed vector must survive the extrinsic incubation period of the disease before it can begin to cause secondary infections. This in fact follows directly from the equations for the fixed points of the PGFs. The equation $f_3(q_1,q_2,q_3,q_4)=q_3$ results in $q_3>q_4$ when $q_3<1$ and ${\mu }_M>0$. For our baseline parameter values, where host species are identical, our findings are also in agreement that the highest overall probability of extinction occurs when the disease is introduced by an exposed vector and the lowest probability of extinction occurs when the disease is introduced by an infectious vector. Our model assumes that the vector biting rate does not depend on the number of hosts present. If the disease is introduced by an infectious vector, that vector will, on average, take a bite on some host every $\frac{1}{a}\approx 3$ days provided it does not die first. That bite will always result in transmission of the disease since avian susceptibility (${\beta }_B$ and ${\beta }_{\hat{B}}$) is assumed to be 1 for both species [5]. When the disease is introduced by an infectious host, however, the chance of that host being bitten by a susceptible vector depends on the number of susceptible vectors, as well as the other hosts present that those vectors may bite instead. The infectivity of that host then determines the probability that the bite results in disease transmission to the vector, with the baseline infectivity (${\delta }_B$ and ${\delta }_{\hat{B}}$) set to 0.36 for both species, based on laboratory estimates for the American robin [5]. Other potential hosts may have higher infectivity, such as the house sparrow, or lower infectivity, such as the mourning dove [5].

As vector abundance increases, we see the lowest probability of extinction change from disease introduction by an infectious vector to introduction by an avian host. Increasing vector abundance does not change the biting rate of each individual vector; if the disease is introduced by an infectious vector, the vector has the same chance of dying before transmitting the disease regardless of how many hosts or other vectors there are. When the disease is introduced by an infectious host, on the other hand, increasing vector abundance increases the expected number of bites on each host, along with the number of potential secondary infections. Thus we see the probability of disease extinction saturate with vector abundance when the disease is introduced by an exposed or infectious vector, but continue to decrease with vector abundance when introduced by an infectious host.

Lloyd et al. [18] developed a stochastic CTMC model for vector-borne disease with multiple host species, with results focused on the case of two host species differing only in relative abundance and vector biting rates. They found the effect of vector feeding preferences on the probability of extinction to depend on whether the disease was introduced by a vector or a randomly chosen host, and the probability of extinction for introduction by an infectious vector to be maximized when hosts are bitten in proportion to abundance. When exploring the effect of vector feeding preferences on equally competent host species in our model, we also see the probability of disease extinction maximized when hosts are equally exposed to vectors, for introduction by both exposed and infectious vectors. When the disease is introduced by a host, we find that as the vector biting rate on one host species increases, the probability of extinction decreases if the disease is introduced by the preferred species, but increases if it is introduced by the non-preferred species, regardless of differences in abundance. Considering only the probability of extinction for introduction by a randomly chosen host can mask the differences in the probability of extinction for introduction by individuals of any one species. Lloyd et al. [18] also found vector feeding preferences to have a minimal effect on the probability of disease extinction when introduced by a vector at high vector abundance. Similarly, we find the probability of extinction when introduced by an exposed or infectious vector to approach the same limiting values as vector abundance increases regardless of differences in host exposure.

Our theoretical and numerical results show that the probability an enzootic WNV outbreak occurs in an area depends strongly on both the avian host and vector community in that area as well as the means by which the disease is introduced into the population. WNV might be introduced into new communities at any point during the transmission season by local movement of either birds or mosquitoes, or from longer distances by migratory birds [12]. The presence and abundance of avian and vector species in any given location depend on many factors, including climate, the type of habitat available (i.e. urban, forest, wetlands), and land use [37–40]. The community composition of hosts and vectors can also vary temporally in a given location, causing the probability of disease extinction to vary depending on when the disease is introduced. In areas where vectors are not present year round, vector abundance may be lower during periods of initial growth in the spring and at the end of the transmission season in fall, compared to summer when temperatures are warmer. The composition of avian host species, and biting rates on those species, may vary throughout the year due to reproduction as well as the arrival and departure of migratory species, with the late summer migration of preferred avian species such as the American robin potentially resulting in increased bites on alternative hosts such as humans and other mammals [41].

Acknowledgments

High Performance Computing resources provided by the High Performance Research Computing (HPRC) core facility at Virginia Commonwealth University (https://hprc.vcu.edu) were used for conducting the research reported in this work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by a grant from the Simons Foundation (426126, SR).