An ABM of biting midge dynamics to understand Bluetongue outbreaks

Bluetongue (BT) is a well-known vector-borne disease that infects ruminants such as sheep, cattle, and deer with high mortality rates. Recent outbreaks in Europe highlight the importance of understanding vector-host dynamics and potential courses of action to mitigate the damage that can done by BT. We present an agent-based model (ABM), entitled MidgePy, that considers the actions of individual Culicoides spp. biting midges and attempts to understand their role as vectors in BT outbreaks. Sensitivity analysis is performed and results indicate that midge survival rate has a significant impact on the probability of a BTV outbreak as well as its severity. This suggests that future methods to control BT spread could combine large-scale vaccination programs with biting midge population control measures such as the use of pesticides.


Introduction
In the late 2000s and early 2010s, multiple Bluetongue (BT) disease outbreaks occurred across Northern Europe, a region which did not previously witness BT spread. The resulting epizootic, fueled by an abundance of immunologically naive livestock, lasted multiple years until it was finally brought under control through mass vaccination programs and the restriction of livestock movement across country borders (Gethmann et al. [2020]). The epizootic had a significant economic impact on the livestock industry, and its persistence represents a continued threat to the European agricultural economy (Conraths et al. [2009]).
Occasional epizootics occur in any given year now and result in losses amounting to hundreds of millions of dollars due to necessary vaccination programs or lost revenue (Mayo et al. [2020]). Furthermore, livestock infected with BT experience high morbidity and mortality rates (Saegerman et al. [2008]). Understanding how BT arrived in Northern Europe is of interest to help prevent or mitigate potential outbreaks in new regions. In fact, BT is now considered a prime example of a disease whose range is changing due to climate change bringing warmer temperatures farther north (Brand and Keeling [2017], Samy and Peterson [2016]).
BT is a vector-borne disease caused by Bluetongue Virus (BTV) that is transmitted by Culicoides spp. biting midges (Diptera: Ceratopogonidae) and affects a wide range of ruminants such as sheep, cattle, and deer. Thus, the dispersal of Culicoides spp. biting midges plays an integral role in the spreading of BT between farms which are geographically distant (Pedgley and Brooksby [1983], Alba et al. [2004], Maarouf [1990, 1989], Chapman et al. [2010]).
Previous mathematical models have focused on various aspects of the spread of BT among animal populations. Using historical wind data, Sedda et al. [2012] retroactively simulated the spread of the 2006 BTV outbreak in Northern Europe by considering a spatio-temporal environment. Gourley et al. [2011] used a delay differential equation model to investigate reproduction numbers for BT. A model that incorporated a temperature-dependent incubation period for BT was designed by Li and Li and Zhao [2019] to predict the reproduction ratio and perform sensitivity analysis. Guis et al. [2012] developed a climate driven model, which they used to predict temporal changes in R 0 , and subsequently predict an increase in the future risk of BT across Northern Europe. Additionally, Turner et al. [2012] modeled farm-to-farm BT spread across Eastern England, which also considered the seasonal effect on transmission rate. Gubbins et al. [2008] performed sensitivity analysis on a temperature-dependent model considering a population with two host species to predict the risk of BT in Great Britain. Szmaragd et al. [2009Szmaragd et al. [ , 2010 designed a model that examined the effects of vaccine uptake among farms in Great Britain following the 2006 BTV-8 outbreak.
The use of agent-based models (ABMs) has been well-established for understanding dynamics of disease spread and the navigation of insects. For example, Smith et al. [2018] published a comprehensive review on the use of ABMs to model malaria transmission. These models were each designed specifically to either understand the role of environment heterogeneity in malaria spread, analyze the effectiveness of intervention strategies, or for parameter estimation. The flight patterns of insects such as moths (Bau and Cardé [2015], Liberzon et al. [2018], Stepien et al. [2020], Golov et al. [2021]), locusts (Topaz et al. [2008], Bernoff et al. [2020]), flies (Lin et al. [2015], Alderton et al. [2018], Leitch et al. [2021], Diouf et al. [2022), honeybees (Dorin et al. [2022]), and butterflies (Grant et al. [2018]) have all been studied with ABMs using a variety of modes including simple random walks and directed movement toward a target.
This paper aims to understand how aspects such as the number of initial infected midges, the daily survival rate of the midges, and the extrinsic incubation period of midges inoculated with BTV can affect the probability of an outbreak on a single farm. The outline of this paper is as follows: in Section 2, we provide biological information on BT and its spread. In Section 3, we describe an agent-based model, entitled MidgePy, that examines the spread of BT on a small-scale large mammal farm of approximately one square kilometer. In Section 4, we perform sensitivity analysis of Culicoides spp. survival rates and the extrinsic incubation period of BTV, as well as an analysis on outbreak probability. Finally, in Section 5, we summarize our results including their application to real-world efforts to combat BT outbreaks.

Bluetongue (BT)
Bluetongue virus (BTV) is a virus from the family Reoviridae, genus Orbivirus (Rivera et al. [2021]). BTV causes Bluetongue (BT), a hemorrhagic disease in ruminants that is transmitted by many species of Culicoides spp. biting midges (Mellor [1990(Mellor [ , 2000, ). There are more than 20 known serotypes of BTV that circulate between different regions including the Middle East, Europe, and North America (Gerbier et al. [2008], Maclachlan et al. [2015], Saegerman et al. [2008]). While BTV is transmitted in many parts of North America, there exists only one confirmed vector: C. sonorensis (Gerry et al. [2001]). This poses a conundrum as the primary range for C. sonorensis does not extend to the Southeastern United States, yet there have been confirmed cases of BT in that region (Smith and Stallknecht Smith and Stallknecht [1996]). Entomologists have identified potential BTV vector species in the Southeastern United States including C. stellifer, C. debilipalpis, and C. pallidicornis, however transmission of BT by these species has yet to be directly observed (McGregor et al. [2019]).
While vaccines exist that prevent BT, the virus is not immunologically simple and thus the vaccines must be serotype specific. With the large number of serotypes, each with varying virulence, it is difficult or impossible to vaccinate livestock against all strains (Noad and Roy [2009]). In addition, recent studies have shown that, due to the segmented genome of BTV, the use of live attenuated viruses to vaccinate ruminants has the possibility of introducing new genetic material to environments, increasing the risk of creating new BTV serotypes (Rojas et al. [2021]).
Transmission begins when a competent Culicoides spp. vector bites a viremic ruminant, ingesting blood that contains BTV. The virus then undergoes an incubation period inside the midge, during which it goes through multiple stages in the incubation cycle before eventually replicating within the midge and reaching its salivary glands (Mellor [2000], ). The length of this process is known as the Extrinsic Incubation Period (EIP) (Carpenter et al. [2011]). In Culicoides spp. vectors, this process typically lasts 14 days, however, it is significantly affected by the ambient temperature, which is accompanied by increasing Culicoides spp. activity (Tsutsui et al. [2011], Mayo et al. [2020], Tugwell et al. [2021]). Additionally, the biting rate in Culicoides spp. has been shown to be positively correlated with disease transmission rates (Elbers and Meiswinkel [2014]). Though it has recently been shown that horizontal transmission of BT between ruminants is possible (Maclachlan et al. [2019]), the goal of this paper is to focus solely on the vector-host transmission cycle, therefore, horizontal transmission of BT will not be considered here.

Methods
We present a spatially and temporally explicit agent-based model (ABM), entitled MidgePy, that characterizes the actions of ruminants and Culicoides spp. biting midges and the spread of Bluetongue (BT) among these populations. The spatial landscape of our study is based on a cervid farm in North Florida that has been the focus of previous BT  Erram and Burkett-Cadena [2018]). However, the spatial landscape in MidgePy is customizable so that outbreak likelihood given the introduction of BTV-infected midges to any region may be examined.
The MidgePy algorithm is depicted as a flow chart in Fig. 1, which we subsequently describe.

Model Description
MidgePy is comprised of two main classes of agents: a midge swarm class, named MidgeSwarm, and a ruminant swarm class, named HostSwarm. Each class holds arrays containing information on each agent such as their location and whether or not they are infected with BT. The length of the array indicates the population size of the respective agent type. In particular, the MidgeSwarm and HostSwarm classes contain arrays of lengths γP d and P d , respectively, where γ : 1 is the ratio of midges to ruminants in a simulation. See Table 1 for a list of model parameters and their values.
In our simulations, we set the spatial domain to be a continuous, homogeneous square of size  At the beginning of each day in the simulation, each ruminant is given a random position within the farm domain. To reduce computational complexity, it was assumed that the ruminant population's locations remained static.
At each time step, of length ∆t, which midges have recently fed is determined. We consider a midge to be not hungry if it had a blood meal from a ruminant within the last T full days, which we set to be 2 days. The midges that have recently fed (time since last meal < T full ) move in a simple random walk with roaming flight velocity v r . If a midge has not recently fed (time since last meal ≥ T full ), it will check if there is a ruminant within a radius of distance d det from itself, which we set to be 300 m. If there are no ruminants within that distance, the midge moves in a simple random walk with roaming flight velocity v r . If ruminants are located within a midge's detection distance d det , the midge will fly directly in the direction of the closest ruminant with active flight velocity v f until it is within a minimum biting distance d bite , after which the ruminant will be bitten. If the host is viremic with BTV, then BTV is transmitted to the midge with probability P HtoV . If the midge is infectious with BTV, then BTV is transmitted to the host with probability P VtoH .
After inoculation by a midge, a ruminant that had been infected by BTV would undergo an extrinsic incubation period for ρ H days, after which it would become viremic and could potentially infect other midges. Due to a lack of data providing a more precise value, we set ρ H to be 2 days. Similarly, a midge that had been infected by BTV would undergo an extrinsic incubation period, ρ. The model did not consider the ability of ruminants or midges to recover from BT, and as such they would remain viremic for the duration of the simulation.
Studies have shown that ruminants are often asymptomatic when infected with BT, and can remain viremic for several weeks (Singer et al. [2001]). Therefore due to the short simulation period and the long ruminant lifespan, even with BT, the model did not consider ruminant death. When assuming a constant daily survival rate for the midges, regardless of age, the lifespan of the midges followed an exponential distribution, consistent with laboratory findings (Lysyk and Danyk [2007]). Therefore the mean expected survival time of the midges was 1 1−α days, and at the end of each day simulated, it was determined whether each midge would survive to the next day with probability α. If the midge was determined to have died, it was replaced by a new midge in a random location that was not infected with BTV so that the population size of midges remained constant.
Each day was simulated to be T d = 5 hours long, which is approximately the amount of time during dawn and dusk where Culicoides spp. midges are most active (Fall et al. [2015], González et al. [2017], Blackwell [1997]). At the start of the simulation, an initial number of midges, I 0 , are randomly selected to be infected with BTV.
We note that while Culicoides spp. activity has been observed to be strongly temperature dependent (Tugwell et al. [2021]), in particular, that midges prefer to fly in warmer temperatures, and as such their flight activity changes throughout the year in temperate regions, we do not take this varying behavior into account in the model since North Florida is subtropical and experiences less seasonality.

Simulation Parameters
An increase in the midges-per-ruminant ratio γ significantly increased the computational cost when running the model, and so we ran simulations to determine a value of γ that would allow for a reasonable computational time yet be representative of the full population dynamics. Initial simulations were run setting γ = 30, and then increased to γ = 100. Simulations run with γ > 100 showed no significant differences. It was decided that γ = 100 was a suitable parameter value for this analysis.
Unless otherwise stated, the simulation period of the model was 60 days. This was chosen as the primary interest for this model was understanding outbreak dynamics, and so longer periods were not necessary.

Results
In this section, we analyze the sensitivity of the model parameters to the model output (Section 4.1), determine the dependence of the number of infected ruminants on Culicoides spp. survival rate and extrinsic incubation period (Section 4.2), and evaluate the dependence of outbreak probability on the initial number of infected midges and their survival rate (Section 4.3).

Sensitivity Analysis
We determine how sensitive the model is to changes in the Culicoides spp. daily survival rate, α, extrinsic incubation period, ρ, and initial number of infected midges, I 0 , using Sobol' sensitivity analysis (Sobol' [2001]) via Saltelli's extension of the Sobol' sequence (Saltelli [2002], Saltelli et al. [2010]) as implemented in the SALib Sensitivity Analysis Library (Herman and Usher [2017]). This is a variance-based method which allows for calculation of model output changes with respect to variation with a single parameter (first-order), and combinations of all parameters (total-order). By comparing the first-order and total-order indices, the presence of higher-order interactions can be inferred.
The parameter values selected for α and ρ were α ∈ [0.6, 0.9] and ρ ∈ [10, 20], which was a range of the typical values observed for these parameters (Table 1). The total number of infected ruminants after 60 days was used as the model output, and the resulting first-and total-order indices were calculated. Due to computational constraints and the stochastic nature of the model, a total of 96 simulations, with parameters given by the Saltelli distribution, were run for each trial.
To clarify terminology used in this section, we define a simulation to be a single execution of the model under a single set of parameters. A trial is a group of simulations. We averaged the results from many trials to find the expected values for the analysis in this section.
Both the first-order indices and total-order indices, shown in Figure 3, indicated that the model was much more sensitive to perturbations in the Culicoides spp. survival rate α than it was to perturbations in the extrinsic incubation period ρ by the taller bars. Thus, this analysis implies that focusing on the survival of midges is more important than the incubation period in reducing the number of infected ruminants.
Additionally, it was found that the model sensitivity to α decreased as I 0 increased, which was to be expected as more BTV-infected Culicoides spp. midges would reduce the need for high survival rates. Initial analysis was conducted using 10 trials for each I 0 , and then was increased to 15 trials. Little change was noticed, and so the analysis shown was conducted using 15 trials.

Heat Map
A heat map, shown in Figure 4, was generated to visualize the total percentage of infected ruminants over a 60 day period given changes in the Culicoides spp. survival rate α and extrinsic incubation rate ρ. Parameter values were chosen to be α ∈ [0.6, 0.9] in increments of 0.015, and ρ ∈ [10, 20] in increments of 0.5. For each combination of parameters α and ρ, 10 simulations of each case were run due to the stochastic nature of the model and the tendency for low values of α or high values of ρ to result in no outbreak at all. The average number of infected ruminants was averaged over all the simulations for each (α, ρ)-pair. with BTV. One interesting note of Figure 4 is that there appears to be a linear bifurcation line with respect to α and ρ that determines the likelihood of an outbreak. Low values of α combined with high values of ρ would result in few to no outbreaks, while high α and low ρ result in large-scale outbreaks. This follows intuition, as a high survival rate and low extrinsic incubation period would allow many more Culicoides spp. midges to become infected with BTV and subsequently pass the virus on to another ruminant.

Outbreak Probability
For the purposes of this model, an outbreak was defined to be the occurrence of a single BTV-infected ruminant. We fixed the initial number of infected midges, I 0 , and daily survival rate, α, for each simulation and would allow the simulation to run until either all BTV-infected midges died without a single ruminant becoming infected or a single α 3UREDELOLW\ I 0 Figure 5: The probability of outbreak given different initial numbers of infected midges, I 0 , depending on the daily survival rate, α. An outbreak was defined to be the infection of a single ruminant.
ruminant had become infected with BTV. The parameter α ranged from 0 to 1 in increments of 0.02, and 500 simulations were run for each value of α. From this, the estimated outbreak probability was calculated by averaging over all the simulations. While it is unrealistic to expect values of α < 0.5 or α > 0.9 in nature (Gerry and Mullens [2000]), they were still calculated as part of this analysis for completeness. Figure 5 indicates that for low values of I 0 , the relationship between α and the probability of outbreak appears to be linear. As I 0 increases, especially past I 0 = 15, the distribution of the probability of outbreak approaches a step function.

Discussion
The goal of this paper was to develop a model that represents the spread of BTV in a controlled environment. This model was investigated through methods that included sensitivity analysis and an examination of the probability of outbreak.
Results indicated that the model was more sensitive to variations in Culicoides spp. survival rate than to the extrinsic incubation period of BTV (Figure 3), which was reinforced by analysis of the probability of outbreak ( Figure 5). Additional visual analysis of the heat map in Figure 4 suggests that there exists a bifurcation when considering the relationship between survival rate and extrinsic incubation period, which specifies parameter regimes in which outbreaks are expected. As such, we conclude from this model that there exists a strong relationship between survival rate and the spread of BTV.
It is first worth noting the limitations of the model. As it is extremely difficult to accurately count the population of biting midges in a given system, the population of midges was set to be 100 times the number of ruminants. This ratio was chosen due to computational constraints. It is very likely that this is an underestimation of the total number of midges, and so the true spread of BTV would occur among a much larger population. Additionally, the model does not consider ruminants that ultimately recover or die from BTV infection, thereby resulting in a non-decreasing number of viremic hosts that will indefinitely contribute to the spread of BTV. Therefore it is possible that the model is overestimating the number of infected ruminants over long simulation periods (greater than one or two months).
Simulations of the model suggest that early in the outbreak, BTV could be in the system while not being detected. This would occur during a short window where there are no viremic ruminants or midges capable of spreading BTV, and instead, the virus is within a host animal or midge after inoculation has occurred by the virus has not yet completed its replication cycle. This suggests that it is much harder for ecologists and epidemiologists to detect and isolate early outbreaks of BTV in a region. We suggest that it would be possible to counter this issue by repeated testing of susceptible hosts over multiple days in at-risk areas.
A key focus of this model was analyzing possible outbreak scenarios. As such, the analysis has shown that survival rate of the midges plays the most important factor in determining the probability of an outbreak. Therefore it would be prudent that in regions at risk of an outbreak, the main focus should be to reduce the survival rate of vector species through means such as pest control. Other options include temporary removal or isolation of BTV host species, which most commonly are livestock used in agriculture.
A very interesting result of the model is that the outbreak probability increases with I 0 . While the trend appears linear for low values of I 0 , it was found that this trend becomes increasingly concave. It appears that in the limit as I 0 increases, the probability of outbreak would trend towards a step function. While this follows intuition, and helps to further validate the model, the results further suggest that high numbers of infected Culicoides spp. midges virtually guarantee the possibility of an outbreak. This implies that population control measures such as pesticides may not be enough to prevent epizootics, as long as there are enough vectors already infected with BTV to continue the cycle. Instead, alternative methods such as mass livestock vaccination would be more effective in stopping BTV spread.
Multiple studies have shown a strong seasonality with the growth and decay of BTV (Gerry et al. [2001], Carpenter et al. [2011]). It is known that BTV can lie dormant for multiple months during colder, more unfavorable conditions, then return during spring and summer where its prevalence reaches its peak. This model did not attempt to simulate outbreaks on large time scales, and so seasonality was not considered in the design. Instead, the model assumes constant parameters and consistent environmental variables. Additionally, one of the key goals of this model was to understand BTV spread in the southeast United States, particularly North Florida, where BTV remains endemic and does not spread through outbreaks like in Europe. In the southeast and North Florida in particular, the subtropical environment offers less seasonality than in more northern regions. As such, the lack of seasonality in this model is not as critical of a limitation as it would be in a model for more temperate regions.
Future work will focus on the development of a continuous model based on differential equations to identify a functional relationship between the Culoicoides spp. daily survival rate α and extrinsic incubation period ρ. While the sensitivity analysis indicates α to be more influential than ρ in determining the size of an outbreak, ρ still plays a significant role. In fact, there appears to be a linear relationship between α and ρ that shows how high survival rates may not be enough to guarantee an outbreak when combined with long incubation periods.
This model was designed with simplicity in mind. For that reason, it would be straightforward to modify the model in order to investigate the spread of other vector-borne diseases. Other future avenues of research using this model as a base include modeling malaria, which is spread through humans and other animal hosts. This would add a new level of complexity as humans move significant distances more frequently than livestock. Diseases such as malaria could be implemented into the model with relative ease, and sensitivity analysis could be performed to understand better methods for controlling spread. Other possibilities for future research include sensitivity analysis on other parameters such as vaccination rate, ruminant intrinsic incubation period, and probability of infection. Such analysis could supplement biological research on BTV-host interactions.

Supplementary information
The source code for MidgePy used to generate the results for this article is available through GitHub at https: //github.com/stepien-lab/MidgePy [v1.0.0]. The code is platform independent and written in Python. The data used for this article is available through OSF at https://www.osf.io/fhven.