The stochastic mechanistic model is first used to examine possible establishment and spread of V. velutina based on the discovery of an active nest near Tetbury, Gloucestershire. The simplest scenario is to examine the potential uncontrolled spread from this location, assuming that queens had already dispersed into the landscape before the nest was destroyed. Figure 1 provides a spatial representation of a single (successful) establishment. From a single nest (Fig. 1A) this realization produces 4 nests in the following year (Fig. 1B) which is slightly above the expected value, 16 nests by year 3 (Fig. 1C) and 129 nests by year 5 (Fig. 1D). After 10 years the invasion is predicted to be widespread, with over 50,000 nests in total and a high density of nests predicted in some regions, exceeding 5 nests per km2 (Fig. 1E). These simulations highlight the long-distance dispersal of V. velutina foundresses to form new nests, such that much of England and Wales can be colonized within a decade and carrying capacity is reached within 20 years – although parameter uncertainty and stochastic variability lead to huge variation between replicates (Fig. 1G). At carrying capacity (Fig. 1F) the impact of latitude and environment on the local abundance of nests become clearer; the early dynamics are, as expected, governed by dispersal and the expected number of successful foundresses per nest, whereas the long-term dynamics are governed by latitude and the suitability on the local habitat. The effects of latitude and the associated reproductive success is now explored in more detail in the Supplementary Material, and in the assessment of early invasion dynamics given below.

Figure 1 Dynamics of V. velutina nests in England and Wales showing its spatial spread and rapid increase in numbers. Maps A to E represent a single realisation of the spatial model with ecological parameters at the median of the MCMC posterior distributions based on the data from Andernos-les-Bains18. In maps A to D (which represent years 1, 2, 3 and 5) individual nest sites are marked with crosses, in map E (at 10 years) we show the density of nests in 1 km squares. Map F displays the long-term state (20 years following invasion), averaged across stochastic replicates and across parameter uncertainty. Graph G shows the temporal dynamics highlighting the mean (black), median (red) and 50% and 95% prediction intervals (shaded). (Map F and Graph G are from a thousand stochastic replicates with random draws from the posterior parameter distributions. Maps are generated from EEA Corine Land Cover data28 and displayed with bespoke software using Matlab29.) Full size image

The announcement that a second V. velutina had previously been found some significant distance from the Tetbury nest, gives extra weight to the assumption that V. velutina may first have invaded in 2015, yet remained undetected. To investigate this possibility further, we perform three pieces of simple spatial analysis to determine: (1) the most likely location of an initial invasion in 2015; (2) over what range should we expect to observe nests in the year of discovery (2016), and how many remain undiscovered; and finally, (3) over what range we should expect to observe nests in the year following discovery (2017). Unsurprisingly, the most likely position of any invasion in 2015 is directly between the two sightings (Fig. 2A) although with some degree of spatial uncertainty. Given that we assume a Poisson distribution of new successful foundresses per nest, simple probabilistic rules mean that the detection and destruction of V. velutina nests does not change the predicted distribution of undiscovered nests. If we assume an initial nest in 2015, then our model suggests that on average an additional 3.3 daughter nests should be dispersed into the wider environment, with a high (96%) probability that at least one nest remains undiscovered. The spatial distribution of these potential nest locations is over a very wide area, although the focus remains around the Bristol area (Fig. 2B) – our high-risk zone is defined such that all but one nest are expected to lie within the red contour. Taking these predictions into 2017 (Fig. 2C), we find that, while the focal location remains, the probability of finding a nest increases in both spatial extent and in magnitude. By 2017 the mean predicted number of nests in the absence of further control is approximately 17, and from this large number the chance of stochastic extinction is extremely low.

Figure 2 The likely spatial locations of nests in the years 2015, 2016 and 2017, and the sensitivity of results to the effect of latitude-on the numbers of nests. Given the two locations in which a nest and a hornet were found in 2016, map (A) shows the most likely positions of a common ancestor nest in 2015. Extrapolating forward from this probable location, maps (B and C) show the likely density of undiscovered nests in 2016 and new nests in 2017 respectively. The red contour is defined such that outside this (high-risk) region we would expect to find less than one nest. The inset histograms show the likely distribution of nests in each year, accounting for parameter uncertainty and stochastic variability. Graphs (D–F) explore the sensitivity of these findings to our assumption about the effects of latitude; in particular, we vary the linear function of latitude and plot on the x-axis the realised reproductive ratio at the likely location of the founder nest in 2015 relative to that inferred from the Andernos-les-Bains data (Supplementary Material). We show the expected number of nests (D), the chance that the invasion naturally dies out (E) and the area within the high-risk red contour (F) – blue lines correspond to Year 1 (2016), red lines to Year 2 (2017). (All results are calculated explicitly from the probabilistic rates in the Supplementary Material without having to use stochastic simulations. Maps are generated from EEA Corine Land Cover data28 and displayed with bespoke software using Matlab29.) Full size image

Such predictions are naturally affected by the assumed reproductive potential at the point of invasion, which in turn depends upon our assumptions about the impact of latitude. We have parsimoniously assumed the reproductive potential drops linearly from the values estimated from Andernos-les-Bains to zero in the north of England; Fig. 2D–F explore this assumption in more detail. Unsurprisingly, both the number of expected nests and the size of the high-risk zone increase with the assumed reproductive ratio, while the chance of natural extinction reduces. All have a non-linear dependence. Our default assumption is that the UK founder nest has approximately 38% of reproductive potential of a nest in Andernos-les-Bains due to its more northerly latitude and therefore colder climate. Any scaling that reduces the reproductive potential below 8% of the value in Andernos-les-Bains leads to the rapid extinction of V. velutina from Great Britain.

Given the huge damage that V. velutina can do to honey bee colonies and other important wild pollinators, an uncontrolled expansion would be disastrous, both in terms of the impact on the apicultural industry and the related decline in pollinator services. We therefore utilise our model to consider the impact of detection and destruction of nests in controlling a rapidly expanding invasion. We investigate the probability of eventual eradication for different probabilities of detection, and for different assumptions about the initial timing when probabilistic detection begins: in the first year (blue) or second year (red dashed) following invasion (Fig. 3). We assume that there is an independent probability of detecting (and subsequently destroying) each nest before it produces the next generation of founder queens. In practice this probability is the product of two independent probabilities: first detecting hornets feeding at an apiary, and second, discovering the nest by back-tracing the hornet flight paths. As anticipated, early initial awareness of nests (blue vs red) and/or high detection probabilities are required to achieve eradication, with detection probabilities in excess of 90% required to be confident of elimination. We note that this threshold represents an average over both stochastic realisations of the model and parameter uncertainty. The potential for the initial invasion to have occurred in 2015 (and hence for there to have been one year without control) would place the current invasion on the lower (red dashed) curves suggesting eradication would be more difficult.

Figure 3 The probability of V. velutina eradication following successive years of nest detection and destruction. All results represent the average over (ten thousand) stochastic replicates and capture full parameter uncertainty. Solid blue lines represent simulations where detection (and destruction) begin in the year of invasion – when there is assumed to be only one active nest; dashed red lines correspond to simulations where the invading nest remained undiscovered in the first year, and probabilistic detection only begins in the second year. Crosses are for simulations without additional radial detection, while open symbols represent differing radii of detection. In (A) we assume that once a nest is probabilistically detected a thorough search of the area will find 99% of all nests within a given radius; in (B) we assume these radial searches only find 48% of nests (comparable to the detection probability in Andernos-les-Bains (Supplementary Material)). In (C) and (D) we investigate the expected time to successful colonisation when the UK is subjected to a constant rate of invasion; detection rates in (C) and (D) correspond to the vertical lines in (A) and (B). Search radii of 8 km or more in (C) and 32 km in (D) lead to permanent exclusion of V. velutina. Full size image

As an additional element to this random detection of nests, we assume that once the presence of V. velutina is identified at a location, then localized radial searches will be conducted to find other nests. Assuming this can be done with great efficiency, such that 99% of all nests within the search radius are discovered, leads to substantial improvements in the chance of eradication (Fig. 3A). However, even with this intensive local searching, large radii are needed to ensure control if the primary chance of detection is low; this is because even in the early stages of an invasion nests can be widely dispersed. The findings from the Andernos-les-Bains region data (Supplementary Material) indicate that even when V. velutina are known to be in an area, only 48% (CI 33–61%) of nests are discovered while still active, such that their destruction halts the spread of V. velutina. When radial searches are reduced in efficacy to these values (Fig. 3B), their impact in controlling an outbreak is marginal except at the largest radii.

Invasion is unlikely to be an isolated incident. As the spread of V. velutina increases throughout mainland Europe, new invasions are set to become more likely. Mathematically, the control of these invasions can be treated independently, although in practice multiple invasions within a year may place additional stress on limited resources or lead to complacency once one nest is eradicated. Keeping with the independence assumption we extend the results to consider a stochastic rate of invasion and the expected time before a successful colonization. In both examples (Fig. 3C and D), limited local searching means that there is a finite (and often short) time until control efforts fail; however, with larger search radii (those not shown) permanent exclusion is possible.