An advection-deposition-survival model to assess the risk of introduction of vector-borne diseases through the wind: Application to bluetongue outbreaks in Spain

Abstract

This work develops a methodology for estimating risk of wind-borne introduction of flying insects into a country, identifying areas and periods of high risk of vector-borne diseases incursion. This risk can be characterized by the role of suitable temperatures and wind currents in small insects’ survival and movements, respectively. The model predicts the number density of introduced insects over space and time based on three processes: the advection due to wind currents, the deposition on the ground and the survival due to climatic conditions. Spanish livestock has suffered many bluetongue outbreaks since 2004 and numerous experts point to Culicoides transported by wind from affected areas in North Africa as a possible cause. This work implements numerical experiments simulating the introduction of Culicoides in 2004. The model identified southern and eastern Spain, particularly between June and November, as being at greatest risk of wind-borne Culicoides introduction, which matches field data on bluetongue outbreaks in Spain this year. This validation suggests that this model may be useful for predicting introduction of airborne pathogens of significance to animal productivity.

About the Authors:

Eduardo Fernández-Carrión

Roles Conceptualization, Formal analysis, Methodology, Resources, Software, Writing – original draft

* E-mail: [email protected]

Affiliation: VISAVET Center and Animal Health Department, Veterinary School, Universidad Complutense de Madrid, Madrid, Spain

ORCID http://orcid.org/0000-0003-4384-1961

Benjamin Ivorra

Roles Conceptualization, Formal analysis, Methodology, Software, Supervision, Validation, Writing – review & editing

Affiliation: MOMAT Research group, IMI-Institute and Applied Mathematics Department, Universidad Complutense de Madrid, Madrid, Spain

Ángel Manuel Ramos

Roles Conceptualization, Formal analysis, Methodology, Project administration, Supervision, Validation, Writing – review & editing

Affiliation: MOMAT Research group, IMI-Institute and Applied Mathematics Department, Universidad Complutense de Madrid, Madrid, Spain

Beatriz Martínez-López

Roles Conceptualization, Data curation, Resources, Writing – original draft, Writing – review & editing

Affiliation: CADMS Center for Animal Disease Modeling and Surveillance, School of Veterinary Medicine, UC Davis, Davis, California, United States of America

Cecilia Aguilar-Vega

Roles Writing – review & editing

Affiliation: VISAVET Center and Animal Health Department, Veterinary School, Universidad Complutense de Madrid, Madrid, Spain

ORCID http://orcid.org/0000-0001-5981-7291

José Manuel Sánchez-Vizcaíno

Roles Conceptualization, Data curation, Project administration, Resources, Writing – review & editing

Affiliation: VISAVET Center and Animal Health Department, Veterinary School, Universidad Complutense de Madrid, Madrid, Spain

Introduction

Introduction of vector-borne diseases can substantially harm public and animal health, causing significant sanitary and financial loss. Long-range wind-borne transportation of infected flying insects has previously been linked to the introduction of viruses affecting humans and/or animals such as West Nile virus, dengue or Rift Valley fever. Therefore, modeling the movement of wind-borne midges or mosquitoes may help identify geographic regions and seasons at higher risk of incursion of several arboviral diseases, improving surveillance efforts.

Bluetongue is a viral disease of ruminants transmitted by biting midges of the genus Culicoides. Culicoides imicola is considered the main vector of Bluetongue in South Europe and Africa [1]. This virus has traditionally been considered endemic to tropical and subtropical regions, between parallels 35°S and 40°N; however, recently several outbreaks have occurred further North and it is considered to be a re-emerging disease in new latitudes [2]. So far, the disease has been reported in Australia, USA, Africa, Middle East, Asia and Europe, although the geographic distribution of the 27 virus serotypes differs by region [3].

Previous studies suggest that Culicoides could be transported by the wind to distances of up to 170 km from their original location [4]. Wind-borne movement of bluetongue virus-infected Culicoides is thought to explain the northward spread of bluetongue in Europe, primarily through routes from Turkey to Greece-Bulgaria, from Algeria-Tunisia to Italy and from Morocco-Algeria to Spain [5]. Specifically, wind-borne introduction of these insects to Spain from North Africa has been suggested as the cause of the numerous bluetongue outbreaks in Spain in 2004–05 [6–8] and in other countries and periods [9, 10]. In October 2004, 121 bluetongue outbreaks were reported in Spain; this number increased to 322 in December 2004 and to 328 at the start of 2005 [11–13].

The goal of the present work was to develop a model that could identify geographic areas and time periods at highest risk of introduction of small insects potentially infected with human and/or animal pathogens. The model takes into account three processes: (1) advection, in which insects are considered as particles in the air, and move according to wind currents; (2) deposition, in which wind-borne insects are deposited on the ground; and (3) survival, reflecting insects’ survival or mortality as a function of climatic factors. This model was numerically implemented to simulate wind-borne Culicoides introduction from North Africa into Spain and was validated using field data on bluetongue outbreaks reported in Spain in 2004.

Materials and methods

Model’s scheme

Culicoides are 1–3 mm in size, they weigh less than 1 mg and, under favorable conditions, they can be found as high as 200 m above the ground [14]. Previous studies consider that the long-range transportation of Culicoides is likely caused by wind currents at high altitudes [1, 5, 15]. Therefore the model assumes that the movement of small insects in the air shows a similar behaviour to wind-borne small particles [16]. Then, the model takes into consideration the following processes:

1. (1). Wind advection. Particle movement in the air due to wind currents is computed using an advection partial differential equation (PDE). The same kind of PDE has been used to predict the movements of dust [17, 18], sediments [19, 20], oil [21, 22] or gaseous substances [23, 24].

2. (2). Vertical deposition. The model assumes that Culicoides deposited onto the ground from the air in a manner similar to other airborne particles. Thus, deposition is computed using particles’ sedimentation and flotation theory based on settling principles through the air, adapted here to the specific characteristics of Culicoides.

3. (3). Survival. Dead insects do not pose an introduction risk. Consequently, the model estimates a mortality rate during wind-borne transport on the basis of climatic variables.

Mathematical model

We consider the spatial domain and time interval [0, T]. We denote by Γ_g and Γ_s(t) the subsets of the boundary of Ω corresponding, respectively, to the ground and to the source of Culicoides entering into Ω at time t ∈ [0, T]. C(x, y, z, t) is the spatial and temporal distribution of the number density of insects (i.e. the number of insects per unit volume) at point (x, y, z) ∈ Ω at time t ∈ [0, T]. Since the number density of insects varies according to wind advection, vertical deposition and survival, the evolution of C is governed by the following equations(1)where C₀ is the initial distribution of the number density of Culicoides in Ω, C_s(x, y, z, t) is the number density of insects at point (x, y, z) ∈ Γ_s(t) at time t ∈ [0, T], w = (w_x, w_y, w_z) is the velocity field related to the wind and deposition effects, which satisfies w_z = 0 on Γ_g × (0, T], and σ(f, C) is the mortality function (i.e. the number of Culicoides that die per unit volume and unit time), which depends on C and the spatio-temporal distribution of climatic factors f associated with survival.

Model implementation for the case of Culicoides

Here, the model formulated above is tailored for studying the case of the introduction of Culicoides from North Africa to Iberian Peninsula and Balearic Islands. The model has been fully developed in Matlab.

Domain’s discretization.

We discretized Γ_g into a two-dimensional regular meshgrid of n_x × n_y grid points, where . Here we considered n_x = n_y = 50, covering the area extending from latitude 34°N to 44°N and from longitude 10°W to 4.25°E (around 22 and 26 km between vertical and horizontal grid points, respectively), comprising parts of Spain, Portugal, northern Morocco and northern Algeria (Fig 1).

[Figure omitted. See PDF.]

Fig 1. Domain Γ_g in numerical experiments, extending from latitude 34°N to 44°N and from longitude 10°W to 4.25°E.

Black lines demarcate country boundaries and black dots locate the fixed 2,500 grid points.

https://doi.org/10.1371/journal.pone.0194573.g001

We also discretized the temporal domain [0, T] into steps, using a fixed time step Δt > 0, given by Δt = T/N. Here, we considered Δt = 1/24 days (1 hour).

Climatic data source and interpolation.

Wind and temperature data on Γ_g were sourced from Forecast.io website (www.forecast.io). This data consisted of hourly measurements of mean temperature (θ, in Celsius degrees °C), mean horizontal wind direction (from 0° to 359°) and mean horizontal wind speed (in m/s) taken along 2004 at each grid point on the surface (2,500 source points). Climatic data at intermediate points on the ground within the domain, Γ_g, were estimated through linear interpolation at every time step.

At high altitudes, z > 10 m, temperature values were estimated through the standard free atmospheric lapse rate of 6.5°C decrease per kilometer of increasing elevation estimated for Mediterranean atmosphere [25]:(2)where θ(z₀) is the temperature on the ground z₀, and γ = 0.0065 °C/m.

Horizontal wind speed values at high altitudes, z > 10 m, were estimated using the following polynomial variation [26]:(3)where w(z₀) is the wind speed on the ground z₀, and α is the Hellman exponent for particles above sea and land surfaces (in this work, α = 0.10 and α = 0.16, respectively).

Deposition.

Vertical deposition of Culicoides from high altitudes was assumed to be governed by processes similar to those affecting the deposition velocity of dust particles. This deposition is due to the combined action of gravitational settling, wind drag capacity and the aerodynamic properties of the particles (here, the Culicoides).

In general, flying insects use integrated systems consisting of wings to generate aerodynamic forces, muscles to move the wings, and sensing and control systems to guide and manoeuvre [27]. The Reynolds number, Re, and the drag coefficient, C_d, are two dimensionless quantities used in fluid mechanics to predict the relative internal movement and drag of particles in different fluid velocities [28]. In particular, as the size of a flier is reduced, the wing-to-body mass ratio tends to decrease paired with the Reynolds number as well [27] (in this work, Re = 100, assumed for small flying insects with slow flights [27, 29]). Following Haider et. al. [30], the drag coefficient was computed as a function of the Reynolds number adjusted for small spherical particles:(4)

The settling velocity for Culicoides was assumed to be governed by the same processes affecting the settling velocity of spherical particles, V_g(z) [31]:(5)where R is the radius of the Culicoides, here assumed to be R = 2 ⋅ 10⁻³ m [32]; M is the weight of the Culicoides, here assumed to be M = 0.5 ⋅ 10⁻³ kg [33]; ρ is the density of the Culicoides, here approximated to an spherical particle ρ = M/(4/3πR³); g is the gravity acceleration, assumed constant, g = 9.81 m/s²; and ρ_air(θ) is the mass air density at temperature θ, see Table 1 [34].

[Figure omitted. See PDF.]

Table 1. Effect of temperature on air density.

https://doi.org/10.1371/journal.pone.0194573.t001

However, when wind becomes sufficiently strong, particles in the air are susceptible to enter into long-term suspension and are therefore transported long-range [35]. Since one of the main variables for modeling sediment transport is the force exerted on the soil by the wind, a fundamental variable used to predict aeolian particle transport is the shear velocity, , whose module is calculated as [36, 37]; where is the surface shear stress, i.e. the vertical flux of horizontal momentum measured near the surface; which can be computed as [38, 39]. Thus,(6)

We may now use the dimensionless parameter to distinguish three types of flight for Culicoides [36, 37]: i) a long term suspension, when λ < 0.1; ii) a short term suspension, when 0.1 < λ < 0.7; iii) and vertical settling, when λ > 0.7. Concretely, the settling velocity of Culicoides, w_z, was computed in terms of V_g (see Eq 5) through the following continuous function:(7)

Since the number density of wind-borne particles is generally greater at lower altitudes than at higher ones [40, 41], we assumed the center of mass of wind-borne living midges to lie at an altitude of 1 km. On the other hand, according to the literature, while long-range transportation of Culicoides is caused by wind currents at high altitudes, favorable temperatures and wind speeds for natural Culicoides outdoor mobility are often located less than 200 m above the ground [1, 5, 14, 15]. Thus, the model assumes that below this altitude, the Culicoides do not move through wind currents (i.e. w = (0, 0, 0) if z < 200 m).

Mortality function.

The mortality function predicts the number density of insects based on temperature, since this is one of the most important determinants of Culicoides life cycle and lifespan [42, 43]. Although other climatic parameters, such as relative humidity and precipitation, can also affect insect survival, neglecting these factors only moderately affects the accuracy of model predictions [44]. In order to simplify the model, we only considered the impact of temperature in the Culicoides mortality function σ; i.e. f = θ, with θ(x, y, z, t) being the temperature value (here in °C) at point (x, y, z) ∈ Ω and at time t ∈ [0, T]. Therefore the mortality function is defined by(8)where μ is the Culicoides mortality rate (i.e. the inverse of lifespan) as a function of temperature. This function was derived based on studies reporting total Culicoides mortality after 2 h at 42°C or −6°C [45] and Culicoides lifespan was assumed to be shorter than 24 h after exposure to 40.5°C or −4°C [45]. Culicoides lifespan of 27.5, 18.8, 13.4 and 10.2 days were observed at temperature of 10°C, 15°C, 20°C and 30°C, respectively [46]. These values were used to estimate intermediate values for arbitrary temperatures through the Hermite cubic interpolation computed with the Matlab function pchip and denoted by 1/μ(.). Thus the inverse function μ(.) returns the estimated mortality rate of the Culicoides depending on the value of the temperature (Table 2).

[Figure omitted. See PDF.]

Table 2. Estimated mean lifespan of Culicoides imicola at different temperatures, based on field data [45, 46].

https://doi.org/10.1371/journal.pone.0194573.t002

Culicoides sampling data and activity function.

Data on the number of Culicoides imicola trapped per day along the Spanish territory between 2004 and 2013 and average temperatures associated with each capture were collected from the Ministry of Agriculture, Fisheries, Food and Environment (www.mapama.gob.es) within the Spanish national surveillance programme.

Field data on the mean number of Culicoides trapped in Spain at temperatures between 0°C and 35°C were fitted to a polynomial function computed with the polyfit command in Matlab, which was defined as the activity function and denoted by δ(.). We used this function to estimate the outdoor number density of Culicoides in source locations, Γ_s(t).

Numerical implementation.

The system of Eq (1) was solved numerically using a Lagrangian particle approach. The source of the Lagrangian particles remained constant (i.e. Γ_s(t) = Γ_s, ∀t ∈ [0, T]) and was defined as lying in fixed cells at an altitude of 1 km over northern Morocco and/or northern Algeria. At each time step, M new particles started from Γ_s, such that the model included (k + 1) × M particles at time step k ∈ {0, …, N − 1}, and the total amount of particles per simulation was N × M. All particles were assigned with an index i ∈ {1, …, N × M}. Therefore, at time step k, M particles with indexes {k × M + 1, ⋯, (k + 1) × M} started from Γ_s.

Particle trajectories were estimated by considering the velocity field w through the explicit scheme:(9)where is the position of the Lagrangian particle with index i = {1, …, (n + 1) × M} at time step n = n(i), ‥, N, is the wind velocity vector at position at time step n and (i.e. n(i) is the greatest integer less than or equal to ). Then, the number density of Culicoides transported at step n by particle i (which was located at ) was denoted as .

Given the difficulty of estimating the number of Culicoides placed in a specific location, several epidemiological works assume an initial ratio of mosquitoes or midges instead [47, 48]. In our case, the number density of Culicoides transported by particle i situated in Γ_s at step n was assumed to be proportional to the Culicoides activity function, δ(.), as follows(10)

Next, the evolution of this particle number density function was governed by the following scheme:(11)where i = {1, …, (n + 1) × M} and n ∈ {n(i), ‥, N − 1}.

According to the features of the Culicoides, the number density of Culicoides deposited on Γ_g was calculated at each time step n as the sum of all particles below 200 m altitude; consequently, Eq (11) was not applied to such particles.

Simulations.

Three numerical experiments were performed to evaluate the risk of Culicoides introduction in Spain from North Africa and thereby allow validation of our model against field data of bluetongue outbreaks reported in 2004. The first experiment compared the introduction of Culicoides along an entire year from the cells located at 1 km above the ground of Morocco alone, Algeria alone, and both countries (i.e. the choice of three different sources Γ_s). Each simulation started on the first day of 2004 and finished after T = 365 days. The goal of this experiment was to identify areas at high risk of midge introduction during one year from different African countries.

The second experiment examined the risk of wind-borne Culicoides introduction from both Morocco and Algeria for each month of a year; thus, 12 simulations were performed, each starting on the first day of every month. The goal of this experiment was to identify months associated with higher risk of Culicoides introduction.

A third experiment aimed to study the influence that variation in the values of input parameters affect in the model outcomes, especially the climatic factors and the physical features of the Culicoides on the final depositions in the target territory. Thus, a sensitivity analysis was carried out by increasing ±1 to ±5% and ±5 to ±25% the value of the following climatic and physical variables, respectively, in an entire year of simulation from initial Culicoides positions located at 1 km above the ground of Morocco and Algeria: temperature, θ, wind speed, w, Hellman exponent, α, temperature lapse rate, γ, and insect’s weight, M, radius, R, Reynolds number, Re. The results of the model is also associated to the initial altitude at source locations, Γ_s, here located at 1 km high. Thus, a sensitivity analysis was carried out by increasing ±5 to ±25% the initial altitude.

The results of the numerical experiments were compared with field data from the bluetongue outbreaks in Spain in 2004. This allowed us to validate the model.

Results

Mortality and activity functions

Fig 2 shows the lifespan associated to different values of temperature (blue dots) displayed in Table 2 and the cubic Hermite interpolation computed for intermediate values (red line). Thus, the mean lifespan of Culicoides is(12)where the value of a_k(θ), for each k ∈ {0, …, 4}, and b(θ) are given in Table 3, depending on the value of θ.

[Figure omitted. See PDF.]

Fig 2. Culicoides lifespan used in the present model according to temperature variations based on field data (blue dots) and cubic Hermite approximation (red line), 1/μ(.).

https://doi.org/10.1371/journal.pone.0194573.g002

[Figure omitted. See PDF.]

Table 3. Piecewise cubic Hermite interpolating polynomial parameters for the mean lifespan of Culicoides according to temperature values from −6°C to 42°C, see Eq (12).

https://doi.org/10.1371/journal.pone.0194573.t003

Fig 3 shows the mean number of Culicoides captured (blue dots) and the approximation polynomial fitting (red line). Thus, considering a midges abundance higher or equal to zero, the Culicoides activity function is defined as(13)where c₀ = −0.226, c₁ = 0.344 °C⁻¹, c₂ = −0.178 °C⁻², c₃ = 4.24 ⋅ 10⁻² °C⁻³, c₄ = −4.85 ⋅ 10⁻³ °C⁻⁴, c₅ = 2.89 ⋅ 10⁻⁴ °C⁻⁵, c₆ = −9.04 ⋅ 10⁻⁶ °C⁻⁶, c₇ = 1.38 ⋅ 10⁻⁷ °C⁻⁷ and c₈ = −7.90 ⋅ 10⁻¹⁰ °C⁻⁸, if θ ∈ [1.25, 33.45] °C; and c_k = 0 °C^−k, for each k ∈ {1, …, 8}, otherwise.

[Figure omitted. See PDF.]

Fig 3. Averaged Culicoides trapped per capture (blue dots) and polynomial approximation curve (red line), δ(.).

https://doi.org/10.1371/journal.pone.0194573.g003

Annual risk of introduction

Fig 4 shows the number density of Culicoides deposited (in logarithmic scale) across the entire year on target territory via wind-borne transport from the skies over Morocco, Algeria or both. The most exposed areas were the south and east of the Iberian Peninsula and Balearic Islands. Midge deposition patterns varied considerably depending on whether the source was Morocco or Algeria.

[Figure omitted. See PDF.]

Fig 4. Simulated spatial distribution of the number density of Culicoides deposited on the target territory throughout 2004, when the midge source was both Morocco and Algeria (left), Morocco alone (center) or Algeria alone (right).

Number densities were plotted in logarithmic scale, log(C + 1) (log1). Blue dots represent the source of Culicoides, Γ_s.

https://doi.org/10.1371/journal.pone.0194573.g004

When Morocco was the source, midge depositions constituted the 50.99% of the total depositions, and they occurred mainly in the south cost of the Iberian Peninsula, and especially close to the Strait of Gibraltar. However, the Culicoides spread might reach the inland west areas of the Iberian Peninsula, reaching even the north of Portugal. On the other hand, a low number density of Culicoides reached the Balearic Islands, discarding the introduction of Culicoides in the center of the Iberian Peninsula.

When Algeria was the source, midge depositions constituted the 49.01% of the total depositions, and they occurred mainly in the region from southeastern Spain to the eastern coast and Balearic Islands. In some cases, the Culicoides spread may reach the center and the northeast of the Iberian Peninsula, but in very low quantity.

Monthly risk of introduction

Since changes in wind currents and temperature during the year can affect the dynamics of midge transport from North Africa, we simulated the spatial density of deposited Culicoides over all 12 months of the year. Table 4 shows the percentage of the number densities of Culicoides deposited on the target territory during each month of 2004. Despite having strong wind currents between December and February, the low temperatures do not favor the Culicoides activity and it represented a low risk for introduction, reaching 12.65% of the total depositions in a year. Despite favorable temperatures between March to May, the wind currents were not sufficiently strong to transport high amount of Culicoides, only reaching 19.40%. The period of high risk of introduction is concentrated in the remaining half of the year, reaching 29.43% between June and August; and 38.52% between September and November, when strong wind streams and suitable temperatures for Culicoides activity in African territories are optimal.

[Figure omitted. See PDF.]

Table 4. Simulated percentage of number density of Culicoides deposited on the ground during each month of 2004.

https://doi.org/10.1371/journal.pone.0194573.t004

Next we examined the spatial distribution of depositions over the 12 months of the year. Fig 5 shows the monthly distribution of Culicoides number density deposited on target territory from the source in both Morocco and Algeria. The largest densities of deposition occurred in the southern and eastern parts of the Iberian Peninsula and Balearic Islands, consistent with the annual-level simulation. Risk of wind-borne incursions were high throughout the year in the southern part of the Iberian Peninsula, specifically close to the Strait of Gibraltar. Central areas of the Iberian Peninsula and the north of Portugal were scarcely reached in specific periods with strong wind streams. However, the Balearic Islands were reached frequently with a high amount of depositions, as mentioned above, mainly from Algeria but also could be reached from Morocco.

[Figure omitted. See PDF.]

Fig 5. Simulated spatial distribution of the number density of Culicoides deposited on Spanish territory for each month of 2004, when the midge source was both Morocco and Algeria.

Number densities were plotted in logarithmic scale, log(C + 1) (log1). Blue dots represent the source of Culicoides, Γ_s.

https://doi.org/10.1371/journal.pone.0194573.g005

Sensitivity analysis

Table 5 shows the changes in the output number density deposited over target territory, under the assumption of slight variations in climatic variables and in physical features of the Culicoides. Among all climatic variables computed, the wind speed and the Hellmann exponent, both associated directly with wind speed at higher altitudes, represented the most relevant variables in terms of particle transportation. Indeed, an increase/decrease of 5% in the wind speed (e.g. from 20 m/s to 21 m/s) may increase/decrease the final number density around ±10%. In contrast, variables affecting temperature slightly varied the final number density. Indeed, a temperature increase/decrease of 5% (e.g. from 20°C to 21°C) may increase/decrease the final number density around ±3%.

[Figure omitted. See PDF.]

Table 5. Percentage of output variation of final number density of Culicoides deposited over target territory, C, after a whole year simulation for input variation of {±5%, ±4%, ±3%, ±2%, ±1%} for climatic variables and {±25%, ±20%, ±15%, ±10%, ±5%} for insect’s variables.

https://doi.org/10.1371/journal.pone.0194573.t005

A reduction in the weight of the Culicoides, or an increment in the radius decreases the density of the wind-borne particle, ρ, which reduces the settling velocity accordingly to Eq 5. On the contrary, a reduction in Reynolds number increases the drag coefficient, C_d, which decreases the settling velocity, see Eqs 4 and 5. These cases implied a reduction in settling velocity, which extends long-suspension in the air reaching longer distances.

Model validation

Since October 2004 to January 2005, a total of 328 bluetongue outbreaks were reported in the Iberian Peninsula, with no previous outbreaks reported since 1960 [6, 11–13]. Among them, 82% were located in the south of Spain, and only 18% were located in the south-west. Since September 2004 to December 2004, a total of 230 outbreaks were reported in Morocco, with no previous outbreaks reported it that year. Concretely, 28 new outbreaks were reported in September in north-northeastern Morocco. Other works consider that a period of 4 weeks is the time necessary to detect an infected animal once the infected vector lands in the surroundings [15, 49]. According to this, the virus could have been introduced in September from infected areas of Morocco, to the Spanish territory, concretely to areas close to the Strait of Gibraltar (see Fig 6).

[Figure omitted. See PDF.]

Fig 6. Bluetongue outbreaks reported in Spain in between October 2004 and January 2005 (green dots).

First outbreaks reported between October 12th and 19th (white diamonds). Simulated spatial distribution of the number density of Culicoides deposited on the Iberian Peninsula throughout 2004 when the midge source was the sky over Morocco (see Fig 5), plotted in logarithmic scale, log(C + 1) (log1).

https://doi.org/10.1371/journal.pone.0194573.g006

Discussion

In this work, we developed a methodology for estimating risk of long-range wind-borne transport and deposition of small insects, while taking into account the likelihood of insect survival given certain climatic conditions. The methodology relies on a mathematical advection model that is solved numerically using Lagrangian particle tracking, based on the wind velocity field, and particles’ sedimentation and flotation theory for deposition, that takes into account temperature-based insect lifespan and activity. In particular, we focused on assessing the risk of introduction of Culicoides from North Africa to Spain, in order to explain the bluetongue outbreaks that occurred in this territory during 2004.

Two numerical experiments were carried out to simulate the incursion of Culicoides through wind currents in 2004 in Spain, one yearly and one monthly. The results showed that the highest number density of Culicoides deposited on regions agreed with the real primary bluetongue reported outbreaks in the period of study. Concretely, since September 2004, several bluetongue outbreaks were reported in Morocco. In this period, the model showed that main incursions of Culicoides from Moroccan territory were deposited close to the Strait of Gibraltar, showing a high rate of matching between model results and bluetongue notifications.

For the remaining reported outbreaks and for future reported ones in 2005 [13], Calvete et al. [50] identified suitable areas for Culicoides based on climatic factors and livestock distribution in Spain. Those regions were distributed mainly in the south and south-west of Spain. The model developed in this work aimed to assess midges’ incursions from North Africa due to wind currents, not the spread within the Iberian Peninsula. In future works, this model should link continuous updated bluetongue status in North African areas, with the host population distribution and with spatio-temporal suitable conditions for vector survival in target, as proposed by Calvete et al. [50], in order to estimate not only the risk of introduction, but also the risk of spread.

Our results suggest that the Culicoides source can substantially affect midge deposition in Spain: when midges came from over northern Morocco, they usually deposited close to the Strait of Gibraltar; when they came from over northern Algeria, they usually deposited along southeastern Spain and the Balearic Islands. Our results further showed risk variation during the year: risk of introduction was relatively high between June and November, and moderate during other months. These findings illustrate the power of our approach for helping to focus surveillance efforts in space and time.

This work was motivated because, first, different authors foretold the Strait of Gibraltar as one of the main bluetongue introduction pathways in Europe [2, 5] and, second, the presence of Culicoides and bluetongue outbreaks were reported along the region and period of study, which allowed the validation of our predictions [11, 12]. The use of the methodology presented here may serve to assess regions and periods at high risk of incursion of several vectors such as midges or mosquitoes. Indeed, Culicoides may also be a vector for African horse sickness, and many species of mosquitoes are vectors for diseases such as West Nile virus, dengue, Rift Valley fever, Zika, etc. Currently, the system developed in this work embeds a wind and climatic databases automatically, which allows surveillance of real-time (hourly) wind patterns and Culicoides deposition under friendly configurable features, which are easily adaptable to other insects. However, an important quantitative validation of this model should be carried out with more bluetongue outbreaks in territories with no previous disease notifications in a long period; e.g. the bluetongue outbreaks in other European territories or possible infections of African territories from Spanish source locations. The application of this system is to launch risk alerts when suspicious initial infected areas that may spread the vector to long-range distances, especially where livestock population may be exposed. This development could contribute to an effective early warning system for preventing and controlling vector-borne disease incursions.

Supporting information

[Figure omitted. See PDF.]

S1 File. Video simulation.

Video file of particle movement in the air and ground concentration on Spanish territory during September 2004.

https://doi.org/10.1371/journal.pone.0194573.s001

(MP4)

Acknowledgments

This study was funded by the EU project ‘Understanding pathogen, livestock, environment interactions involving bluetongue’ (PALE-BLU, H2020-SFS-2016-2); the Spanish Ministry of Economy and Competitiveness (project MTM2015-64865-P); the Junta de Andalucía and European Regional Development Fund (project P12-TIC301); and the MOMAT research group (Ref. 910480), supported by the Universidad Complutense de Madrid. This work was carried out using data from the Spanish Meteorological Service (AEMET) and the Ministry of Agriculture, Fisheries, Food and Environment (MAPAMA). The authors thank Javier Lucientes (Universidad de Zaragoza) for providing support and data on Culicoides captures and Pablo Ayala (Universidad Pontificia Comillas) for support in fluid mechanics theory.

Citation: Fernández-Carrión E, Ivorra B, Ramos ÁM, Martínez-López B, Aguilar-Vega C, Sánchez-Vizcaíno JM (2018) An advection-deposition-survival model to assess the risk of introduction of vector-borne diseases through the wind: Application to bluetongue outbreaks in Spain. PLoS ONE 13(3): e0194573. https://doi.org/10.1371/journal.pone.0194573

References

1. Mellor PS, Boorman J, Baylis M. Culicoides biting midges: their role as arbovirus vectors. Annual review of entomology. 2000;45(1):307–340. pmid:10761580

2. Wilson AJ, Mellor PS. Bluetongue in Europe: past, present and future. Philosophical Transactions of the Royal Society of London B: Biological Sciences. 2009;364(1530):2669–2681. pmid:19687037

3. Maan S, Maan NS, Nomikou K, Veronesi E, Bachanek-Bankowska K, Belaganahalli MN, et al. Complete genome characterisation of a novel 26th bluetongue virus serotype from Kuwait. PLoS One. 2011;6(10):e26147. pmid:22031822

4. Gloster J, Mellor PS, Burgin L, Sanders C, Carpenter S. Will bluetongue come on the wind to the United Kingdom in 2007? The Veterinary Record. 2007;160(13):422–426. pmid:17400899

5. Wilson A, Mellor P. Bluetongue in Europe: vectors, epidemiology and climate change. Parasitology Research. 2008;103(1):69–77.

6. Diego AC, Sánchez-Cordón PJ, Sánchez-Vizcaíno JM. Bluetongue in Spain: from the first outbreak to 2012. Transboundary and Emerging Diseases. 2014;61(6):e1–e11. pmid:23480789

7. Martínez-López B, Pérez C, Baldasano JM, Sánchez-Vizcaíno JM. BTV: Risk of introduction into Spain associated with wind streams. In: 12th Symposium of the International Society for Veterinary Epidemiology and Economics; 2009. p. 94–96.

8. Allepuz A, García-Bocanegra I, Napp S, Casal J, Arenas A, Sáez M, et al. Monitoring bluetongue disease (BTV-1) epidemic in southern Spain during 2007. Preventive Veterinary Medicine. 2010;96(3):263–271. pmid:20663576

9. Gubbins S, Szmaragd C, Burgin L, Wilson A, Volkova V, Gloster J, Gunn GJ. Assessing the consequences of an incursion of a vector-borne disease: I. identifying feasible incursion scenarios for bluetongue in Scotland Epidemics. 2010;2(3):148–154. pmid:21352785

10. Eagles D, Melville L, Weir R, Davis S, Bellis G, Zalucki MP, Walker PJ, Durr PA. Long-distance aerial dispersal modelling of Culicoides biting midges: case studies of incursions into Australia. BMC veterinary research. 2014;10(1):135. pmid:24943652

11. Informe de la situación de la lengua azul en España. Dirección general de la ganadería. Subdirección general de sanidad animal; 2004 (29 de octubre).

12. Informe de la situación de la lengua azul en España. Dirección general de la ganadería. Subdirección general de sanidad animal; 2004 (7 de diciembre).

13. Informe de la situación de la lengua azul en España. Dirección general de la ganadería. Subdirección general de sanidad animal; 2005 (14 de enero).

14. Sanders CJ, Selby R, Carpenter S, Reynolds DR. High-altitude flight of Culicoides biting midges. Bulletin of Entomological Research. 2011;94:123–136.

15. Hendrickx G, Gilbert M, Staubach C, Elbers A, Mintiens K, Gerbier G, et al. A wind density model to quantify the airborne spread of Culicoides species during north-western Europe bluetongue epidemic, 2006. Preventive Veterinary Medicine. 2008;87(1):162–181. pmid:18639355

16. Lillie TH, Marquardt WC, Jones RH. The flight range of Culicoides variipennis (Diptera: Ceratopogonidae). The Canadian Entomologist. 1981;113(05):419–426.

17. In HJ, Park SU. A simulation of long-range transport of Yellow Sand observed in April 1998 in Korea. Atmospheric Environment. 2002;36(26):4173–4187.

18. In HJ, Park SU. The soil particle size dependent emission parameterization for an Asian dust (Yellow Sand) observed in Korea in April 2002. Atmospheric Environment. 2003;37(33):4625–4636.

19. Harris CK, Wiberg PL. A two-dimensional, time-dependent model of suspended sediment transport and bed reworking for continental shelves. Computers & Geosciences. 2001;27(6):675–690.

20. Warner JC, Sherwood CR, Signell RP, Harris CK, Arango HG. Development of a three-dimensional, regional, coupled wave, current, and sediment-transport model. Computers & Geosciences. 2008;34(10):1284–1306.

21. Alavani C, Glowinski R, Gómez S, Ivorra B, Joshi P, Ramos AM. Modelling and simulation of a polluted water pumping process. Mathematical and Computer Modelling. 2010;51(5):461–472.

22. Ivorra B, Gómez S, Glowinski R, Ramos AM. Nonlinear Advection–Diffusion–Reaction Phenomena Involved in the Evolution and Pumping of Oil in Open Sea: Modeling, Numerical Simulation and Validation Considering the Prestige and Oleg Naydenov Oil Spill Cases. Journal of Scientific Computing. 2015;70(3):1–27.

23. Wortmann S, Vilhena MT, Moreira DM, Buske D. A new analytical approach to simulate the pollutant dispersion in the PBL. Atmospheric Environment. 2005;39(12):2171–2178.

24. Bruse M, Fleer H. Simulating surface–plant–air interactions inside urban environments with a three dimensional numerical model. Environmental Modelling & Software. 1998;13(3):373–384.

25. Kuhlemann J, Rohling EJ, Krumrei I, Kubik P, Ivy-Ochs S, Kucera M. Regional synthesis of Mediterranean atmospheric circulation during the last glacial maximum. Science. 2008;321(5894):1338–1340. pmid:18669823

26. Heier S. Grid integration of wind energy: onshore and offshore conversion systems. John Wiley & Sons. 2014.

27. Shyy W, Kang C, Chirarattananon P, Ravi S, Liu H. Aerodynamics, sensing and control of insect-scale flapping-wing flight. Proc. R. Soc. A. 2016;472:20150712. pmid:27118897

28. Gabitto J, Tsouris C. Drag coefficient and settling velocity for particles of cylindrical shape. Powder Technology. 2008;183(2):314–322.

29. Shyy W, Liu H. Flapping wings and aerodynamic lift: the role of leading-edge vortices. AIAA Journal. 2007;45(12):December.

30. Haider A, Levenspiel O. Drag coefficient and terminal velocity of spherical and nonspherical particles. Powder technology. 1989;58(1):63–70.

31. Gregory R, Zabel TF, Edzwald JK. Sedimentation and flotation. Water Quality and Treatment. 1999;7:1–87.

32. Blanton FS, Wirth WW. The sand flies (Culicoides) of Florida (Diptera: Ceratopogonidae). Arthropods of Florida and Neighboring Land Areas. 1979;10.

33. Leprince DJ, Higgins JA, Church GE, Issel CJ, McManus JM, Foil LD. Body size of Culicoides variipennis (Diptera: Ceratopogonidae) in relation to bloodmeal size estimates and the ingestion of Onchocerca cervicalis (Nematoda: Filarioidea) microfiliariae. Journal of the American Mosquito Control Association. 1989;.

34. Davis RS. Equation for the determination of the density of moist air (1981/91). Metrologia. 1992;29(1):67.

35. Kok JF, Parteli EJR, Michaels TI, Karam DB. The physics of wind-blown sand and dust. Reports on Progress in Physics. 2012;75(10). pmid:22982806

36. Tsoar H, Pye K. Dust transport and the question of desert loess formation. Sedimentology. 1987;34(1):139–153.

37. Brenig L, Offer Z. Airbone particles dynamics: towards a theoretical approach Environmental Modeling and Assessment. 2001;6(1):1–5.

38. Powell MD, Ginis I. Drag coefficient distribution and wind speed dependence in tropical cyclones. Final Report to the National Oceanic and Atmospheric Administration (NOAA) Joint Hurricane Testbed (JHT) Program. 2006.

39. Kraus EB, Businger JA. Atmosphere-ocean interaction. Oxford University Press. 1994;27.

40. Solomos S, Kallos G, Mavromatidis E, Kushta J. Density currents as a desert dust mobilization mechanism. Atmospheric Chemistry and Physics. 2012;12(22).

41. Tsidulko M, Krichak SO, Alpert P, Kakaliagou O, Kallos G, Papadopoulos A. Numerical study of a very intensive eastern Mediterranean dust storm, 13–16 March 1998. Journal of Geophysical Research: Atmospheres. 2002;107(D21).

42. Gloster J, Burgin L, Witham C, Athanassiadou M, Mellor PS. Bluetongue in the United Kingdom and northern Europe in 2007 and key issues for 2008. Veterinary Record. 2008;162(10):298. pmid:18326840

43. Sellers RF, Mellor PS. Temperature and the persistence of viruses in Culicoides spp. during adverse conditions. Revue scientifique et technique (International Office of Epizootics). 1993;12(3):733–755.

44. Braverman Y, Rechtman S, Frish A, Braverman R. Dynamics of biting activity of C. imicola Kieffer (Diptera: Ceratopogonidae) during the year. Israel Journal of Veterinary Medicine. 2003;58(2/3):46–56.

45. Verhoef FA, Venter GJ, Weldon CW. Thermal limits of two biting midges, Culicoides imicola Kieffer and C. bolitinos Meiswinkel (Diptera: Ceratopogonidae). Parasites & vectors. 2014;7(1):384.

46. Wittmann EJ, Mellor PS, Baylis M. Effect of temperature on the transmission of orbiviruses by the biting midge, Culicoides sonorensis. Medical and veterinary entomology. 2002;16(2):147–156. pmid:12109708

47. Backer JA, Nodelijk G. Transmission and control of African horse sickness in the Netherlands: a model analysis. PLoS One. 2011;6(8):e23066. pmid:21850252

48. Ritchie-Dunham JL, Méndez Galván JF. Evaluating epidemic intervention policies with systems thinking: a case study of dengue fever in Mexico. System Dynamics Review. 1999;15(2):119–138.

49. García-Lastra R, Leginagoikoa I, Plazaola JM, Ocabo B, Aduriz G, Nunes T, Juste RA. Bluetongue virus serotype 1 outbreak in the Basque Country (Northern Spain) 2007–2008. Data support a primary vector windborne transport PLoS One. 2012;7(3). pmid:22479628

50. Calvete C, Estrada R, Miranda MA, Borrás D, Calvo JH, Lucientes J. Modelling the distributions and spatial coincidence of bluetongue vectors Culicoides imicola and the Culicoides obsoletus group throughout the Iberian peninsula. Medical and Veterinary Entomology. 2008;22(2):124–134. pmid:18498611