Most of the current biophysical models designed to address the large-scale distribution of malaria assume that transmission of the disease is independent of the vector involved. Another common assumption in these type of model is that the mortality rate of mosquitoes is constant over their life span and that their dispersion is negligible. Mosquito models are important in the prediction of malaria and hence there is a need for a realistic representation of the vectors involved.
We construct a biophysical model including two competing species, Anopheles gambiae s.s. and Anopheles arabiensis. Sensitivity analysis highlight the importance of relative humidity and mosquito size, the initial conditions and dispersion, and a rarely used parameter, the probability of finding blood. We also show that the assumption of exponential mortality of adult mosquitoes does not match the observed data, and suggest that an age dimension can overcome this problem.
This study highlights some of the assumptions commonly used when constructing mosquito-malaria models and presents a realistic model of An. gambiae s.s. and An. arabiensis and their interaction. This new mosquito model, OMaWa, can improve our understanding of the dynamics of these vectors, which in turn can be used to understand the dynamics of malaria.
Anopheles gambiae complexModelMalaria
This is the first of two papers describing a dynamic model (Open Malaria Warning; OMaWa) of Anopheles arabiensis and Anopheles gambiae s.s. Our aims in this article are 1) to formulate recent research on the Anopheles gambiae complex in a mathematical framework, and 2) to show how the new formulations influence the dynamics of malaria and mosquito populations.
In this paper, we describe a model of the dynamics of the two species and then show how parameters can influence the success of the two species, and how temperature, humidity and mosquito size can influence malaria transmission.
Climate and malaria
Most of the 149-274 million cases and 537,000-907,000 deaths from malaria occur in sub-Saharan Africa [1, 2]. Climate has been one of the main drivers of this disease , governing the spatial extent and year-to-year variations. The pathway from climate to malaria goes through the parasite and the mosquito. Although it is well established  how parasite development is influenced by temperature , the vector’s response to weather and climate is more complex. Mosquito density depends not only on temperature but also on the abundance of breeding sites (rainfall and evaporation) , desiccation (humidity) , and competition between mosquitoes . In the past 20 years, a shift in the distribution of An. arabiensis and An. gambiae s.s. has been observed in Kenya , showing that the species composition is not static over time. In the context of climate change , variability in vector populations is a factor that has not been considered so far.
Malaria and mosquito models
At the turn of the 20th century the work of several researchers, including Battista Grassi and Ronald Ross, resulted in the discovery that mosquitoes of the Anopheles genus transmit malaria [11, 12]. Over the next 20 years, Ross, and later Lotka and Waite, developed mathematical models that became central in malaria control [13–19]. In the 1950s, George MacDonald refined these models and showed that DDT could be used to interrupt malaria transmission . Since then, several modelers have followed in the footprints of Ross, Lotka, and MacDonald [21–30]. Some have designed models to show how temperature alone influences malaria transmission , while others have focused on the theoretical effect of bed nets , multiple interventions  or climate change [34–36]. There is also a growing number of models that address the dynamics of immunity within individuals  and in communities [21, 38].
In 2011, The malERA Consultative Group on Modeling  provided a review of the current state of mathematical models and pointed to the importance of good mosquito models for assessing the impact of climate change on malaria.
Many traditional models rely on a threshold principle. The idea has been to find thresholds for longevity, number of bites or days to recovery that must be reduced to interrupt the transmission. With increased computational power it is now possible to make more complex models and hence explore a wider range for the dynamics of malaria and mosquito survival. By integrating the knowledge from simpler models into a complex system, it is possible to test if the assumptions are true over a wider geographical range. In addition, these complex models can make quantitative predictions about strategies for control .
Model summary and motivation
A model is mental copy that describes one possible representation of a system. We present an alternative formulation of the dynamics of An. gambiae s.s. and An. arabiensis. The model is a system of ordinary differential equations (ODEs) with three compartments: eggs, first to fourth instar larvae, and pupae; an age-structured formulation of adult mosquitoes; and size prediction for adult mosquitoes (measured as wing length in mm). This can be considered the skeleton of the model. As demonstrated later, the model structure can be simplified when mosquito size can be neglected or when we assume no births. The model can be run with a spatial structure in which we include or exclude mosquito dispersion, or as an idealized model in which the model is evaluated at a single point.
The ODEs parametrize daily mortality rates, which are size-dependent for adult mosquitoes; development rates in the aquatic stages; biting rates; fecundity; the probability of finding a blood meal; and mortality related to flushing of eggs, larva and pupa out of oviposition sites. These parametrization schemes are driven by air temperature, relative humidity, relative soil moisture, water temperature, and runoff. As already mentioned, the model can be applied in a spatial domain. In this case, temperature and other environmental data are taken from a regional climate model, the Weather Research and Forecasting Model (WRF) . In the examples shown later, we run the model at a resolution of approximately 50 km and a temporal resolution of 5-20 years in steps of 3 h. In addition to weather data, human  and cattle [43, 44] densities are introduced to estimate the probability of feeding.
At this spatial resolution, the model should potentially be able to define larger foci of mosquito productivity, while the ability to identify hotspots will be limited . However, 50 km is the standard for regional climate models addressing long-term changes in climate . In addition, the true accuracy of historical cattle and human population density estimates for Africa in general is not likely to be greater than 50 km.
The mosquito model described here is designed to capture the spatial distribution and the time-dependent density of An. gambiae s.s. and An. arabiensis. If the model can capture the current distribution and density of the two species and how they are related to malaria, a future version of this model, including infections, could be used to explore the long-term impact of current interventions under a changing climate. To have confidence that the model has these abilities, several aspects not considered here should be evaluated (papers under preparation). In addition, if malaria modelers move towards the ensemble thinking widely adopted in the climate community, this model could be one representation of historical and future changes for malaria. The aim of such an ensemble would be to deal with uncertainties in the system. Ultimately, the goal would be to produce policy-relevant information including uncertainty.
We have chosen to represent the non-exponential mortality of An. gambiae s.s. and An. arabiensis as observed in laboratory settings , semi-field conditions , and in the field . A common assumption is that in the field, mortality rates are constant with age because of predation . To date, few studies have confirmed this, while there is field-based evidence of age-dependent Ae. aegypti mortality , which has implications for malaria transmission . In the model, we also describe how mosquito size changes over the season. This might seem to be an overcomplication of the model. The motivation, however, is that we have observed substantial improvements for arid regions such the Sahel when we included mosquito size prediction. Fouet et al. reported that mosquito size is an important adaptation strategy in arid environments .
We do not claim that the additional complexity adds any value. Stating this before the model has been fully evaluated and compared to simpler models would be dangerous. The model is thus one possible way of describing the dynamics of An. gambiae s.s. and An. arabiensis. It is under continuous development, and we expect to add and alter components as new data become available.
To highlight some of the components that contribute to the dynamics of An. gambiae s.s. and An. arabiensis in the model, five sensitivity experiments focus on the effect of temperature, relative humidity and mosquito size on malaria transmission. We also show how An. gambiae s.s. and An. arabiensis respond to changes in the probability of finding blood, carrying capacity, initial conditions, and dispersion.
Material and methods: model description
Summary of the model
Figure 1 provides an overview of the model. In the following sections we present the ideas behind the model and its general structure, how a climate model is used to drive the mosquito model, and the parametrization schemes used in the model. It should be possible to read each part independently; for example, data from a climate model can be used to drive any malaria model; the parametrization scheme can be used in any malaria model; and the malaria model described here can be used with different parametrization schemes, with or without data from a climate model.
As mentioned above, the model comprises a system of ODEs for eggs, first to fourth instar larvae, and pupae; an age-structured formulation for adult mosquitoes; and size prediction for adult mosquitoes (measured as wing length in mm). The first limitation in the aquatic stage is the availability of ovipositing sites, which is parametrized in terms of relative soil moisture and the potential for puddle formation in a specific location. Once ovipositing sites have been formed, adult female mosquitoes are allowed to deposit eggs until the site is full, defined as the biomass relative to the carrying capacity for the location. To account for density-dependent mortality, first instar larvae can be preyed on by fourth instar larvae , and an extra density-dependent mortality term is added to account for prey-independent mortality . The numbers of eggs, larvae and pupae are reduced when the precipitation rate exceeds the infiltration rate. The larval density in the aquatic habitat influences the size of adult mosquitoes . We account for this by predicting mosquito size at emergence as a function of larval density. In addition to temperature and relative humidity , mosquito size influences the daily adult survival probability [7, 51, 54, 55] (, Aedes aegypti). We therefore describe an adult survival model that takes temperature, relative humidity and mosquito size into consideration. In addition, adult mortality and fecundity can increase if there are no or few sources of blood. This follows the idea that a mosquito living in an environment where much energy has to be used to find blood will do this at the cost of survival.
We adopt these general ideas for two species, An. gambiae s.s. and An. arabiensis. It should be noted that we have less confidence in the model for the An. gambiae s.s. M form, since aestivation (as documented by Lehmann et al.  and Adamou et al. ) is not included. In addition, there are some indications that the M form breeds in larger pools  and hence the puddle parametrization might have limited validity for this form.
In addition to time, the model can include two (three, since space is two-dimensional) additional dimensions, namely age and space. The space dimension allows dispersion of mosquitoes, meaning that (re)establishment through migration to areas that were previously free of An. gambiae s.l. is possible. The gradual invasion of Brazil by An. arabiensis in the 1930s  is one example of dispersion.
The ODEs were solved using the ODE solver lsoda [61–63]. The relative and absolute error tolerances were not modified from the original lsoda implementation (1e-6). The model can be run either as a spatial model (with or without mosquito dispersion) or evaluated at a single point at which movement is neglected. A detailed overview of the possible model parameters can be found in Table 1.
Differential equations for the aquatic compartment
The aquatic compartment consists of six stages: eggs (E), four larval stages (L1,L2,L3,L4), and pupae (P). Transitions between the different compartments can be expressed in terms of delayed equations. To simplify the solution and avoid numerical instabilities, we approximate the model as ODEs . Lunde et al. reported on the errors introduced by this approximation .
New eggs added to the population depend on the number of adult mosquitoes (m), the size of adult mosquitoes (msize), the inverse length of the gonotrophic cycle (G(T)), how much water is available (SMr, dimensionless) and the larval biomass already present in puddles (BL):
where represents potential new eggs from each age group, G(T) is either constant or dependent on temperature T, SMris a function of the relative soil moisture and the potential puddle formation area, K is the maximum larval biomass a grid cell can hold, βN,E(T) is natural mortality rate for eggs [Eqs. (16) and (18)], βI,Eis the induced mortality rate for eggs (not specified) and τEis the inverse of development time from eggs to first instar larvae.
The term 1-BL/K is used as a scaling factor to modify the growth rate. When the population is low compared to the breeding sites available, its growth is high. As the population grows, there is more competition for food, predators become more abundant, and the growth slows. In the egg compartment this represents the idea that the mosquitoes will lay fewer eggs when breeding sites are already occupied .
First instar larvae (L1) are added as eggs develop into larvae. Additional mortality is added in the transition stage in relation to how much biomass there already is in a given location . This approximation of increased (density-dependent) mortality arises because of competition and predators; if a puddle already is full, the number of eggs developing to first instar larvae is reduced, whereas if a puddle is empty (1-BL/K = 1), no extra mortality occurs. Similar terms could have been added to the second, third and fourth instar larvae, but we assume that earlier life stages will be affected more by density-dependent competition and predation.
Shoukry looked at how fourth instar larvae of An. pharoensis prey on first instar larvae during a 24-h experiment . Using these data, we add additional mortality for first instar larvae according to the density of fourth over first instar larvae. The constant Cpredis tunable to both limit the predation on L1 and make it more specific to species in the future. At most temperatures, this constant does not influence the density of mosquitoes (Additional file 1).
The number of first instar larva is given by:
Second (L2), third (L3) and fourth instar larvae (L4) and pupae (P) are controlled by the development rate τ and mortality β:
where β is the daily mortality rate, with the first subscript denoting natural (N) or induced (I) mortality and the second subscript denoting the aquatic stage. The subscript for the development rate, τ, corresponds to the aquatic stage. The parametrization schemes and data sources used to estimate the rate at which eggs are laid (G(T) and ε), mortality (β) and the development rate (τ) are discussed later.
Differential equations for adult mosquitoes
The life history and mortality rate vary over the lifespan of a mosquito population. We formulated a model to account for this variation. Adult mosquitoes are denoted by mn, where n indicates the age group; n = 1 is the youngest group and n = 9 refers to the oldest mosquitoes. The age groups in the model are m1 =[0,1], m2 = (2,4], m3 = (5,8], m4 = (9,13], m5 = (14,19], m6 = (20,26], m7 = (27,34], m8 = (35,43] and m9 = (44,∞] days, with ageing coefficients anof 1.000, 0.500, 0.333, 0.250, 0.200, 0.167, 0.143, 0.125 and 0.067 for , respectively. Mosquito ageing is represented by Ψn, where n denotes the age group. Ageing is time-invariant and is thus not related to the number of gonotrophic cycles.
Although there is no ageing from age group 9, the term Ψ9 is included to limit the concentration of old mosquitoes. This is a user-specified variable and in the model results shown here we set this to for An. arabiensis and An. gambiae s.s.; this value should be set to ensure that mosquito populations can survive during dry periods [66, 67], but still hinder accumulation of old mosquitoes. This can be particularly useful if the mortality model described later is replaced with a model in which mortality is independent of age.
When m is written with subscripts ı and ȷ in addition to n, this denotes inclusion of mosquitoes from neighboring areas. For example, subscript ı-1 indicates that mosquitoes to the west of the point of interest are interacting with the point of interest. The formulation presented here includes movement of mosquitoes, and where appropriate we denote mosquitoes by mn,ı,ȷ.
Again, β denotes mortality, with the first subscript denoting natural (N) or induced (I) mortality and the second subscript denoting the age group (mn) of the mosquitoes. Φ represents the mosquito flux (transport) and subscripts ı and ȷ define which boundaries are evaluated. This is discussed in the section “Movement of mosquitoes”.
The number of adult mosquitoes of a specific age in a grid point is controlled by new mosquitoes from mn-1, as well as the flux to and from the point of interest , natural mortality , induced mortality , ageing to mn+1, and mortality due to lack of food (P(B)). Parametrization schemes related to mortality are discussed later.
This results in the following equation for the first age group:
The equations for age groups n = [2,9] are
Differential equations predicting mosquito size
Mosquito size (msize) is important for the efficiency of mosquito multiplication. There are also some indications that increased body size is a strategy for survival in arid environments . In general, high larval density leads to a smaller body size as adults, and vice verse . Where only one species is competing for a resource, such as in a small puddle, mosquito size, and hence the number of eggs laid by each mosquito, will be of less importance. If two species are competing for the same resource (e.g. An. arabiensis and An. gambiae s.s.), the trade off between development time and size can be important in competition for breeding sites. An. gambiae s.s. generally develop faster than An. arabiensis, but end up with a smaller body size. An. arabiensis spends more time in the aquatic stages and develops larger bodies, and can thus produce more eggs. Since our model includes competition between those species, we describe mosquito size as a function of competition for breeding sites. In theory this should improve our ability to separate geographical and seasonal distributions of An. arabiensis and An. gambiae s.s.
Since the size of An. arabiensis and An. gambiae s.s. stabilizes after approximately 4 days  and ovoposition does not start before this, it is not necessary to differentiate the maximum and minimum size depending on age to mimic changes in the number of eggs per mosquito with age. However, this may be required if mortality based on desiccation [7, 69] is used. Although mosquito size at a given time can be approximated using finite differences, we develop a different approach that is more efficient in terms of computational time in our model framework. Mosquito size for the first age group depends on larval size. Since the pupation time is short, this assumption is justified, although it might introduce minor errors. In a future version of the model, we plan to predict larval size dynamically. The limitations set on mosquito size (described in “Parametrization schemes in the aquatic stages”) in this model might lead to An. arabiensis that are slightly too small compared size in the field study of Ye-Ebiyo et al. , but the size is in line with studies by Huestis et al.  and Kirby et al. . Kirby et al. also noted that mixed populations of An. arabiensis and An. gambiae s.s. had a negative effect on mosquito size at some temperatures. This mechanism is not included in the current work. However, the most important aspect of modelling of mosquito size is to capture seasonal and spatial variations.
For size prediction we use the symbol , where n is the age group as described above.
The size (wing length in mm)of newly emerged mosquitoes is approximated according to the linear relationship
where larva size Lsize(in mg) is approximated as:
The constants asppand bsppare 0.45 and 0.12 for An. arabiensis and 0.383 and 0.147 for An. gambiae s.s., respectively .
The size of mosquitoes in the first age group at any time is given by
Therefore, the size of newly emerged mosquitoes () depends on the number of newly emerged pupae and the relative density of larva at the breeding site.
For the remaining age groups, size is estimated as
Therefore, the size in age groups 2-9 only depends on the number of mosquitoes surviving from one age group to the next (mn-1) and the size of mosquitoes in the younger age group ().
To drive a dynamic malaria model it is necessary to have boundary conditions that are consistent over time and space. Temperature, relative humidity, and rainfall data from weather stations are point measures. Hence, they might not be representative of larger areas over shorter time scales. This is especially true in areas with varying topography or where convective rainfall is dominant [73–75]. Despite the limitations of rainfall stations, they can provide a robust estimate of large-scale events. By pooling data from several stations, the error for a single station is reduced and the data can provide a good estimate for dry and wet years, for example. Hence, weather stations are useful tools for validating climate models.
The problems of point measurements are described later, and represent one of the reasons why OMaWa is tightly linked to a climate model. As shown in sensitivity experiments, the model can also be run with constant forcing (e.g. temperature) or with data from weather stations.
Where we present results for Africa as a whole, OMaWa is driven by data from WRF 3.3.1. This realization (TC50), described in part two of this paper, has a tropical channel set-up in which set-up, the domain consists of boundaries above and below a certain latitude and no side boundaries. The model was run at 50-km resolution from January 1, 1989 to January 1, 2009. At the northern (45°N) and southern (-45°N) boundaries the model was driven by Era Interim. The Kain Frisch cumulus parametrization scheme was used [76, 77]. This experiment was not designed to reproduce observed year-to-year weather variability, but to assess the mean mosquito density and distribution. The driving experiment is described in the section on model validation.
Climate and weather models
Currently, our best guess of (future) climate at multidecadal time scales comes from general circulation models (GCMs). These models are designed to close the energy budget of the Earth and include an interactive representation of the atmosphere, ocean, land, and sea ice. A set of scenarios with different emissions describes how sensitive the climate is to atmospheric constituents (greenhouse gasses) . While climate is the average weather over time and space, weather can change over minutes, hours, days and seasons. The same equations used to predict climate are used to predict weather. However, weather forecasts are more dependent on current observations of the atmosphere. Hence, weather predictions are initial value problems, whereas climate simulations are rather boundary value problems.
Both climate and weather models are mostly structured on a grid, with coordinates from west to east (x), north to south (y) and bottom to top (z). In the grid, one square (or polygon) represents the weather within that square. While climate models often have a horizontal resolution of more than 10000 km2, operational weather models such as the European Centre for Medium-Range Weather Forecast (ECMWF) model are run at approximately 160 km2. If the state of the atmosphere is observed correctly, higher resolution can lead to better local skill in predicting the weather. A hybrid between a weather model and a climate model is a limited-area model (LAM), which relies on initial and boundary conditions from a weather or climate model. Given these conditions (weather), the LAM can be run at a higher resolution over a limited area, which potentially improves the spatial accuracy of the coarse model . The WRF model is a widely used LAM .
In tropical regions, most rainfall comes from convective clouds. This type of rainfall is generally intense and of short duration. The geographical extent of such rainfall episodes may be limited. Therefore, rainfall measurements in regions where convective rainfall is dominant should be handled with care [74, 75, 80, 81], especially when extrapolating station data to areas with no data. While station data are accurate at a specific point, climate models and satellite estimates give a more general description of the weather within a certain area; Chen and Knutson reviewed how models compare to observations at varying scales . Since future climate is projected using climate models and considering the limitations of weather stations, construction of a mosquito/malaria model around a LAM is a good choice. The LAM will have higher resolution than most climate models, with higher-resolution orography, coastlines, and land use, but will still give a general description of the weather within a certain area.
Parametrization schemes in the aquatic stages
To relate a variable such as mortality to the physical environment, we need simplified equations that describe this relationship. An equation in which temperature influences mortality only states that there is a relationship between the two, but does not explain why temperature modifies mortality. In this paper we use parametrization schemes to represent the influence of the environment on mosquitoes. This section describes the aquatic parametrization schemes used, excluding water availability, which is discussed later.
The aquatic stages comprise eggs, four instar stages, and pupae. The number of eggs in a location at any time is controlled by the number of potential new eggs laid (ε), available water (K), natural and induced mortality (βN/I,E/L/P) and movement from the E to the L1 compartment. In addition, 20% instant mortality is introduced when rainfall exceeds the infiltration rate. This is in line with observations by Paaijmans et al. . The number of new eggs is simplified to a function of the number of gravid mosquitoes in each age group and their size (measured as wing length) based on observations [55, 84–86]. The critical size is set to a wing length of 2.6 mm, which is less than that observed by Lyimo and Takken  but greater than observations by Yaro et al. . Maximum wing length is set to 3.3 mm for An. gambiae s.s.[88, 89] and 3.7 mm for An. arabiensis. The relationship between the number of eggs (ε) and wing length is then approximated according to the linear relationship
where mn is the number of mosquitoes in age group n. Note that this limits the number of eggs laid by a single mosquito per gonotrophic cycle to approximately 184, which is somewhat less than the number observed by Yaro et al. , but in line with that reported by Howard et al. .
Estimation of water temperature
Using the 0-10-cm soil temperature (Tsoil) from the NOAH land surface model [91–94] to approximate the mean water temperature (Twater) in larval habitats, we assume that evaporative cooling and heat fluxes at the water boundaries are negligible. Hence, the water temperature is equal to the top soil temperature. Paaijmans et al. showed that the 5-cm soil temperature represents the water temperature in small ponds reasonably well . Therefore, the model will have limited validity in areas where larger puddles are the main breeding sites. There is also a chance that diurnal fluctuations will be slightly over- or underestimated. When a grid cell covers several km2, this effect should be negligible, although we do not have data to support this. We hope to improve the prediction of water temperature in the future, either by modelling this explicitly or using a parametrized version based on data from Huang et al. .
Parametrization of mortality
We used two approaches to calculate mortality in the aquatic stages. In the simpler approach, we assume that mortality and development time in the aquatic stages are independent of the species. We also assume that the relationship between the mortality rate and temperature is the same for eggs, instars and pupae. In this method we do not consider competition effects as described by Paaijmans et al. . This type of parametrization is suitable when the model is used for one species only (e.g. if the model represents an area where only one of the two species is present).
Species-independent mortality (BLL)
Data provided by Bayoh and Lindsay  were used to describe the mortality rate according to Eq. (14) (p < 0.01, R2 = 0.81). We call this the BLL method. Mortality rate data are plotted in Figure 2b.
where βN,L(Twater) = βN,E(Twater) = βN,P(Twater) is the aquatic mortality rate per day and Twateris the water temperature (°C). The constants knare given in Table 2.
Constants for equation 14 and 33
Species-dependent mortality (KBLL)
Kirby et al. reported that the mortality rate of An. gambiae s.s. and An. arabiensis is modulated by the presence of each other in the temperature range 25-35°C . To account for this we developed two mortality models, one for An. gambiae s.s. and one for An. arabiensis. We call this parametrization scheme KBLL. The mortality rates are based on data from Bayoh and Lindsay  and from Kirby et al. . Although Holstein also reported larval mortality for (An. gambiae s.s.) when exposed to extreme low and high temperatures , we did not include these data when estimating the mortality curves. However, the data are plotted in Figure 2 for comparison. According to our curves, the An. arabiensis mortality rate will increase in the range 25-35°C as the relative presence of An. gambiae s.s. increases. Conversely, the mortality rate of An. gambiae s.s. will decrease as the proportion of An. arabiensis increases. The mortality rate βN,Lis given by
fgamband farabare the ratio of An. gambiae s.s. to An. arabiensis larvae and An. arabiensis to An. gambiae s.s. larvae, respectively. At each time step, Lsizeis estimated as a function of BLand K. As the density increases, there will be more competition and hence less food for each larva, which leads to smaller larvae.
Parametrization of the development rate
The rate of development between the different aquatic stages follows the corrected version of Bayoh and Lindsay . Since these data are only valid for An. gambiae s.s., we made a small modification to prolong the development times for An. arabiensis. Data from Kirby et al.  and Paaijmans et al.  suggest that time for development from a larva to an adult is approximately 5.5% longer for An. arabiensis than for An. gambiae s.s. Hence, we increased the development time for An. arabiensis by 5.5%. The reason for this longer development time is that An. arabiensis takes longer to develop a larger body. Curves of the development rate are shown in Figure 3.
The two previous studies also suggest that the development rate  and mortality  of the two species are modulated by the presence of each other, so we take account of this in out model. The development time for An. arabiensis is prolonged in the presence of An. gambiae s.s., while the time is shortened for An. gambiae s.s. as the relative proportion of An. arabiensis increases. Using data from Paaijmans et al. , the development rate τ is modified according to
for An. gambiae s.s. and
for An. arabiensis. farab and fgambis the fraction of An. arabiensis and An. gambiae s.s., respectively.
Parametrization of breeding sites
The formation of puddles can be described as a balance of runoff, infiltration, evaporation, and rainfall entering the puddle. The formulation of an idealized puddle can be found in Additional file 2.
Modelling of every single breeding site requires high enough resolution to resolve the puddle. In practice this is not possible and the problem has to be simplified.
Mushinzimana et al. described typical breeding sites in a Kenyan highland area . Most of the puddles were located at less than 100 m from rivers, which means we can assume that semi-permanent puddles will mostly form in the proximity of rivers and lakes. They also found that the number of breeding sites was close to threefold higher in the rainy season compared to the dry season, and grouped breeding sites by surface area.
If we assume that breeding mainly occurs in the vicinity of potential rivers and lakes, the availability of breeding sites can be expressed as a function of potential river length and soil saturation. At high resolution this might not always be true , but since the model is designed to be applied to coarser grids, we believe the assumption is as reasonable as or more reasonable than the common assumption that puddle formation is only dependent on rainfall . The newest version of the NOAH land surface model in WRF 3.4 also includes groundwater and dynamic vegetation, and future versions might change the way in which puddles are parametrized. In OMaWa we introduce a simple parametrization scheme to represent breeding sites.
The Hydrological Data and Maps based on SHuttle Elevation Derivatives at Multiple Scales (HydroSHEDS) 15s river data set from the US Geological Survey (USGS)  was used to derive the total potential river length within a grid cell. Since the algorithm used to develop this data set describes where water would collect if it were available within the catchment, it also represents a general description of the potential for water aggregation within an area. However, the validity might decrease on moving to finer scales .
Here we divide rivers into three different classes: perennial, intermittent and ephemeral streams. For each class, potential river length (Rp, km) within a grid is defined as
where Ξ is the equally spaced river data-set resolution in degrees, where Δlon = Δlat, ER is the radius of the Earth (6371.22 km) and φ is latitude in radians.
In a simplified model we estimate puddle volume as a function of river length and relative soil moisture. Although this is a very crude estimate, we compared this simple model with data from Mushinzimana et al.  and derived a simple expression for the carrying capacity in a grid cell:
where is the maximum larval biomass per km of river (2400 mg, estimated from data collected by Munga et al. ) and SMris the relative soil moisture content (fraction).
In the current implementation we do not distinguish between fast- and slow-flowing rivers. It should be noted that this way of approximating breeding sites has limited validity in areas with irrigation or around rivers where breeding sites could form as rivers recede [66, 67, 102]. Some special cases, such as along the River Nile in Sudan, where breeding sites form as a result of rainfall hundreds of kilometers away, will not be captured at all .
Parametrization of the gonotrophic cycle
The gonotrophic cycle depends on temperature and is important for the vectorial capacity of mosquitoes. Lardeux et al. studied the gonotrophic cycle for An. pseudopunctipennis. We combine their data with other published studies on anophelines to estimate the length of the gonotropic cycle. There are few studies on An. gambiae s.l., and hence we have to assume that other anophelines share the same physiology and strategy with respect to the gonotropic cycle. Ruiz et al. showed there are some differences , but until further evidence of the reproductive strategies of different members of Anopheles genera becomes available, we will not consider this effect. Studies used to develop the formula include those by Guillermo et al. (, An. albimanus), Afrane et al. (, An. gambiae s.l.), and Maharaj (, An. arabiensis). We also include the formula given by Hoshen and Morse . Their model is based on degree days and is included according to Eq. (26). The gonotropic rate (day-1) and data used to develop the formula are shown in Figure 4.
where Tairis the air temperature (°C), Ddis degree days, and Tcis the critical temperature from Hoshen and Morse , with Dd= 37, and Tc= 7.7.
Parametrization of the age-dependent mortality of adult mosquitoes
The mortality of adult anophelines differs according to age and species [7, 107, 109]. This has often been overlooked in mosquito models [23, 110]. To show how this assumption can influence the stability of mosquito populations and malaria transmission, we use the mortality model of Martens  as a reference. We also plot Eq. 7 from Ermert  in Figure 5 to highlight the differences between this model and established models. For convenience, we repeated Marten’s equation, as follows:
Our new survival curves are based on unpublished data from Bayoh and Lindsay . The validity ranges from 5 to 40°C by 5°C and 40-100% by 20% relative humidity. We name the scheme BLLad (Bayoh-Lindsay-Lunde adult mortality). The data set and the curves are valid for An. gambiae s.s. The lowest agreement between the model and the data is at 40% relative humidity and 40°C. While the data suggest that all An. gambiae s.s. would be dead after approximately 2 days, the survival curve would result in no mosquitoes after approximately 4 days at 40% relative humidity and 40°C. To correct for this error, we include data from Kirby and Lindsay , who described the responses of An. gambiae s.s. and An. arabiensis to high temperatures. By assuming that maximum survival is 480 min for An. gambiae s.s. and 1440 min for An. arabiensis at temperatures greater than 40°C, we can set the mortality rate to 3day-1 and 1day-1, independent of age group. However, there are uncertainties at relative humidity below 40%. The lack of studies in this range is a limitation of this survival model, and could make the model less accurate for An. gambiae s.l. in some regions. The basic principle of these survival curves is that mortality will be low in the first few days after emergence. In addition, mosquitoes that survive up to a certain age have a higher survival probability (depending on Tairand relative humidity). In Figure 5, survival at 60% relative humidity and 0, 10, 20, 30, and 40°C is plotted.
Size affects the survival of adult mosquitoes [7, 51, 54, 55] (, Aedes aegypti). If we assume that the major differences in mortality between An. gambiae s.s. and An. arabiensis can be attributed to mosquito size, we can modify α as a linear function of mosquito size. Here we subjectively choose reasonable constants for h(msize). Tairmay be completely or partly replaced by indoor temperature (Tindoor, described later), depending on the proportion of mosquitoes indoors. In experiments covering the African domain, we assumed that 80% of An. gambiae s.s. and 20% of An. arabiensis are located indoors.
where ζ = 6, g is a function of mosquito size, and RH is relative humidity. The mortality rate for each age interval can then be approximated as
If we assume that differences in adult mortality for An. gambiae s.s. and An. arabiensis can be explained by differences in body size, these BLLad curves can be used for both species. We explore this mortality model in .
AL adult mortality
A similar approach can be used for An. arabiensis. Using survival curves reported by Afrane et al. (, Figure two) (copyedited with g3data ), we can estimate mortality based on the daily maximum temperature. Because of the few data points, this approach is much more uncertain and should be considered experimental. The advantage of this mortality model is that the data are not estimated from a laboratory setting. The maximum temperature reflects some aspects, such as radiation, albedo, and humidity, of the environment in which mosquitoes live. In some of the results presented in part two, we use this model for adult survival.
Constants c1,…,11 are listed in Table 2. By setting ζ = 2 we can simplify the survival curve for An. arabiensis to
Paaijmans et al. discussed the importance of using indoor rather than outdoor temperature, to describe the environment for mosquitoes and parasites . They included two studies that showed the relationship between indoor and outdoor temperature in Kenya  and Tanzania . Here we add two additional studies, one from Kenya  and one describing the temperature in traditional and low-cost modern housing in the Eastern Cape, South Africa . The data used to parametrize equation 36 came from; 1, Afrane et al. ; 2, Makaka and Meyer ; and 3, Paaijmans et al. [114–116] (R2 = 0.89). It is clear that temperatures inside a house are more stable than outdoor temperatures. House type greatly influences daily temperature fluctuations [117, 118], and the model used here might not be valid for all house types. While some studies have assumed that houses are always hotter than the surroundings , we approximate the indoor temperature as
Since the data are based on maximum and minimum temperatures, the timing of the indoor temperature might be offset by a couple of hours. This is evident in a study by Makaka and Meyer , who delayed the maximum indoor temperature by a couple of hours compared to the environmental temperature. At present we do not account for this delay, since the diurnal temperature ranges will be correct even if we do not. The data and regression line are shown in Figure 7. Further studies on indoor compared to outdoor temperatures are needed to make this correction more accurate.
Hence, Taircan be partly or fully replaced by Tindoor, depending on the proportion of mosquitoes indoors.
It should be noted that we still do not include temperatures in resting places described by Holstein, such as holes in rocks and cracks in soil, covered pigsties, rabbit hutches, hen coops and dry wells , and by de Meillon (, under stones).
Approximation of mosquito movement
The role of diffusion and advection in vector borne diseases have been explored in several papers [102, 121–127]. Considering the gradual invasion of Brazil in the 1930s by An. arabiensisit can be argued that movement of mosquitoes is important over decades. Here we include the active and passive transport of mosquitoes as fluxes across grid boundaries. Passive transport is movement of mosquitoes caused by wind, while active transport is movement due to flying. On shorter time scales the role of such movement will be limited. However, on long time scales it is necessary to allow mosquitoes to travel to allow them to establish in new locations.
Transport of mosquitoes is defined by fluxes (s-1) at the grid boundaries. In the model we allow fluxes from the eight neighboring grid points. A special case is implemented when a neighbouring cell is water. In this case, fluxes to water are reduced to 0.1% of the original flux to avoid large losses of mosquitoes along the coastline. Given strong winds from land to the ocean, such an assumption could lead to accumulation of mosquitoes along the coast. Conversely, allowing free movement to the ocean could lead to undesired loss of mosquitoes.
Since the movement of mosquitoes has a high computational cost, the spatial fluxes do not change the size calculations. This will introduce some minor errors when the movement of mosquitoes is low compared to their density, with larger errors if many mosquitoes are moved relative to their density. When a cell free of mosquitoes is colonized, the size is set to 3.05 mm.
The possible flight range of anophelines varies with food availability . We do not include vegetation types in the model and hence it is hard to justify differences in flight performance based on, for example, land use. The dispersion coefficient describes how far mosquitoes can move in a day. We assume that the dispersion coefficient D is constant, independent of geographical location. For An. gambiae s.s. and An. arabiensis, real flight performance outside the laboratory of only a few hundred meters per day (approx. 300-700 m) has been reported [102, 129, 130]. In this experiment we subjectively chose D=30mday-1 independent of age group. Anophelinae also travel with humans , which adds to the transport equation and makes the dispersion coefficient uncertain. Gillies noted that wind direction mostly has a minor effect on dispersal , while de Meillon  and Adams  reported distances of 2-4.5 miles (3-7 km) in the direction of the prevailing wind. Thus, it cannot be ruled out that wind plays a role on longer time scales. Hence, we express movement caused by wind as a function of 10-m zonal (u) and meridional (v) wind components (ms-1). This can be understood by considering the following example. For a constant u-wind of 10ms-1 and v-wind set to 0, mosquitoes will be moved a distance related to a scale factor Sf, which is equal to the distance travelled at 20ms-1 to the east. For example, with Sf= 750mday-1, the eastward distance traveled will be m in 1 day, but since each mosquito is not modelled individually, it would be more natural to describe this as a fraction moving a certain distance. Different wind directions and speeds will result in other distances/fractions and directions. D and Sfare unknown tunable constants.
Since the species considered here are most active at night , movement will be suppressed between 06:00 and 18:00 h (local time) and amplified at night according to
where LT is local time, and
Transport of mosquitoes and mosquito sizes inside and outside a grid are defined by
More specifically, during a time Δt, movement can be calculated as follows. On a day with no wind, transport is equal in all directions, D = 30mday-1, and the flux at a boundary is defined as
and transport ηı,ȷnis then equal to
In the presence of wind, we obtain additional transport as a function of zonal and meridional wind components.
Mortality related to feeding
One factor that is often overlooked in malaria (mosquito) models is survival related to food availability (P(B)). Ye-Ebiyo et al. reported that maize pollen availability has a positive effect on larval (and hence mosquito) fitness [70, 134]. Creating maps of plant types is beyond the scope of this study, and hence we chose not to account for mortality related to crops. However, we performed initial tests in which we included GlobCover Land Cover version V2.2 (European Space Agency ) to give a rough estimate of regions where increased fitness could be expected. The other source of food for female anophelines is blood. Compared to a starved mosquito, a mosquito that has had access to blood on days 1-3 has a theoretical flight distance that is increased by a factor of 6-7 . Therefore, it is plausible that the higher (lower) the probability of finding a blood meal (P(B)), the higher (lower) is survival in the early life stages of adult mosquitoes. Bouma and Rowland reported higher parasite prevalence among children of families who kept cattle compared to those who did not , which can indicate either higher survival (older mosquitoes) or simply that some anophelines are attracted to cattle. If we assume that a newly emerged mosquito has a flight range of Frm= 0.5km2day-1, the daily probability of finding a blood meal can be calculated as
where ρhumanand ρbovineis the probability of finding a human and bovine source, respectively. ρhumansis defined as the human population density per km2 multiplied by 0.1 (since a smaller area on a human is accessible) and ρbovineis defined as the bovine density per km2, each with a user-defined threshold at which the density is so low that P(B) is virtually zero. Since P(B) is a conceptual parameter, it can be tuned.
Since blood meals, besides sugar meals, are important for the mobility  and survival of female anophelines , the success of a species is likely to be linked to the presence of the preferred host. The dominant blood source for An. arabiensis is bovine and human blood, while it is human blood for An. gambiae s.s. . In reality there are strong indications that the human blood index is a dynamic quantity rather than a constant [139–142]. In the current implementation, HBI is a static number and hence there are probably errors related to this term. To find the probability of feeding on humans at each time step, we combine two data sets. Between 2000 and 2010 we use population densities from the Gridded Population of the World (GPW) , and for before 2000 and after 2010 we use growth rates from the Population Division of the Department of Economic, and Social Affairs of the United Nations Secretariat . Since there are no projections of cattle densities, this quantity is time-invariant and based on Food and Agriculture Organization (FAO) 2005 estimates . We are currently working to include time-varying cattle densities.
In the model, mortality caused by food limitations is a function of how many humans or cattle are available per mosquito and the human blood index. We assume that HBI is time- and space-invariant, and only depends on the species. For simplicity we chose available humans to be humans who are not sleeping under a bed net. In the simulations presented here, we set bed net usage to zero, and hence the results represent mosquito distribution without interventions. Bayoh et al. hypothesized that the survival of the different species is related to the availability of the preferred host . The daily mortality rate caused by limited human blood is expressed as
The functional form of of equation 42 can be seen in Additional file 3.
Figure 8 shows the probability of finding a blood meal for the sibling species on January 1, 1999.
Results and discussion
Sensitivity experiments are useful in understanding which parameters are important for the success of An. arabiensis and An. gambiae s.s. and which are important for malaria transmission. Classical sensitivity analysis investigates the robustness of a study when parameters are estimated from statistical modelling. Our model uses parametrization schemes to represent the influence of the environment on the two species. We show how the model responds to changing temperature, humidity, mosquito size, dispersion and the probability of finding blood. This approach does not allow us to directly measure the robustness of each parametrization scheme, but gives us an insight into which external factors influence the model and where it is of importance to have improved parametrization schemes. We use the term sensitivity experiments for this analysis.
To demonstrate some of the capabilities of the model, we set up a series of experiments. Some aspects are best visualized as a one-dimensional model (time and age), while other features are shown using a spatial domain (time, age, and space). For the one-dimensional experiments, the water temperature is set to the air temperature, except for temperature greater than 33°C, for which we set temperature to 33°C. This modification is required since pupae and fourth instar larvae will not develop below 18°C or above 34°C . The results are therefore less robust when temperature is greater than 33°C. Unless otherwise stated, we use size-dependent mortality, correction for indoor temperature, the KBLL method to estimate mortality in the aquatic stages, correction for the development rate in the aquatic stages depending on the ratio of each species, and movement of mosquitoes (in the spatial cases).
Sensitivity to temperature, relative humidity and mosquito size (TempHumSize)
The age-dependent mortality is influenced by temperature, relative humidity and mosquito size [Eq. (32)]. This experiment explores how the dynamics of malaria is sensitive to temperature, relative humidity and mosquito size (measured as mm). We assume that no births occur to isolate the effect of the transmission process, and consequently constant mosquito body size in the course of integration, but include mortality and the biting rate. In this experiment we assume that only one species is present (since the main competition occurs in the aquatic stages). This sensitivity test is designed to observe how the proportion of mosquitoes becomes infected as a function of temperature, relative humidity and mosquito size, given that we start with 1000 newly emerged mosquitoes, with m1 = 1000 and m2-9 = 0 as the initial conditions. In this experiment, 1% of the human population is infectious for Plasmodium falciparum. Mosquitoes are infected with an efficiency of 100%, meaning that biting an infectious human results in gametocyte transmission to the mosquito. In practice, this would be the same as saying that 10% of humans were infectious and gametocyte transmission had an efficiency of 10%. We also neglect the effect of heterogeneous biting. This is the only experiment in which we model the proportion of infectious mosquitoes explicitly. The modified equations describing the transmission process are described in .
The rate of sporozoite development within mosquitoes is expressed as 
where a = 9.5907, b = 0.0051029, c = 0.7349, and d = 17.0325. This expression was derived from the figure in MacDonald page 119  using g3data , and fitted using non-linear least-squares .
The gonotrophic cycle and biting rate are defined in Eq. (26).
The integrations are repeated with different combinations of temperature and relative humidity. This is a simple representation of gametocyte transmission to mosquitoes and is an idealized approach for exploring the proportion of mosquitoes (of the original 1000) that would become infected under different temperature, RH and mosquito size. Figure 9 shows how the percentage of infectious mosquitoes changes with temperature, RH and mosquito size. Lyimo and Koella reported that the largest mosquitoes were less likely to have sporozoites, but had more oocysts than smaller mosquitoes . They attributed this to increased mortality in the presence of many oocysts, an effect that is not included in our model. Figure 9 shows that the potential percentage of infected mosquitoes is sensitive to all three parameters in the model. Although higher survival has been attributed to body size in dry [7, 51] and semi-arid environments , the advantage or disadvantage of a larger body has been poorly described in saturated environments. Therefore, the sensitivity to body size at 80% RH should be interpreted with care. According to the model, temperature is not the only factor that governs the transmission of malaria (in areas with no interventions); humidity and how mosquitoes adapt to dehydration stress are also important factors. The most efficient transmission, expressed as the integral, with respect to days, occurs at 25°C at 40% and 80% RH, and at 24.5°C at 10% RH, independent of mosquito size.
These results should be viewed in light of recent findings by Paaijmans et al. that optimal transmission occurs at lower temperatures .
Sensitivity to temperature and carrying capacity (TempCar)
The aim of this sensitivity test was to investigate how carrying capacity and temperature determine the relative proportion of An. arabiensis and An. gambiae s.s. We set the relative humidity to 80% and the probability of getting a blood meal to one. We assumed that the soil was saturated and we varied the temperature between 16 and 38°C (with corrections over 33°C for water temperature) and the carrying capacity between 0.0625 and 125 mgkm-2.
Carrying capacity in the aquatic stages influences larval growth and adult survival. While An. arabiensis invests more time in growth than An. gambiae s.s., the former develops a larger body, and consequently has the potential to oviposit more eggs than the latter. If the two species experience the same mortality rate in the aquatic stages, more An. gambiae s.s. will emerge, but over time An. arabiensis can face this challenge by outnumbering the eggs of An. gambiae s.s. in the habitat. Thus, we are interested in testing how the carrying capacity in the aquatic stages alters the relative proportion of each of the adult species. In this model we only consider the competition between these two species, and hence neglect other competing species .
As observed in Figure 10, An. gambiae s.s. dominates between 27 and 30°C. This is the effect of the development rate modifications described by Kirby et al.  and Paaijmans et al.  (Figure 2 and “Species-dependent mortality (KBLL)”). Interestingly, the dominance of An. arabiensis is most pronounced in the drier simulations, meaning that high competition, compared to adult survival, is favourable for this species. This can be attributed to the strategy of larger body size and higher egg production. Lehmann et al. found that An. arabiensis dominated during the dry season, while An. gambiae s.s. dominated in the rainy season . The advantage of An. arabiensis in crowded breeding places might be one factor contributing to the shift in species composition as the surface area of puddles starts to shrink.
Sensitivity to temperature and the probability of finding blood (pBlood1D)
This experiment shows how the model responds to changes in the probability of finding a blood meal, which influences the rate at which mosquitoes can oviposit and increases energy consumption if hosts are hard to locate. If, for example, cattle are easier to find compared to humans, An. arabiensis will potentially use less energy per batch of eggs and will also be able to utilize breeding sites at a higher rate than An. gambiae s.s. It is also possible that An. arabiensis uses cattle for navigation . Over time, such differences might lead to dominance by one species. In this experiment, we varied the probability of finding blood, P(B), for An. arabiensis from zero to one, as well as varying the temperature as described for TempCar.
We set the probability of finding blood to one for An. gambiae s.s., independent of the probability of An. arabiensis finding a blood meal. This is a purely theoretical experiment designed to demonstrate a concept. The probability of finding blood is varied between zero and one for An. arabiensis. The scenario in which P(B) = 1 for An. gambiae s.s. and zero for An. arabiensis is not a realistic scenario, but the difference in P(B) is grounded in differences in their feeding behaviour, whereby An. arabiensis can utilize cattle more efficiently than An. gambiae s.s., for example.
Figure 11 shows the relative fraction of An. arabiensis. In addition to the pattern observed in Figure 10, it is also evident that if P(B) is low for An. arabiensis, An. gambiae s.s. dominates. P(B) can be interpreted as a parameter that describes both the probability of finding blood for reproduction and survival, and the energy spent in the search for a blood meal. For example, easy access to cattle might give An. arabiensis an advantage in exploiting breeding sites, which could lead to suppression of the number of An. gambiae s.s. if increased use of bed nets reduces the effective human population or causes higher mortality of anthropophilic species. This mechanism might help to explain the decline in An. gambiae s.s. observed by Bayoh et al. .
Sensitivity to the probability of finding blood in a spatial domain (pBlood2D)
This experiment is similar to pBlood1D, but this time we integrate the model for 5 years over the African domain. The experiment consists of two runs, for which the first has P(B) similar to Figure 8 and the second has P(B)=1 over all land areas for both species. The population density is space-invariant at 400 humans/km2 (remember that the number of mosquitoes is limited by the number of hosts). Thus, the only limitation in this experiment is the physical environment (air and water temperatures, relative humidity, wind and run-off), which is updated every 3 h. The initial conditions for the mosquito populations were the same for the two runs.
Even though we have stated that the probability of finding blood P(B) is an expression of the cost of finding a host, it might well be that P(B) also includes a component that describes the environment shaped by cattle and humans. Therefore, it should be noted that it is difficult to distinguish between the true probability of finding blood and the environmental changes caused by the presence of humans or cattle.
Under the scenario of equal probability of finding blood for the two species, An. gambiae s.s. loses the competition after 5 years (Figures 12 and 13), probably because of the greater reproductive potential of An. arabiensis. The only strongholds left for this species are DRC, Congo, and Gabon. Hence, the strategy of An. arabiensis to develop a larger body, produce more eggs, and possibly reduce adult mortality at the cost of spending more time in the aquatic stages is successful when access to blood is unlimited. An. arabiensis has extended its distribution as far north as the southern tip of Western Sahara. While the original set-up of the model (P0) limits the distribution of An. gambiae s.l. to approximately 17°N in the Sahel, the experiment with P(B)=1 (P1) has a distribution up to 22°N in Mali, Niger, Chad and Sudan. This is in line with observations of the northerly limit of An. gambiae s.l. [148–150]. The lack of An. gambiae s.l. north of 17° in the original set-up (P0) might be a result of the way the model is formulated. The population density is calculated within a box of approximately 50 km×50 km. It might well be the case that pockets of higher population/cattle densities within this box could sustain a mosquito population. This is not resolved in the model. It is also worth mentioning the study by De Meillon  of the anophelines of Namibia, which revealed that An. gambiae s.l. is present in large parts of the country. The original set-up (P0) allows sustainable mosquito populations in Namibia, while the density of An. gambiae s.l. in P1 is more comparable to the observations of De Meillon. The problems of capturing the distribution of An. gambiae s.l. in Namibia may originate from the problems of resolving pockets of high host density or changes in cattle density and distribution at the time of the study compared to the present day [43, 44, 152].
It is also worth mentioning that the density of An. gambiae s.l. in South Africa is not very sensitive to the probability of finding a blood meal. Hence, the distribution of An. gambiae s.l. is mainly restricted by climate according to the model.
Figures 12 and 13 show the distribution and density of An. arabiensis and An. gambiae s.s. under realistic (P0) and space-invariant (P1) P(B) after 6, 12, and 18 months. The integration was started on January 1 and the model was run for 5 years.
Mosquito transport (mosqTran)
The purpose of this experiment was to demonstrate how the initial conditions and competition influence the distribution of An. gambiae s.s. and An. arabiensis. To explore the theoretical dispersion distance and the influence of the initial conditions, we set up a simple experiment. In mosqTran(a) the model was initialized with An. arabiensis at -4.494381°E, 14.0154°N (Sahel), and An. gambiae s.s. at -4.494381°E, 6.502846°N (Cte d’Ivore, Ivory Coast) on January 1, 1989. The second experiment, mosqTran(b), had the same setup, but without An. arabiensis.
The purpose of this demonstration was to show the importance of mosquito movement and how new areas can or cannot be colonized. In a model in which movement is restricted, the vector range would also be restricted by the initial model conditions. For example, if only one point was specified for mosquitoes at the beginning of the integration, only the same point would have mosquitoes after 100 years. With dynamic movement the mosquitoes could colonize new areas if the environmental conditions, or the probability of finding blood, change over time.
Figure 14 shows the relative difference in An. gambiae s.s. distribution in the two experiments. It is evident that in the presence of An. arabiensis, An. gambiae s.s. fails to colonize large parts of Mali and Burkina Faso. It can be argued that this is not a result of the initial conditions, but of competition. Additional file 4 illustrates why this is indeed a result of the initial conditions, although the initial conditions would not play a role in the absence of competition.
Figure 15 shows the number of months required to reach a density of 20 mosquitoes/km2. It is interesting to note that dispersal occurs in pulses. The dispersal of An. arabiensis is slower than that of An. gambiae s.s., probably because of the drier conditions in the Sahel and An. gambiae s.s. reached the area before An. arabiensis (Figure 15). The simulations show that establishment in an already occupied area is a much slower process compared to the case of no competition. From the simulations we can also speculate on whether the dominance of one species can act as a barrier to genetic flow, like a mountain range or dessert. This also raises some questions regarding whether hibernation or dispersal is the mechanism behind the dominance of the An. gambiae s.s. M form in parts of Mali. Although there are strong indications that the An. gambiae s.s. M undergoes aestivation during the dry season [57, 58, 71], it is also possible that the persistence of the An. gambiae s.s. M form in the Niono district in Mali can serve as a refuge during the dry season . In both cases the M form receives a kick-start at the beginning of the rainy season, and might slow down the dispersal of An. arabiensis and the S form of An. gambiae s.s. A similar mechanism could contribute to the dominance of An. arabiensis in Ethiopia in the Turkana district, where the presence of An. arabiensis prevents rapid invasion by An. gambiae s.s.
We developed a model to predict the presence and abundance of An. arabiensis and An. gambiae s.s. The model is age-structured and includes mosquito dispersal.
Sensitivity tests showed that as well as temperature, relative humidity and mosquito size are important factors in malaria transmission. The result for body size is in line with several studies [7, 51, 54, 55, 88, 154] and thus the model captures some of the aspects related to higher survival among larger individuals. Note that we have not accounted for the higher metabolism in large mosquitoes , which might reduce survival under warm and dry conditions. There are also contrasting results with respect to body size and egg production . It is likely that there is an optimum size that depends on the environment and is a function of temperature and humidity. Currently there are few results to back up this statement. However, Sanford et al. found significant differences in Anopheles gambiae s.s. wing length between Mali and Guinea-Bissau .
We show that relative humidity can be important for malaria transmission. Several models have neglected the role of (relative) humidity [29, 157] and it is true that desiccation might not be a driver of mortality at moderate humidity (>70%?). The main argument for leaving out this parameter is the corresponding reduction in model complexity. As long as rainfall drive the carrying capacity, mosquito numbers will be restricted at lower humidity (no rain), and as a consequence the resulting number of mosquitoes can be limited for the wrong reasons, but with the correct result. For example, Ermert et al.  handle this deficiency by reducing vector survival during dry atmospheric conditions, defined as a function of 10-day accumulated rainfall. More studies on the survival of An. gambiae s.l. in relation to size and relative humidity in the range 5-40% are needed for more confidence in the role of humidity in the survival of An. gambiae s.l.
Assumption of exponential mortality has several advantages (see Figure 5 for examples of models in which exponential mortality is used). The model becomes fast to solve and it is easier to analyse the equations analytically. However, several studies have shown that mortality of An. gambiae s.l. is not exponential, and that inclusion of an age dimension alters the expected outcome of interventions targeted to reduce the vector population . Therefore, we believe that models in which age-dependent mortality is assumed should be further explored. The sensitivity tests also suggest that carrying capacity within a restricted area plays a role in the distribution of An. arabiensis and An. gambiae s.s. The true carrying capacity is hard to estimate on a continental scale and thus relies on qualified guesswork taking into account rainfall, groundwater and soil saturation, for example. Carrying capacity influences not only the relative distribution of the two species but also the total number of mosquitoes. To correctly estimate the biting rate, a correct estimate of carrying capacity is required, and thus more work is needed to parametrize puddle formation. It should also be noted that no current large-scale models can describe the formation of puddles as rivers retreat, as described by Animut et al. .
Experiment pBlood2D showed how the model responds to the parameter P(B), the probability of finding a blood meal. P(B) is important in describing a realistic distribution of An. arabiensis and An. gambiae s.s. Thus, we hypothesize that the large-scale distribution of bovines is key to the success of An. arabiensis. Likewise, large-scale human density favours the presence of An. gambiae s.s.
Finally, experiment mosqTran showed how the initial conditions influence the dispersal of An. gambiae s.s. (and An. arabiensis). The distribution of An. gambiae s.s. changes dramatically with the presence of An. arabiensis, and thus the initial model conditions are highly relevant for correct description of the distribution of the two species. When rainfall is highly seasonal, the first come, first served principle seems to be important for the success of a species in drier conditions. Whether or not this plays a role in the evolution of aestivation in An. gambiae s.s. M form  is a question that should be further investigated.
The strong influence of initial conditions on dispersal of the An. gambiae complex is not irrelevant when assessing the impact of climate change, since vectorial capacity varies between species.
The availability of mosquito models allows researchers to build on and improve our understanding of the role of the An. gambiae complex in malaria transmission. We hope to refine the model as new data on mosquito biology become available, and to incorporate the effects of interventions.
We are grateful to the National Center for Atmospheric Research (NCAR) for making their WRF model available in the public domain. We also thank the Bergen Centre for Computational Science for computational and other resources provided during this study. This work was made possible by grants from The Norwegian Programme for Development, Research and Education (NUFU) and the University of Bergen. We thank Steve Lindsay for providing data on the survival of An. gambiae s.s. under different temperatures and relative humidities. We thank two anonymous reviewers for their constructive comments, which helped us to improve the manuscript.
Centre for International Health, University of Bergen
Bjerknes Centre for Climate Research, University of Bergen/Uni Research
Lotka AJ: Contribution to the analysis of malaria epidemiology. II. General part (continued). Comparison of two formulae given by sir Ronald Ross.Am J Epidemiol 1923, 3:38–54. [http://aje.oxfordjournals.org]
Mordecai EA, Paaijmans KP, Johnson LR, Balzer C, Ben-Horin T, de Moor E, McNally A, Pawar S, Ryan SJ, Smith TC, Lafferty KD, Thrall P: Optimal temperature for malaria transmission is dramatically lower than previously predicted. 2012. [http://www.ncbi.nlm.nih.gov/pubmed/23050931]
Skamarock WC, Klemp JB, Dudhia J, Gill DO, Barker DM, Wang W, Powers JG: A Description of the Advanced Research WRF Version 2. Tech. rep., The National Center for Atmospheric Research 2005
Center for International Earth Science Information Network (CIESIN) - Columbia University; and Centro Internacional de Agricultura Tropical (CIAT): Gridded Population of the World Version 3 (GPWv3): population density grids. [http://sedac.ciesin.columbia.edu/gpw]
Robinson T, Fao’s Animal Production and Health Division: Observed livestock densities. 2011.
Wint W, Robinson T: Gridded livestock of the world - 2007. Tech. rep Food and agriculture organization of the United Nations, Rome 2007
Bayoh N: Studies on the development and survival of Anopheles gambiae sensu stricto at various temperatures and relative humidities. PhD thesis, University of Durham 2001
Afrane YA, Zhou G, Lawson BW, Githeko AK, Yan G: Effects of microclimatic changes caused by deforestation on the survivorship and reproductive fitness of Anopheles gambiae in western Kenya highlands.Am J Trop Med Hyg 2006, 74:772–778. [http://www.ncbi.nlm.nih.gov/pubmed/16687679]PubMed
Harrington LC, Vermeylen F, Jones JJ, Kitthawee S, Sithiprasasna R, Edman JD, Scott TW: Age-dependent survival of the dengue vector Aedes aegypti (Diptera: Culicidae) demonstrated by simultaneous release-recapture of different age cohorts.J Med Entomol 2008, 45:307–313. [http://www.ncbi.nlm.nih.gov/pubmed/18402147]PubMedView Article
Gimonneau G, Pombi M, Choisy M, Morand S, Dabiré RK, Simard F: Larval habitat segregation between the molecular forms of the mosquito Anopheles gambiae in a rice field area of Burkina Faso, West Africa. 2011.
Howard A, Adongo E, Vulule J, Githure J: Effects of a botanical larvicide derived from Azadirachta indica (the neem tree) on oviposition behaviour in Anopheles gambiae ss mosquitoes.J Med Plants Res 2011, 5:1948–1954.
Chen F, Mitchell K, Schaake J, Xue Y, Pan H, Koren V, Duan Q, Ek M, Betts A: Modeling of land surface evaporation by four schemes and comparison with FIFE observations.J Geophysical Res-Atmospheres 1996, 101:7251–7268.View Article
Koren V, Schaake J, Mitchell K, Duan Q, Chen F, Baker J: A parameterization of snowpack and frozen ground intended for NCEP weather and climate models.J Geophysical Res-Atmospheres 1999, 104:19569–19585.View Article
Betts A, Chen F, Mitchell K, Janjic Z: Assessment of the land surface and boundary layer models in two operational versions of the NCEP Eta Model using FIFE data.Mon Weather Rev 1997, 125:2896–2916.View Article
Ek M, Mitchell K, Lin Y, Rogers E, Grunmann P, Koren V, Gayno G, Tarpley J: Implementation of Noah land surface model advances in the national centers for environmental prediction operational mesoscale Eta model.J Geophysical Res-Atmospheres 2003, 108:GCP 12 1–16.View Article
Paaijmans KP, Heusinkveld BG, Jacobs AFG: A simplified model to predict diurnal water temperature dynamics in a shallow tropical water pool.Int J Biometeorol 2008, 52:797–803.PubMedView Article
Hoshen M, Morse A: A model structure for estimating malaria risk. In Environmental Change and Malaria Risk: Global and Local Implications, Volume 9. Edited by: Takken W, Martens P, Bogers RJ. : ; 2005:10–10.
Makaka G, Meyer E: Temperature stability of traditional and low-cost modern housing in the Eastern Cape, South Africa.J Build Phys 2006, 30:71.View Article
Okech BA, Gouagna LC, Knols BGJ, Kabiru EW, Killeen GF, Beier JC, Yan G, Githure JI: Influence of indoor microclimate and diet on survival of Anopheles gambiae s.s. (Diptera: Culicidae) in village house conditions in western Kenya.Int J Trop Insect Sci 2004,24(3):207–212.View Article
Okech BA, Gouagna LC, Killeen GF, Knols BGJ, Kabiru EW, Beier JC, Yan G, Githure JI: Influence of sugar availability and indoor microclimate on survival of Anopheles gambiae (Diptera: Culicidae) under semifield conditions in western Kenya.J Med Entomol 2003, 40:657–663. [http://www.ncbi.nlm.nih.gov/pubmed/14596279]PubMedView Article
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.