Modelling Large Scale Invasion of Aedes aegypti and Aedes albopictus Mosquitoes

Show more

1. 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] by adding a predation term and introducing temperature dependent reproduction rate. Moreover, by importing Geo-referenced model parameters into a finite element tool, invasion of Aedes albopictus was simulated at a scale of central Europe. Besides deterministic models mentioned above, stochastic models have also been applied to simulate the population dynamics of mosquitoes. Otero et al. [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.

2. 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.

3. Modelling the spatial-Temporal Dynamics of Aedes aegypti and Aedes albopictus in one Dimension

3.1. Model Equations

The system of reaction diffusion equations is applied in the present work. The general form of reaction diffusion equations is as follows

$\frac{\partial {u}_{i}}{\partial t}=L\left[u\right]+{f}_{i}\left({u}_{1},\cdots ,{u}_{n}\right)$ (1)

In the present study, the spatial operator has the simple form

$L\left[u\right]=\nabla \left(D\nabla u\right)$ (2)

The system of equations is presented as follows:

$\begin{array}{l}\frac{\partial A}{\partial t}=f\left(T\right)\text{\hspace{0.05em}}\varphi \text{\hspace{0.05em}}\left(1-\frac{A}{C\left(P\right)}\right)M-\left(\gamma \left(T\right)+{\mu}_{A}\left(T\right)\right)A-\frac{\beta \text{\hspace{0.17em}}A}{A+{K}_{s}}\text{\hspace{0.05em}}f\left(T\right)\\ \frac{\partial M}{\partial t}=\nabla \cdot \left(D\text{\hspace{0.05em}}\text{\hspace{0.05em}}\nabla M-vM\right)+r\gamma \left(T\right)A-{\mu}_{M}\left(T\right)M\end{array}$ (3)

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 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:

$C\left(P\right)={C}_{0}+k\text{\hspace{0.05em}}P$ (4)

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] .

3.2. 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.

$f\left(T\right)=\text{\hspace{0.17em}}{\left(\frac{{T}_{\mathrm{max}}-T}{{T}_{\mathrm{max}}-{T}_{\text{opt}}}\right)}^{p}\text{exp}\left(\frac{p\left(T-{T}_{\text{opt}}\right)}{{T}_{\mathrm{max}}-{T}_{\text{opt}}}\right)$ (5)

with

$p=\frac{1}{400}\text{\hspace{0.17em}}{{\displaystyle W}}^{2}{\left[1+\sqrt{1+\frac{40}{W}}\right]}^{2}$ (6)

and

$W=\left({Q}_{10}-1\right)\left({T}_{\mathrm{max}}-{T}_{\text{opt}}\right)$ (7)

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.

$\gamma \left(T\right)={n}_{0}\text{\hspace{0.05em}}\left[1-\text{exp}\left(-{\left(\frac{T}{{T}_{1}}\right)}^{r1}\right)\right]\text{\hspace{0.05em}}\text{exp}\left[-{\left(\frac{T}{{T}_{2}}\right)}^{r2}\right]$ (8)

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.

$\mu \left(T\right)=1-\text{\hspace{0.05em}}\left[1-\text{exp}\left(-{\left(\frac{T}{{T}_{1}}\right)}^{r1}\right)\right]\text{\hspace{0.05em}}\text{exp}\left[-{\left(\frac{T}{{T}_{2}}\right)}^{r2}\right]$ (9)

Figure 1. The temperature response curve of oviposition rate of Aedes aegypti. The measurement data are courtesy of Yang et al. [12] . Estimated model parameter values: T_{max}: 39˚C (fixed), T_{opt}: 31˚C (fixed), Φ: 9.35 (day^{−1}), Q_{10}: 2.20.

Figure 2. The temperature response curve of transition rate of Aedes aegypti. The measurement data are courtesy of Yang et al. [12] . Estimated model parameter values: n_{0} = 0.27 (day^{−1}), T_{1} = 33.25 (˚C), r_{1} = 2.86, T_{2} = 37.16 (˚C), r_{2} = 63.80.

Figure 3. The temperature response curves of aquatic and winged mortality rate of Aedes aegypti. The measurement data are courtesy of Yang et al. [12] . (a) fitting curve of aquatic mortality rate, estimated parameter values: T_{a}_{1} = 7.44 (˚C), r_{a}_{1} = 1.58, T_{a2} = 48.19 (˚C), r_{a}_{2} = 6.45; (b) fitting curve of winged mortality rate, estimated parameter values: T_{w}_{1} = 1.16 (˚C), r_{w}_{1} = 0.42, T_{w}_{2} = 42 (˚C) (fixed), r_{w}_{2} = 14.05.

3.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.

$M=\frac{r\gamma A}{{\mu}_{M}}$ (10)

By 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.

$G\left(A\right)=f\left(T\right)\varphi \text{\hspace{0.05em}}\left(1-\frac{A}{C\left(P\right)}\right)\text{\hspace{0.17em}}\frac{r\gamma \left(T\right)A}{{\mu}_{M}\left(T\right)}-\left(\gamma \left(T\right)+{\mu}_{A}\left(T\right)\right)A-\frac{\beta \text{\hspace{0.17em}}A}{A+{K}_{s}}f\left(T\right)=0$ (11)

where G(A) is the growth rate of aquatic phase density.

Equation (11) has three stationary solutions:

${A}_{S1}=0$ (12)

the minimal viable aquatic phase density:

$\begin{array}{l}{A}_{s2}=-\frac{1}{2\text{\hspace{0.17em}}r\text{\hspace{0.05em}}\gamma \text{\hspace{0.05em}}\left(T\right)\varphi f\text{\hspace{0.05em}}\left(T\right)}\cdot [C\left(P\right))\gamma \text{\hspace{0.05em}}\left(T\right){\mu}_{M}\text{\hspace{0.05em}}\left(T\right)+C\left(P\right)\text{\hspace{0.05em}}{\mu}_{A}\text{\hspace{0.05em}}\left(T\right){\mu}_{M}\text{\hspace{0.05em}}\left(T\right)\\ \text{}-C\left(P\right)\text{\hspace{0.05em}}r\text{\hspace{0.05em}}\gamma \text{\hspace{0.05em}}\left(T\right)\varphi f\text{\hspace{0.05em}}\left(T\right)+{K}_{s}\text{\hspace{0.05em}}r\gamma \text{\hspace{0.05em}}\left(T\right)\text{\hspace{0.05em}}\varphi f\text{\hspace{0.05em}}\left(T\right)-W]\end{array}$ (13)

and the maximum aquatic phase density:

$\begin{array}{l}{A}_{s3}=-\frac{1}{2\text{\hspace{0.17em}}r\text{\hspace{0.05em}}\gamma \left(T\right)\varphi f\left(T\right)}\cdot [C\left(P\right)\gamma \left(T\right){\mu}_{M}\left(T\right)+C\left(P\right){\mu}_{A}\left(T\right){\mu}_{M}\left(T\right)\\ \text{}-C\left(P\right)\text{\hspace{0.05em}}r\text{\hspace{0.05em}}\gamma \left(T\right)\text{\hspace{0.05em}}\varphi f\left(T\right)+{K}_{s}\text{\hspace{0.05em}}r\gamma \left(T\right)\text{\hspace{0.05em}}\text{\hspace{0.05em}}\varphi f\left(T\right)+W]\end{array}$ (14)

with

$W=\sqrt{{C}^{2}{\left(\gamma +{\mu}_{A}\right)}^{2}{\mu}_{M}{}^{2}+r\gamma \varphi f\left(T\right)(-2C\text{\hspace{0.17em}}\left(2\beta \text{\hspace{0.05em}}f\left(T\right)+\left(C+{K}_{s}\right)\left(\gamma +{\mu}_{A}\right)\right){\mu}_{M}+{\left(C+{K}_{s}\right)}^{2}r\gamma \varphi f(T))}$

(15)
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).

The criterion of stability
$\lambda =\frac{\partial G\left(A\right)}{\partial A}|A={A}_{\text{stationary}}$ [19] as a function of temperature is shown in Figure 5. The states A = 0 and A = A_{s}_{3} are stable, whereas state A = A_{s}_{2} 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. 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 4. Analysis of stationary solutions without the consideration of human population (modified after Richter et al. [13] ). (a) Growth rate as a function of aquatic phase density (T = 20˚C, β = 0.7 C_{0} (day^{−1}), K_{s} = 0.01 C_{0}); (b) Bifurcation diagram with predation pressure as bifurcation parameter (T = 20˚C, β = 0.7 C_{0} (day^{−1}), K_{s} = 0.01 C_{0}); (c) Bifurcation diagram with temperature as bifurcation parameter; (d) Bifurcation surface with both temperature and predation pressure as bifurcation parameters.

Figure 5. Stability analysis of the three stationary states as a function of temperature. C_{0} = 1000 (larvae/habitat). (a) β = 0.7 C_{0} (day^{−1}), (b) β = 0.3 C_{0} (day^{−1}).

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.

3.4. 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).

$T(t)=Tm+Tr\text{\hspace{0.05em}}\text{sin}\text{\hspace{0.05em}}\left(2\text{\pi}\text{\hspace{0.05em}}\frac{t-{t}_{0}}{365}\right)$ (16)

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

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.

mean temperature range of 6˚C. In Figure 6(b) and Figure 6(c) the propagation of Aedes aegypti as well as seasonal variation of aquatic phase density and winged phase density are illustrated with a three dimensional plot, respectively. In Figure 6(d) densities of both phases at origin as a function of time are shown. As demonstrated in Figure 6(d), the patterns of population dynamics of both phases coincide with each other very well. In both phases, the maximum population density is reached in mid-August, whereas the minimum density comes at the end of February.

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).

$WP\left(T\right)=\frac{1}{1+{\text{e}}^{-{k}_{0}\left(T-{T}_{0}\right)}}$ (17)

$\begin{array}{l}\frac{\partial A}{\partial t}=f\left(T\right)\varphi \left(1-\frac{A}{C\left(P\right)}\right)M-WP\left(T\right)\text{\hspace{0.05em}}\text{\hspace{0.05em}}\left(\gamma \left(T\right)+{\mu}_{A}\left(T\right)\right)A\text{\hspace{0.05em}}\text{\hspace{0.05em}}-\frac{\beta \text{\hspace{0.17em}}A}{A+{K}_{s}}\text{\hspace{0.05em}}f\left(T\right)WP\left(T\right)\\ \frac{\partial M}{\partial t}=\nabla \cdot \left(D\text{\hspace{0.05em}}\text{\hspace{0.05em}}\nabla M-vM\right)+WP\left(T\right)\text{\hspace{0.05em}}r\gamma \left(T\right)A\text{\hspace{0.05em}}\text{\hspace{0.05em}}-{\mu}_{M}\left(T\right)M\end{array}$ (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.

Figure 7. Winter diapause function of Aedes albopictus. T_{0} = 11˚C, k_{0} = 5.

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.

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.

4. Modelling large Scale invasion and control of Aedes aegypti and Aedes albopictus in Central Europe

4.1. 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] .

$\beta \left(x,y\right)={\beta}_{\mathrm{max}}\mathrm{exp}\left(-{\left(\frac{P\left(x,y\right)}{{P}_{tr}}\right)}^{n}\right)$ (19)

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 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

Table 1. Summary of model parameter values.

Figure 9. Population dynamics of Aedes aegypti (aquatic phase). (a) middle June (year 1); (b) middle August (year 1); (c) middle October (year 1); (d) middle November (year 1).

Figure 10. Population dynamics of Aedes aegypti (winged phase). (a) middle June (year 1); (b) middle August (year 1); (c) middle October (year 1); (d) middle November (year 1).

Figure 11. Population dynamics of Aedes albopictus (winged phase). (a) middle November (year 1); (b) end May (year 2); (c) middle August (year 2); (d) end August(year 2).

Figure 12. Population dynamics of Aedes albopictus (aquatic phase). (a) middle November (year 1); (b) end May (year 2); (c) middle August (year 2); (d) end August (year 2).

Figure 13. Population dynamics of Aedes aegypti and Aedes albopictus in Frankfurt. (a) winged phase; (b) aquatic phase.

regions. The most abundant population densities of both phases are found in middle August (see Figure 9(b) and Figure 10(b)). After that, this species spreads although further, its population density decreases as a result of the reduction of ambient temperature. Till the middle of November, both phases of this species have almost disappeared from the whole study area.

Similar as adults Aedes aegypti, adults Aedes albopictus have almost died out since the middle of November (see Figure 11(a)). However, this species can survive the cold winter in the form of diapausing eggs, which is demonstrated by Figure 12(a). At the end of May in the following year adult mosquitoes begin to be active after winter diapause (see Figure 11(b)). However, the population density of aquatic phase reduces (see Figure 12(b)). The winter eggs begin to hatch, but only a small portion of eggs can be reproduced since the population density of adult mosquitoes is still rather small during this period. This species establishes itself gradually faster and spreads further and further as the ambient temperature rises, colonizes new urban areas and reaches their highest population density in the middle of August once more (see Figure 11(c) and Figure 12(c)). At the end of August, the population densities of both phases begin to reduce rapidly (see Figure 11(d) and Figure 12(d)), at the mean time winter diapausing eggs start being produced.

In Figure 13 the population dynamics of both mosquito species in Frankfurt during the whole simulation period are illustrated. The population dynamics of both species are quite similar with each other until the end of September of the first year since the same biological parameters were applied. Aedes aegypti cannot survive the cold winter in the Upper Rhine Valley. In contrast, Aedes albopictus can survive in this region and thrill in urban areas during summer. In Rome, where the annual mean temperature is 17.2˚C during 2000 (see section 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.

4.2. Effect of temperature on thedispersal of Aedes albopictus

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 are available in the IPCC’s Fourth Assessment Report (AR4), which are based on variety of simulations under different emissions scenarios. Under the A1B scenario, an annual mean warming for 2080 to 2090 (with respect to 1980 to 1990) varies from 2.3˚C to 5.3˚C in the northern Europe and 2.2˚C to 5.1˚C in the southern Europe and Mediterranean region [22] .

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 temperature (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. 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.

5. 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

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

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.

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.

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.

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 species 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.

Acknowledgements

We thank Prof. Dr. Hyun Mo Yang for his kind permission to use the entomological data.

References

[1] WHO, Dengue and Sever Dengue. Fact Sheet No 117.

http://www.who.int/mediacentre/factsheets/fs117/en/index.html

[2] Reiter, P. (2010) Yellow Fever and Dengue: A Threat to Europe? Eurosurveillance, 15, 19509.

[3] Invasive Species Specialist Group, Global Invasive Species Database—Aedes albopictus (Insect).

http://www.issg.org/database/species/ecology.asp?si=109&fr=1&sts=sss&lang=EN

[4] Vezzani, D. and Carbajo, A.E. (2008) Aedes aegypti, Aedes albopictus and Dengue in Argentina: Current Knowledge and Future Directions. Memórias do Instituto Oswaldo Cruz, 103, 66-74.

https://doi.org/10.1590/S0074-02762008005000003

[5] Scholte, E., Den Hartog, W., Dik, M., Schoelitsz, B., Brooks, M., Schaffner, F., Foussadier, R., Braks, M. and Beeuwkes, J. (2010) Introduction and Control of Three Invasive Mosquito Species in the Netherlands, July-October 2010. Eurosurveillance, 15, 19710.

[6] Werner, D., Kronefeld, M., Schaffner, F. and Kampen, H. (2012) Two Invasive Mosquito Species, Aedes albopictus and Aedes japonicus japonicus, Trappd in South-West Germany, July to August 2011. Eurosurveillance, 17, 20067.

https://doi.org/10.2807/ese.17.04.20067-en

[7] Patz, J.A., Martens, W.J.M., Focks, D.A. and Jetten T.H. (1998) Dengue Fever Epidemic Potential as Projected by General Circulation Models of Global Climate Change. Environmental Health Perspectives, 106, 147-153.

https://doi.org/10.1289/ehp.98106147

[8] Kearney, M. Porter, W.P., Williams, C., Ritchie, S. and Hoffmann, A.A. (2009) Integrating Biophysical Models and Evolutionary Theory to Predict Climatic Impacts on Species’ Ranges: The Dengue Mosquito Aedes aegypti in Australia. Functional Ecology, 23, 528-538.

https://doi.org/10.1111/j.1365-2435.2008.01538.x

[9] Fisher, D., Thomas, S.M., Niemitz, F., Reineking, B. and Beierkuhnlein, C. (2011) Projection of Climate Suitability for Aedes albopictus Skuse (Culicidae) in Europe under Climate Change Conditions. Global and Planetary Change, 78, 54-65.

https://doi.org/10.1016/j.gloplacha.2011.05.008

[10] Caminade, C., Medlock, J. M., Ducheyne, E., Mclntyre, K.M., Leach, S., Baylis, M. and Morse, A.P. (2012) Suitability of European Climate for the Asian Tiger Mosquito Aedes albopictus: Recent Trends and Future Scenarios. Journal of the Royal Society Interface, 9, 2708-2717.

https://doi.org/10.1098/rsif.2012.0138

[11] Takahashi, L.T., Maidana, N.A., Ferreira Jr., W.C., Pulino, P. and Yang, H.M. (2005) Mathematical Models for the Aedes aegypti Dispersal Dynamics: Travelling Waves by Wing and Wind. Bulletin of Mathematical Biology, 67, 509-528.

https://doi.org/10.1016/j.bulm.2004.08.005

[12] Yang, H.M., Marcoris, M.L.G., Galvani, K.C., Andrighetti, M.T.M. and Wanderely, D.M.V. (2009) Assessing the Effects of Temperature on the Population of Aedes aegypti, the Vector of Dengue. Epidemiological Infection, 137, 1188-1202.

https://doi.org/10.1017/S0950268809002040

[13] Richter, O. and He, W. (2012) Modelling Large Scale Invasion of New Species under Temperature Change by Reaction-Diffusion Equations. International Environmental Modelling and Software Society (iEMSs) 2012 International Congress on Environmental Modelling and Software Managing Resources of a Limited Planet, Sixth Biennial Meeting, Leipzig, Germany, 1-5 July 2012, 2187-2193.

[14] Otero, M., Schweigmann, N. and Solari, H.G. (2008) A Stochastic Spatial Dynamical Model for Aedes aegypti. Bulletin of Mathematical Biology, 70, 1297-1325.

https://doi.org/10.1007/s11538-008-9300-y

[15] Delatte, H., Gimonneau, G., Triboire, A. and Fontenille, D. (2009) Influence of Temperature on Immature Development, Survival, Longevity, Fecundity, and Gonotrophic Cycles of Aedes albopictus, Vector of Chikungunya and Dengue in the Indian Ocean. Journal of Medical Entomology, 46, 33-41.

https://doi.org/10.1603/033.046.0105

[16] Maidana, A. and Yang, H.M. (2008) Describing the Geographic Spread of Dengue Disease by Traveling Waves. Mathematical Biosciences, 215, 64-77.

https://doi.org/10.1016/j.mbs.2008.05.008

[17] Anguelov, R., Dumont, Y. and Lubama, J. (2011) Mathematical Modeling of Sterile Insect Technology for Control of Anopheles Mosquito. AIP Conference Proceedings, 1404, 155-161.

https://doi.org/10.1063/1.3659915

[18] Dumont, Y., Chiroleu, F. and Domerg, C. (2008) On a Temporal Model for the Chikungungya Disease: Modeling Theory and Numeric. Mathematical Biosciences, 213, 80-91.

https://doi.org/10.1016/j.mbs.2008.02.008

[19] Richter, O., Moenickes, S. and Suhling, F. (2011) Modelling the Effect of Temperature on the Range Expansion of Species by Reaction-Diffusion Equations. Mathematical Biosciences, 235, 171-181.

https://doi.org/10.1016/j.mbs.2011.12.001

[20] Toma, L., Severini, F., Luca, M.D., Bella, A. and Romi, R. (2003) Seasonal Patterns of Oviposition and Egg Hatching Rate of Aedes albopictus in Rome. Journal of the American Mosquito Control Association, 19, 19-22.

[21] UCEA (2000) Osservazioni metereologiche dell’anno 2000. In: Mangianti, I.L. and Perini, L., Eds., Osservatorio meteorologico Torre Calandrelli, Collegio Romano, UCEA, Rome, Italy, 20 p.

[22] IPCC (2007) Climate Change 2007: The Physical Science Basis. In: Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., Averyt, K.B., Tignor, M. and Miller, H.L., Eds., Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK/New York, USA, 996 p.

[23] Hawley, W.A., Reiter, P., Copeland, R.S., Pumpuni, C.B. and Craig Jr., G.B. (1987) Aedes albopictus in North America: Probable Introduction in Used Tires from Northern Asia. Science, 236, 1114-1116.

https://doi.org/10.1126/science.3576225

[24] Brown, J.E., Scholte, E.J., Dik, M., Den Hartog, W., Beeuwkes, J. and Powell, J.R. (2011) Aedes aegypti Mosquitoes Imported into the Netherlands, 2010. Emerging Infectious Diseases, 17, 2335-2337.

https://doi.org/10.3201/eid1712.110992

[25] Medlock, J.M., Hansford, K.M., Schaffner, F., Versteirt, V., Hendrickx, G., Zeller, H. and Bortel, W.V. (2012) A Review of the Invasive Mosquitoes in Europe: Ecology, Public Health Risks, and Control Options. Vector-Borne and Zoonotic Diseases, 12, 435-447.

https://doi.org/10.1089/vbz.2011.0814

[26] Rozeboom, L.E. (1940) The overwintering of Aedes aegypti L in Stillwater, Oklahoma. Proceedings of the Oklahoma Academy of Science, 19, 81-82.