MODELING SPATIAL SPREAD OF WEST NILE VIRUS AND IMPACT OF DIRECTIONAL DISPERSAL OF BIRDS

A patchy model for the spatial spread of West Nile virus is formulated and analyzed. The basic reproduction number is calculated and com- pared for different long-range dispersal patterns of birds, and simulations are carried out to demonstrate discontinuous or jump spatial spread of the virus when the birds' long-range dispersal dominates the nearest neighborhood interaction and diffusion of mosquitoes and birds.

1. Introduction.Although West Nile virus (WNv) was isolated in the West Nile district of Uganda in 1937 [18], and WNv in the Eastern Hemisphere has been maintained in an enzootic cycle involving culicine mosquitoes and birds [9,10], WNv activities in North America were not recorded until August 1999, in the borough of Queens, New York [3,14].Despite this long delay, the virus has spread rapidly within the subsequent five years and evidence seems to indicate that WNv has become a permanent fixture of the North American medical landscape [16].
There has been increasing interests in modelling the transmission dynamics of WNv, and several models have been developed to address various aspects of the disease transmission; see, for example, [4,12,13,19,20,23].Here, we focus on the spatial spread patterns and speed of transmission, which are obviously of paramount importance in terms of prevention and control.One of the goals of this research is to understand the jump or discontinuous spatial spread patterns in the establishment phase of WNv, as shown in the 2000-03 Health Canada map of dead birds submitted for WNv diagnosis by health region (Figure 1).This discontinuous spatial spread seems to be the consequence of the combination of the local interaction and spatial diffusion of birds and mosquitoes and long-range dispersal of birds, and this motivates the use of patchy models instead of the reaction-diffusion model.Our second goal is to see how the interaction of the ecology of birds and mosquitoes, the epidemiology of bird-mosquito cycles, and the diffusion and immigration patterns of birds affect the long-term and transient transmission of the diseases within the whole region consisting of multiple patches.We will do so by calculating the basic reproduction number of the region as a function of the basic reproduction number of each patch, the spatial dispersal rates and patterns of birds, and the spatial scale of the birds' flying range in comparison with the mosquitoes' flying range.
In section 2, we formulate a patchy model, based on the spatially homogeneous single season ordinary differential equations model in [23] for the local interaction of birds and mosquitoes within a patch, and we use linear dispersal among patches.We use the the average distance a female mosquito can traverse during its lifetime as a measure for the partition of the region under consideration, and hence the number of patches that a bird can fly during the life span of a female mosquito becomes a natural and important parameter in the modeling and for the analysis and simulations.In section 3, we calculate the reproduction number R 0 and study how this number is affected by the direction-selective dispersal pattern of birds.Our analysis in section 3 and the simulations carried out in section 4 show that this direction-selective dispersal of birds decreases the reproduction number and slows the spatial spread of WNv.Our simulations seem to indicate that the jump spatial spread of WNv arises if the birds' long-range dispersal dominates the nearest neighborhood interaction and diffusion of mosquitoes and birds.
The current work focuses on the one-dimensional patch model, which can only be regarded as a theoretical approximation of the West Nile virus landscape in Canada by connecting Ontario with British Columbia with a straight line and ignoring the transmission heterogeneities along other directions.Better understanding of the West Nile virus spread in the Canadian medical landscape can be achieved only by extending this work to a two-dimensional model and by incorporating more spatiotemporal heterogeneities.
2. The Model and Disease-free Equilibrium.To formulate the patchy model for the spatial spread of the West Nile virus, we assume that there are N patches under consideration where, depending on the purpose of modeling, availability of data, and implementation of surveillance, control, and prevention measures, the partition of the whole region into non-overlapping patches changes.In the modeling and simulations below, we use the average distance a female mosquito can fly during its lifetime as a measuring unit for the partition.Therefore, if we assume the region is one-dimension and we use 1, . . ., N to denote the corresponding patches, then mosquitoes belonging to the ith patch can fly only to their nearest neighboring patches i − 1 and i + 1, while birds belonging to the i-th patch fly to their mth neighbor patches i − m, . . ., i − 1, i + 1, . . ., i + m, with m ≥ 1 (see Figure 2).
We also make the following assumptions: (H1) The virus does not have any adverse effect on mosquitoes, and vertical transmission in mosquitoes can be ignored.(H2) Most birds will recover from the virus and become immune to further infection, and new-born birds have no immunity.(H3) Birds and mosquitoes have fixed recruitment rates in each patch.We use the model set up in [4] for the dynamics between birds and mosquitoes within a patch and linear spatial dispersal among patches.We denote the number of individuals of birds and mosquitoes on the ith patch at time t respectively by B Si : the susceptible birds in patch i; B Ii : the infectious birds in patch i; M Si : the susceptible mosquitoes in patch i; M Ii : the infectious mosquitoes in patch i.
Following the flow chart in Figure 3, we have the following model equations (2.1) with 1 ≤ i ≤ N , N being the number of patches.The total number of birds in patch i is N Bi = B Si +B Ii .All parameters are defined in Table 2. Based on the biological fact that the death rate of infected birds is greater than that of susceptible birds, we assume d b2i ≥ d bi for all i.
We assume the dispersion rates of birds depend on the distance from the starting patch to their destination, but these rates may depend on the direction; namely, we have where The dispersal rates of mosquitoes are given by where We now calculate the disease-free equilibrium (DFE).Setting B Ii = 0 and Diffusion rate of mosquitoes from the ith patch to the jth patch with 1 ≤ i ≤ N .For convenience, we rewrite (2.4) in vector form as where • M m is a tridiagonal matrix with Observing that the ith column of the matrix we have, from the Gershgorin circle theorem [1], that all eigenvalues of (2.6) 3. Basic Reproduction Number and Impact of Directional Dispersal.We now calculate the basic reproduction number, denoted by R 0 , that is "the expected number of secondary cases produced, in a completely susceptible population, by one typical infectious individual" [6].If R 0 < 1, then on average an infected individual produces less than one new infected individual over the course of its infectious period, and the infection cannot grow.Conversely, if R 0 > 1, then each infected individual produces, on average, more than one new infection, and the disease can invade and spread in the population.Here, we use the formula in [7] to calculate the reproduction number R 0 .
We use the vector notation to rewrite the equations in which the infectious appear in terms of the difference between f j , the rate of appearance of new infectious in compartment j, and v j , the rate of transfer of individuals into and out of compartment j by all other processes: The corresponding Jacobian matrices, F and V, describing the linearization of this reduced system about the DFE, are given by In this and in what follows, an empty element or block in a matrix means zero (number or matrix), and by the dispersal condition (2.2), and by the dispersal condition of (2.3).
Theorem 3.1.For model (2.1) Proof.Since the sum of the ith column of B is positive, by Gershgorin circle theorem, all eigenvalues of B have positive real parts.At the same time, B has the Z sign pattern.By [2], B is a nonsingular M-matrix.Similarly, M is also a nonsingular M-matrix.According to the definition of V, it is a nonsingular M-matrix.Then V −1 is nonnegative.Since F and V −1 are nonnegative, we have FV −1 is nonnegative.Using the formula in [7], we have R 0 = ρ(FV −1 ).
The basic reproduction number R 0 is well defined, but it is a nontrivial task to find the explicit form of R 0 in the general case.In what follows, we consider two special cases to illustrate the impact of symmetric spatial dispersal of birds.
3.1.Identical Patches and Symmetric Dispersal.Here, we consider the special case where all patches are identical from the aspect of ecology and epidemiology.In other words, we assume that In this case, a straightforward calculation gives the coordinates for the DFE as In the following theorem, we give an explicit expression for the reproduction number R 0 under the above assumptions.
Theorem 3.2.If all patches are identical and the dispersal of birds is symmetric, then the basic reproduction number is given by Before proving this theorem, we need the following lemma.
Lemma 3.3.Let A 1 , A 2 be two positive definite matrices, and γ be the spectrum radius of the matrix 0 A 1 A 2 0 .Then where || • || is the operator norm.
Proof.Let λ be an eigenvalue of 0 A 1 A 2 0 .We show that λ 2 is an eigenvalue We have By the second equation of (3.3), we have Plugging the first equation of (3.3) into (3.4),we get Therefore, λ 2 is an eigenvalue of A 1 A 2 .It follows that where ρ(A 1 A 2 ) is the spectrum radius of A 1 A 2 .On the other hand, we have Therefore, Remark 3.4.In the case where A is a positive definite matrix, ||A|| is the spectrum radius of A.

Now we can assert the following.
Proof of Theorem 3.
So, the vector Denoted by ρ(FV −1 ), the spectrum radius of FV −1 , we have By Lemma 3.3 and Remark 3.4, we obtain Thus, Finally, by the definition of R 0 , we have as desired.
Note that if all patches are isolated from each other (or the whole region is homogeneous), one can calculate the basic reproduction number R 0H to obtain In other words, a region consisting of identical patches coupled by symmetric dispersal of birds has the same reproduction number as if each patch is isolated from the others.This conclusion is not true anymore, however, if the dispersal of birds is not symmetric, as shown in next subsection.
3.2.Nonsymmetric Dispersal of Birds.We now consider the case where the dispersal rates of birds depend on the direction.We shall use the perturbation theory to calculate the basic reproduction number of model (2.1) in the special case of three identical patches.We denote the diffusion rate of birds to the left by D bl and to the right by D br .Other notations remain the same as above.
After some straightforward calculations, we obtain the DFE as and We denote the corresponding F and V matrices by F 3 and V 3 .They are given by To address the impact of the direction-selective dispersal of birds on R 0 , we write D br = D bl + , where is a small positive number.Note that = 0 implies the symmetric dispersal of birds.
Let A = V 3 F −1 3 .As shown above, if diffusion rates of birds are symmetric; that is = 0, (D bl = D br ), the smallest eigenvalue of A =0 = A 0 is Let p denote the eigenvalue of A with a corresponding eigenvector − → v so that where Here, after some straightforward calculations, we have that , and Substituting (3.7) to (3.6) and comparing the coefficients in and 2 terms, we have where I is the identity matrix.Notice that det(A 0 − p 0 I) = 0, we conclude that equation (3.8) in − → v 1 has solutions if and only if the rank of coefficient matrix equals that of the augmented matrix; that is (3.10) Since p 0 is an eigenvalue of A with multiplicity 1 if both D bl and d m are positive, we must have p 1 = 0.
We now proceed to find p 2 .Substituting p 1 = 0 to (3.8), we can solve for − → v 1 .Substituting − → v 1 into (3.9),we can solve (3.9) in the same manner as we did for (3.8).Namely, (3.9) has solutions if and only if (3.11) From (3.11), we get Then, we have Note that if we fix D bl > 0, p 2 is always positive.Therefore the breaking of symmetry in spatial dispersal of birds always decreases R 0 .
4. Discontinuous Spread and Surveillance Data.In this section, we report some numerical simulation results to demonstrate the effect of different dispersal patterns of birds on the spatial spread of WNv, and to illustrate possible discrepancy between surveillance data and the model-based simulations in different time scales.Our focus is on the time when a particular patch has recorded WNv activities, namely, when at least one bird has died of WNv infection.This allows us to compare the simulation results with surveillance data, since Health Canada uses dead birds with WNv as an indicator for determining whether a region has WNv activities.0.00001 estimated * *Since most mosquitoes will stay where they are born, we take a small number for mosquitoes' diffusion rates.
We first discuss ranges of parameters involved.We assume all patches are identical from the aspect of ecology and epidemiology as discussed before, and we fix all parameters as in Table 2.The dispersal rates of birds are a decreasing function of the distance from the origin, but the spatial dispersal may be nonsymmetric in terms of spatial direction (left vs. right, in the case of one-dimensional space).For the sake of simplicity, we assume the dispersal rate of birds g b (k), with k = i − j, from patch i to j, is a piecewise linear function, given below: where m ± i are the furthest patches that a bird can fly during the average life span of female mosquitoes and h 1 measures the diffusivity rate of birds to the left, while h 2 measures the diffusivity rate of birds to the right.The net rate at which a bird flies out of a given patch should be less than 1; therefore, 0 ≤ (h 1 + h 2 )m/2 < 1.Notice that h 1 = h 2 corresponds to the bidirectional dispersal symmetric with respect to the spatial direction, while h 1 = h 2 corresponds to the nonsymmetric spatial direction selective dispersal that seems to be more closer to the ecological reality of birds in Canada within the time scale under consideration.
On average, birds can fly 13.4km per day or 1000km per year [15].During the average life span of 30 days, most female mosquitoes remain within 1.6 km of their breeding site.A few species may range up to 10 miles or more.Thus in the average lifespan of female mosquitoes, the flying range of a bird is about 40 times than that of mosquitoes.Hence, we shall assume m = 40.
The distance from British Columbia Province to South Ontario is about 3000km, and hence we assume, in what follows, the total number of patches is N = 300.In the simulation, the time unit is one day.Since WNv is a seasonal disease, we consider the period from late April to early October to be a total of about 180 days.We assume each patch has 10,000 birds and 200,000 female mosquitoes initially and only the first patch has 5 infected birds and 20 infected mosquitoes.
In Fig. 4, we use the continuous dispersal function (4.1) of birds with m = 40.Figure 4a deals with the symmetric dispersal case with h 1 = h 2 = 0.005, In the symmetric case, we can see the spread speed of WNv is about 1000 km/year, which coincides with the observed spread rate in North American [13].Increasing h 2 slightly while keeping h 1 unchanged will yield nonsymmetric dispersal of birds, but this minor breaking of symmetry has limited impact on the number of infected birds and their spatial spread as shown in Figure 4b.However, if we continue to increase h 2 to h 2 = 0.01, we obtain the graph of Figure 4c.In this case the spread speed of the disease is much faster and the magnitude of outbreaks is higher compared to the cases shown in Figure 4(a, b).Naturally, we notice that the spatial spread is continuous in the sense that there is no patch i escaping from WNv if patch j > i has WNv activities (i.e., had infected birds).This is due to the continuous spread of birds.
In reality, birds may skip some patches during their long-range dispersal.To model this special dispersal pattern, we consider a dispersal function of the following form In other words, the birds in patch i jump to patches i ± 4, . . ., i ± 4J, where J is the integer part of m/4.We shall remark that the model (2.1) is still a continuous model even the dispersal function of birds is not continuous.Figure 5 provides simulation results using the above jump dispersal function and in the case h 1 = h 2 = 0.01.We observe obvious jumps in the transmission of WNv and the disease spread speed is about 1000km/year.In this case, some patches avoid the disease because of the discontinuous dispersal of birds. 5. Discussion.In this paper, we use a mathematical model to understand the spatial spread patterns in the establishment phase of WNv in a region consisting of multiple patches.In the case where patches are identical (from the aspect of ecology and epidemiology) and where the spatial dispersal of birds and mosquitoes are symmetric with respect to all possible spatial directions, we show that the basic reproduction number coincides with that when each patch is isolated from others.This number changes when the spatial dispersal of birds depends on the spatial direction, and our calculation based on a perturbation argument gives an analytic formula to show how this number is changed by the spatial direction selective spatial dispersal.
Figure 1, based on the surveillance date of Health Canada, indicates the jump or discontinuous spread of WNv activities; namely, some patches were free of recorded birds dead from the virus, while their neighboring patches do have recorded dead birds.This phenomena cannot be explained via continuous (partial differential   equations) models, and suggests that the spatial spread of WNv is largely determined by the long-range dispersal of birds.Here we use a special function (4.2) to describe the spatial dispersal of birds, and our simulations illustrate that the discontinuous spatial spread arises very naturally as the consequence of the local interaction between birds and mosquitoes, the short distance diffusion of mosquitoes, and the long-range jump dispersal of birds.
More biological realities need to be incorporated into the model to have a full simulation and comparison with surveillance data.The geographic characters, such as lake areas, or forest and bush areas, which have huge populations of mosquitoes and birds compared to those of patches consisting of mainly cities, should be considered.Also important for model simulation is the dimensionality: more realistic models should be at least two-dimensional.The biological characters of birds must be taken into consideration: different birds have different dispersal patterns and different epidemiological features for WNv.Seasonal factors must also be explored.

Figure 1 .
Figure 1.The maps of dead birds submitted for West Nile virus diagnosis by health region in Canada 2000-03 [5].

Figure 2 .
Figure 2. The mosquitoes can fly to their nearest neighboring patches only, but birds can reach as far as the mth neighboring patches.In other words, partition of the region into patches is based on the flying range of mosquitoes.The pictures of birds and mosquitoes are taken from http://westnilemaps.usgs.gov/.

Figure 3 .
Figure 3. Flow chart of the West Nile virus in a patchy environment.
B b have positive real parts.At the same time, B b has the Z sign pattern.By [2], B b is a nonsingular M-matrix.Then B −1 b is a positive matrix (all elements of this matrix are positive).Similarly, M m is also a nonsingular M-matrix and the inverse matrix of M m is positive.Since − → b and − → m are nonnegative and nonzero vectors, (2.5) has exactly one positive solution, denoted by − → B S = (B * S1 , . . ., B * SN ) T , − − → M S = (M * S1 , . . ., M * SN ) T .

2 .
Under our assumptions, B and M are strictly positive-definite.Now, we claim that λ 1 = d b2 is the smallest eigenvalue of B. In fact, all eigenvalues of B satisfy det|B − λI| = 0. Adding all rows of B − λI to its first row, the new first row becomes (d b2 − λ)[1, . . ., 1].That means, d b2 is the eigenvalue of B. We claim this is also the smallest eigenvalue.This is because, by the Gershgorin circle theorem, we know all eigenvalues of B are greater than or equal to d b2 .Now we look for the corresponding eigenvector of λ 1 .Solving the linear system (B − λ 1 I)ζ = 0, where I is the identity, we obtain ζ = (1, . . ., 1)T .Similarly, we know λ 2 = d m is the smallest eigenvalue of M with the same eigenvector ζ.Noting that F 1 = C mb I and F 2 = C bm M * S B * S I in this case, we have the maximum eigenvalues µ 1 = C mb /λ 1 and µ 2 = C bm M * S /B * S λ 2 of F 1 M −1 and F 2 B −1 , respectively, sharing the same eigenvector ζ.Note that

Figure 4 .
Figure 4.The numbers of infected birds on each patch versus time.To simulate the effects of different birds' diffusion patterns, h 1 , the diffusivity to left, is fixed at 0.005, while h 2 , the diffusivity to right, is varied.

Figure 5 .
Figure 5.The numbers of infected birds on each patch versus time when birds skip some patches.

Table 1 .
Definitions for parameters in the model.
bij Diffusion rate of birds from the ith patch to the jth patch D mij

Table 2 .
Values of parameters in the simulations.