# Modelling the population dynamics of Rift Valley fever virus mosquito vectors in the western Mediterranean Basin

Alex Drouin, Thomas Balenghien, Benoit Durand, Carles Aranda, Amal Bennouna, Ali Bouattour, Said C Boubidi, Annamaria Conte, Sarah Delacour, Maria Goffredo, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Alex Drouin, Thomas Balenghien, Benoit Durand, Carles Aranda, Amal Bennouna, et al. Modelling the population dynamics of Rift Valley fever virus mosquito vectors in the western Mediterranean Basin. Ecological Modelling, 2025, 502, pp.111013. 10.1016/j.ecolmodel.2024.111013. pasteur-05024977

HAL Id: pasteur-05024977

https://riip.hal.science/pasteur-05024977v1

Submitted on 8 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Check for
updates

# Modelling the population dynamics of Rift Valley fever virus mosquito vectors in the western Mediterranean Basin

Alex Drouin<sup>a,b,c</sup>, Thomas Balenghien<sup>a,b</sup>, Benoit Durand<sup>c</sup>, Carles Aranda<sup>d,e</sup>, Amal Bennouna<sup>f</sup>, Ali Bouattour<sup>g</sup>, Said C Boubidi<sup>h</sup>, Annamaria Conte<sup>i</sup>, Sarah Delacour<sup>j</sup>, Maria Goffredo<sup>i</sup>, Oumnia Himmi<sup>k</sup>, Grégory L'Ambert<sup>l</sup>, Francis Schaffner<sup>m</sup>, Véronique Chevalier<sup>a,b,n,o,*</sup>

## ARTICLE INFO

**Keywords:**

Rift valley fever

Population dynamics

Modelling

Mediterranean basin

Mosquitoes

## ABSTRACT

Rift Valley fever (RVF) is a zoonotic vector-borne disease mainly transmitted by mosquitoes, and present in Africa, the Arabian Peninsula, and the Indian Ocean. The endemic situation in Mauritania, and the recent outbreaks in Libya have raised concerns about the potential spread of the virus in the western Mediterranean Basin, where competent mosquitoes are present. However, given the large diversity of climates and landscapes in this region, the areas and periods at risk of RVF virus (RVFV) transmission remain unknown. Vector abundance is one of the drivers of arboviruses transmission, therefore knowledge on mosquito species distributions and population dynamics is needed to implement surveillance and to assess the risk of RVFV circulation. Here, we adapted a published modelling framework of mosquito population dynamics to five potential RVFV vectors in the western Mediterranean Basin (*Aedes caspius*, *Aedes detritus*, *Aedes vexans*, *Culex pipiens* and *Culex theileri*). The mechanistic model was designed with a daily time step and a 0.1° x 0.1° spatial resolution and takes temperature and precipitations data as inputs, along with published vector distribution maps. We used mosquito trapping data from Spain, France, Italy and Morocco to calibrate the model, and we produced monthly maps of abundance of the five vectors for the whole studied area. We then evaluated the model performances by assessing the correlation between field data and model predictions. Finally, we performed a sensitivity analysis to identify the main influential parameters. The model was able to reproduce most of the abundance peaks for the five mosquito species. Goodness-of-fit was high for *Aedes* species, especially for *Ae. caspius*, a highly competent mosquito for RVFV transmission, but lower for *Culex* species, with potential overpredictions in some regions. More knowledge is required about the presence and abundance of potential RVFV vectors in the Mediterranean Basin to improve predictions. However, this first model allows to identify seasons and areas with high vectors abundances that could be used in the future for surveillance of the disease.

<footnote><sup>a</sup> CIRAD, UMR ASTRE, F-34398 Montpellier, France</footnote>

<footnote><sup>b</sup> ASTRE, Univ Montpellier, CIRAD, INRAE, Montpellier, France</footnote>

<footnote><sup>c</sup> EPIMIM, Laboratoire de Santé Animale, Anses, Ecole Nationale Vétérinaire d'Alfort, 94700 Maisons-Alfort, France</footnote>

<footnote><sup>d</sup> OIRTA, Centre de Recerca en Sanitat Animal (CRESA, IRTA-UAB), Campus de la Universitat Autònoma de Barcelona, Bellaterra, Cerdanyola del Vallès, Spain</footnote>

<footnote><sup>e</sup> Servei de Control de Mosquits del Consell Comarcal del Baix Llobregat, Barcelona, Spain</footnote>

<footnote><sup>f</sup> Department of Virology. Pathogen Discovery Laboratory. Institut Pasteur, Paris, France</footnote>

<footnote><sup>g</sup> Laboratoire Des Virus, Vecteurs et Hôtes (LR20IPT02), Institut Pasteur de Tunis, Université Tunis El Manar, Tunis 1002, Tunisia</footnote>

<footnote><sup>h</sup> Laboratoire d'Eco-épidémiologie Parasitaire et Génétique des Populations, Route du Petit Staoueli, Institut Pasteur d'Algérie, Dely-Brahim 16047, Algeria</footnote>

<footnote><sup>i</sup> Istituto Zooprofilattico Sperimentale dell'Abruzzo e del Molise 'G. Caporale', Campo Boario, 64100 Teramo, Italy</footnote>

<footnote><sup>j</sup> Animal Health Department, Faculty of Veterinary Medicine of Zaragoza, University of Zaragoza, C/Miguel Servet 177, 50013 Zaragoza, Spain</footnote>

<footnote><sup>k</sup> Geo-Biodiversity and Natural Patrimony Laboratory (GEOBIOL), Scientific Institute, Mohammed V University, Rabat, Morocco</footnote>

<footnote><sup>l</sup> EID Méditerranée, Direction Technique, Montpellier, France</footnote>

<footnote><sup>m</sup> Francis Schaffner Consultancy (FSC), Lorracherstrasse 50, 4125 Riehen, Switzerland</footnote>

<footnote><sup>n</sup> Epidemiology and Clinical Research Unit, Institut Pasteur de Madagascar, Antananarivo, Madagascar</footnote>

<footnote><sup>o</sup> CIRAD, UMR ASTRE, Antananarivo, Madagascar</footnote>

<footnote>* Corresponding author.
E-mail address: veronique.chevalier@cirad.fr (V. Chevalier).</footnote>
# 1. Introduction

Rift Valley fever (RVF) is a vector-borne zoonosis caused by a *Phlebovirus*, transmitted by mosquitoes, and affecting livestock, including dromedary camels (*Linthicum* et al., 2016). The disease causes high abortion rates in pregnant females and severe mortality in young animals. In humans, RVF is often asymptomatic or induces a flu-like syndrome, but complications can occur with hepatic, ocular or neurologic consequences (Javelle et al., 2020). There is no treatment or vaccine for humans. For livestock, live and killed vaccines are available, but few countries are implementing vaccination campaigns (Dungu et al., 2018). RVF thus impacts public health and leads to important economic losses in affected countries (Peyre et al., 2015): the cost of the 2006-2007 outbreaks observed in East Africa has been estimated to 6.7 million dollars in Tanzania, 66 million dollars in Kenya and 471 million dollars in Somalia. Since its first observation in 1930, the RVF virus (RVFV) has been reported in almost all countries of sub-Saharan Africa, in the Indian Ocean and in the Arabian Peninsula (Nielsen et al., 2020).

RVFV is known to mainly spread through the movements of infected livestock. Indeed, importation of infected ruminants has been suggested as the probable introduction route during the first documented incursion of RVF in the Mediterranean Basin, in the Nile delta (1977-1978), that led to 200,000 human cases and 600 deaths (Kenawy et al., 2018). The disease has since been regularly reported in Egypt (Kenawy et al., 2018), in Mauritania (resulting in 47 confirmed human cases and 23 deaths during the last outbreak in 2022) (WHO, 2022), and recently emerged in Libya, where 3 cases have been reported in small ruminants on the Mediterranean coast near Tripoli in 2021 (WOAH, 2021). Virus emergence in Libya was again suggested to be linked to illegal introduction of livestock. This situation, associated with the existence of unmonitored live animals trade routes from the sub-Saharan region to the northern part of Africa (EFSA Panel on Animal Health and Welfare (AHAW), 2013; Nielsen et al., 2020), raise concerns about the further risk of introduction of the virus in Maghreb, where some serological evidences of RVFV exposure has already been found in livestock (Di Nardo et al., 2014; Hellal et al., 2021). Moreover, among the 47 mosquito species theoretically able to transmit the RVFV after exposure (Lumley et al., 2017), 10 are present in the Mediterranean Basin, and *Aedes caspius*, *Aedes detritus*, *Aedes vexans*, *Culex pipiens* and *Culex theileri* are the most likely to eventually transmit the virus due to their ecological characteristics (Drouin et al., 2022).

RVFV transmission mainly depends on the abundance of competent vectors and on the density of ruminant hosts. Very few studies have assessed the risk of RVFV circulation in the Mediterranean Basin, using either MCDA methodology (Arsevska et al., 2016; Sánchez-Vizcaíno et al., 2013; Tran et al., 2013a) or modelling studies (Nielsen et al., 2020). In absence of sufficient data, the authors used temperature, rainfall, proximity to aquatic areas, elevation, or land cover data, coupled with expert knowledge, to approximate mosquito presence or abundance. As mosquito biology is mainly driven by meteorological factors, especially temperature (Bellone and Failloux, 2020) and rainfall (Arsevska et al., 2016; Balenghien et al., 2010), incorporating these parameters into population dynamics or transmission modelling studies is crucial to improve the predictions, both geographically and temporally. However, these elements are often lacking in models of mosquito-borne diseases (Reiner et al., 2013). This is also the case for RVF, for which only 22 % of the models reviewed by Cecilia and colleagues took into account these abiotic factors (Cecilia et al., 2022), leading to limitations in both identifying areas at risk and accounting for the seasonal variability of transmission. Understanding and predicting RVFV vector population dynamics is thus an essential step in evaluating the risk of circulation of the virus in a given area. To overcome the lack of sufficient field data, mechanistic models of vector population dynamics have been elaborated, but at a local scale only in the Mediterranean Basin (Balenghien et al., 2010; Cailly et al., 2012; Ezanno et al., 2015; Groen et al., 2017). Moreover, these models were restricted to European locations and, even if *Cx. pipiens* and *Ae. caspius*

were included, the studies did not focus on RVF and therefore did not consider all the main potential RVFV vector species. Lastly, other authors have developed wild-ranging distribution models (Amdouni et al., 2022; Ducheyne et al., 2013; Gangoso et al., 2020; Mughini-Gras et al., 2014; Outammassine et al., 2022, 2021; Schaffner et al., 2016; Versteirt et al., 2013; Wint et al., 2020) but neglected seasonal variation of abundances.

A generic mechanistic model for mosquito population dynamics, using temperature as the main driver of development and mortality, was first developed by Cailly and colleagues (Cailly et al., 2012) and adapted to *Culex* and *Aedes* genus (Ezanno et al., 2015). Then, this model was successfully implemented under several climatic conditions. In northern Senegal, it has been adapted to *Ae. vexans*, *Cx. poicilipes* and *Cx. tritaeniorhynchus* at a regional scale, using satellite meteorological data and water body detection, but the best fit was obtained using *in situ* meteorological data (Tran et al., 2019). In Botswana, the same model was adapted to *Cx. pipiens*, in four different sites experiencing RVFV circulation (Hammami et al., 2016). The authors compared four scenarios in order to reproduce the observed mosquito dynamics, taking into account *in situ* temperature data, precipitations and/or surface of flooded areas obtained through satellite imagery. The best model fit was obtained by integrating temperature and precipitations data, with or without considering surfaces of flooded areas, depending on the study site.

Here, as a first step to identify areas and periods at risk of RVFV transmission in the Mediterranean Basin, we adapted and spatialized the model of Cailly and colleagues (Cailly et al., 2012; Ezanno et al., 2015), and used it to simulate the population dynamics of the five main potential RVFV vectors in the western part of the Mediterranean Basin, namely *Ae. caspius*, *Ae. detritus*, *Ae. vexans*, *Cx. pipiens* and *Cx. theileri*. We parametrized and calibrated the model for each mosquito species using a literature review and mosquito trapping data from Spain, France, Italy and Morocco. We then used this model to produce monthly maps of abundance for the five species and assessed the quality of its predictions by computing the temporal correlation between the model outputs and the observed entomological data.

# 2. Materials and methods

## 2.1. Data

### 2.1.1. Study area

It covered the western part of the Mediterranean Basin, including territories from Portugal, Spain, France, Italy, Morocco, Algeria and Tunisia, from longitudes -13.23° to 18.67°, and from latitudes 27.62° to 47.12° (Supplementary Figs. S1-S3). This area encompassed different ecoregions belonging to the following types: temperate broadleaf and mixed forests, temperate conifer forests, Mediterranean forests woodlands and scrubs, flooded grasslands, deserts and xeric shrublands, montane grasslands and shrublands (Dinerstein et al., 2017).

### 2.1.2. Meteorological data

They were extracted from the ERA5-Land reanalysis dataset (Muñoz Sabater, 2019), containing hourly estimates of temperature and precipitations with a 0.1° x 0.1° spatial resolution and covering the period from 1981 to present. We obtained these data from the Copernicus Climate Change Service (C3S) Climate Data Store, for the whole zone of interest and between 01/01/1990 and 31/12/2021. The selection of the data for the study area resulted in a 35,789-pixel grid, for which we computed the daily mean temperature ('2m temperature' variable in the ERA5-Land dataset, Supplementary Fig. S1) and extracted the total daily precipitations by selecting the last cumulative value of the day ('total precipitation' variable, Supplementary Fig. S2).

### 2.1.3. Probability of vectors presence data

We used the data provided by Wint and colleagues (Wint et al., 2020). In this study, the authors used Boosted Regression Trees and Random Forest to extrapolate data of mosquito presence, based on
landcover information of each pixel. They provided the probabilities of presence (between 0 and 1) of *Ae. albopictus*, *Ae. caspius*, *Ae. japonicus*, *Ae. vexans*, *Cx. pipiens* and *Cx. theileri* for Europe and North Africa with a 0.01° spatial resolution. Because we needed the probabilities of presence at a 0.1° resolution to match the ERA5-Land dataset, we aggregated the data. First, we isolated the land pixels by cropping the probabilities of presence rasters using the ESRI World Countries layers (ESRI, 2022). Then we computed the mean values of the 0.01° pixels included in the 0.1° pixels. Because the probability of presence data for *Ae. detritus* is not provided by Wint and colleagues (Wint et al., 2020), we used the data available for *Ae. caspius* for both species, considering that these two halophilic mosquitoes share a large proportion of their coastal breeding sites (Trari, 2017).

### 2.1.4. Entomological data

We used geolocated daily mosquito trapping data from 10 different sources: 2 from Spain (SP1 and SP2), 4 from France (FR1 to FR4), 1 from Italy (IT1) and 3 from Morocco (MO1 to MO3), sampled from 1994 to 2021 (see the 'Entomological data' section in Supplementary Information and Supplementary Fig. S3). Mosquitoes were collected using bird- and horse-baited traps, traps with chemical lure, CO2-baited traps (with or without a light source), human landing catches, light traps, and aquatic nets for larvae (Supplementary Table S1).

Given the diversity of mosquito trapping methods and the potential variability in implementation of the same method among the different collection protocols (material, duration of the trapping session, etc.), we did not pool the results from different sources. Hereinafter, data

collected using a given trapping method in a given source will be referred to as a 'dataset', and longitudinal data from a given trap as a 'time series'. Overall, the entomological data available included 15 datasets for *Ae. caspius*, 14 for *Ae. detritus*, 13 for *Ae. vexans*, 16 for *Cx. pipiens* and 15 for *Cx. theileri*, representing 686, 685, 682, 620 and 615 time series, respectively (Supplementary Table S1). The length of the times series ranged from 1 mosquito count to 143 (Supplementary Fig. S4). All these data were geolocated in 570 pixels, with at least 505 different pixels per species (Supplementary Fig. S3).

## 2.2. Mosquito population dynamics model

We tailored the model of mosquito population dynamics of Cailly and colleagues (Cailly et al., 2012; Ezanno et al., 2015) to the five species of interest. This generic model is mechanistic and compartmental, and is based on ordinary differential equations (ODE). It describes the life cycle of mosquitoes using 10 compartments: the aquatic stages, i.e., eggs (E), larvae (L) and pupae (P), and the aerial stages, i.e., emerging adults ($A_{em}$), host-seeking adults ($A_h$), adults during the digestion of the blood meal and egg maturation ($A_g$) and ovipositing adults ($A_o$) (Fig. 1). Only adult females are considered, and the three last stages are split between nulliparous females (subscript 1) and parous females (subscript 2). Because mosquitoes are poikilotherm organisms, their dynamics are tightly linked to environmental conditions. Thus, several transition processes of the model depend on meteorological factors. Definitions and values of the parameters and functions used in the model are given in Table 1.

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Flow diagram of the mosquito population dynamics model in a given pixel, adapted from Ezanno et al. (2015). The blue background refers to the aquatic stages (E: eggs, L: larvae, P: pupae) and the white background to the aerial stages ($A_{em}$: emerging adults, $A_h$: host-seeking adults, $A_g$: adults during the digestion of the blood meal and egg maturation, $A_o$: ovipositing adults). $A_h$, $A_g$, and $A_o$ are split between nulliparous females (subscript 1) and multiparous females (subscript 2). Greek letters refer to constant parameters ($\beta_1$: number of eggs laid/ovipositing nulliparous female, $\beta_2$: number of eggs laid per ovipositing parous female, $\sigma_{em}$: sex ratio at emergence, $\gamma_{A_{em}}$: development rate of emerging adults, $\gamma_{A_h}$: transition rate from host-seeking to engorged adults, $\gamma_{A_g}$: transition rate from ovipositing to host-seeking adults, $\mu_E$: egg mortality rate, $\mu_L$: minimum larval mortality rate, $\mu_P$: minimum pupal mortality rate, $\mu_A$: minimum adult mortality rate, $\mu_{em}$: mortality rate during emergence, $\mu_r$: mortality rate related to seeking behavior, $\kappa_L$: environment carrying capacity for larvae and $\kappa_P$: environment carrying capacity for pupae) and Latin letters to meteorological-dependent functions taking rainfall and temperature as inputs ($f_E$: transition rate from eggs to larvae, $f_L$: transition rate from larvae to pupae, $f_P$: transition rate from pupae to emerging adults, $f_A_g$: transition rate from engorged to ovipositing adults, $m_L$: larval mortality rate, $m_P$: pupal mortality rate and $m_A$: adult mortality rate), as defined in Table 1.</figcaption></figure>
<table><caption>**Table 1**

Constant parameters and meteorological-dependent functions of the model, for each of the five species of interest.</caption><thead><tr><th>Parameter</th><th>Definition</th><th>Ae. caspius</th><th>Ae. detritus</th><th>Ae. vexans</th><th>Cx. pipiens</th><th>Cx. theileri</th><th>References</th></tr></thead><tbody><tr><td>β<sub>1</sub></td><td>Number of eggs laid/ovipositing nulliparous female</td><td>160</td><td>120</td><td>160</td><td>141</td><td>92</td><td>Berchi et al., 2012, Bogojević et al., 2011, Schaeffer et al., 2008, Cailly et al., 2012, Cecilia et al., 2022</td></tr><tr><td>β<sub>2</sub></td><td>Number of eggs laid/ovipositing parous female</td><td>80</td><td>60*</td><td>80*</td><td>80</td><td>78</td><td>Berchi et al., 2012, Bogojević et al., 2011, Schaeffer et al., 2008, Cecilia et al., 2022</td></tr><tr><td>f<sub>E</sub><br>(day<sup>-1</sup>)</td><td>Transition rate from eggs to larvae</td><td>Boolean<br>{ 1 if θ > 15 & R > 15<br>0 otherwise </td><td>Boolean<br>{ 1 if R > 15<br>0 otherwise </td><td>Boolean<br>{ 1 if θ > 10 & R > 15<br>0 otherwise </td><td>Logan model with<br>A = 0.16<br>B = 0.11<br>C = 5.01<br>θ<sub>min</sub> = 10<br>θ<sub>max</sub> = 35</td><td>Logan model with<br>A = 0.06<br>B = 0.17<br>C = 3.62<br>θ<sub>min</sub> = 12<br>θ<sub>max</sub> = 39</td><td>Bogojević et al., 2011, Cecilia et al., 2020, Clarkson and Enevoldson, 2020, De Ascentis et al., 2022</td></tr><tr><td>f<sub>L</sub><br>(day<sup>-1</sup>)</td><td>Transition rate from larvae to pupae</td><td>f<sub>P</sub><br>1.65</td><td>f<sub>P</sub><sup>*</sup><br>1.65</td><td>Logan model with<br>A = 0.02<br>B = 0.15<br>C = 3.56<br>θ<sub>min</sub> = 10<br>θ<sub>max</sub> = 35</td><td>Logan model with<br>A = 0.017<br>B = 0.13<br>C = 3.98<br>θ<sub>min</sub> = 7<br>θ<sub>max</sub> = 35</td><td>Logan model with<br>A = 0.02<br>B = 0.12<br>C = 3.28<br>θ<sub>min</sub> = 12<br>θ<sub>max</sub> = 36</td><td>Bogojević et al., 2011, Di Nardo et al., 2014, Dinerstein et al., 2017, Drouin et al., 2022</td></tr><tr><td>f<sub>P</sub><br>(day<sup>-1</sup>)</td><td>Transition rate from pupae to emerging adults</td><td>Logan model with<br>A = 0.14<br>B = 0.16<br>C = 5.01<br>θ<sub>min</sub> = 10<br>θ<sub>max</sub> = 35</td><td>Logan model with *<br>A = 0.14<br>B = 0.16<br>C = 5.01<br>θ<sub>min</sub> = 10<br>θ<sub>max</sub> = 35</td><td>Logan model with<br>A = 0.21<br>B = 0.11<br>C = 3.41<br>θ<sub>min</sub> = 15<br>θ<sub>max</sub> = 35</td><td>Logan model with<br>A = 0.13<br>B = 0.23<br>C = 4.09<br>θ<sub>min</sub> = 10<br>θ<sub>max</sub> = 35</td><td>Logan model with<br>A = 0.15<br>B = 0.14<br>C = 4.67<br>θ<sub>min</sub> = 12<br>θ<sub>max</sub> = 36</td><td>Bogojević et al., 2011, Cecilia et al., 2020, Di Nardo et al., 2014, Dinerstein et al., 2017, Drouin et al., 2022</td></tr><tr><td>σ<sub>em</sub></td><td>Sex ratio at emergence</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td><td>Bogojević et al., 2011, Schaeffer et al., 2008, Cecilia et al., 2020, Ducheyne et al., 2013, Dungu et al., 2018</td></tr><tr><td>γ<sub>Aem</sub><br>(day<sup>-1</sup>)</td><td>Development rate of emerging adults</td><td>0.4</td><td>0.4*</td><td>0.5</td><td>1.143</td><td>1.143*</td><td>Bogojević et al., 2011, Cecilia et al., 2020, Durand et al., 2020, EFSA Panel on Animal Health and Welfare (AHAW) 2013, El Ouali Lalami et al., 2009</td></tr><tr><td>γ<sub>Ah</sub><br>(day<sup>-1</sup>)</td><td>Transition rate from host-seeking to engorged adults</td><td>0.222</td><td>0.222*</td><td>0.33</td><td>0.885</td><td>0.885*</td><td>Bogojević et al., 2011, Schaeffer et al., 2008, Ducheyne et al., 2013, ESRI, 2022</td></tr><tr><td>f<sub>Ag</sub><br>(day<sup>-1</sup>)</td><td>Transition rate from engorged to ovipositing adults</td><td>θ - T<sub>Ag</sub>/TDD<sub>Ag</sub><br>T<sub>Ag</sub> = 10<br>TDD<sub>Ag</sub> = 78<br>f<sub>Ag</sub> > 0</td><td>θ - T<sub>Ag</sub><sup>*</sup>/TDD<sub>Ag</sub><br>T<sub>Ag</sub> = 10<br>TDD<sub>Ag</sub> = 78<br>f<sub>Ag</sub> > 0</td><td>θ - T<sub>Ag</sub><sup>*</sup>/TDD<sub>Ag</sub><br>T<sub>Ag</sub> = 10<br>TDD<sub>Ag</sub> = 78<br>f<sub>Ag</sub> > 0</td><td>θ - T<sub>Ag</sub>/TDD<sub>Ag</sub><br>T<sub>Ag</sub> = 9.8<br>TDD<sub>Ag</sub> = 64.4<br>f<sub>Ag</sub> > 0</td><td>θ - T<sub>Ag</sub>/TDD<sub>Ag</sub><br>T<sub>Ag</sub> = 9.8<br>TDD<sub>Ag</sub> = 73.09769<br>f<sub>Ag</sub> > 0</td><td>Bogojević et al., 2011, El Ouali Lalami et al., 2009, Ezanno et al., 2015</td></tr><tr><td>γ<sub>Ah</sub><br>(day<sup>-1</sup>)</td><td>Transition rate from ovipositing to host-seeking adults</td><td>0.222</td><td>0.222*</td><td>0.33</td><td>2</td><td>2*</td><td>Berchi et al., 2012, Bogojević et al., 2011, Ducheyne et al., 2013, Gabinaud, 1975</td></tr><tr><td>μ<sub>E</sub><br>(day<sup>-1</sup>)</td><td>Egg mortality rate</td><td>0.001*</td><td>0.001*</td><td>0.001</td><td>0.262</td><td>0.0290</td><td>Bogojević et al., 2011, Schaeffer et al., 2008, De Ascentis et al., 2022, Ducheyne et al., 2013</td></tr><tr><td>m<sub>L</sub><br>(day<sup>-1</sup>)</td><td>Larval mortality rate</td><td>e<sup>-θ/2</sup> + μ<sub>L</sub></td><td>e<sup>-θ/2</sup> + μ<sub>L</sub></td><td>e<sup>-θ/2</sup> + μ<sub>L</sub></td><td>e<sup>-θ/2</sup> + μ<sub>L</sub></td><td>e<sup>-θ/2</sup> + μ<sub>L</sub></td><td>Bogojević et al., 2011, Schaeffer et al., 2008, Ducheyne et al., 2013</td></tr><tr><td>μ<sub>L</sub><br>(day<sup>-1</sup>)</td><td>Minimum larval mortality rate</td><td>0.0367</td><td>0.003113426</td><td>0.04671658</td><td>0.0304</td><td>0.05392209</td><td>Bogojević et al., 2011, Drouin et al., 2022, Gangoso et al., 2020, Groen et al., 2017</td></tr></tbody></table>
<table><caption>**Table 1 (continued)**</caption><thead><tr><th>Parameter</th><th>Definition</th><th>Ae. caspius</th><th>Ae. detritus</th><th>Ae. vexans</th><th>Cx. pipiens</th><th>Cx. theileri</th><th>References</th></tr></thead><tbody><tr><td>m<sub>p</sub><br>(day<sup>-1</sup>)</td><td>Pupal mortality rate</td><td>e<sup>-θ/2</sup> + μ<sub>p</sub></td><td>e<sup>-θ/2</sup> + μ<sub>p</sub></td><td>e<sup>-θ/2</sup> + μ<sub>p</sub></td><td>e<sup>-θ/2</sup> + μ<sub>p</sub></td><td>e<sup>-θ/2</sup> + μ<sub>p</sub></td><td>Bogojević et al., 2011,<br>Schaeffer et al., 2008,<br>Ducheyne et al., 2013</td></tr><tr><td>μ<sub>p</sub><br>(day<sup>-1</sup>)</td><td>Minimum pupal<br>mortality rate</td><td>0.12</td><td>0.001137951</td><td>0.02200325</td><td>0.0146</td><td>0.06593466</td><td>Bogojević et al., 2011,<br>Drouin et al., 2022,<br>Gangoso et al., 2020,<br>Groen et al., 2017</td></tr><tr><td>m<sub>A</sub><br>(day<sup>-1</sup>)</td><td>Adult mortality rate</td><td>-0.005941 +<br>0.002965 × θ<br>for m<sub>A</sub> > μ<sub>A</sub></td><td>-0.005941 +<br>0.002965 × θ<sup>*</sup><br>for m<sub>A</sub> > μ<sub>A</sub></td><td>0.02320402 +<br>0.002075247 × θ<br>for m<sub>A</sub> > μ<sub>A</sub></td><td>-0.06578724 +<br>0.004460362 × θ<br>for m<sub>A</sub> > μ<sub>A</sub></td><td>0.01841533 +<br>0.002431452 × θ<br>for m<sub>A</sub> > μ<sub>A</sub></td><td>Bogojević et al., 2011,<br>Hammadi et al., 2009,<br>Hammami et al., 2016,<br>Hartvigsen, 2013</td></tr><tr><td>μ<sub>A</sub><br>(day<sup>-1</sup>)</td><td>Minimum adult<br>mortality rate</td><td>0.03773585</td><td>0.03773585<sup>*</sup></td><td>0.02673797</td><td>0.01083424</td><td>0.01083424<sup>*</sup></td><td>Hammadi et al., 2009,<br>Hammami et al., 2016,<br>Hellal et al., 2021, Himmi,<br>2007</td></tr><tr><td>m<sup>dilia</sup><sub>A</sub><br>(day<sup>-1</sup>)</td><td>Adult mortality rate<br>during diapause<br>stage</td><td>m<sub>A</sub></td><td>m<sub>A</sub><sup>*</sup></td><td>m<sub>A</sub><sup>*</sup></td><td>μ<sub>A</sub></td><td>μ<sub>A</sub><sup>*</sup></td><td>Bogojević et al., 2011</td></tr><tr><td>μ<sub>em</sub><br>(day<sup>-1</sup>)</td><td>Mortality rate during<br>emergence</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>Bogojević et al., 2011,<br>Ducheyne et al., 2013</td></tr><tr><td>μ<sub>r</sub><br>(day<sup>-1</sup>)</td><td>Mortality rate<br>related to seeking<br>behavior</td><td>0.08</td><td>0.08</td><td>0.08</td><td>0.08</td><td>0.08</td><td>Bogojević et al., 2011,<br>Ducheyne et al., 2013</td></tr><tr><td>κ<sub>L,max</sub><br>(larvae.ha<sup>-1</sup>)</td><td>Maximal<br>environment<br>carrying capacity for<br>larvae</td><td>10<sup>10</sup></td><td>10<sup>10</sup></td><td>10<sup>10</sup></td><td>8 × 10<sup>8</sup></td><td>8 × 10<sup>8</sup></td><td>Bogojević et al., 2011</td></tr><tr><td>κ<sub>P,max</sub><br>(pupae.ha<sup>-1</sup>)</td><td>Maximal<br>environment<br>carrying capacity for<br>pupae</td><td>10<sup>8</sup></td><td>10<sup>8</sup></td><td>10<sup>8</sup></td><td>10<sup>7</sup></td><td>10<sup>7</sup></td><td>Bogojević et al., 2011</td></tr><tr><td>φ</td><td>Boolean variable<br>accounting for the<br>differences in<br>resistant stage<br>during diapause</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>Bogojević et al., 2011,<br>Ducheyne et al., 2013</td></tr><tr><td>start</td><td>First day of the<br>favorable season</td><td>15<sup>th</sup> March</td><td>No diapause</td><td>15<sup>th</sup> March</td><td>15<sup>th</sup> March</td><td>15<sup>th</sup> March</td><td>Bogojević et al., 2011,<br>Cecilia et al., 2020</td></tr><tr><td>end</td><td>Last day of the<br>favorable season</td><td>14<sup>th</sup> October</td><td></td><td>14<sup>th</sup> October</td><td>21<sup>st</sup> September</td><td>21<sup>st</sup> September</td><td>Bogojević et al., 2011,<br>Cecilia et al., 2020</td></tr></tbody></table>

Greek letters refer to constant parameters and Latin to meteorological-dependent functions. θ refers to the temperature (Celsius degrees) and R to the precipitations (in millimeters). *Parameter extrapolated from another species of the same genus, in absence of specific data for the species considered.

Regarding the aquatic stages (Eq. 1a and Fig. 1), the egg compartment (E) is renewed by the ovipositing females (A<sub>1o</sub> and A<sub>2o</sub>), with a constant ovipositing rate, but with a different number of eggs laid by nulliparous females (A<sub>1o</sub>) and multiparous females (A<sub>2o</sub>). Egg hatching rates depend on temperature for *Culex* species, but hatching is triggered by precipitations for *Aedes* species. The transition rates from larvae (L) to pupae (P) and the emergence of adults (A<sub>em</sub>) are dependent on temperature for all species. Egg mortality rates are constant, but larvae and pupae mortality rates depend on temperature. Larvae mortality rates and emergence success of pupae are also dependent on the population density in the breeding site (Cailly et al., 2012), represented by the ratios L/κ<sub>L</sub> and P/κ<sub>P</sub>, with κ<sub>L</sub> and κ<sub>P</sub> being the carrying capacities of the environment allowing to limit the population growth.

$$ 
\left\{
\begin{array}{l}
\text{Aquatic stages} \\
\dot{E} = \gamma_{A_o}(\beta_1 A_{1_o} + \beta_2 A_{2_o}) - (\mu_E + f_E)E \\
\dot{L} = f_E E - \left(m_L \left(1 + \frac{L}{\kappa_L}\right) + f_L\right)L \\
\dot{P} = f_L L - (m_P + f_P)P
\end{array}
\right.
\qquad (1a) $$

Only females enter the aerial stages modelled by Eq. 1b, with a defined sex ratio σ<sub>em</sub>. Transition rates between adult stages are all constant, except the transition rates from engorged to ovipositing adults, representing the duration of egg maturation, and which depend on temperature (Cailly et al., 2012). Adult mortality rates are also temperature-dependent, with increased rates for the most mobile stages (host-seeking females A<sub>h</sub> and ovipositing females A<sub>o</sub>), represented by μ<sub>r</sub>.
$$
\left\{
\begin{aligned}
\text{Aerial stages} \\
\dot{A}_{em} &= f_p P \sigma_{em} \exp \left( - \mu_{em} \left( 1 + \frac{P}{\kappa_P} \right) \right) - (m_A + \gamma_{Aem}) A_{em} \\
\dot{A}_{1h} &= \gamma_{Aem} A_{em} - (m_A + \mu_r + \gamma_{A_h}) A_{1h} \\
\dot{A}_{1g} &= \gamma_{Ah} A_{1h} - (m_A + f_{A_g}) A_{1g} \\
\dot{A}_{1o} &= f_{A_g} A_{1g} - (m_A + \mu_r + \gamma_{A_o}) A_{1o} \\
\dot{A}_{2h} &= \gamma_{A_o} (A_{1o} + A_{2o}) - (m_A + \mu_r + \gamma_{A_h}) A_{2h} \\
\dot{A}_{2g} &= \gamma_{A_h} A_{2h} - (m_A + f_{A_g}) A_{2g} \\
\dot{A}_{2o} &= f_{A_g} A_{2g} - (m_A + \mu_r + \gamma_{A_o}) A_{2o}
\end{aligned}
\right.
\qquad (1b)
$$

Under temperate climates with a marked winter season, mosquitoes go through an inactive period, commonly called diapause. *Aedes* mosquitoes survive through this unfavorable season in the egg stage (E), and *Culex* mosquitoes in the adult stage (A<sub>em</sub>). Therefore, we used two ODE systems to represent the life cycle of mosquitoes in the European region of the study area, one for the favorable season (Eqs. 1a and 1b), i.e., the season of mosquito activity (defined by the parameters *start* and *end*), and one for the unfavorable season (Eq. 2, figuring only the states for which the equation differed from the favorable season). In Eq. 2, Φ is a Boolean accounting for differences in the diapause stage (Φ = 0 for *Aedes* species and Φ = 1 for *Culex* species) and *m*<sub>A</sub><sup>dia</sup> represent the adult mortality rate during the diapause period, decreasing to a minimal value for the *Culex* species.

$$
\left\{
\begin{aligned}
\text{Aquatic stages} \\
\dot{E} &= \gamma_{A_o} (\beta_1 A_{1o} + \beta_2 A_{2o}) - \Phi f_E E - \mu_E E \\
\dot{L} &= \Phi f_E E - \left(m_L \left(1 + \frac{L}{\kappa_L}\right) + f_L\right) L \\
\text{Aerial stages} \\
\dot{A}_{em} &= f_p P \sigma_{em} \exp \left( - \mu_{em} \left( 1 + \frac{P}{\kappa_P} \right) \right) - (m_A^{dia} + (1 - \Phi) \gamma_{A_{em}}) A_{em} \\
\dot{A}_{1h} &= (1 - \Phi) \gamma_{A_{em}} A_{em} - (m_A + \mu_r + \gamma_{A_h}) A_{1h}
\end{aligned}
\right.
\qquad (2)
$$

Considering the absence of an effective diapause period for the five
mosquito species in the Maghreb region, and the yearlong presence of
larvae in habitats (Berchi et al., 2012; El Ouali Lalami et al., 2009;
Himmi, 2007), only one set of equations (Eqs. 1a and 1b) was used for
this area. The model was implemented in R 4.1.2 (R Core Team, 2020),
in discrete time, using a daily time step.

2.3. Model parametrization

2.3.1. *Constant parameters and meteorological-dependent functions*

All constant parameters and meteorological-dependent functions
were obtained from a literature review and are provided in Table 1.
When information was missing for a given species, parameters or
functions were extrapolated from another species of the same genus.
Several models have been developed to represent temperature-
dependent development rates in insects, but the Logan model, as
described by Eq. 3, is one of the most used (Logan et al., 1976; Rebaudo
and Rabhi, 2018). This nonlinear model allows to describe the increase
in development rate between the minimal temperature and the opti-
mum, and then the rapid decrease between this optimum and the
maximal temperature.

$$
r(\theta) = A \left( e^{B(\theta - \theta_{\min})} - e^{B(\theta_{\max} - \theta_{\min})} \frac{\theta_{\max} - \theta}{C} \right) \quad (3)
$$

In Eq. 3, *r* is the temperature-dependent rate, *θ* the temperature and *θ*<sub>min</sub>-*θ*<sub>max</sub> the range of temperature of possible development. *A*, *B* and *C* are shape parameters. In our study, we attributed values to the parameters of the Logan model for each temperature-dependent rate, based on experimental data when they were available, using the range of temperature for which a development was observed for *θ*<sub>min</sub> and *θ*<sub>max</sub>, and calibrating *A*, *B* and *C* using a nonlinear least-squares method (Table 1). Parameters of Logan models were estimated for *f*<sub>*E*</sub> (Cx. *theileri*), *f*<sub>*L*</sub> (Ae. *vexans*, Cx. *pipiens* and Cx. *theileri*) and *f*<sub>*P*</sub> (Ae. *vexans*, Cx. *pipiens* and Cx. *theileri*), and were extracted from previous published parametrizations of the generic model for the other functions.

2.3.2. Parametrization of the carrying capacities

In each pixel, the model was run for all five species, using the pixel-
specific meteorological data as inputs. Considering the size of each pixel
(around 9 km²) and the active flight capacities of both *Aedes* and *Culex*
mosquitoes (Becker et al., 2020; Bogojević et al., 2011), we ignored their
diffusion between spatial units. Each model run thus described the dy-
namics of a closed mosquito population for a given location within a
given pixel.

The carrying capacity can be defined as the ‘maximum population
size that can be supported indefinitely by a given environment’ (Hixon,
2008). This parameter is used in ecological models to limit the popu-
lation size. Mathematically, it allows the population growth rate to
decrease as the population size approaches the carrying capacity, and
the population size to decline if it exceeds it (Hartvigsen, 2013). In a
given pixel of our model, this carrying capacity was represented for
larvae and pupae using a density-dependent mortality in the larval stage
(controlled by parameter κ<sub>L</sub> in Eqs. 1 and 2), and a density-dependent
transition success from pupae to emerging adults (controlled by
parameter κ<sub>P</sub> in Eqs. 1 and 2). We assumed that these carrying capacities
could vary between spatial units, and we defined two pixel-specific
parameters: κ<sub>Lxy</sub> and κ<sub>Pxy</sub> for a given pixel (x,y), for larvae and pupae,
respectively. For the aquatic stages, the carrying capacity reflects the
availability of breeding sites. However, it was not possible to assess the
presence of suitable breeding sites for the five mosquito species in the
whole region of interest, and we thus used the probability of presence
data from Wint and colleagues (Wint et al., 2020) as a proxy. To define
the mathematical relation between the pixel carrying capacities κ<sub>Lxy</sub> and
κ<sub>Pxy</sub> and the probability of presence, we compared three methods
(M1-3):

- M1: the carrying capacities were identical for all pixels, i.e., $\kappa_{Lxy} = \kappa_{Lmax}$ and $\kappa_{Pxy} = \kappa_{Pmax}$, $\kappa_{Lxy}$ and $\kappa_{Pxy}$ being parameters of the model (Table 1). Meteorological data were thus the sole source of geographic variation of mosquito abundance.

- M2: for a given pixel (x,y), the carrying capacity was proportional to the probability of presence $p_{x,y}$ for the considered mosquito species, provided by Wint and colleagues (Wint et al., 2020), i.e., $\kappa_{Lxy} = p_{x,y} \times \kappa_{Lmax}$ and $\kappa_{Pxy} = p_{x,y} \times \kappa_{Pmax}$.

- M3: for a given pixel (x,y), the carrying capacity was proportional to a power function of $p_{x,y}$, i.e., $\kappa_{Lxy} = \kappa_{Lmax}^{(p_{x,y})^a}$ and $\kappa_{Pxy} = \kappa_{Pmax}^{(p_{x,y})^a}$, with $a$ being a shape parameter. $a$ was calibrated using entomological data, by minimizing the squares of distances between the observed data and rescaled predicted abundances (i.e., rescaled between zero and the maximal observed value) at the same dates and locations (see the ‘Parametrization of the carrying capacities’ section in Supplementary Information). Because the biology of *Aedes* and *Culex* mosquitoes are different, we hypothesized that the relationship between the carrying capacities and the probability of presence was different between these genera. We thus computed one value of $a$ per genus.

We compared M1, M2 and M3 using the datasets for which data were
observed in pixels with at least 10 different values of probability of
presence. We studied the distribution of distances between the rescaled predicted abundance and observed data for each unique value of the probability of presence, and selected the method for which these distances were the smallest.

## 2.4. Model fit assessment

Using the selected method to parametrize the larvae and pupae carrying capacities, and after initialization (see the ‘Model initialization’ section in Supplementary information), we ran the model in all pixels during the whole period, i.e., from 01/01/1990 to 31/12/2021, with a daily time step, and created monthly maps by computing the mean. To assess the model fit, we first computed the percentage of agreement between observed data and model predictions for the low versus high abundances of mosquitoes, based on an arbitrary threshold of the abundance peak size of 10 mosquitoes (Supplementary Fig. S5). We then evaluated the correlation between the observed and predicted entomological abundances at daily, weekly, biweekly and monthly time steps, by computing the corresponding means of observed and predicted abundances and their Spearman coefficients. To avoid evaluating correlation for observed time series with no mosquitoes trapped, or for these with few counts, we arbitrary selected the time series with at least 10 trapping results and 5 nonzero values for this analysis (Supplementary Table S1).

## 2.5. Sensitivity analysis

To assess the influence of parameters and functions on the model outputs, we used the ‘one-factor-at-a-time’ Morris screening method (Morris, 1991), implemented in the dedicated function of the R package ‘sensitivity’ (Iooss et al., 2024). This method requires few simulations but does not allow to evaluate parameters interaction. The Morris method is based on the definition of a variation space for each

parameter, and proceeds by sampling trajectories in this space, each point of a trajectory defining a set of parameters used for the model simulations. Variations of the outputs (‘elementary effects’) are recorded for each simulation, and two sensitivity indicators are computed: $\mu^*$ which is the mean of the absolute values of the elementary effects, and $\sigma$ their standard deviation. High values of $\mu^*$ indicate an important effect of the parameter on the output variable, whereas high values of $\sigma$ indicate a nonlinear effect on the output. We used the relative variation of the maximum mosquito abundance during the last year of simulation, i.e., 2021, as the model output. We discretized the variation of each parameter using 6 levels within a range of $\pm 10\%$ around the baseline value and generated 10 trajectories in the parameter space thus defined. Meteorological-dependent functions were assessed by introducing a multiplicative parameter taking a baseline value of one and performing the sensitivity analysis on this parameter. Because model equations differed between Europe and Africa, we conducted the sensitivity analysis for each mosquito species separately for the two continents. Due to the high computational time needed (one run of the model was up to 10 hours long), we randomly selected 30 pixels in Europe and 30 in Africa, for each mosquito species, among pixels with a high probability of presence ($> 0.8$) and a predicted peak of at least 1000 mosquitoes during 2021. The sensitivity analysis was performed using the outputs obtained in these 60 pixels.

## 3. Results

### 3.1. Parametrization of the carrying capacities

M1, M2 and M3 were evaluated using SP1 and IT1 datasets only, because in both cases data were available for each of the five species and were observed in pixels with a large range of probabilities of presence (for instance, time series were available in pixels with 42 different

<figure><img src="image_6.png" /><figcaption>Fig. 2. Predictions of mosquito abundances for 2021 (monthly means of host-seeking adults ($A_{1h} + A_{2h}$))</figcaption></figure>
values of the probability of presence for *Ae. caspius* in IT1 data, ranging from 0.07 to 0.69, and in pixels with 506 different values in SP1 data, ranging from 0.08 to 0.97; other values are presented in Supplementary Figs. S8 to S17). Overall, the differences between the rescaled predicted data and the observed data were lower for M3 than for M1 and M2, and this method using a power link and a calibrated $\alpha$ value was selected to define the pixel-specific carrying capacities for the five mosquito species. The differences between the rescaled predicted data and the observed data were higher for *Culex* species, and they increased for pixels with a higher probability of presence for the five mosquito species.

## 3.2. Model predictions

We used the model to produce daily predictions of abundance for the five mosquito species, in each pixel of the studied area, from 1990 to 2021. According to predictions, *Ae. caspius* and *Ae. detritus* distributions were restricted to the coasts, the valleys of large rivers (e.g., the Tagus in Portugal, the Guadalquivir in Spain, the Po Valley in Italy or the Sebou in Morocco), the Gulf of Lion in France, the gulfs of Hammamet and Tunis in Tunisia, Chott Melrhir and Chott el Djerid in Algeria and Tunisia, respectively, and in sabkha areas in Tunisia (Fig. 2 and Supplementary videos S1 to S5). Predicted *Ae. caspius* abundance was high in these areas from May to November in Europe, from May to June and from October to February in North Africa. According to predictions, *Ae. vexans* mosquitoes were mainly located in the Po Valley region. In contrast, predicted *Culex* populations had a broad distribution, in the whole area of interest for *Cx. pipiens* and mostly in the Iberian Peninsula and in the Maghreb region for *Cx. theileri*. Predicted *Culex* abundance was particularly high from April to October. As expected, and for all five species, the model predicted fewer mosquitoes during winter than during the rest of the year.

## 3.3. Model fit assessment

The model was able to predict a low abundance of *Aedes* mosquitoes in most of the pixels with low observed trapping values, and conversely to predict high abundances in pixels with high observed trapping values. The proportion of agreement between predicted and observed time series reached 79 %, 75 % and 92 % for *Ae. caspius*, *Ae. detritus* and *Ae. vexans*, respectively, taking an arbitrary threshold of 10 individuals at peak to define low and high abundances (Table 2). However, for the *Culex* species, the predictions reproduced the observed values with less precision, showing an agreement of 33 % and 4 % for *Cx. pipiens* and *Cx. theileri*, respectively. Overall, the model over-predicted the abundance of *Culex* species.

<caption>**Table 2**</caption>

Number of time series for which there was a concordance between the peak of observed and predicted abundances (weekly time step).

<table><thead><tr><th rowspan="2"></th><th rowspan="2">Maximum observed value</th><th colspan="2">Maximum predicted abundance</th></tr><tr><th>&lt;10</th><th>&ge;10</th></tr></thead><tbody><tr><td rowspan="2"><em>Ae. caspius</em></td><td>&lt;10</td><td>502</td><td>138</td></tr><tr><td>&ge;10</td><td>5</td><td>41</td></tr><tr><td rowspan="2"><em>Ae. detritus</em></td><td>&lt;10</td><td>497</td><td>168</td></tr><tr><td>&ge;10</td><td>0</td><td>20</td></tr><tr><td rowspan="2"><em>Ae. vexans</em></td><td>&lt;10</td><td>615</td><td>41</td></tr><tr><td>&ge;10</td><td>12</td><td>14</td></tr><tr><td rowspan="2"><em>Cx. pipiens</em></td><td>&lt;10</td><td>0</td><td>418</td></tr><tr><td>&ge;10</td><td>0</td><td>202</td></tr><tr><td rowspan="2"><em>Cx. theileri</em></td><td>&lt;10</td><td>0</td><td>591</td></tr><tr><td>&ge;10</td><td>0</td><td>24</td></tr></tbody></table>

The values represent the number of time series meeting the conditions of the cell. Bold cells highlight an agreement between observed and predicted time series, i.e., the time series for which the mosquito abundance was low or high for both the predicted and observed data.

Considering only the observed time series with at least 10 observations and at least 5 nonzero values of mosquito counts to assess the correlation with predicted abundance (weekly mean values), the median value of the Spearman correlation coefficients was of 0.15, 0.32, 0.13, 0.31 and 0.05, for *Ae. caspius*, *Ae. detritus*, *Ae. vexans*, *Cx. pipiens* and *Cx. theileri*, respectively (Fig. 3). The percentage of time series for which the correlation coefficient was greater than 0.2 reached 78 % and 69 % for *Ae. detritus* and *Cx. pipiens*, respectively. The median value of the correlation coefficient was similar when using daily predicted abundance, as well as when using biweekly and monthly mean values (Supplementary Table S2 and Supplementary Figs. S18 to S20).

Fig. 4 provides three selected examples of weekly observed and predicted time series from the IT1, FR1 and FR4 datasets. In these examples, the model was able to reproduce *Cx. pipiens* activity periods, as well as abundance peaks of *Ae. caspius* and *Ae. detritus* when they followed precipitations periods of the meteorological data.

## 3.4. Sensitivity analysis

Due to the low probability of presence for *Ae. vexans* and *Cx. pipiens* in North Africa (the maximum values being respectively 0.66 and 0.79), the threshold used to randomly select the 30 pixels was set to 0.5 and 0.6 for these two species, respectively, instead of the 0.8 threshold used for Europe and for the other species (Supplementary Figs. S21 to S24). Moreover, we set the threshold of peak of abundance to 1 mosquito for *Ae. vexans* in North Africa, as it was not possible to identify pixels with an annual peak greater than 1000 mosquitoes in 2021. Results of the sensitivity analysis for the 30 pixels in Europe and 30 in Africa (Supplementary Figs. S25 to S29) are graphically represented in Fig. 5 for *Ae. caspius* and *Cx. pipiens* ($\mu^*$ values) and in Supplementary Figs. S30 to S32 for the three other species. In these figures, a red cell indicates that the parameter was influential on the predictions of the corresponding pixel.

Table 3 shows that 11 of the 22 parameters included in the sensitivity analysis were influential ($\mu^* > 0.1$) for at least one mosquito species, for the majority of pixels of at least one of the two regions. The sex ratio at emergence $\sigma_{em}$ was influential for all species except for *Cx. pipiens*, in both regions. The mortality rate during emergence $\mu_{em}$ and the maximum environment carrying capacity for pupae $\kappa_{Pmax}$ were also influential for all species, in both regions, except for *Ae. vexans* in North Africa. Results for nonlinear effects are presented in Supplementary Figs. S33 to S37. These effects were mainly observed for the favorable season start date for *Ae. caspius*.

# 4. Discussion

## 4.1. Summary of the results

The model was able to simulate population dynamics of five potential RVFV vector species in the western Mediterranean Basin, between 1990 and 2021. Predicted abundances of *Aedes* mosquitoes were restricted to specific regions suitable for their development, whereas *Culex* populations were widespread in the Mediterranean Basin. When comparing the predictions with mosquito trap data from Spain, France, Italy and Morocco, the model was able to predict low abundances of *Aedes* mosquitoes in the regions where few mosquitoes were trapped, but tended to overpredict *Culex* populations. The model was also able to reproduce most of the peaks of mosquito abundance.

The model used here was adapted from the work of Cailly and colleagues (Cailly et al., 2012) and based on a generic framework that had already been used with satisfying results in Botswana (Hammami et al., 2016), France (mainland and Reunion Island) (Ezanno et al., 2015; Tran et al., 2020, 2013b), Mauritius (Iyaloo et al., 2021) and Senegal (Tran et al., 2019) for several species belonging to *Aedes*, *Anopheles* and *Culex* mosquito genera.
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Distribution of the Spearman correlation coefficient between observed and predicted abundances (temporal correlation — weekly mean abundances). N represents the number of time series included (having at least 10 observations and 5 values greater than 0). The dotted lines indicate the median values.</figcaption></figure>

## 4.2. Availability of entomological data

To spatialize the model at the western Mediterranean Basin scale, we used the probabilities of presence computed by Wint and colleagues (Wint et al., 2020) for each of the five mosquito species. Based on entomological datasets from four different countries, we then calibrated the model and assessed its outputs. Datasets from Spain (SP1) and Italy (IT1) provided time series observed in pixels with a large number of different probabilities of presence values (for example, 42 and 506 for *Ae. caspius* and *Cx. pipiens*, respectively), ranging in large intervals (0.07 to 0.69 and 0.08 to 0.97, respectively) for the five mosquito species. However, in the Maghreb region, data were available from Morocco only, with a reduced number of associated probabilities of presence (6 values for *Cx. pipiens* in MO1, 3 in MO2 and 2 in MO3). Moreover, these available time series mainly concerned larvae (MO1). The inclusion of additional data from this part of the Mediterranean Basin would improve the goodness-of-fit of the model, with a particular need for trapping data for *Cx. theileri*, an abundant species in North Africa.

## 4.3. Quality of temporal predictions of the model

Predictions of the model were able to reproduce most of the periods of mosquito activity (Fig. 4), but some discrepancies remained. *Ae. detritus* field populations often decrease during summer months, as observed in most of the datasets used here and in the literature (Clarkson and Enevoldson, 2020; Veronesi et al., 2012). However, this phenomenon was not well reproduced by the model, at least for some areas such as the Gulf of Lion region, the coast of the Balearic Sea in Spain or the Naples region in Italy (Fig. 2 and Supplementary videos). In addition, our assumption of an absence of diapause in Maghreb may be subject to local variations, as such a period of inactivity has been observed for *Ae.*

*caspius* in the Algerian Sahara (Boubidi, unpublished data). On the other hand, some *Aedes* peaks of abundance were not predicted (Figs. 4b and 4c). This is likely due to the hatching of these species being only triggered by rainfall in the model, as a proxy of the increase in water levels and the subsequent flooding of eggs. However, other processes may induce an increase in water levels, whether natural (as the rise of the water table level or marine flooding due to storms occurring during high tide) or artificial (as flooding for cultivation in rice fields, hunting or fauna protection in marshes) (Balenghien et al., 2006; Ben Ayed et al., 2019; Hammadi et al., 2009). Due to the large size of the study area and to its environmental diversity, it was not possible to account for these specific and local mechanisms, as it has been done in the Camargue region (France) for instance (Ezanno et al., 2015). Therefore, *Aedes* abundance peaks resulting from flooding that were not linked to rainfall could not be predicted by our model. Moreover, egg hatching was triggered by rainfall with a threshold set to 15 mm of precipitations per day in our model. In larval developmental sites, the majority of *Aedes* eggs seems to be laid at the same height above the water, at an average level between the highest and lowest water levels (Gabinaud, 1975), and should therefore hatch at the same time. However, it is likely that the emergence of mosquito populations can be triggered by lower rainfall, depending on local conditions. The use of *in situ* meteorological data, reported by local meteorological stations, may foster a stronger correlation between observed and predicted abundance, as demonstrated by Tran and colleagues (Tran et al., 2019), and would reduce the uncertainty linked to the input data, but such information was not homogeneously available for the entire zone of interest. The next step to improve our model could be to use a hydrological model simulating the variations in water levels, as used in Senegal at the pond scale (Soti et al., 2012) or more broadly (Tran et al., 2019). However, this would be difficult to achieve across the entire region of interest, and it would not
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Examples of observed and predicted abundances of host-seeking adults (*A*<sub>1h</sub> + *A*<sub>2h</sub>), at a weekly time step. **a.** *Culex pipiens*, Campli, Teramo (time series from IT1 data, from 2019 to 2021). **b.** *Aedes caspius*, Tour Carbonnière — Camargue wetland (time series from FR1 data, 2006). **c.** *Aedes vexans*, Tour du Valat — Camargue wetland (time series from FR4 data, 2004). Mean weekly precipitations are displayed for *Aedes* mosquitoes.</figcaption></figure>

solve the issue of water level increases that are not linked to rainfall.

4.4. Quality of spatial predictions of the model

Spatial predictions of abundance were good for *Aedes* species.
However, the model tended to over-predict the abundance of *Culex*
mosquitoes, in particular of *Cx. theileri*, as in Tunisia (Bouattour, un-
published data), or in Italy where our model predicted a high abundance
of this species whereas a few individuals were recorded in the IT1

datasets at the same time and locations (De Ascentis et al., 2022). The
time series with the highest observed abundance of *Culex* mosquitoes
were located in pixels with low to medium probabilities of presence,
leading to the computation of a low value of *α*. Therefore, the carrying
capacities were relatively high even in pixels with a low probability of
presence and the model predicted important *Culex* populations in these
pixels. This resulted in less contrasted predicted abundances than for the
*Aedes* mosquitoes, and to a lower dependency of predictions on the
probability of presence data computed by Wint and colleagues (Wint
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Graphic representation of the results of the sensitivity analysis for *Ae. caspius* and *Cx. pipiens*, performed on a random sample of pixels for both Europe and North Africa. The y-axis represents the parameters included in the sensitivity analysis, and the x-axis the 60 pixels in which the analysis was conducted (the left panels represent the 30 European pixels and the right panels the 30 North African pixels). Each cell is colored according to the corresponding value of μ* (a darker color indicates a higher value, and therefore a higher influence of the corresponding parameter on the outputs, for this pixel).</figcaption></figure>

<table><caption>**Table 3**
Influential parameters (μ* > 0.1) for the five mosquito species.</caption><thead><tr><th></th><th><em>Ae. caspius</em></th><th><em>Ae. detritus</em></th><th><em>Ae. vexans</em></th><th><em>Cx. pipiens</em></th><th><em>Cx. theileri</em></th></tr></thead><tbody><tr><td>Europe</td><td>σ<sub>em</sub>, μ<sub>em</sub>, κ<sub>Pmax</sub>,<br>p<sub>fj</sub>, start*</td><td>σ<sub>em</sub>, μ<sub>em</sub>, κ<sub>Pmax</sub>,<br>p<sub>fj</sub></td><td>σ<sub>em</sub>, γ<sub>Ah</sub>, μ<sub>em</sub>,<br>κ<sub>Pmax</sub>, P<sub>fj</sub></td><td>γ<sub>Ah</sub>, μ<sub>em</sub>, κ<sub>Pmax</sub>,<br>p<sub>fj</sub>, P<sub>fj</sub></td><td>σ<sub>em</sub>, μ<sub>A</sub>, μ<sub>em</sub>,<br>κ<sub>Pmax</sub>, p<sub>fj</sub>, p<sub>mA</sub>,<br>end, start</td></tr><tr><td>North Africa</td><td>σ<sub>em</sub>, γ<sub>Ah</sub>, μ<sub>em</sub>,<br>κ<sub>Pmax</sub>, P<sub>fj</sub></td><td>σ<sub>em</sub>, μ<sub>em</sub>, κ<sub>Pmax</sub>,<br>p<sub>fj</sub></td><td>σ<sub>em</sub>, p<sub>fj</sub>, p<sub>mL</sub>, p<sub>mA</sub></td><td>γ<sub>Ah</sub>, μ<sub>em</sub>, κ<sub>Pmax</sub>, p<sub>fj</sub>, p<sub>fj</sub></td><td>σ<sub>em</sub>, γ<sub>Ah</sub>, μ<sub>em</sub>,<br>κ<sub>Pmax</sub>, p<sub>fj</sub>, p<sub>fj</sub></td></tr></tbody></table>

p<sub>fj</sub>: multiplicative parameter of the transition rate from larvae to pupae, p<sub>fj</sub>: multiplicative parameter of the transition rate from pupae to emerging adults, σ<sub>em</sub>: sex ratio at emergence, γ<sub>Ah</sub>: transition rate from host-seeking to engorged adults, p<sub>m</sub>: multiplicative parameter of the larval mortality rate, p<sub>mA</sub>: multiplicative parameter of the adult mortality rate, μ<sub>A</sub>: minimum adult mortality rate, μ<sub>em</sub>: mortality rate during emergence, κ<sub>Pmax</sub>: maximum environment carrying capacity for pupae, start: first day of the favorable season, end: last day of the favorable season. *highly influential in some pixels only.

et al., 2020). However, these results were expected to a certain extent, as *Cx. pipiens* and *Cx. theileri* are widespread and ubiquitous species, in the whole Mediterranean Basin and in southern Spain and North Africa, respectively (Amara Korba et al., 2016; Becker et al., 2020; Trari, 2017). Finally, very few data were available for *Cx. theileri*, explaining the poor quality of predictions. Using a different functional link between the carrying capacity and the probability of presence, rather than the power function selected here, may potentially improve the model predictions for these species, and further investigations are needed. In Senegal, the availability of larval developmental sites has been assessed using

satellite images in Senegal (Tran et al., 2019), but this approach would be difficult to achieve in our study, considering the size of the region of interest and the diversity of biotopes of the five mosquito species considered. Computation of the correlation between predicted abundance and observed data at a given date and for different pixels would have been an alternative, appropriate way to evaluate the model's ability to reproduce the geographic distribution of the RVFV vectors. However, the trapping conditions have a major impact on the number of mosquitoes caught, and the spatial resolution used in this study (around 9 km²) may have led to high variability of settings for the trapping area.
Thus, data may not be comparable between different traps for a given date. Correlation between traps would not have been properly assessed in this case, because we would have not been able to separate the effects of spatial correlation from trapping conditions.

Finally, in the absence of sufficient data, we parametrized the carrying capacity of Ae. detritus using the probability of presence of Ae. caspius computed by Wint and colleagues (section ‘Probability of presence data’ of the Supplementary Information), assuming that these two species shared the same habitats. However, although both species are found together in some larval developmental sites, immature stages of Ae. caspius tolerate waters with lower saline levels than Ae. detritus, and thus have a higher broader geographic distribution (Ben Ayed et al., 2019; Gabinaud, 1975; Trari, 2017). Specific distribution data on Ae. detritus would allow to refine the predictions.

### 4.5. Sources of uncertainty in the model

Despite the efforts that have already been made to parametrize the model (Cailly et al., 2012; Ezanno et al., 2015; Hammami et al., 2016; Tran et al., 2019) and supplemented here, uncertainties remain in our study, linked to the input data, or to the model structure itself. Indeed, the spatialized data we used, whether meteorological (Munoz Sabater, 2019) or providing the probability of mosquito presence in each pixel (Wint et al., 2020), were extracted from single sources, guaranteeing the same quality of data for the whole study area, and thus the comparability of predictions across pixels. However, they were produced using a combination of field observation and modelling tools, and thus include uncertainty. Our model was run independently in each pixel, without considering diffusion between spatial units, and the differences in predicted dynamics were thus only linked to the meteorological data used as inputs, and to the carrying capacities of the aquatic stages. However, specific biological information is still lacking for some parameters and functions (Table 1) especially for Ae. detritus and Cx. theileri. Inclusion of such data would improve the model performances for these species. As highlighted by the sensitivity analysis we conducted, 11 of the 22 parameters exerted measurable influence on the model outputs. A sensitivity analysis performed on the whole results could have led to a more precise identification of the influent parameters. On the other side the introduction of stochasticity could have been a way of reducing the dependence of outputs on parameter values. However, these two solutions have been discarded due to computational time needed to perform multiple simulations. Multi-input-multi-output modelling approaches recently used in epidemiology could also be explored in this context (Tutsoy, 2023).

### 4.6. Implication of our results and future avenues

A previous study showed that Ae. caspius would be the most competent vector for the RVFV should it spread through the Mediterranean Basin (Drouin et al., 2022). According to our model results, the risk of circulation of the virus linked to the predicted presence of this vector could be high, but mainly circumscribed to particular regions, and limited to certain periods of mosquito activity (from spring to fall in Europe and during spring and from fall to winter in North Africa). Regarding the moderately competent species Ae. vexans, we observed a large zone of potential risk of RVFV transmission in the Po Valley in Italy. In the same study (Drouin et al., 2022), Cx. pipiens was also identified as a moderately competent species but, according to our predictions, this vector is widespread in the whole area of interest, and particularly abundant from spring to fall. Therefore, the risk of transmission due to this species remains to be evaluated. Finally, the transmission of RVFV in North Africa may occur all year long, mainly due to the permanent presence of Cx. theileri.

Modelling approaches have proven to be able to predict vector population dynamics for various mosquito species (Ahumada et al., 2004; Balenghien et al., 2010; Schaeffer et al., 2008; Cailly et al., 2012). Entomological model outputs can be integrated in transmission models,

as already applied for RVF in Senegal to create maps of epidemic potential of the virus, and to disentangle the mechanisms underlying inter-epidemic viral recurrence (Cecilia et al., 2020; Durand et al., 2020). These models are also useful for optimizing surveillance and control actions, as demonstrated by the development of the ‘ALBORUN’ mapping tool (Tran et al., 2020). This forecasting model, based on meteorological data, is used by the vector control service of the French Regional Health Agency on Reunion Island to monitor the abundance of *Aedes albopictus* mosquitoes and the associated risk of chikungunya and dengue outbreaks. However, vector presence and their competence are only two of the several parameters influencing the risk of circulation of RVFV in a given area. To improve these predictions, other aspects of mosquito biology, such as the trophic preferences of these five species, and the presence and densities of susceptible livestock species, should be considered.

The evolving epidemiological situation in North Africa calls for efforts to prepare for the introduction of the virus into the western part of the Mediterranean Basin, including Southern Europe, where the risk of introduction has been considered as low but not null by the European Food Safety Authority (EFSA). Overall, the model presented here is a first step to assess the risk of RVFV transmission in the region, and results can be used to target surveillance activities, especially regions of high abundance of Ae. caspius, which has a more circumscribed geographical distribution, and which has been identified as the most competent mosquito species for RVFV transmission in the Mediterranean Basin (Drouin et al., 2022).

## Code accessibility

Entomological model R code is available on the following GitHub repository: https://github.com/drouinalex/spatial_modelling_rvf{vectors}

## CRediT authorship contribution statement

**Alex Drouin:** Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Thomas Balenghien:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Conceptualization. **Benoit Durand:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Conceptualization. **Carles Aranda:** Writing – review & editing, Investigation, Data curation. **Amal Bennouna:** Writing – review & editing, Investigation, Data curation. **Ali Bouattour:** Writing – review & editing, Investigation, Data curation. **Said C Boubidi:** Writing – review & editing, Validation, Investigation, Data curation. **Annamaria Conte:** Writing – review & editing, Investigation, Data curation. **Sarah Delacour:** Writing – review & editing, Investigation, Data curation. **Maria Goffredo:** Writing – review & editing, Investigation, Data curation. **Oumnia Himmi:** Writing – review & editing, Investigation, Data curation. **Grégory L’Ambert:** Writing – review & editing, Investigation, Data curation. **Francis Schaffner:** Writing – review & editing, Investigation, Data curation. **Véronique Chevalier:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Conceptualization.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Alex Drouin reports financial support was provided by French Ministry of Agriculture. Other authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
## Acknowledgement

Muñoz Sabater, (2019) was downloaded from the Copernicus Climate Change Service (C3S) Climate Data Store. The results contain modified Copernicus Climate Change Service information 2020. Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data it contains.

We are thankful to the Spanish Ministry of Agriculture, Fishery and Food (MAPA) for allowing access to data from the Spanish National Surveillance control and eradication of bluetongue program and to the French Ministry of Agriculture and Food, which funded this research. We are grateful to Pauline Ezanno, Annelise Tran, William Wint and François Rebaudo for their help in parametrizing the model.

## Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ecolmodel.2024.111013.

## Data availability

Published datasets analysed during the current study are available from authors of the associated articles, and other data are available from C.A., A.C., S.D. and M.G., upon reasonable request.

## References

* Ahumada, J.A., Laoointe, D., Samuel, M.D., 2004. Modeling the population dynamics of *Culex quinquefasciatus* (Diptera: Culicidae), along an elevational gradient in Hawaii. J. Med. Entomol. 41, 1157-1170. https://doi.org/10.1603/0022-2585-41.6.1157.

* Amara Korba, R., Alayat, M.S., Bouiba, L., Boudrissa, A., Bouslama, Z., Boukraa, S., Francis, F., Failloux, A.-B., Boubidi, S.C., 2016. Ecological differentiation of members of the *Culex pipiens* complex, potential vectors of West Nile virus and Rift Valley fever virus in Algeria. Parasit. Vectors 9, 455. https://doi.org/10.1186/s13071-016-1725-9.

* Amdouni, J., Conte, A., Ippoliti, C., Candeloro, L., Tora, S., Sghaier, S., Hassine, T.B., Fakhfekh, E.A., Savini, G., Hammami, S., 2022. *Culex pipiens* distribution in Tunisia: identification of suitable areas through Random Forest and MaxEnt approaches. Vet. Med. Sci. 8, 2703-2715. https://doi.org/10.1002/vms3.897.

* Arsevska, E., Hellal, J., Mejri, S., Hammami, S., Marianneau, P., Calavas, D., Hénaux, V., 2016. Identifying areas suitable for the occurrence of Rift Valley fever in North Africa: implications for surveillance. Transbound. Emerg. Dis. 63, 658-674. https://doi.org/10.1111/tbed.12331.

* Balenghien, T., Carron, A., Sinègre, G., Bicout, D.J., 2010. Mosquito density forecast from flooding: population dynamics model for *Aedes caspius* (Pallas). Bull. Entomol. Res. 100, 247-254. https://doi.org/10.1017/S0007485309990745.

* Balenghien, T., Fouque, F., Sabatier, P., Bicout, D.J., 2006. Horse-, bird-, and human-seeking behavior and seasonal abundance of mosquitoes in a West Nile virus focus of southern France. J. Med. Entomol. 43, 936-946. https://doi.org/10.1093/jmedent/43.5.936.

* Becker, N., Petrić, D., Zgomba, M., Boase, C., Madon, M.B., Dahl, C., Kaiser, A., 2020. Mosquitoes: identification, ecology and control. Fascinating Life Sciences, 3rd ed. Springer International Publishing, Switzerland. https://doi.org/10.1007/978-3-030-11623-1.

* Bellone, R., Failloux, A.-B., 2020. The role of temperature in shaping mosquito-borne viruses transmission. Front. Microbiol. 11. https://doi.org/10.3389/fmicb.2020.584846, 584846:1-584846:15.

* Ben Ayed, W., Amraoui, F., M'ghirbi, Y., Schaffner, F., Rhaim, A., Failloux, A.-B., Bouattour, A., 2019. A survey of *Aedes* (Diptera: Culicidae) mosquitoes in Tunisia and the potential role of *Aedes detritus* and *Aedes caspius* in the transmission of Zika virus. J. Med. Entomol. 56, 1377-1383. https://doi.org/10.1093/jme/tjz067.

* Berchi, S., Aouati, A., Louadi, K., 2012. Typologie des gîtes propices au développement larvaire de *Culex pipiens* L. 1758 (Diptera-Culicidae), source de nuisance à Constantine (Algérie). Ecol. Mediterr. 38, 5-16. https://doi.org/10.3406/ecmed.2012.1311.

* Bogojević, M.S., Merdić, E., Bogdanović, T., 2011. The flight distances of floodwater mosquitoes (*Aedes vexans*, *Ochlerotatus sticticus* and *Ochlerotatus caspius*. Osijek, Eastern Croatia. Biologia (Bratisl.) 66, 678-683. https://doi.org/10.2478/s11756-011-0073-7.

* Cailly, P., Tran, A., Balenghien, T., L'Ambert, G., Toty, C., Ezanno, P., 2012. A climate-driven abundance model to assess mosquito control strategies. Ecol. Model. 227, 7-17. https://doi.org/10.1016/j.ecolmodel.2011.10.027.

* Cecilia, H., Drouin, A., Métras, R., Balenghien, T., Durand, B., Chevalier, V., Ezanno, P., 2022. Mechanistic models of Rift Valley fever virus transmission: a systematic review. PLoS Negl. Trop. Dis. 16, e0010339. https://doi.org/10.1371/journal.pntd.00010339.

* Cecilia, H., Métras, R., Fall, A.G., Lo, M.M., Lancelot, R., Ezanno, P., 2020. It's risky to wander in September: modelling the epidemic potential of Rift Valley fever in a Sahelian setting. Epidemics 33, 2020. https://doi.org/10.1016/j.epidem.2020.100409, 02.25.20027821.

* Clarkson, M., Enevoldson, T., 2020. The factors which influence the breeding and number of *Aedes detritus* in the Neston area of Cheshire, UK, the production of a local mosquito forecast and public bite reporting. J. Eur. Mosq. Control Assoc. 38, 17-32.

* De Ascentis, M., Quaglia, M., D'Alessio, S.G., Iapaolo, F., Pizzurro, F., Ruggeri, F., Rossi, N., Blardi, L., Ippoliti, C., Gioci, D., Portanti, O., Pisciella, M., Di Lorenzo, A., Ciarrocchi, E., Irelli, R., Conte, A., Morelli, D., Monaco, F., Savini, G., Goffredo, M., 2022. Species of mosquitoes present in Abruzzo and Molise and their possible role as vector of Usutu and West Nile viruses. Vet. Ital. 58. https://doi.org/10.12834/Vetlt.3046.20276.1.

* Di Nardo, A., Rossi, D., Saleh, S.M.L., Lejlifa, S.M., Hamdi, S.J., Di Gennaro, A., Savini, G., Thrusfield, M.V., 2014. Evidence of rift valley fever seroprevalence in the Sahrawi semi-nomadic pastoralist system, Western Sahara. BMC Vet. Res. 10 (92), 1-92. https://doi.org/10.1186/1746-6148-10-92, 9.

* Dinerstein, E., Olson, D., Joshi, A., Vynne, C., Burgess, N.D., Wikramanayake, E., Hahn, N., Palminteri, S., Hedao, P., Noss, R., Hansen, M., Locke, H., Ellis, E.C., Jones, B., Barber, C.V., Hayes, R., Kormos, C., Martin, V., Crist, E., Sechrest, W., Price, L., Baillie, J.E.M., Weeden, D., Suckling, K., Davis, C., Sizer, N., Moore, R., Thau, D., Birch, T., Potapov, P., Turubanova, S., Tyukavina, A., de Souza, N., Pintea, L., Brito, J.C., Llewellyn, O.A., Miller, A.G., Patzelt, A., Ghazanfar, S.A., Timberlake, J., Klöser, H., Shennan-Farpón, Y., Kindt, R., Lillesø, J.-P.B., van Breugel, P., Graudal, L., Voge, M., Al-Shammari, K.F., Saleem, M., 2017. An Ecoregion-based approach to protecting half the terrestrial realm. BioScience 67, 534-545. https://doi.org/10.1093/biosci/bix014.

* Drouin, A., Chevalier, V., Durand, B., Balenghien, T., 2022. Vector competence of Mediterranean mosquitoes for Rift Valley fever virus: a meta-analysis. Pathogens 11, 503. https://doi.org/10.3390/pathogens11050503.

* Ducheyne, E., Versteirt, V., Hendrickx, G., 2013. Abundance of Rift Valley fever vectors in Europe and the Mediterranean Basin. EFSA Support. Publ. 10, 420E:1-420E:25. https://doi.org/10.2903/sp.efsa.2013.EN-420.

* Dungu, B., Lubisi, B.A., Ikegami, T., 2018. Rift Valley fever vaccines: current and future needs. Curr. Opin. Virol. 29, 8-15. https://doi.org/10.1016/j.coviro.2018.02.001.

* Durand, B., Lo Modou, M., Tran, A., Ba, A., Sow, B., Belkhiria, J., Fall, A.G., Bitéye, B., Grosbois, V., Chevalier, V., 2020. Rift Valley fever in northern Senegal: a modelling approach to analyses the processes underlying virus circulation recurrence. PLoS Negl. Trop. Dis. 14, e0008009. https://doi.org/10.1371/journal.pntd.0008009.

* EFSA Panel on Animal Health and Welfare (AHAW), 2013. Scientific opinion on Rift Valley fever. EFSA J. 11 (3180), 1-3180. https://doi.org/10.2903/j.efsa.2013.3180, 48.

* El Ouali Lalami, A., Hindi, T., Azzouzi, A., Elghadraoui, L., Maniar, S., Faraj, C., Adlaoui, E.B., Ameur, I., Ibsnouda Koraichi, S., 2009. Inventaire et répartition saisonnière des Culicidae dans le centre du Maroc. Entomol. Faun. 62, 131-138.

* ESRI, 2022. World countries [WWW Document]. URL https://www.egris.com/home/item.html?id=ac80670eb213440ea5899bbf92a04998#overview (accessed 1.24.23).

* Ezanno, P., Aubry-Kientz, M., Arnoux, S., Cailly, P., L'Ambert, G., Toty, C., Balenghien, T., Tran, A., 2015. A generic weather-driven model to predict mosquito population dynamics applied to species of *Anopheles* *Culex* and *Aedes* genera of southern France. Prev. Vet. Med. 120, 39-50. https://doi.org/10.1016/j.prevetmed.2014.12.018.

* Gabinaud, A., 1975. Écologie de deux *Aedes* halophiles du littoral méditerranéen français: *Aedes (Ochlerotatus) caspius* (Pallas, 1771): *Aedes (Ochlerotatus) detritus* (Haliday, 1833), Nematocera-Culicidae: Utilisation de la vegetation comme indicateur biotique pour l'établissement d'une carte écologique: application en dynamique des populations (Thèse de doctorat). Université des Sciences et Techniques de Montpellier 2, 1970-2014.

* Gangoso, L., Aragonés, D., Martínez-de la Puente, J., Lucientes, J., Delacour-Estrella, S., Estrada Pena, R., Montalvo, T., Bueno-Marí, R., Bravo-Barriga, D., Frontera, E., Marqués, E., Ruiz-Arrondo, I., Muñoz, A., Oteo, J.A., Miranda, M.A., Barceló, C., Árias Vázquez, M.S., Silva-Torres, M.I., Ferraguti, M., Magallanes, S., Muriel, J., Marzal, A., Aranda, C., Ruiz, S., González, M.A., Morchón, R., Gómez-Barroso, D., Figueroa, J., 2020. Determinants of the current and future distribution of the West Nile virus mosquito vector *Culex pipiens* in Spain. Environ. Res. 188, 109837. https://doi.org/10.1016/j.envres.2020.109837.

* Groen, T.A., L'Ambert, G., Bellini, R., Chaskopoulou, A., Petric, D., Zgomba, M., Marrama, L., Bicout, D.J., 2017. Ecology of West Nile virus across four European countries: empirical modelling of the *Culex pipiens* abundance dynamics as a function of weather. Parasit. Vect. 10, 524. https://doi.org/10.1186/s13071-017-2484-y.

* Hammadi, D., Boubidi, S., Chaib, S.E., Saber, A., Khechache, Y., Gasmi, M., Harrat, Z., 2009. Le paludisme au Sahara algérien. Bull. Société Pathol. Exot. 102, 185-192. https://doi.org/10.3185/pathex03356.

* Hammami, P., Tran, A., Kemp, A., Tshikae, P., Kgori, P., Chevalier, V., Paweska, J., Jori, F., 2016. Rift Valley fever vector diversity and impact of meteorological and environmental factors on *Culex pipiens* dynamics in the Okavango Delta, Botswana. Parasit. Vectors 9, 434. https://doi.org/10.1186/s13071-016-1712-1.

* Hartvigsen, G., 2013. Carrying capacity, concept of. In: Levin, S.A. (Ed.), Encyclopedia of Biodiversity, Second Edition. Academic Press, Waltham, pp. 695-701. https://doi.org/10.1016/B978-0-12-384719-5.00020-4.

* Hellal, J., Mejri, S., Lacote, S., Sghaier, S., Dkhil, A., Arsevska, E., Calavas, D., Hénaux, V., Marianeau, P., Hammami, S., 2021. Rosological evidence of Rift Valley fever in domestic ruminants in Tunisia underlines the need for effective surveillance. Open Vet. J. 11, 337-341. https://doi.org/10.5455/OVJ.2021.v11.i3.1.

* Himmi, O., 2007. Les culicides (Insectes, Diptères) au Maroc : Systématique écologique et études épidémiologiques pilotes. Université Mohammed V, Faculté des Sciences, Rabat, Maroc.
Hixon, M.A., 2008. Carrying capacity. In: Jørgensen, S.E., Fath, B.D. (Eds.), Encyclopedia of Ecology. Academic Press, Oxford, pp. 528–530. https://doi.org/10.1016/B978-008045405-4.00468-7.

Iooss, B., Veiga, S.D., Janon, A., Pujol, G., 2024. sensitivity: Global Sensitivity Analysis of Model Outputs and Importance Measures. https://cran.r-project.org/web/packges/sensitivity/index.html.

Iyaloo, D.P., Degenne, P., Elahee, K.B., Lo Seen, D., Bheecarry, A., Tran, A., 2021. ALBOMAURICE: a predictive model for mapping *Aedes albopictus* mosquito populations in Mauritius. SoftwareX 13, 100638. https://doi.org/10.1016/j.softx.2020.100638.

Javelle, E., Lesueur, A., Pommier de Santi, V., de Laval, F., Lefebvre, T., Holweck, G., Durand, G.A., Leparc-Goffart, I., Texier, G., Simon, F., 2020. The challenging management of Rift Valley Fever in humans: literature review of the clinical disease and algorithm proposal. Ann. Clin. Microbiol. Antimicrob. 19 (4), 1–4. https://doi.org/10.1186/s12941-020-0346-5, 18.

Kenawy, M.A., Abdel-Hamid, Y.M., Beier, J.C., 2018. Rift Valley Fever in Egypt and other African countries: historical review, recent outbreaks and possibility of disease occurrence in Egypt. Acta Trop 181, 40–49. https://doi.org/10.1016/j.actatropica.2018.01.015.

Linthicum, K.J., Britch, S.C., Anyamba, A., 2016. Rift Valley fever: an emerging mosquito-borne disease. Annu. Rev. Entomol. 61, 395–415. https://doi.org/10.1146/annurev-ento-010715-023819.

Logan, J.A., Wollkind, D.J., Hoyt, S.C., Tanigoshi, L.K., 1976. An analytic model for description of temperature dependent rate phenomena in arthropods. Environ. Entomol. 5, 1133–1140. https://doi.org/10.1093/ee/5.6.1133.

Lumley, S., Horton, D.L., Hernandez-Triana, L.L.M., Johnson, N., Fooks, A.R., Hewson, R., 2017. Rift Valley fever virus: strategies for maintenance, survival and vertical transmission in mosquitoes. J. Gen. Virol. 98, 875–887. https://doi.org/10.1099/jgv.0.000765.

Morris, M.D., 1991. Factorial sampling plans for preliminary computational experiments. Technometrics 33, 161–174. https://doi.org/10.2307/1269043.

Mughini-Gras, L., Mulatti, P., Severini, F., Bocolini, D., Romi, R., Bongiorno, G., Khoury, C., Bianchi, R., Montarsi, F., Patregnani, T., Bonfanti, L., Rezza, G., Capelli, G., Busani, L., 2014. Ecological niche modelling of potential West Nile virus vector mosquito species and their geographical association with equine epizootics in Italy. EcoHealth 11, 120–132. https://doi.org/10.1007/s10393-013-0878-7.

Muñoz Sabater, J., 2019. ERA5-Land hourly data from 1981 to present [WWW Document]. Copernic. Clim. Change Serv. Clim. Data Store. URL https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=overview. accessed 1.18.23.

Nielsen, S.S., Alvarez, J., Bicout, D.J., Calistri, P., Depner, K., Drewe, J.A., Garin-Bastuji, B., Rojas, J.L.G., Schmidt, C.G., Michel, V., Chueca, M.A.M., Roberts, H.C., Sihvonen, L.H., Stahl, K., Calvo, A.V., Viltrop, A., Winckler, C., Bett, B., Cetre-Sossah, C., Chevalier, V., Devos, C., Gubbins, S., Monaco, F., Sotiria-Eleni, A., Broglia, A., Abrahantes, J.C., Dhollander, S., Stede, Y.V.D., Zancanaro, G., 2020. Rift Valley fever — epidemiological update and risk of introduction into Europe. EFSA J. 18, e06041. https://doi.org/10.2903/j.efsa.2020.6041.

Outammassine, A., Boussaa, S., Zouhair, S., Loqman, S., 2021. Modeling and mapping the habitat suitability and the potential distribution of Arboviruses vectors in Morocco. Parasite 28, 37. https://doi.org/10.1051/parasite/2021030.

Outammassine, A., Zouhair, S., Loqman, S., 2022. Global potential distribution of three underappreciated arboviruses vectors (*Aedes japonicus*, *Aedes vexans* and *Aedes vittatus*) under current and future climate conditions. Transbound. Emerg. Dis. 69. https://doi.org/10.1111/tbed.14404.

Peyre, M., Chevalier, V., Abdo-Salem, S., Velthuis, A., Antoine-Moussiaux, N., Thiry, E., Roger, F., 2015. A systematic scoping study of the socio-economic impact of Rift Valley fever: research gaps and needs. Zoonoses Public Health 62, 309–325. https://doi.org/10.1111/zph.12153.

R Core Team, 2020. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing.

Rebaudo, F., Rabhi, V.-B., 2018. Modeling temperature-dependent development rate and phenology in insects: review of major developments, challenges, and future directions. Entomol. Exp. Appl. 166, 607–617. https://doi.org/10.1111/eea.12693.

Reiner, R.C., Perkins, T.A., Barker, C.M., Niu, T., Chaves, L.F., Ellis, A.M., George, D.B., Le Menach, A., Pulliam, J.R.C., Bisanzio, D., Buckee, C., Chiyaka, C., Cummings, D.A. T., Garcia, A.J., Gatton, M.L., Gething, P.W., Hartley, D.M., Johnston, G., Klein, E.Y., Michael, E., Lindsay, S.W., Lloyd, A.L., Pigott, D.M., Reisen, W.K., Ruktanonchai, N., Singh, B.K., Tatem, A.J., Kitron, U., Hay, S.I., Scott, T.W., Smith, D.L., 2013. A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970–2010. J. R. Soc. Interface 10, 20120921. https://doi.org/10.1098/rsif.2012.0921.

Sánchez-Vizcaíno, F., Martínez-López, B., Sánchez-Vizcaíno, J.M., 2013. Identification of suitable areas for the occurrence of Rift Valley fever outbreaks in Spain using a multiple criteria decision framework. Vet. Microbiol. 165, 71–78. https://doi.org/10.1016/j.vetmic.2013.03.016.

Schaeffer, B., Mondet, B., Touzeau, S., 2008. Using a climate-dependent model to predict mosquito abundance: application to *Aedes (Stegomyia) africanus* and *Aedes (Diceromyia) furcifer* (Diptera: Culicidae). Infect. Genet. Evol. J. Mol. Epidemiol. Evol. Genet. Infect. Dis. 8, 422–432. https://doi.org/10.1016/j.meegid.2007.07.002.

Schaffner, F., Versteirt, V., Bortel, W.V., Zeller, H., Wint, W., Alexander, N., 2016. VBORNET gap analysis: mosquito vector distribution models utilised to identify areas of potential species distribution in areas lacking records. Open Health Data 4 (e6), 8. https://doi.org/10.5334/ohd.27, 1-e6.

Soti, V., Tran, A., Degenne, P., Chevalier, V., Seen, D.L., Thiongane, Y., Diallo, M., Guéán, J.-F., Fontenille, D., 2012. Combining hydrology and mosquito population models to identify the drivers of Rift Valley fever emergence in semi-arid regions of West Africa. PLoS Negl. Trop. Dis. 6, e1795. https://doi.org/10.1371/journal.pntd.0001795.

Tran, A., Fall, A.G., Biteye, B., Ciss, M., Gimonneau, G., Castets, M., Seck, M.T., Chevalier, V., 2019. Spatial modeling of mosquito vectors for Rift Valley fever virus in Northern Senegal: integrating satellite-derived meteorological estimates in population dynamics models. Remote Sens. 11, 1024. https://doi.org/10.3390/rs11091024.

Tran, A., Ippoliti, C., Balenghien, T., Conte, A., Gely, M., Calistri, P., Goffredo, M., Baldet, T., Chevalier, V., 2013a. A geographical information system-based multicriteria evaluation to map areas at risk for Rift Valley fever vector-borne transmission in Italy. Transbound. Emerg. Dis. 60 (Suppl 2), 14–23. https://doi.org/10.1111/tbed.12156.

Tran, A., L'Ambert, G., Lacour, G., Benoît, R., Demarchi, M., Cros, M., Cailly, P., Aubry-Kientz, M., Balenghien, T., Ezanno, P., 2013b. A rainfall- and temperature-driven abundance model for *Aedes albopictus* populations. Int. J. Environ. Res. Public Health 10, 1698–1719. https://doi.org/10.3390/ijeph10051698.

Tran, A., Mangeas, M., Demarchi, M., Roux, E., Degenne, P., Haramboure, M., Le Goff, G., Damiens, D., Gouagna, L.-C., Herbreteau, V., Dehecq, J.-S., 2020. Complementarity of empirical and process-based approaches to modelling mosquito population dynamics with *Aedes albopictus* as an example—application to the development of an operational mapping tool of vector populations. PLoS ONE 15, e0227407. https://doi.org/10.1371/journal.pone.0227407.

Trari, B., 2017. Les moustiques (Insectes, Diptères) du Maroc : atlas de répartition et études épidémiologiques (Thèse de doctorat d'Etat és Sciences Biologiques). Faculté des Sciences, Université Mohammed V, Rabat, Maroc.

Tutsoy, O., 2023. Graph theory based large-scale machine learning with multi-dimensional constrained optimization approaches for exact epidemiological modeling of pandemic diseases. IEEE Trans. Pattern Anal. Mach. Intell. 45, 9836–9845. https://doi.org/10.1109/TPAMI.2023.3256421.

Veronesi, R., Gentile, G., Carrieri, M., Maccagnani, B., Stermieri, L., Bellini, R., 2012. Seasonal pattern of daily activity of *Aedes caspius*, *Aedes detritus*, *Culex modestus*, and *Culex pipiens* in the Po Delta of northern Italy and significance for vector-borne disease risk assessment. J. Vector Ecol. 37, 49–61. https://doi.org/10.1111/j.1948-7134.2012.00199.x.

Versteirt, V., Ducheyne, E., Schaffner, F., Hendrickx, G., 2013. Systematic literature review on the geographic distribution of Rift Valley fever vectors in Europe and the neighbouring countries of the Mediterranean Basin. EFSA Support. Publ. 10 (412E), 59. https://doi.org/10.2903/sp.efsa.2013.EN-412, 1-412E.

Wint, W., Van Bortel, W., Schaffner, F., 2020. RVF vector spatial distribution models: probability of presence. EFSA Support. Publ. 17 (1800E), 30. https://doi.org/10.2903/sp.efsa.2020.EN-1800, 1-1800E.

WOAH, 2021. Libya — Rift Valley fever virus (Inf. with) — Follow up report 9 [WWW Document]. WOAH — WAHIS. URL https://wahis.woah.org/#/in-review/3099?reportId=151770&fromPage=event-dashboard-url (accessed 5.2.23).

WHO, 2022. Rift Valley fever — Mauritania [WWW Document]. Dis. Outbreak News. URL https://www.who.int/emergencies/disease-outbreak-news/item/2022-DON417 (accessed 2.20.23).