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

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.


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

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 longrange 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) 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) 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) 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 O & R 3 and time interval [0, T]. We denote by Γ g and Γ s (t) the subsets of the boundary of O corresponding, respectively, to the ground and to the source of Culicoides entering into O at time t 2 [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) 2 O at time t 2 [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 C ¼ C s on fðx; y; z; tÞ : ðx; y; zÞ 2 G s ðtÞg; where C 0 is the initial distribution of the number density of Culicoides in O, C s (x, y, z, t) is the number density of insects at point (x, y, z) 2 Γ s (t) at time t 2 [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 n x ; n y 2 N. 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).
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]: where θ(z 0 ) is the temperature on the ground z 0 , and γ = 0.0065˚C/m. Horizontal wind speed values at high altitudes, z > 10 m, were estimated using the following polynomial variation [26]: where w(z 0 ) is the wind speed on the ground z 0 , 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-tobody 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: 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]: where R is the radius of the Culicoides, here assumed to be R = 2 Á 10 −3 m [32]; M is the weight of the Culicoides, here assumed to be M = 0.5 Á 10 −3 kg [33]; ρ is the density of the Culicoides, here approximated to an spherical particle ρ = M/(4/3πR 3 ); g is the gravity acceleration, assumed constant, g = 9.81 m/s 2 ; and ρ air (θ) is the mass air density at temperature θ, see Table 1 [34]. 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 u Ã ¼ jũ Ã j ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi jtj=r air p [36,37]; wheret is the surface shear stress, i.e. the vertical flux of horizontal momentum measured near the surface; which can be computed ast ¼ r air C d wðz 0 Þ 2 [38,39]. Thus, We may now use the dimensionless parameter lðzÞ ¼ V g ðzÞ u Ã 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: lðzÞ À 0:1 0:6 V g ðzÞ if 0:1 < lðzÞ < 0:7; V g ðzÞ if lðzÞ ! 0: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) 2 O and at time t 2 [0, T]. Therefore the mortality function is defined by 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].  (Table 2). 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 , 8t 2 [0, T]) and was defined as lying in M 2 N 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 2 {0, . . ., N − 1}, and the total amount of particles per simulation was N × M. All particles were assigned with an index Particle trajectories were estimated by considering the velocity field w through the explicit scheme: where X i n is the position of the Lagrangian particle with index i = {1, . . ., (n + 1) × M} at time step n = n(i), ‥, N, w i n is the wind velocity vector at position X i n at time step n and nðiÞ ¼ b iÀ 1 M c (i.e. n(i) is the greatest integer less than or equal to iÀ 1 M ). Then, the number density of Culicoides transported at step n by particle i (which was located at X i n ) was denoted as C i n . 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 Next, the evolution of this particle number density function was governed by the following scheme: where i = {1, . . ., (n + 1) × M} and n 2 {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.

Mortality and activity functions
where the value of a k (θ), for each k 2 {0, . . ., 4}, and b(θ) are given in Table 3, depending on the value of θ. 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    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.

Annual risk of introduction
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.
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. 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%.

Sensitivity analysis
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][12][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 Advection-deposition-survival model to assess the risk of introduction of bluetongue in Spain 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).

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  5), plotted in logarithmic scale, log(C + 1) (log1). https://doi.org/10.1371/journal.pone.0194573.g006 Advection-deposition-survival model to assess the risk of introduction of bluetongue in Spain 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 S1 File. Video simulation. Video file of particle movement in the air and ground concentration on Spanish territory during September 2004. (MP4)