Modelling Wolbachia infection in a sex-structured mosquito population carrying West Nile virus

Wolbachia is possibly the most studied reproductive parasite of arthropod species. It appears to be a promising candidate for biocontrol of some mosquito borne diseases. We begin by developing a sex-structured model for a Wolbachia infected mosquito population. Our model incorporates the key effects of Wolbachia infection including cytoplasmic incompatibility and male killing. We also allow the possibility of reduced reproductive output, incomplete maternal transmission, and different mortality rates for uninfected/infected male/female individuals. We study the existence and local stability of equilibria, including the biologically relevant and interesting boundary equilibria. For some biologically relevant parameter regimes there may be multiple coexistence steady states including, very importantly, a coexistence steady state in which Wolbachia infected individuals dominate. We also extend the model to incorporate West Nile virus (WNv) dynamics, using an SEI modelling approach. Recent evidence suggests that a particular strain of Wolbachia infection significantly reduces WNv replication in Aedes aegypti. We model this via increased time spent in the WNv-exposed compartment for Wolbachia infected female mosquitoes. A basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}R0 is computed for the WNv infection. Our results suggest that, if the mosquito population consists mainly of Wolbachia infected individuals, WNv eradication is likely if WNv replication in Wolbachia infected individuals is sufficiently reduced.


Introduction
Wolbachia is a maternally transmitted intracellular symbiont, and it is the most common reproductive parasite infecting a significant proportion of insect species, see e.g. O'Neill et al. (1997), Werren (1997). Wolbachia typically inhibits testes and ovaries of its host, and it is also present in its host's eggs. It interferes with its host's reproductive mechanism in a remarkable fashion. This allows Wolbachia to successfully establish itself in a number of arthropod species. Well-known effects of Wolbachia infections include cytoplasmic incompatibility (CI for short) and feminization of genetic males also known as male killing (MK for short), see e.g. Caspari and Watson (1959), Hoffmann and Turelli (1997), Telschow et al. (2005a, b). Another important well-known effect of Wolbachia infections is the inducement of parthenogenesis, see e.g. Engelstädter et al. (2004), Stouthamer (1997. All of these contribute to the fact that the mathematical modelling of Wolbachia infection dynamics is both interesting and challenging. In recent decades a substantial number of mathematical modelling approaches have been applied to model different types of Wolbachia infections in a variety of arthropod species. Perhaps most frequently researchers have been focusing on the development of mathematical models for Wolbachia infections in mosquito species. Many of the earlier models took the form of discrete time matrix models, written for population frequencies, see e.g. Turelli (1994), Vautrin (2007), and the references therein. Using frequency-type models a number of researchers investigated the possibility of coexistence of multiple Wolbachia strains, each of which exhibits different types of the reproductive mechanisms mentioned earlier, see e.g. Engelstädter et al. (2004), Farkas and Hinow (2010), Keeling et al. (2003), Vautrin (2007). Among others, Wolbachia strains have been investigated as a potential biological control tool to eradicate mosquito borne diseases. Originally the focus has been on Wolbachia strains that induce life-shortening of their hosts. This is because for many vector borne diseases only older mosquitoes are of interest from the point of view of disease transmission. Therefore the use of (discrete) age-structured population models has become increasingly prevalent, see e.g. Rasgon and Scott (2004) and the references therein. Fairly recently, in McMeniman (2009) the results of laboratory experiments were reported envisaging a successful introduction of a life-shortening Wolbachia strain in the mosquito species Aedes aegypti. In McMeniman (2009) three key fac-tors, namely, strong CI, low fitness cost and high maternal transmission rate, were identified as drivers of a successful introduction of the new Wolbachia strain into an Aedes population. To this end researchers have developed and analysed continuous age-structured population models for Wolbachia infection dynamics, which take the form of partial differential equations, see Farkas and Hinow (2010); which can often be recast as delay equations, see e.g. Hancock et al. (2011a, b).
In recent years there have been substantial modelling efforts to theoretically investigate the potential of biological control tools for limiting the impact of mosquito borne diseases. It is now widely recognised that biological control represents a viable alternative to established methods such as the use of insecticides and bed nets. Among others, the sterile insect technique has been investigated in the recent papers (Dufourd and Dumont 2012;Li 2011;Li and Yuan 2015). More recently, it was reported that particular strains of Wolbachia (completely or almost completely) block dengue virus replication inside the mosquito hosts, see for example (Blagrove 2012;Hoffmann 2011;Walker 2011). To this end Hughes and Britton (2013) developed a mathematical model for Wolbachia infection as a potential control tool for dengue fever. Their work suggests that Wolbachia may be effective as such a control measure in areas where the basic reproduction number R 0 is not too large. These recent results underpin the possibility that Wolbachia may be a promising candidate for biocontrol of mosquito borne diseases, in general. Besides dengue, West Nile virus (WNv) is another wellknown mosquito borne disease of current interest. WNv infection cycles between mosquitoes (especially Culex species) and a number of species, particularly birds. Some infected birds develop high levels of virus in their bloodstream and mosquitoes can become infected by biting these infectious birds. After about a week, infected mosquitoes can transmit the virus to susceptible birds. Mosquitoes infected with West Nile virus also bite and infect people, horses, and other mammals. However, humans, horses, and other mammals are 'dead end' hosts. This virus was first isolated in the West Nile region of Uganda, and since then has spread rapidly, for example in North America during the past 12 years. Since there is no vaccine available, the emphasis has been mainly on controlling the vector mosquito species. Some recent experiments, see Hussain (2013), have confirmed that replication of the virus in orally fed mosquitoes was largely inhibited in the wMelPop strain of Wolbachia. Interestingly, in a recent paper, Dodson et al. (2014) demonstrated in laboratory experiments that the wAlbB Wolbachia strain in fact enhances WNv infection rates in the mosquito species Culex transalis. However, in Dodson et al. (2014) the Wolbachia was not a stable maternally inherited infection, but rather they infected transiently somatic mosquito tissues, and hence the wAlbB infection did not induce significant immune response in the mosquitoes. This is probably key to their findings. Here we will focus on modelling a maternally inherited Wolbachia infection in a population model, which hypothesizes a large number of successive generations. Nevertheless, the findings in Dodson et al. (2014) underpin the importance of Wolbachia research in general and highlight the importance of contrasting the findings of new theoretical, laboratory and field investigations.
In this work we introduce sex-structured models for Wolbachia infection dynamics in a mosquito population. This will allow us to incorporate and study the well-known effects of CI and MK of particular Wolbachia infections, simultaneously. First we will treat a model which only involves the mosquito population itself. Then we will use this model as a basis for a much more complex scenario incorporating WNv dynamics in a Wolbachia infected mosquito population. The full WNv model will naturally include the bird population, too.

Model derivation
We start by introducing a model for a Wolbachia infection in a sex-structured mosquito population, incorporating sex-structure using a well established approach originally due to Kendall (1949). More recent papers of Hadeler (2012) and Hadeler et al. (1988) derive and discuss sex-structured pair formation models in depth. We only model (explicitly) the adult population of mosquitoes. Our model allows us to take into account the well-known effects of cytoplasmic incompatibility (CI), incomplete maternal transmission, fertility cost of the Wolbachia infection to reproductive output, and male killing (MK), at the same time. We note that it was shown in Engelstädter et al. (2004) that a stable coexistence of MK and CI inducing Wolbachia strains is possible, in principle. Introduction of male killing Wolbachia strains in vector populations may have a significant effect on the disease dynamics, as typically only female mosquitoes are transmitting the disease. Also note that according to Walker (2011), those Wolbachia strains which cause greater disruption, as in the case of dengue transmission, confer greater fitness costs to the mosquitoes. This may well be the case for West Nile virus, hence we account for the reduced reproductive output in our model.
We deduce our starting model from basic principles. In particular, first we deduce mating rules arising at the individual level. Starting with the adult population of size N , we construct a random mating graph. This is a bipartite graph, not necessarily complete, in which each vertex has degree at most one. The vertices represent male and female individuals and edges represent realized matings. Let us denote by M, M w , F, F w the numbers of un/infected males and females, respectively. For every adult mating pair, offspring is created according to the following rules. Below, m, f and m w , f w denote uninfected/infected male/female individual, respectively. The parameter β models the reduced reproductive output of Wolbachia infected females, τ measures maternal transmission in the sense that it is the probability that a Wolbachia infected mother passes on the infection to its offspring, q measures CI in the sense that when a Wolbachia infected male mates with an uninfected female, q is the probability that there is no viable offspring. Finally, γ measures MK in the sense that it is the probability that a Wolbachia infected male larva dies during its development. A complete list of parameter values will be given later on. With this notation the mating rules are described below.
(1) m × f : create one pair of the same type (m, f ).
(2) m × f w : with probability β, create no offspring. This reflects the fecundity reduction due to the Wolbachia infection. In the complementary case, with probability (4) m w × f : with probability q, create no offspring. This is the effect of cytoplasmic incompatibility (CI). With probability 1 − q, create a new pair (m, f ).
Notice that, in contrast to Farkas and Hinow (2010), the sex ratio at birth will not be 1 : 1, it is distorted due to male killing. Also we allow different mortality rates for males and females, in general. Therefore, even in the case when there is no male killing, the sex ratio would be distorted, in general. Also, we assume that any offspring resulting from CI crossing is uninfected. We apply the mating rules described above to construct the birth function in our model. If the population sizes in the four compartments are denoted by M, M w , F, F w , respectively, then the total number of possible matings is (M + M w )(F + F w ). The total number of matings for example between uninfected males and infected females is M F w . Hence the probability that a given mating of type m × f w takes place is (F+F w ) . To compute the total number of matings per unit time we follow the harmonic mean birth function approach from Keyfitz (1972). Accordingly, the total number of matings is proportional to (2.1) Hence the birth rate of offspring arising for example from the mating between an uninfected male and an infected female is proportional to We also naturally assume that there is competition between female individuals for finding an appropriate water reservoir to lay eggs. This is taken into account via a function λ(F total ) which we assume (at least in the first instance) to be a monotonically decreasing function of the total number of females F total , to allow for this competition for nesting places. Though λ(F total ) is decreasing, it may approach a positive limit as F total → ∞. This is to allow for the fact that gravid females that cannot find a place to lay their eggs may destroy eggs previously laid by others, and lay theirs instead. Thus the overall egg-laying rate should approach a positive limit as F total → ∞, and therefore we assume that λ(F total ) is a decreasing function such that λ(∞) > 0. Based on the individual mating rules explained earlier, our model reads as follows. (2.5) A complete list of the variables, parameters and coefficient functions appearing in model (2.2)-(2.5) is given below.
• M: number of uninfected male mosquitoes.
• F: number of uninfected female mosquitoes.
• M w : number of Wolbachia infected male mosquitoes.
• F w : number of Wolbachia infected female mosquitoes.
• M total = M + M w , total number of male mosquitoes.
• F total = F + F w , total number of female mosquitoes.
• N = M total + F total , total number of mosquitoes.
• β: reduction in reproductive output of Wolbachia infected females.
• λ(F total ): average egg laying rate, which depends on the total number of female mosquitoes. • μ m : per-capita mortality rate of uninfected male mosquitoes.
• μ f : per-capita mortality rate of uninfected female mosquitoes.
• μ mw : per-capita mortality rate of Wolbachia infected male mosquitoes.
• μ f w : per-capita mortality rate of Wolbachia infected female mosquitoes.
Model (2.2)-(2.5) is our starting point for a study of the Wolbachia infection dynamics in a sex-structured mosquito population. Later, we will expand this model by introducing WNv infection.

Positivity and boundedness
First we begin by establishing positivity and boundedness of solutions of model (2.2)-(2.5).

Proposition 2.1 Assume that λ is a monotone decreasing function such that
(2.6) hold. Then, the variables (M, F, M w , F w ) satisfying equations (2.2)-(2.5) remain non-negative if they are non-negative initially, and they remain bounded for all times.
Since F total remains bounded it follows that M total is bounded as well, because adding (2.2) and (2.4) we have where B is an upper bound for λ(F total (t))F total (t), which exists since F total (t) is bounded and therefore so is λ(F total (t))F total (t). From the differential inequality (2.8), we can conclude that M total is bounded, too.

Boundary equilibria and their stability
It is straightforward to see that model (2.2)-(2.5) has only one non-trivial Wolbachia free boundary equilibrium and under the assumptions that λ(0) > μ f + μ m , λ is a decreasing non-negative function, and λ(F total ) → λ min (with λ min sufficiently small) as F total → ∞.
Note that there is no Wolbachia infected boundary equilibrium unless τ = 1, a case that we shall treat separately later.
Theorem 2.1 Suppose that λ is a monotone decreasing non-negative function such that λ(F total ) → λ min as F total → ∞, with λ min sufficiently small, λ(0) > μ f + μ m and Then, the Wolbachia free boundary equilibrium E * = (M * , F * , 0, 0) of model (2.2)-(2.5) is locally asymptotically stable. Proof Linearisation of system (2.2)-(2.5) at the equilibrium E * yields the following partially decoupled systems. The first system below is just the linearisation of equations (2.4)-(2.5) at the steady state, which we shall use to show that (2.13) From the second equation of (2.12), it is clear that if follows from the first equation of (2.12). Since M * and F * are given by (2.9) and (2.10), inequality (2.14) is equivalent to assumption (2.11). It remains to prove the local stability of (M, F) = (M * , F * ) as a solution of system (2.13). The Jacobian matrix of system (2.13) evaluated at (M * , F * ) is given by The eigenvalues of J (M * , F * ) satisfy the characteristic equation If τ = 1, i.e. we have complete maternal transmission of Wolbachia, then a boundary equilibrium of the form (0, 0, M * w , F * w ) may exist. The components of such an equilibrium solution must satisfy Next we study the linear stability of such equilibrium, showing that it is linearly stable under condition (2.16) below. Inequality (2.16) does not depend on γ , the male killing rate, but the steady state components M * w and F * w do depend on γ in the manner expected (for example, M * w = 0 when γ = 1). Although Theorem 2.2 only apples if τ = 1, we will be interested later on in the case when τ is just slightly less than 1. Then, maternal transmission is imperfect and Wolbachia infected females produce small numbers of uninfected offspring. We anticipate that as τ decreases from 1 to a value just less than 1, the equilibrium (0, 0, M * w , F * w ) shifts to another nearby position with small numbers of Wolbachia uninfected individuals and large numbers of infected ones; with no change of stability for τ close enough to 1. The existence and stability of such an equilibrium will be important later on when we introduce West Nile virus (WNv) disease dynamics because, at a WNv-free equilibrium with large numbers of Wolbachia infected mosquitoes, the basic reproduction number R 0 for WNv is likely to be less than 1. The implication is that Wolbachia infection in mosquitoes has the potential to control WNv infection. It does so by disrupting WNv virus replication causing WNv infected mosquitoes effectively to remain permanently (or for a very long time) in the latent stage of WNv.
Theorem 2.2 Assume that τ = 1, λ is monotone decreasing with lim F total →∞ λ(F total ) = λ min , (λ min sufficiently small) and The proof is similar to that of Theorem 2.1. The linearisation around (0, 0, M * w , F * w ) yields a system of linear equations for (M, F) that (when τ = 1) are decoupled from the rest of the system. Moreover, it may be shown that holds, which becomes inequality (2.16) when the equilibrium equations (2.15) are invoked. Then, the F w and M w equations are considered, in the case when τ = 1 and F = M ≡ 0. Tedious computations yield that the linearisation of that system around the steady state (M * w , F * w ) has the Jacobian matrix equal to and it may be further shown that its eigenvalues both have negative real parts. Thus we conclude that the steady state (0, 0, M * w , F * w ) is locally asymptotically stable.

Theorem 2.3 Suppose that λ is monotone decreasing and that
(2.17) we are reduced to system (2.13) and we now show that (M(t), F(t)) → (0, 0) as t → ∞, though this result is local, i.e. for small introductions of F and M. Note that (M, F) = (0, 0) is not an equilibrium of (2.13), due to the singularity, but we can remove that singularity by introducing the new variable ξ = F/M. In terms of the variables ξ and M, system (2.13) becomes , (2.20) We now show that (ξ, M) = (ξ * , 0) is a locally stable steady state of (2.20), where provided ξ * > 0. The latter is not automatic. However, note that, from the first equation of (2.13), By similar reasoning we arrive at the same conclusion if λ(0) < μ f . Therefore, we may assume henceforth that λ(0) ≥ max(μ f , μ m ) and, under these circumstances, ξ * > 0. The linearisation of the second equation of (2.20) near the equilibrium (ξ, M) = (ξ * , 0) reads and therefore, since λ(0) < μ m + μ f , we have M(t) → 0. In this limit the ξ equation becomes and to show that ξ * is locally stable as a solution of this equation, it suffices to show that F (ξ * ) < 0. But, after some algebra, we have To show that F (ξ * ) < 0 holds, it suffices to show that λ(0) > |μ m − μ f |, i.e. that both λ(0) > μ m − μ f and λ(0) > μ f − μ m hold. But this follows from the fact that we are now restricting to the case when λ(0) ≥ max(μ f , μ m ). Therefore, the proof of the theorem is complete.

Existence of strictly positive steady states
In this section we examine the possible existence of coexistence steady states (M, F, M w , F w ) of model (2.2)-(2.5), i.e. steady states in which each component is strictly positive. It turns out that in some parameter regimes multiple coexistence steady states may exist while, in others, there is just one or none at all. An understanding of these properties helps us to understand how one might exploit Wolbachia infection in mosquitoes to effectively control WNv. In this section we simplify by assuming that γ = 0, i.e. that there is no male killing. At the steady state, dividing (2.2) by (2.3), and (2.4) by (2.5), we obtain From (2.3) and (2.5) we then obtain Note that if at the strictly positive steady state, we have κ 1 (λ * ) = 0, then this necessarily implies that κ 2 (λ * ) = 0. This is only possible if holds. This case is excluded from Theorem 2.4 but is treated in the next subsection. If (2.28) does not hold then, using (2.27), from (2.23) we obtain (2.29) The right hand side of (2.29) is, in general, a cubic polynomial in λ * . If there exists a positive root λ 1 * , then since λ is a strictly monotone function, a corresponding unique F 1 + F 1 w value may be found. From (2.23) we may then determine a unique solution (F 1 , F 1 w ). Further analytic progress is possible in certain particular cases of interest, which we now investigate. The first concerns the case when τ = 1, or when τ is very close to 1, meaning that maternal transmission of Wolbachia is complete or nearly complete. This is in fact the biologically relevant case for a number of CI inducing Wolbachia strains in mosquito species treated in the literature, see e.g. Engelstädter et al. (2004). This leads us to expect the existence of a steady state of (2.2)-(2.5) with large numbers of Wolbachia infected mosquitoes and few, or no, uninfected ones. The stability of such a steady state still depends on the other parameter values, and it will be stable if there is a high probability of mating between infected males and uninfected females resulting in no offspring (the effect of CI), i.e. q is close to 1; see also inequality (2.16) for the case τ = 1. The existence of a stable steady state of (2.2)-(2.5) with the above mentioned properties is important because the low number of Wolbachia uninfected females implies that the quantity F * s , featuring in the first term of the parameter R 0 defined later in (3.47), is small. The likelihood of R 0 being less than 1 (the condition for WNv-eradication) depends mostly on that first term involving F * s , since the second term in (3.47) involves a small parameter ε and is automatically small. For these reasons, we are interested in stable steady states of model (2.2)-(2.5) of the form (M * , F * , M * w , F * w ), with M * and F * small compared to M * w and F * w (and, ideally, M * = F * = 0). Therefore, referring to Theorems 2.1 and 2.2, and restricting for now to the case τ = 1, we ideally would like the boundary equilibrium (M * , F * , 0, 0) of Theorem 2.1 to be unstable, and the boundary equilibrium (0, 0, M * w , F * w ) of Theorem 2.2 to be linearly stable. The conditions for this, when τ = 1, are that the inequalities μ f (1 − β) > μ f w and (1 − q)μ f w < (1 − β)μ f should hold simultaneously, but note that the second of these follows from the first. For τ = 1, guided by elementary competition theory, instability of one boundary equilibrium and stability of the other suggests that there will be no coexistence equilibrium, and this is what we prove in Theorem 2.4 below. If τ is decreased from 1 to a value slightly below 1, there is no longer a boundary equilibrium with only Wolbachia infected mosquitoes present. What happens is that the equilibrium (0, 0, M * w , F * w ), which exists when τ = 1, moves to another nearby point in R 4 + , so that Wolbachia infected mosquitoes now coexist with uninfected ones, the former being dominant. Very importantly, this will be the only coexistence steady state if τ is sufficiently close to 1, and it is a desirable steady state for WNv eradication because particular Wolbachia strains can significantly reduce WNv virus replication in mosquitoes, see Hussain (2013).

Theorem 2.4 Suppose that
and that λ is a strictly positive decreasing function with λ(∞) = λ min (and λ min sufficiently small). Then system (2.2)-(2.5) has no coexistence equilibrium with M * , F * , M * w , F * w > 0. If the foregoing hypotheses hold, except that τ is slightly less than 1, then system (2.2)-(2.5) has precisely one coexistence equilibrium in which Wolbachia uninfected mosquitoes exist in very small numbers relative to Wolbachia infected ones.
Proof Since we assume τ = 1, the form of (2.29) simplifies and in fact we may cancel κ 1 (λ * ) since we seek equilibria in which M * , F * , M * w , F * w > 0. After some further algebra, we find that there is just one value for λ * , given by (2.30) If λ * ≤ 0 then the equation λ * = λ(F + F w ) cannot be solved for F + F w , so we may restrict to the case that λ * > 0. Recalling that κ 1 and κ 2 are defined by (2.27), we find, with λ * given by (2.30) and with τ = 1, that and and, since we assume μ m μ f w = μ f μ mw , it follows that κ 1 (λ * )κ 2 (λ * ) < 0. This makes it impossible to find F > 0 and F w > 0 satisfying (2.26), and so there is no coexistence equilibrium. Next we prove the second assertion of the theorem. Let τ = 1 − ε. Since we expect the equilibrium (0, 0, M * w , F * w ), which exists when τ = 1, to move to another nearby point when τ = 1 − , for sufficiently small, we seek an equilibrium of (2.2)-(2.5) of the form Coefficients of ε yield and (2.34) reads (2.35) Note that F (1) > 0 because of the assumption μ f (1 − β) > μ f w , hence from (2.33) we conclude that M (1) > 0 holds, and (2.36) In conclusion, if τ is reduced from 1 to the value 1 − ε then the equilibrium , with M (1) and F (1) given by (2.36) and (2.35), respectively.
Note that from the proof of Theorem 2.4 we can see that at the equilibrium the male/female ratio for Wolbachia uninfected mosquitoes is given approximately by Theorem 2.4 excludes the case when μ m μ f w = μ f μ mw , which we treat now. In particular, we assume that γ = 0, and In this situation, we have M = μF, and M w = μF w . From equations (2.3) and (2.5) we obtain (2.38) From equation (2.38) we find that (2.41) For the existence of the positive steady state we need to guarantee that the quadratic equation above has (at least one) positive (real) solution, and that the solution is less than c.
From equation (2.41) it is clear that in case of complete maternal transmission, i.e. for τ = 1 there cannot be more than one coexistence steady state. In this case, from equation (2.41) we obtain holds, then we have 0 < F < c, and if λ(0) > (1+μ)μ f w μ(1−β)τ also holds, then a unique strictly positive steady state exists.
Circumstances under which (2.43) is likely to hold include that q is sufficiently close to 1 and, at the same time, β is sufficiently close to 1 or μ f μ f w is less than 1. The biological interpretation is clear: for the existence of a coexistence steady state, the fertility cost (or mortality increase) due to Wolbachia infection should be sufficiently large.
It is clear from (2.41) that we can never have more than two coexistence steady states. On the other hand it is interesting to show that for some realistic parameter values it is possible to have two coexistence steady states.
To this end we consider the case q = 1, μ f μ f w = 1, and we assume that τ = 1, β = 1. In this case, from (2.41), we have That is, for 0 < F 1/2 < c to hold, we need to assume that hold simultaneously. It is easy to verify that this can be achieved for any c, for example with τ = 0.99, β = 0.5. Note that the condition μ f μ f w = 1 can be relaxed, too. Also, by continuity arguments, it follows that for parameter values close enough we still have two coexistence steady states. We summarise our findings in the following proposition.

Proposition 2.2
In the case when γ = 0, q = 1 and μ f μ m = μ f w μ mw , there exists a set of values for the remaining parameters such that system (2.2)-(2.5) admits two coexistence steady states.
In summary, we have shown that our Wolbachia model (2.2)-(2.5) may exhibit all of the three qualitatively different possible scenarios, i.e. when there are 0, 1 or 2 coexistence steady states.

Model incorporating West Nile virus (WNv)
There has been considerable recent interest in West Nile virus (WNv), with the great majority of mathematical papers on the topic having appeared in the last 15 years. Numerous types of models have appeared, some including spatial effects and others giving consideration to issues such as age-structure in hosts, optimal control or backward bifurcation. See, for example, Blayneh et al. (2010), Bowman et al. (2005), Gourley et al. (2006/07), Lewis et al. (2006) and Wonham and Lewis (2008).
We introduce the following model as an extension of model (2.2)-(2.5) to include WNv dynamics. Our model for WNv dynamics has similarities to that in Bergsman et al. (2016) with one resident bird population. Hence, as in Bergsman et al. (2016) and in some references therein, we compartmentalise the vector and bird population into SEI and SEIR classes, respectively. Taking into account our earlier model the complete WNv-Wolbachia mosquito-bird population model takes the following form.
Some of the parameters of system (3.46) are defined after (2.2)-(2.5); the rest are defined as follows: • α f : biting rate of female Wolbachia uninfected mosquitoes; • α f w : biting rate of female Wolbachia infected mosquitoes; • p b f : transmission probability of WNv from infectious birds to WNv-susceptible female mosquitoes; • p f b : transmission probability of WNv from WNv-infectious female mosquitoes to susceptible birds; • ν f : per-capita transition rate of WNv-exposed female Wolbachia uninfected mosquitoes to the infectious stage of WNv; • ε ∈ [0, 1]: small parameter modelling increased time that Wolbachia infected mosquitoes spend in the latent stage of WNv, due to the tendency of Wolbachia infection to hamper the replication of WNv in mosquitoes; • ν b : per-capita transition rate of WNv-exposed birds to the infectious stage of WNv; • ν i : per-capita rate at which infectious birds recover; • μ b : per-capita natural death rate for birds; • μ bi : per-capita WNv-induced death rate for infectious birds.
In this section, it must be emphasized that we are considering two different kinds of infection. Mosquitoes may be infected by either Wolbachia or WNv, or both. WNv infection is assumed possible only for female mosquitoes (since it is females that bite) and is modelled using an SEI (susceptible-exposed-infectious) approach with subscripts s, e and i. The variables F s , F e and F i denote the numbers of Wolbachia uninfected mosquitoes that have susceptible, exposed and infectious status with respect to WNv. A subscript w indicates Wolbachia infection, so that F ws , F we and F wi denote the numbers of Wolbachia infected mosquitoes that have susceptible, exposed and infectious status with respect to WNv. The variables M and M w are the numbers of Wolbachia uninfected and Wolbachia infected male mosquitoes, none of which have WNv. Birds are only susceptible to WNv and their numbers are given by the variables B s , B e , B i and B r denoting susceptible, exposed, infectious and recovered birds. Many of the terms of system (3.46) are also present in (2.2)-(2.5) without change, here we just discuss the extra terms that model the addition of WNv dynamics. WNvsusceptible mosquitoes, whether Wolbachia infected or not, acquire WNv infection by biting infectious birds; this is modelled via the last term in the first and fourth equations of (3.46) using the idea of mass action normalised by total host density, the biting rates (the α parameters defined above) having been separated out, rather than being absorbed into the transmission coefficients p b f and p f b as is often customary. Having acquired WNv infection from a bird, a mosquito enters the latent phase of WNv and is classed as an exposed mosquito. Exposed mosquitoes become WNv-infectious at rates ν f F e and εν f F we for Wolbachia uninfected and Wolbachia infected mosquitoes, respectively. In the latter, the presence of ε ∈ [0, 1] models the tendency of Wolbachia infected mosquitoes that have contracted WNv infection to spend a greater amount of time in the latent stage of WNv, since Wolbachia infection tends to block WNv replication making it less likely that such a mosquito would ever become WNv-infectious. Of course, we are at liberty to take ε very small indeed, with the implication that the Wolbachia infected mosquito spends so long in the latent stage of WNv that it probably dies in that stage. This is our approach to modelling the blocking of WNv replication by Wolbachia.
Birds acquire WNv from bites by WNv-infectious mosquitoes, which may or may not have Wolbachia infection as well. Thus there are two infection rates for birds, these can be found in the right hand side of the seventh equation of (3.46), and also in the eighth equation since birds initially enter the exposed stage of WNv. This has a mean duration of 1/ν b for birds, after which they become WNv-infectious. Birds may recover from WNv, at a per-capita rate ν i . Note that, for birds, death due to WNv is modelled using a separate parameter μ bi to distinguish from natural death, accounted for by μ b . The function (B total ) is the birth rate function for birds.
The approach we use here to model the latency stage of WNv (in either birds or mosquitoes) is not the only possible approach. Our approach permits individuals to spend different amounts of time in the latency stage, and we may only speak of the mean time spent in that stage. There are other approaches in which all individuals of a particular status (for example, all Wolbachia uninfected mosquitoes) spend the same amount of time in the latent stage of WNv. The time could be different for Wolbachia infected mosquitoes. These approaches result in models with time delays.

Local stability of the WNv-free equilibria
Equilibria of system (3.46) may exist in which WNv is absent. Such WNv-free equilibria include the equilibrium (M * , F * , 0, 0) considered in Theorem 2.1, in which both WNv and Wolbachia are absent, and equilibria in which WNv is absent but Wolbachia are present. We show that multiple WNv-free equilibria may coexist that have both Wolbachia uninfected and Wolbachia infected mosquitoes, we present a necessary and sufficient condition for any particular WNv-free equilibrium to be locally stable, and we show that the most likely scenario for eradication of WNv is to have large numbers of Wolbachia infected mosquitoes, with solutions of system (3.46) evolving to a WNv-free equilibrium that has large numbers of Wolbachia infected mosquitoes and relatively few uninfected ones. Theorem 3.1 applies to any WNv-free equilibrium, of which there may be several. Of course, we may have R 0 < 1 at one WNv-free equilibrium and R 0 > 1 at another. It depends on the values of F * s , F * ws and B * s for the particular WNv-free equilibrium under consideration. For clarity of exposition, we include as a hypothesis that the equilibrium be stable to the subset of perturbations in which WNv is absent (i.e. stable as a solution of the subsystem (2.2)-(2.5)), rather than including explicit conditions for stability of an equilibrium as a solution of that subsystem. The latter stability problem is a tedious one in its own right and is under consideration elsewhere in this paper. Theorem 3.1, in the form presented below, highlights clearly the particular role played by R 0 .
Then, if R 0 < 1, the WNv-free equilibrium under consideration is locally stable as a solution of the full system (3.46), if it is stable to perturbations in which the exposed and infectious variables remain zero.
Proof At any WNv-free equilibrium, the linearisation of system (3.46) decouples to some extent making it sufficient to show that, when R 0 < 1, each component of the solution of the following system: (3.48) tends to zero. It is taken as a hypothesis that the susceptible variables then approach their respective steady state values. Note that system (3.48) has a structure that allows the application of Theorem 5.5.1 in Smith (1995), making it possible to restrict attention to the real roots of the characteristic equation associated with (3.48). That characteristic equation, corresponding to trial solutions with temporal dependence exp(λt), is most easily analysed when written in the form (3.49) As functions of the real variable λ, the right hand side of (3.49) is decreasing, at least for λ ≥ 0, while the left hand side is a quadratic with two real negative roots. A simple graphical argument shows that if the left hand side exceeds the right hand side when λ = 0 (i.e., if R 0 < 1, with R 0 defined by (3.47)), then any real roots of the characteristic equation are negative which, since only the real roots need to be considered, implies that each component of the solution of (3.48) approaches zero as t → ∞. The proof of the theorem is now complete.

Numerical simulations
In the simulations shown in Figs. 1-8, we set λ(F total ) := re −F total /k , where r is the maximum per-capita mosquito egg-laying rate, and k measures intra-specific competition among female mosquitoes. It should be noted that our model assumes that the mosquito population persists annually, as for example in the tropical climates of South-East Asia.  Table 1. If τ > 0.9 the basic reproduction number R 0 < 1, and WNv will die out  Table 1. If q > 0.8 the basic reproduction number R 0 < 1, and WNv will die out

Conclusion
In this paper we have derived a detailed sex-structured model for a mosquito population infected with Wolbachia . The model captures many of the well-known key effects of Wolbachia infection, including cytoplasmic incompatibility, male killing, reduction in reproductive output and incomplete maternal transmission of the Wolbachia infection. Our analysis shows that the mosquito population can stabilise at  Table 1. If γ < 0.65 the basic reproduction number R 0 < 1, and WNv will die out a Wolbachia free equilibrium under certain circumstances, which include situations when inequality (2.11) holds. Such circumstances include, for example, if Wolbachia infection significantly reduces reproductive output, and/or Wolbachia infection significantly lowers female life expectancy. We also showed that if τ = 1, i.e. maternal transmission of Wolbachia is complete, then the mosquito population can stabilise at an equilibrium in which all mosquitoes are infected with Wolbachia. This happens in the case of sufficiently high cytoplasmic incompatibility. In the case of τ close to 1 we have shown that Wolbachia infected mosquitoes can coexist with small numbers of uninfected mosquitoes. We have also shown that under some additional assumptions our model has multiple coexistence steady states.
We extended the sex-structured mosquito population model (2.2)-(2.5) to include West Nile virus, which is spread by birds and mosquitoes, treating WNv as an SEI infection for mosquitoes, and as an SEIR infection for birds. We were motivated by results recently reported in Hussain (2013), which suggest that a particular strain of Wolbachia substantially reduces WNv replication in the mosquito species Aedes aegypti. We modelled this crucial phenomenon by incorporating a small parameter ε, the reciprocal of which is proportional to the time spent in the WNv exposed class for Wolbachia infected mosquitoes. This enabled us to assess the potential of Wolbachia infection to eradicate WNv via its effect on WNv replication in Wolbachia infected mosquitoes. Notably the expression we obtained for the basic reproduction number R 0 suggests that WNv will be eradicated if at the steady state the overwhelming majority of mosquitoes are infected with Wolbachia, and the Wolbachia infection substantially reduces WNv replication in mosquitoes. The first of these hypotheses is in fact shown to hold for a number of Wolbachia strains and mosquito species, see e.g. Engelstädter and Telschow (2009).