Modeling and Simulations of Mosquito Dispersal . The Case of Aedes albopictus

To prevent epidemics of mosquito-transmitted diseases like Chikungunya in Ŕeunion Island, we develop tools to control its principal vector, Aedes albopictus . Biological control tools, like the Sterile Insect Technique (SIT), are of great interest as an alternative to chemical control tools which are very detrimental to environment. The success of SIT is based on a good knowledge of the biology of the insect, but also on an accurate modeling of insects distribution. We model the mosquito dispersal with a system of coupled parabolic PDEs. Considering vector control scenarii, we show that the environment can have a strong influence on mosquito distribution and in the efficiency of vector control tools. Keywords-parabolic equation; existence; mosquito dispersal; modeling; numerical simulation; splitting algorithm; vector control; Sterile Insect Technique.


I. INTRODUCTION
Chikungunya is an unusual vector-borne disease.First isolated in Tanzania in 1953, it is now geographically distributed in Africa, India and South-East Asia.After a huge epidemic in Réunion Island and in India in 2006, it appeared for the first time in Europe, in Italy, in 2007 (see [1], [2] and references therein).The symptoms that characterize Chikungunya are high fever, headache, persistant joint pain that can last several weeks.One way to reduce the risk of infection for the population is to control the vector populations.The principal vector for Chikungunya, is Aedes albopictus mosquito, commonly called the "Asian tiger" [3].
Standard vector control tools, like adulticide and lavicide, together with mechanical control are useful but cannot always be used for several reasons.Firstly, because Réunion Island is a hot spot of endemicity.Secondly, because mosquito can develop resistance to insecticides.Thirdly, because only approximately 10% of the island can be treated, due to its chaotic landscape.Therefore, it is necessary to consider new sustainable alternatives or additional tools, like the Sterile Insect Technique (SIT).SIT consists in releasing sterilized male mosquitoes that will mate with wild females which won't be able to have offspring.Consequently, this will lead to the decrease of the vector population [4], [5].The success of SIT is based on a good knowledge of the biology and the behavior of the vector, but also on an accurate modeling of its dispersal to optimize the impact of sterile males.The previous published models were temporal models, and didn't take into account the spatial component.But, it is necessary to consider a spatiotemporal model to obtain realistic simulations and to simulate several vector control strategies.
In a previous paper [6], we have considered a dispersal model with only adult females splitted in two compartments: the blood meal searchers and the breeding site searchers.This led to a system of two partial differential equations.Here, we add the aquatic stage, immature females, resting females, wild males, and finally, sterile males.This leads to a system of coupled (nonlinear) partial differential equations.After some theoretical results and the presentation of the compartmental model, we present briefly the numerical methods based on the operator splitting technique [7].Finally, we end the paper with some numerical simulations with and without chemical or biological control.

II. THE MATHEMATICAL MODEL
Aedes albopictus is found in South-East Asia, the Pacific and the Indian Ocean islands, and up North through China and Japan.In the last twenty years, it has invaded several developed countries in Europe, the USA, Africa and South America (see [8] and references therein).
When mosquitoes are not subject to stimuli, it is possible to assume that they move randomly in any direction [9].This leads to a diffusion equation which can be extended to take into account the landscape heterogeneity or correlated random walk.Therefore, we have to incorporate advection terms or drift terms when mosquitoes, stimulated by attractants, move preferably in certain directions.
For simplicity, we present the generic PDE that models the spread of a mosquito (sub)population.Of course, it is now well recognized that the environment heterogeneity can have an important effect [10], [11].This is taken into account in the model by assuming spatial and temporal variations in the parameters.So, let u represent the number of insects.A possible modeling is to consider the following general advection-diffusion-reaction-like equation: in a bounded domain Ω ⊂ R n (where 1 ≤ n ≤ 3) with a smooth boundary ∂Ω.D (x, t) ≥ 0 is the diffusion (dispersion) coefficient or the diffusivity (D can eventually depend on u).Entomologists usually admit that there is no passive transportation of Aedes albopictus mosquito by the wind.Conversely, the blood-seeking mosquitoes will follow odors and carbon dioxide carried by the wind; this is modeled by the term V (x, t) .∇u.Indeed, it is well known that carbon dioxide (CO 2 ), in interaction with other components, acts as an attractant and induces a direct response to guide the mosquito towards the host.The breeding sites or the blood feeding sites attractions are modeled by the term ∇ • (C (x, t) u).
Thus it is necessary to take this important fact into account: the term ∇. • (C(x, t)u) represents a localized attractants (or a repellings) due to the presence or not of (blood or sugar) meals (like animals, humans or fruits...).C(x, t) represents an advection velocity toward favorable "places".The definition of C takes into account wind's direction and strength.For the sake of simplicity, the effective attraction areas are represented as ellipses.The attractor is set as one of the foci of the ellipse.The other focus point is calculated as a function of wind's direction and strength.Outside the ellipse, there is no attraction and the related advection term is equal to zero.Note that if there is no wind, the attraction area is reduced to a disk of which the center is the attractor.The reaction term g can be nonlinear, and represents the time-evolution of the population (death-birth rates, migration, vector control...), and thus may depend on the mosquito population u, and some environmental parameters (temperature, position in the domain,...).We suppose that u (x, 0) = u 0 (x) for x ∈ Ω, where u 0 can be a continuous (or possibly discontinuous) function.It may be possible to consider generalized boundary conditions, like Robin conditions, − → ∇u • n + αu = ρ(x, t), for all (x, t) ∈ ∂Ω × (0, T ), where n stands for the exterior unit normal to the boundary ∂Ω.For the sake of simplicity we will consider homogeneous Neumann boundary conditions, i.e. α = 0. Finally, we deduce the following (quasi)linear parabolic equation: for all x ∈ ∂Ω, and t > 0, Problem ( 2) is a (quasi)linear parabolic equation, for which it is possible to show the existence of a local solution.Moreover, under mild conditions it may be possible to show that the solution is global [12].

A. The Compartmental Model
In order to obtain some biologically interesting simulations, we take into account some biological facts about Aedes albopictus [3].There are two main stages in the development of mosquitoes: an aquatic stage and an adult stage.The aquatic stage gathers eggs, larvae and pupae.The adult stage can be divided into several compartments: immature females, blood feeding females, breeding females, resting females and males.Since we assume no sex differences in the aquatic stage, mosquitoes, after emergence, are distributed between the immature female compartment and the male compartment, according to r, the ratio of adult females mosquitoes to the total mosquito population.After mating with males, immature females enter the feeding female compartment and seek for blood meals before going into the resting compartment.Afterwards, the females pass into the breeding compartment seeking for a breeding site to deposit eggs.Once egg deposit is done, breeding females need blood and pass into the feeding compartment again.The eggs laid by the breeding females supply the aquatic stage.
In order to simulate SIT control, we add a sterile male compartment.Sterile males are released in specific places.The transmission rate between the immature females and the feeding females is conditioned by the proportion of non-sterile males to the whole male population, near the immature females.Contrary to Anopheles mosquito, Aedes males are looking for young females and thus, in general, they are located near the breeding sites or near the hosts.
We assume that resting females are not subjected to the attraction of blood meals or breeding sites.Resting females diffuse slowly, and their direction can be affected by the wind.The rate of transmission from one compartment to another allows us to take into account the average time spent by mosquitoes in each compartment.According to the previous explanations, we derive model (3).After rescaling, we consider Ω = [−1, 1] 2 , Q T = Ω × (0, T ].We set three indicator functions, 1 b , 1 f and 1 s .Function 1 b defines the area where the breeding females found a breeding site to lay eggs and become feeding females. 1 f defines the area where feeding females found a blood meal and become resting females. 1 s defines the area where sterile males are released.C f represents the attraction due to blood meals, like houses, and C b represents the attraction due to the breeding sites.M A , M Y , M f , M m and M ms are respectively the mortality rates for the aquatic stage, the immature females, the mature females, the wild males and the sterilized males.η A is the egg hatching rate, β Y is the rate at which immature females become bloodfeeding females, µ f r is the rate at which blood-feeding females become resting females, µ rb is the rate at which resting females become breeding females, µ bf is the rate at which breeding females become blood-feeding females, N Egg is the average number of eggs laid per female, K is the carrying capacity of a breeding site, and Λ s is the number of sterile males released periodically each τ days.Following [4], [5], we assume that the probability that an immature female becomes a feeding female depends on the ratio uM uM +uMs , when sterile males are released.System (3) can be rewritten in the following way, with with 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 Problem ( 4) is a (quasi)linear weakly coupled parabolic system, for which it is possible to show the existence of a local solution.The existence of a unique global solution may be proved using, for instance, [13].
To provide numerical simulations of system (3), we need to construct a reliable algorithm, that preserves most of the properties of the system (positivity of the solution, equilibrium, if any, and its (un)stability properties, ...).

B. The Numerical Methods
We consider an operator splitting method [7].This is an interesting method that enables us to solve separately each term of equation (1).So, we will solve successively the convective term, the diffusive term, and the reaction term, using the most efficient numerical method for each process.The full system can also be rewritten as follows where A represents the advective (or convective) terms, D the diffusive terms and R the reaction terms.Briefly, • the advection process is solved using the Corner Transport Upwind scheme (CTU) [14].It is a total variation diminishing scheme, that preserves the monotonicity of the solution providing that we verify a CFL condition between the convective parameters, the space-steps and the time-step.Moreover, this scheme minimizes the numerical diffusion and is even exact when ∆t is chosen appropriately.• the diffusion process is solved using the method of lines (MOL) for which we consider the secondorder finite difference method for the spatial discretization, and the TR-BDF2 (Trapezoide Rule -2nd order Backward Difference Formula) method for the time discretization.• the reaction process is solved using the nonstandard finite difference method [15].More details are developed in a forthcoming paper (see also [6]).
Our scheme permits to provide several simulations with and without constant parameters.In particular, we show that the solution converges, at least numerically, to a steady state.The numerical algorithm is implemented in Scilab [16], while the visualization is obtained with "R" [17].

III. SIMULATIONS AND DISCUSSIONS
We present some simulations in a homogeneous landscape, with time independent parameters, where there are 4 breeding sites and 5 blood-feeding sites as in Fig. 1.We assume that all these attractors have the same attraction force and domains of attraction.For each attractor, the area of attraction is defined as an ellipse depending on the wind's direction and strength.Note that, if there is no wind, the area of attraction is reduced to a disk.When the landscape is homogeneous, without wind, the mature females, i.e., the blood feeding females, the  resting females and the breeding females, tend to gather around the houses and the breeding sites (Fig. 2).This distribution is quite disturbed when wind and vector control get involved.Fig. 3, show the distribution of females when we consider North Wind, the removal of the two East breeding sites and the releases of sterile males.We notice that the majority of the mature females tend to migrate North, following odors or CO 2 carried by the wind.They gather around the Western attractors since the two East breeding sites no longer exist.Moreover, we notice that the number of mature females has decreased.This decrease can be explained by the fact that sterile males have been released, but also because of the North migration that allows mosquitoes to leave the domain.
SIT is respectful towards the environment, and is likely to be used as an alternative to the use of chemical products.On a temporal scale, we compare the impact of a 7-days periodic massive spraying of Deltamethrin around houses over a one-month treatment with 7-days  pulsed sterile males releases near the breeding sites over the same period of time.Fig. 6 shows that with an appropriate choice in the periodicity of the releases, and the number of released sterile males, SIT could be a promising alternative to massive spraying.
However, in order to have the most efficient results using SIT, it is important to have a good knowledge on the environment and take it into account.For instance, compare Fig. 4 and Fig. 5: they show the distribution of mature females controlled with SIT near the North breeding sites, and near the South breeding sites, respectively.We know that, with North wind, females tend to migrate towards the North; this is also the case for the released sterile males.That is why Northern releases have a lower impact on the number of mature females compared to Southern releases that are more efficient.Thus, it is more efficient to release sterile males downwind rather then  upwind (on condition that the wind strength is not to high, otherwise mosquitoes will hide wherever they can rather than fly upwind).This result is even more obvious when we have a look at the temporal evolution of the number of mature females depending on which spatial strategy is chosen for SIT releases (Fig. 7).

IV. CONCLUSION
This work presents promising tool for modeling mosquito dispersal.Thanks to the numerical simulations we could point out interesting results that can be a great use for optimizing SIT.Our work can be helpful to propose new field experiments; in particular it may be possible to consider temperature-varying parameters.This will be presented in a forthcoming paper.We could also add other compartments (sugar feeding compartment, sterile female compartment,... ).As we could see, the environment has a non-negligible influence on mosquito dispersal, and some works need to be done about landscape ecology because, so far, little is known about the interactions between landscape, vegetation and Aedes albopictus dispersal [11].

Fig. 2 .
Fig. 2. Mature Females distribution with no wind, no vector control.

Fig. 3 .
Fig. 3. Mature Females distribution with North wind, mechanical control and releases of 1000 sterile males per week, from the 20th simulation day.

Fig. 4 .
Fig. 4. Mature Females distribution with North wind, and weekly releases of 1000 sterile males, near the North breading sites, from the 20th simulation day.

Fig. 5 .
Fig. 5. Mature Females distribution with North wind, and weekly releases of 1000 sterile males, near the South breading sites, from the 20th simulation day.

Fig. 6 .
Fig. 6.Comparison of the evolution of the number of mature females with a weekly massive spraying control, and SIT with weekly releases of 2000, 4000, 10000 and 20000 sterile males over one month.

Fig. 7 .
Fig. 7. Evolution of the density of mature females with Northern wind with two spatial strategies of SIT ( 1000 sterile males released per week) from the 20th day.