Modelling Large Scale Invasion of Aedes aegypti and Aedes albopictus Mosquitoes

The principle aim of this work is to simulate the invasion of two invasive mosquito species Aedes aegypti and Aedes albopictus in central Europe at a landscape scale. The spatial-temporal dynamics of invasion is investigated in dependence of predation pressure, seasonal variation of ambient temperature as well as human population density. The introduction of temperature dependent entomological parameters enables the simulation of seasonal pattern of population dynamics. The influence of temperature, predation pressure and human population density on invasion is studied in one-dimensional cases. In two dimensions, georeferenced parameters such as annual mean temperature and human population density are prepared by a geographical information system and introduced into the finite element tool COMSOL Multiphysics. The results show that under the current temperature, central Europe cannot become a permanent breeding region for Aedes aegypti. However, southwest Germany especially the regions along the Upper Rhine Valley may provide suitable habitats for the permanent establishment of Aedes albopictus. An annual temperature rise of two degrees would lead to dramatic increase of invasion speed and extension range of Aedes albopictus.


Introduction
Vector-borne disease dengue fever is mainly transmitted by two invasive mosquito species: Aedes (Stegomyia) aegypti (L.) and Aedes (Stegomyia) albopictus (Skuse) [1].Viruses of this disease can be transmitted from human to human through the biting of infected female mosquitoes of both species.
The primary vector of dengue fever Aedes aegypti is also known as the yellow fever mosquito, since it is well-known for the transmission of yellow fever virus.
It came originally from Africa, but now it distributes widely in many tropical and subtropical regions around the world.Its domestic form feeds almost solely on human blood [2] and can use a broad range of artificial containers like cans or flower pots as its breeding sites, which facilitate its thrive in densely populated urban areas.
Another important vector of dengue fever is Aedes albopictus.It came originally from Asia, and is also known as Asian tiger mosquito.It belongs to the top 100 invasive species of the world [3].The domestic form of Aedes albopictus shares similar ecological niches as Aedes aegypti.However, besides artificial containers it can also use a wide range of natural water containers such as tree holes or bamboo stumps as breeding sites [4].Moreover, temperate strains of Aedes albopictus can survive the temperate climate in the form of diapausing eggs.Hence, besides tropical and subtropical regions it can also thrive in many temperate regions such as southern Europe.Recently, it is even found in central and northern European countries like Germany and the Netherlands [5] [6].
The widely geographical distribution of these two mosquito species poses big threats to the public health worldwide.Even more concerns have been aroused nowadays, since these two species may extend their expansion range to less favorable regions as a result of climate change [7] [8] [9] [10].
Since there is no dengue vaccine, and no commercial chikungunya vaccine is available [1], minimizing the population of these two mosquito species is the most effective way to avoid the outbreak of related vector-borne diseases.In order to draw up effective mosquito control plans a good understanding of biology and habitat preference of these two species is necessary.Mechanistic models which consider the life cycles and dispersal of mosquitoes are useful tools to quantitatively analyze the spatial-temporal population dynamic of mosquitoes.
Moreover, the effects of climate change on the range expansion of mosquitoes can be simulated.
Several studies have already been carried out.Takahashi et al. [11] have developed a one-dimensional model to describe the dispersal dynamics of Aedes aegypti.The model system consists of two compartments, by which diffusion and advection of adult females were considered.Travelling wave solutions were obtained numerically.Yang et al. [12] have obtained temperature response functions for different entomological parameters of Aedes aegypti based on temperature-controlled experiments.Then, these functions were inserted into a compartment model to assess the influences of temperature on the population dynamics of Aedes aegypti.Richter et al. [13] have modified the model of Takahashi et al. [11] [14] have adopted a stochastic spatial model to describe the spatial-temporal dynamics of Aedes aegypti in Buenos Aires.The mean daily temperature variation was described by a simple cosinus function.
Most of the above-mentioned mechanistic models were built in one-dimensional case.Variables like temperature were considered as spatially homogeneous.
Moreover, many entomological parameters were assumed as temperature independent.Additionally, many models are not suitable for application in temperate regions, since the seasonal temperature variation and overwintering strategy of mosquitoes were not taken into account, which may play decisive roles for the permanent establishment of mosquito in these regions.
The aim of the present work is to simulate the invasion of Aedes agypti and Aedes albopictus in central Europe.The model system employed here is based on the compartment model of Richter et al. [13].Besides reproduction rate, other entomological parameters such as mortality rate are modified as different temperature response functions based on the experiment data of Yang et al. [12].In the work of Richter et al. [13], besides temperature human population density was also introduced into the model system as an important variable to differentiate the predation pressure between urban areas and natural regions.In the present work, human population density is also applied to determine the environmental capacity of mosquito larvae.In this way, habitat availability of both mosquito species and their close associations with humans can be better described.Furthermore, in order to simulate the seasonal pattern of population dynamic and overwintering of both species, the daily mean temperature variation and winter diapause of Aedes albopictus are taken into account.

Materials and Methods
The invasion of both Aedes aegypti and Aedes albopictus was simulated in one-dimensional cases in Wolfram Mathematica 8.0 at first.The system of reaction-diffusion equations was solved by the numerical methods of lines.Based on the work of Richter et al. [13], the model systems was implemented in COMSOL Multiphysics (version 4.2a) in order to simulate the invasion of both mosquito species at a landscape scale.Here finite element method was applied as the numerical technique.Moreover, temperature and human population density data as well as landscape structure were prepared in a geographical information system and linked to COMSOL Multiphysics.Geo-referenced annual mean temperature datasets of central Europe were imported from the WorldClim Global Climate Data (http://www.worldclim.org/),which provides a series of global climate layers with a spatial resolution of about 1 km 2 .These temperature data were then interpolated into COMSOL environment using the two-dimensional linear interpolation option in the interpolation menu.Geo-referenced human population density data were obtained from DLM250 (digital landscape models at a scale of 1:250,000) of ATKIS (Amtliches Topographisch-Kartographisches Information System) of year 2007.Landscape structures were used to control the sizes of finite element grids.They were converted from shape files to CAD files by using the "Export to CAD" tool from ArcToolbox in ArcGIS.The CAD files were then imported into the COMSOL Multiphysics.

Model Equations
The system of reaction diffusion equations is applied in the present work.The general form of reaction diffusion equations is as follows In the present study, the spatial operator has the simple form [ ] ( ) The system of equations is presented as follows: ( where A is the population density of aquatic phase, M is the population density of winged phase (adult females); Φ is the intrinsic oviposition rate; C(P) is the environmental capacity of aquatic phase, which is defined as a function of human population density P as both species are highly domestic and the available breeding sites through the use of artificial containers are closely related to P; f(T), γ(T), μA(T), μM(T) are temperature response function of oviposition rate, emergence rate, mortality rate of aquatic phase and mortality rate of winged phase, respectively; β is the predation rate, whereas K s is the saturation constant of predators; D and v are the diffusion coefficient and advection coefficient of winged phase respectively; r is the proportion of females.In the model system the life stages of mosquito are divided into an aquatic phase (egg, larva and pupa) and a winged (adult) phase.It is assumed that all the eggs laid by females can survive and hatch to larvae.The sex ratio of mosquitoes is assumed to be 1:1 throughout their whole life stage i.e. r = 0.5, which is plausible as the study of Delatte et al. [15] showed that the sex ratios of Aedes albopictus is not statistically different from 0.5 within a wide range of temperatures.Furthermore, the temperature response function of mortality rate and the ability of dispersal for both males and females are supposed to be identical.Additionally, it is assumed that all females can be always copulated effectively, which was adopted in the work of Yang et al. [12].Hence the male phase can be disregarded in the model system.The diffusion and advection are only considered for the Advances in Pure Mathematics winged phase.Since the movements of larvae and pupae are mainly constrained in small water bodies such as artificial containers, their dispersals are neglected here.However, the aquatic phase is subject to predator pressure.The temperature response of predator pressure is modelled simply by multiplying the predator pressure term with f(T).
In the work of Maidana et al. [16], environmental capacity was introduced not only to the aquatic phase but also to the winged phase in order to regulate the size of the population.In the present work however, only an environmental capacity for the aquatic phase was applied, since in a natural environment the mosquito population is mainly limited by the availability of breeding sites [17].
The environmental capacity of aquatic phase is given by a function of human population density: where coefficient C 0 can be interpreted as the environmental capacity of aquatic phase in an environment without human population.In a peridomestic environment, the increase of breeding sites through the use of artificial containers by adult females for oviposition is considered as a number k of larvae per human.k is closely related to the Breteau Index, i.e. the number of breeding sites per one hundred houses.In the following simulations, k will be set as 3, which corresponds to a high Breteau Index [18].

Temperature Dependent Entomological Parameters
The entomological data applied in the present work come from the temperature-controlled experiment of Yang et al. [12] to access the entomological parameters of Aedes aegypti, in which Aedes aegypti strain was captured from Brazil.
In the work of Yang et al. [12] only polynomials are used as temperature response functions for all the entomological parameters.In the present work, different types of temperature response functions which correspond better to the biological features are applied for different entomological parameters.The parameters of these functions are obtained through nonlinear regression (subroutine "Non Linear Model Fit" in Mathematica, based on Quasi-Newton method and others).
To describe the temperature response of oviposition rate f(T), the O'Neil equation with biological parameters T opt (optimum temperature), T max (lethal temperature) and Q 10 is applied.

( ) (
) with ) According to the experimental data, T max and T opt are fixed as 39˚C and 31˚C, respectively.The data and temperature response curve of oviposition rate is shown in Figure 1.
The data and fitting curve of transition rate from aquatic phase to winged phase are shown in Figure 2.For the temperature response function following equation is applied.
( ) Equation ( 9) is applied as the temperature response function of mortality rate of both aquatic phase and winged phase.The data and fitting curves are shown in Figure 3.

Analysis of the Homogeneous Part
Based on the work of Richter et al. [13], the stationary solution of winged phase M can be obtained as a function of aquatic phase density A easily. ( where G(A) is the growth rate of aquatic phase density.
Equation (11) has three stationary solutions: the minimal viable aquatic phase density: and the maximum aquatic phase density: with The predation term causes a strong Allee effect, i.e. the growth rate is negative when the population density lies below a critical value [19], which is illustrated in Figure 4(a).Advances in Pure Mathematics The criterion of stability ( ) as a function of temperature is shown in Figure 5.The states A = 0 and A = A s3 are stable, whereas state A = A s2 is unstable.
The stationary solutions are governed by temperature T, predation pressure β and human population density P. All these three parameters determine the regions of viability.In the case of Figure 4 the influence of human population density is not considered.higher than β1.If predation pressure is larger than β2, the Allee effect is so strong that it leads to only one stable state i.e. zero state.If the predation pressure lies between the two bifurcation points, there are two stable states which are separated by an unstable state (the dashed line).In Figure 4(c) the stationary solutions as a function of temperature T is represented (β fixed).Two bifurcation points are noted as T 1 and T 2 .A region of viability exists between T 1 and T 2 .
In Figure 4(d) a three dimensional plot of the solution surfaces are obtained by making both temperature and predation pressure as bifurcation parameters.The upper solution surface is locally stable, whereas the lower surface is unstable separating the regions of attraction of the upper surface (maximum aquatic density) from that of the zero plane, which is locally stable in a certain region due to the strong Allee effect.

Seasonal Pattern of Population Dynamic of Aedes aegypti and Aedes albopictus
In a real environment mosquitoes are exposed to seasonal temperature variations.The population dynamics of mosquitoes as well as the transmission rate of Dengue fever fluctuate in different seasons significantly.The winter temperature is a decisive factor determining the permanent breeding zones of Aedes aegypti and Aedes albopictus in a region.Unlike Aedes albopictus, Aedes aegypti cannot produce winter diapausing eggs, hence its main distribution regions are limited in tropical and subtropical regions since it is not able to resist the cold winter temperatures in temperate regions.Aedes albopictus however, can survive the winter temperature in many temperate regions in the form of diapausing eggs.
In Europe, it has already established itself stably in many regions of southern Europe, and poses recently a great threat to central Europe.Therefore, the seasonality is introduced to the model system.At the mean time the mechanism of winter diapause of Aedes albopictus will be modeled as the primary difference to Aedes aegypti.
To describe the seasonal temperature variation a simple sine-function (Equation ( 16)) is applied, where T m and T r are annual mean temperature and daily mean temperature range, respectively (see Figure 6(a)).Coefficient t 0 determines the phase of the sinusoid.If t 0 is set as 107, then the annual minimum temperature comes on day 16 (January 16), whereas the highest temperature occurs on day 198 (July 16).

( ) sin 2π 365 t t T t Tm Tr
Aedes aegypti Figure 6 demonstrates the seasonal population dynamics of Aedes aegypti for both aquatic and winged phase.Aedes aegypti is introduced on day 0 (mid-March) at origin.The human population density in the study area is 2000 (1/km), while predation pressure β is 0.3 C 0 (day −1 ).The temperature variation is shown in Figure 6(a) with an annual mean temperature of 18˚C and a daily

Aedes albopictus
Winter diapause of Aedes albopictus is mainly governed by temperature and photoperiod.According to the study of the seasonal dynamics of hatching and oviposition rate of Aedes albopictus in Rome [20], the onset of hatching in spring happens, when the mean temperature reaches 10˚C to 11˚C and day length surpass over 11 -11.5 h of light.Dramatic decrease of the hatching rate of summer eggs during autumn occurred during mid-September, when the day length declines below 14 h of light.However, about 21% of eggs kept hatching when the temperature was still higher than 10˚C.In the present work, only temperature was considered as the trigger factor for the onset and ending of winter diapause.In order to simplify the model, a sigmoid function i.e.Equation ( 17) is applied to control both the starting of hatching in spring and the entering of diapause in autumn (see Figure 7).The model system is then modified from Equation (3) to Equation (18).

( ) ( )
A compare of the observation data and simulation results is made in order to evaluate the plausibility of the model.According to the study of [20] in Rome during 2000, adult Aedes albopictus mosquitoes were found to be active since late March, and were most abundant at the end of August (week 35).After that, summer egg production decreased rapidly in week 36.A small number of adults could be found until the end of December.Totally, this species was found to be active for about 10 months during a year.
To simulate the population dynamics of Aedes albopictus in Rome, a temperature variation curve was obtained by fitting the temperature function i.e.Equation (16) to the temperature data provided by the Ufficio Centrale di Ecologia Agraria [20] (see Figure 8(a)).The simulation results are shown in Figure 8(b).
Since the unfavorable temperatures for the population development between 10˚C and 15˚C last for a relative long period by the given temperature course, a predation pressure of 0.3 C 0 (day −1 ) leads to the dying out of this species.Hence a relative low predation pressure of 0.1 C 0 (day −1 ) was given.This is reasonable, since in a densely populated city like Rome the predation pressure is relative low.According to the simulation results, the winged phase is active since around day 100 (early April), and its density peaks in week 34.During the winter, a small population of winged phase stays active until the beginning of January.The simulation results described above fit well with the observation data of [19], which verifies that the model for the winter diapause is reliable.However, a rapid decrease of the summer eggs is only found until the end of October (around week 48), which is not consistent with the observation.The reason lies in the fact that photoperiod is the primary trigger for the production of winter diapausing eggs [20], which is however not taken into account in the model.

Dispersal of Aedes aegypti and Aedes albopictus in Central Europe
Based on the simulations in one dimension, the dispersal of both mosquito species are simulated at a landscape scale.Since temporary breeders like Aedes aegypti and Aedes albopictus can adapt well to the peridomestic environment, where predator pressure are much lower than that in a natural habitat, Equation ( 19) is applied to distinguish the significantly different predator pressures between urban areas and natural regions [13].
( ) ( ) max , , exp By using the parameters listed in Table 1 the invasions and spreads of Aedes aegypti and Aedes albopictus in central Europe were simulated for a period of three years.Both mosquito species are introduced in Frankfurt during middle April.For all the four boundaries Dirichlet boundary condition was applied which assumes that the mosquito population density stays zero on the boundaries of the domain.The simulation results are shown in Figure 9 to Figure 13.
In Figure 9        3.4), adults Aedes albopictus begin to be active since late march and can extend its period of activity for about 10 months per year [20].In Frankfurt, where the annual mean temperature is 9.7˚C, the active time reduces to about 6 month.[22].

Effect of Temperature on the Dispersal of Aedes albopictus
In the following scenarios the range expansion of Aedes albopictus will be simulated under current temperature as well as a mean temperature rise of 2˚C.
The population is introduced in Frankfurt during middle April.A totally dispersal of 20 years of this species was simulated.The range expansion of Aedes albopictus under current temperature and under a temperature increase of two degrees were demonstrated in Figure 14 and Figure 15, respectively.
Comparing these two figures one can find that under a temperature rise of 2˚C Aedes albopictus can spread much faster.For example, 10 years after introduction the expansion range of this species under a temperature increase of two degrees is considerably larger than that of the species under the current temper- Aedes albopictus is mainly limited in the west of the Swabian Jura and Franconian Jura (two mountain ranges in southern Germany).However, this species can pass through the small warm areas between mountain regions and even colonize new urban areas in south-eastern Germany under a temperature increase of two degrees, which is demonstrated in Figure 15(d).
Under a higher temperature Aedes albopictus can establish itself more abundantly.Additionally, the increase of temperature can reduce the minimum viable population density for this species.As a consequence, the regions, which are not favorable for this species under the current temperature, can become suitable areas as the mean temperature rises.

Conclusions and Recommendations
In the present work, the spatial-temporal population dynamics of two invasive mosquito species, Aedes aegypti and Aedes albopictus, were qualitatively and quantitatively analyzed by using compartment models based on biology and habitat preference of both species.As specialties, the Allee effect, seasonal variation of ambient temperature, dependency of predation pressure on the human population density as well as winter diapause as the overwintering strategy of  Aedes albopictus were considered in the simulation.Moreover, the range expansions of both mosquito species in central Europe were simulated at a landscape scale.The effect of temperature rise relating to global warming on the range expansion of Aedes albopictus was also simulated.
The main conclusions of the present study are listed as follows.
1) Since both Aedes agypti and Aedes albopictus are highly domestic, human population density can affect the invasion speed as well as expansion range of both species.Densely populated urban areas provide not only abundant breeding sites, but also safe habitats with low predator pressure for both mosquito species.
2) Introduction of seasonal temperature variation is vital for the simulation of permanent establishment of both mosquito species in a region.Minimum winter temperature and predation pressure are two primary limiting factors for the permanent breeding of Aedes aegypti in temperate regions, whereas predation pressure and duration of unfavorable temperature after winter diapause determine the successful overwintering of Aedes albopictus in temperate regions.
3) Central Europe is not likely to become a permanent breeding region for Aedes aegypti, since this species is not able to survive the winter temperatures in this region.However, some densely populated urban areas like Frankfurt may become its "temporary regions", when it is introduced during months with favorable temperature.5) Ambient temperature exerts a significant influence on the propagation of Aedes albopictus in central Europe [13].While this species spreads relatively slowly and its expansion range is largely limited by the mountain areas at the current temperature, with an annual mean temperature rise of two degrees it can propagate much faster and cross some mountainous areas in southwest Germany.
The model employed in the present work is based on the biology as well as ecological niche of both mosquito species.However, related to both aspects there are some limitations in the model as a result of the shortage of available data and insufficient research.
One of the major objectives of this work was to simulate the invasion of Aedes aegypti and Aedes albopictus in central Europe.Since there are abundant variations in genes, morphology and ecology of both species around the world [23] [24], temperature-controlled experiments with both species conducted with the strains that have already been present in Europe or adjacent regions are necessary.For example, Aedes albopictus strains for study should be directly captured from Italy, southern France or Switzerland, since the risk of invasion of this spe-cies from these regions into central Europe is greatest.
Besides studies focusing on entomological parameters, more quantitative studies or surveillance programs on habitat preferences and biological invasions of both mosquito species should be performed.Some parameters in the model system like environmental capacity, predation pressure or diffusion coefficient are dependent on many factors, which makes direct measurements very difficult.So it would be more practical if these parameters were indirectly estimated through model calibration.For example it would be quite meaningful if the present model was first applied to simulate the range expansion of Aedes albopictus in Italy.This species has invaded Italy since 1990, and has now spread to most areas of the country with an altitude below 600 m [25].The model would be more plausible if parameters were adjusted by comparing the simulation results with observation data collected for more than twenty years.
In the present work the influence of microclimate on the invasion and permanent establishment of both species was not considered in the simulation.It is very likely that Aedes aegypti will be able to find necessary protections against winter weather (for example, through breeding in artificial containers) in some cities in North America [2] [26].Other factors which can significantly affect microclimate and hence determine the survivorship of mosquito species should also be taken into account in further work.Additionally, besides temperature, the change of other environmental conditions such as precipitation should also be considered when simulating the effect of climate change on the establishment of both species in central Europe.
inserting Equation(10) into the right hand side of the first equation of Equation (3), the stationary solutions of aquatic phase density A can be obtained.
Figure 4(b) shows the three stationary solutions as functions of predation pressure β (T fixed).The bifurcation points β1 and β2 are specially marked.The Allee effect occurs only when the predation pressure is

Figure 6 .
Figure 6.Population dynamics of Aedes aegypti under the consideration of seasonality.(a) mean daily temperature variation with T m = 18˚C, T r = 6˚C, t 0 = 107; (b) seasonal variation of aquatic phase density; (c) seasonal variation of winged phase density; (d) seasonal variation of both aquatic and winged phase density at origin.

Figure 8 .
Figure 8.The temperature variation and population dynamics of Aedes albopictus in Rome.(a) The temperature variation of Rome in 2000, data from UCEA [21]; (b) Simulated population dynamics of aquatic and winged phase of Aedes albopictus.
and Figure 10 the population dynamic of Aedes aegypti in both aquatic as well as winged phase in different times of year are demonstrated.After the introduction in Frankfurt during late March, this species has established itself well in densely populated urban areas and gradually spread towards adjacent

ature (Figure 14
(b) and Figure 15(b)).As shown in Figure 14(c) and Figure 14(d), under current temperature the range expansion of Aedes albopictus has not changed significantly from year 15 to year 20.The mountainous areas have obviously acted as great barriers for the dispersal of this species from north area towards southeast area in the study area.

Figure 15 .
Figure 15.Range expansion of Aedes albopictus under a temperature increase of two degrees.(a) 5 years after introduction; (b) 10 years after introduction; (c) 15 years after introduction; (d) 20 years after introduction.

4 )
Under the current temperature, southwestern Germany, especially the areas along the Upper Rhine Valley, can provide suitable habitats for the permanent establishment of Aedes albopictus.It can also survive the winter temperatures in this region in the form of diapausing eggs.Hence persistent surveillance and control of Aedes albopictus in this area are necessary.

Table 1 .
Summary of model parameter values.
Belonging to the top 100 invasive species of the world, Aedes albopictus has successfully invaded the southern Europe during the last decades.The distribution of this species in Europe may be altered dramatically in the future as a result of the climate change.Reliable regional climate change projections for Europe