## Abstract

The containment of genetically modified (GM) pollen is an issue of significant concern for many countries. For crops that are bee-pollinated, model predictions of outcrossing rates depend on the movement hypothesis used for the pollinators. Previous work studying pollen spread by honeybees, the most important pollinator worldwide, was based on the assumption that honeybee movement can be well approximated by Brownian motion. A number of recent studies, however, suggest that pollinating insects such as bees perform Lévy flights in their search for food. Such flight patterns yield much larger rates of spread, and so the Brownian motion assumption might significantly underestimate the risk associated with GM pollen outcrossing in conventional crops. In this work, we propose a mechanistic model for pollen dispersal in which the bees perform truncated Lévy flights. This assumption leads to a fractional-order diffusion model for pollen that can be tuned to model motion ranging from pure Brownian to pure Lévy. We parametrize our new model by taking the same pollen dispersal dataset used in Brownian motion modelling studies. By numerically solving the model equations, we show that the isolation distances required to keep outcrossing levels below a certain threshold are substantially increased by comparison with the original predictions, suggesting that isolation distances may need to be much larger than originally thought.

## 1. Introduction

Genetically modified organisms (GMOs) were first introduced in the late twentieth century. Since then, concerns have been expressed regarding possible negative effects of GMOs on human health and/or biodiversity [1,2]. These concerns have prompted regulatory bodies to put in place directives to control the use of genetically modified (GM) products. In the European Union (EU), if the content of GM material in a non-GM product exceeds the established tolerance threshold of 0.9%, the product must be labelled as containing GM material [3], which may affect the product's marketability [4]. In the EU, member states are empowered to take appropriate measures to avoid the unintentional presence of GM material in non-GM crops. The simplest measure for minimizing cross fertilization between GM and non-GM crops is to spatially isolate them. This has led member states to issue mandatory isolation distances, which vary considerably from one state to another. For example, GM and non-GM conventional maize crops can be 25 m apart in The Netherlands while they must be separated by at least 600 m in Luxembourg [5]. While such differences are partly due to political and societal opposition towards agro-food biotechnology, they also suggest a lack of scientific knowledge on the processes driving pollen transport.

The scientific challenge is highly nontrivial, as the speed of spread of a species (e.g. a species carrying a transgene) depend in large part on the frequency and distance of rare long-distance dispersal events [6]. Precisely because they are rare however, few are observed and so it is extremely difficult to estimate their frequency and size. Mathematical modelling is thus used to predict how far transgenic pollen will spread. Current modelling theory offers two options: Brownian motion and pure Lévy movement. These two options, however, have serious flaws: it is known that Brownian motion (the standard diffusion assumption) often seriously underestimates dispersal distances, while pure Lévy movement often overestimates dispersal distances. Policymakers and scientists are thus left with the unenviable task of determining appropriate isolation distances from two flawed models that disagree substantially.

For terrestrial plants, the transport of pollen is mostly driven by wind or insect pollinators. Some plants are exclusively wind-pollinated, while others are exclusively insect-pollinated, and still others can be pollinated via both means. Examples of wind-pollinated plants include grasses and conifers, which release pollen in huge quantities that float on the wind and often pile up in drifts on the ground. Canola (oilseed rape) pollen is moved by both wind and insects. Many plant species, however, do not rely on wind-pollination at all, and so require insect pollinators for reproduction. Plants in this category include fruit- and nut-bearing trees (e.g. apple, pear, almond), sunflowers and sweet potatoes [7]. At least 80% of flowering plants are directly dependent on insect pollinators [8], of which bees (honeybees in particular) are the most prolific. While pollination by wind is well understood and can be modelled using high-resolution atmospheric-dispersion models [9–11], modelling pollination by insects is more difficult [12]. This difficulty is partly due to the complexity of insect displacement patterns. Observations of honeybee flights have shown that up to 80% of these flights are less than 1 m in distance ([13], and references therein). However, bees can also travel much longer distances and some studies have measured flight distances of up to 4 km [14]. These displacements can be represented by a kernel function measuring the probability that a bee flies a certain distance. For oilseed rape pollen, which is dispersed both by wind and insects, Devaux *et al.* [15] have suggested that exponentially decaying kernel functions are not well suited to represent observed long-distance dispersal. Instead, they propose leptokurtic (heavy-tailed) kernels that have power-law type decay in the tails.

The conclusions of Devaux *et al.* [15] are completely in line with the observations of Lévy flight (henceforth, LF) patterns in the foraging movements of many living species [16]. LFs consist of frequently occurring short displacements with more occasional longer displacements. These longer displacements are in turn punctuated by even rarer and even longer displacements, and so on. Such displacement patterns can be represented by a Lévy dispersal kernel. Unlike exponential kernels, Lévy kernels have power-law asymptotic behaviour and hence heavier tails. The probability of a displacement of length *l* is drawn from a distribution that decays as , where 0 < *α* < 2. They are, therefore, well suited to represent movement which includes a non-negligible probability for very large displacements. Lévy-type foraging patterns have been observed for a broad range of organisms, including honeybees [17,18], bumblebees [19], fruit flies [20], large marine predators [21,22], human T cells [23], jellyfishes [24], albatrosses [25] and even human hunter–gatherers [26]. All of these Lévy patterns have been identified through sound statistical techniques, but there is still some debate on the underlying driving mechanism: natural selection to optimize search efficiency [27] or spontaneous emergence from innate behaviours [28]. The goal of this paper is not to elaborate that discussion but to investigate how the bee-mediated spread of pollen is altered when it is assumed that the distribution of foraging honeybee populations is well described by a movement paradigm of Lévy type.

Fractional diffusion equations have not seen widespread use in the ecological literature, but they yield the asymptotic behaviour of both LF and Lévy walk (henceforth, LW) processes (see appendix A for more details), and so are a useful tool in the study of organisms exhibiting Lévy-type dispersal. The main issue with pure Lévy dispersal kernels is that particles can perform arbitrarily large displacements that can go well beyond the physiological limitations of living organisms or landscape constraints. To avoid unrealistic displacements, it is necessary to truncate the tail of the dispersal kernel beyond a certain threshold. In that case, the power-law decay does not extend to arbitrarily large scales but is replaced by a quicker decay or cut off beyond a certain limit. Truncated Lévy processes were first introduced by Mantegna & Stanley [29] and Koponen [30] to eliminate arbitrarily large displacements. Exponentially truncated, also called tempered, Lévy processes were proposed by Rosiński [31] as a smoother alternative. In these processes, the sharp cut-off is replaced by a smooth exponential damping of the tails of the dispersal kernel.

The density of particles following a tempered LF is the solution of a tempered fractional diffusion equation [32–34]. In such an equation, the diffusion process is represented by a non-local integro-differential operator that results in a diffusion rate faster than second-order diffusion but slower than non-tempered fractional-order diffusion [35,36]. The use of a fractional diffusion model allows us to ‘scale-up’ LF foraging patterns of individual organisms and hence study the large-scale effect of long-range displacement on the dispersal and spread of a population. Such an approach has been used to study the large-scale impacts of LF displacements on the spread of an epidemic [37], the dispersal of invasive species [38] and on animal site fidelity [39].

In their review, Beckie & Hall [40] point out that few mechanistical models have been developed to model crop pollen-mediated gene flow. Since that time, statistical models to fit outcrossing data have been developed, but there has been little to no mechanistic modelling of insect-mediated pollen dispersal. Morris [41] paved the way for partial differential equation models by fitting his honeybee movement data to diffusion and advection–diffusion equations. In order to generate a leptokurtic pollen distribution, Tyson *et al.* [42] split the behaviour of the dispersing bee population into two states, corresponding to intensive and extensive searching, with transitions between both states. The authors show that this model fits the available bee dispersal data better than a diffusion or diffusion–advection model, and also that it provides a remarkably good fit to near-distance and long-distance pollen data taken together [42]. The authors do not, however, look closely at the model fit at the largest dispersal distances, and this omission provides us with an opportunity to investigate the relevance of other dispersal models.

In this paper, we derive tempered fractional-order diffusion models to simulate the spatial dynamics of foraging bees (§2), and the distribution of pollen dispersed by bees (§3). We show that the pollen-dispersal model, calibrated with both published parameter values and the data of Tyson *et al.* [42], predicts larger dispersal distances for pollen and provides a better fit to the data in the tails of the distribution (§4). We use this model to provide new estimates of required isolation distances between GM and non-GM crops and hence provide a more comprehensive assessment of the risk associated with GM pollen outcrossing (§5). We compare our results with those of earlier diffusion models, and finish with suggestions for future work (§6).

## 2. Modelling the anomalous dispersal of honeybees

The searching patterns of foraging honeybees have been studied extensively by Reynolds *et al.* [17,18,43]. They have conducted field experiments using harmonic radar to track the movements of displaced honeybees searching for their hive, or for alternative food resources after a known source of food has been removed. These experiments were conducted over a scale of several hundreds of metres. All of these experiments provide strong evidence that honeybees adopt an LF searching strategy to optimize their chances of success. This strategy is particularly favoured when honeybees have no access to waggle dance information from nest-mates, or when their navigational mechanisms are partially disrupted and fail to bring them home. By necessity however, the power-law scaling in the flight data must be truncated at some point. Reynolds *et al.* [43] estimate that this truncation takes place somewhere between 100 and 200 m from the hive. This truncation phenomenon is clearly visible in fig. 2 of [17] which shows that the distribution of flight lengths decays like an inverse square law, until a distance of about 100 m, after which the distribution decays exponentially.

The field experiments of Reynolds *et al.* [17,18,43], therefore, suggest that the displacements of honeybees could be represented by an exponentially truncated Lévy density. With such a density function, the algebraic decay of the tails is damped by an exponential factor. The asymptotic behaviour of the Lévy density is then for 0 < *α* < 2 and λ ≥ 0. The exponent of the distribution, *α*, controls the power-law decay of the tails. As the value of *α* decreases, the thickness of the tails increases and large displacements become more frequent. It is interesting to note that the Lévy density function reduces to the Gaussian density function when *α* = 2. Hence, an uncorrelated Brownian random walk (RW) can be seen as a particular case of a LF.

The truncation parameter *λ* controls the exponential tempering of the distribution tails. It has SI units of m^{−1} and thus allows for the definition of a truncation length , which represents the characteristic length beyond which algebraic power-law decay is overtaken by exponential decay. When *λ* ≠ 0, all the moments of the distribution are finite and the mean displacement length can be defined. It should be noted that the exponential damping of the Lévy density changes the area below the curve and, if the distribution is skewed, the mean is changed also. The exponentially damped Lévy density, therefore, needs to be renormalized and recentred. More details on these issues can be found in [34]. Figure 1 shows three sample LF trajectories whose displacement lengths have been drawn from a truncated Lévy distribution. The exponent of the Lévy distribution is set to *α* = 1.8 and only the truncation parameter *λ* is changed. It can be seen that the probability of large jumps is reduced as *λ* increases. When *λ* = 0, a pure LF is recovered.

At the macroscopic level, replacing an uncorrelated Brownian RW by a truncated LF amounts to replacing the second-order diffusion operator by a truncated fractional-order operator. In one dimension, the resulting diffusion equation for the bee density *B*(*x, t*) can then be expressed as follows:
2.1where *K*_{α} is the fractional diffusivity with SI units of m^{α} s^{−1}. If we assume that the model equation is defined for 0 ≤ *x* ≤ *L*, where *L* is the domain length, the truncated fractional diffusion term is then expressed mathematically as
The left- and right-sided fractional derivatives, and , are defined in the Caputo sense as follows
where and *α* is the largest integer not greater than *α*, i.e. *n* = 2 for , and *Γ* (.) is Euler's gamma function. Unlike ordinary derivatives that are local operators, fractional derivatives are global operators that depend on the entire function. This explains why fractional-order diffusion is faster than normal second-order diffusion. More details on fractional differential equations can be found in Podlubny [44].

Equation (2.1) can generally not be solved analytically and a numerical scheme is thus necessary for computation of an approximate solution. Since fractional derivatives are global operators, their numerical discretization is usually much more complex than for integer-order derivatives. In this study, we have used a finite-element scheme, whose details are given in appendix B. To illustrate the behaviour of the solution for different values of the parameters *α* and *λ*, we solved equation (2.1) with a Gaussian function as the initial condition (figure 2, black curve). When *α* = 2, we observe the classical diffusion pattern and the solution tails decay exponentially (figure 2, blue curve). As soon as the value of *α* is smaller than 2, the tails of the solution start decaying algebraically and hence the solution decays much more slowly than it does in the Gaussian case (figure 2, red curves). Asymptotically, the solution decays as . The spreading of the solution can be controlled by increasing the value of the truncation parameter *λ*. When *λ* = 0, there is no truncation and the solution thus corresponds to a pure Lévy process. By increasing the value of *λ*, we impose an upper threshold on the displacements beyond which large displacements become exponentially rare. By doing so, we can ‘tame’ the tails of the solution and control its spread.

## 3. Modelling bee-mediated pollen transport

Following the approach of Tyson *et al.* [42], we assume that pollen can either be motile or stationary. The former represents the pollen unwittingly taken up by honeybees as they gather nectar. Motile pollen sticks to body hairs on individual bees [45], and is thus transported by the bees as they fly from one blossom to another. As honeybees visit other flowers, some of the motile pollen is transferred to the stigma where it sticks and then becomes stationary pollen. Once deposited, pollen grains remain in the flower and can initiate the formation of seeds. Gene flow from a GM crop to a conventional crop will, therefore, be directly related to the distribution of stationary pollen originating from the GM crop.

As motile pollen is transported by honeybees, we assume that it has exactly the same spatial dynamics as bees and can thus be modelled by a fractional-order diffusion term. The input of motile pollen is assumed to be directly related to the distribution of flowers and hence the location of crop fields or orchards. Finally, the pollen deposition process that represents the transfer of pollen from the motile to stationary states is assumed to depend linearly on the motile pollen concentration. By taking only these three processes into account, we obtain the following model for bee-mediated pollen dispersal:
3.1and
3.2where *P*_{m}(*x, t*) and *P*_{s}(*x, t*) represent the concentrations of motile and stationary pollen, respectively. The diffusion term is exactly the same as the one used in equation (2.1). The function representing the source of motile pollen is denoted *F*(*x, t*), and can account for both the spatial and temporal variability in pollen sources. Time variability can appear, for instance, in different flowering periods. The constant deposition rate of motile pollen is denoted *β*.

It should be noted that equations (3.1)–(3.2) reduce to a classical pollen dispersal model when *α* = 2. Such a model relies on the assumption that bees carrying motile pollen follow a Brownian RW. The same model, therefore, allows us to compare dispersal patterns obtained with the classical Brownian diffusion assumption (*α* = 2) and with the anomalous Lévy diffusion assumption (*α* < 2).

## 4. Model calibration

Model calibration is based both on published parameter values and on outcrossing data collected in apple tree orchards with both GM and non-GM trees. We first consider published data to calibrate the parameters *α*, *λ* and *K*_{α} that govern the spatial dynamics of motile pollen. Reynolds & Frye [20] and Reynolds *et al.* [17] used harmonic radar to track honeybees searching for feeders or attempting to locate their hive when deprived of navigational clues. They found that the honeybee flights are composed of segments whose length follows a truncated power-law distribution with slope corresponding to a value of *α* close to 1. According to their observations, the power-law scaling is truncated somewhere between 100 and 200 m [43]. We have, therefore, taken *α* = 1 and . To estimate the diffusivity, we followed the approach of Morris [41] who expressed the constant in terms of a diffusion length () and time scale (). Since *K*_{α} has SI units of m^{α} s^{−1}, we can define the diffusion length scale as , where is the mean displacement length and is the corresponding time interval. According to the observations of Morris [41], it seems realistic to assume a mean displacement length of 10 m over a time interval of 1 h. In that case, we can define the diffusivity as s^{−1}.

Unlike the parameters *α*, *λ* and *K*_{α}, the deposition parameter *β* cannot easily be measured directly. We, therefore, estimate that parameter value through inverse modelling by using a pollen dispersal dataset collected in an isolated apple orchard [42]. Pollen dispersal was tracked via a transgenic marker. This marker was present in the pollen of one row of transgenic trees in the orchard, and then dispersed by honeybees to produce appleseeds throughout the rest of the orchard. Transgene presence in appleseeds was measured along eight transects in the orchard for two consecutive years. Since the orchard was arranged in four rectangular blocks separated by bare ground, a portion of each transect crossed bare ground (see Tyson *et al.* [42] for details).

The measurements from both years and all transects were joined together in order to obtain a more substantial dataset and hence have a better sampling of the distribution tails. The value of parameter *β* was then estimated by minimizing the least squares error between data and model outputs. This optimization process was performed with the *fmincon* function in Matlab. Since the canopy of orchard trees was approximately 7 m across, the model was run with a pollen input function representing a GM orchard for and conventional orchards for and , with the gap representing the average position of the bare ground portion of each transect. The rest of the domain is assumed to be bare ground where there is no production of pollen. Note that the absolute value of the pollen source term is not important as we are chiefly interested in the proportion of GM pollen among the entire pool of stationary pollen *P*_{s} (*Φ*^{GM}). Following Tyson *et al.* [42], this indicator is defined as
4.1where the superscripts ‘GM’ and ‘non-GM’ stand for transgenic and non-transgenic pollen, respectively. The factor of ‘’ in the numerator reflects the fact that only half of the pollen grains carry the transgene. We assume throughout that the proportion of pollen grains carrying the transgene at any one location is 50% of the pollen originating from the transgenic source, and so in (4.1) we can assume that, of those seeds fertilized with transgenic pollen, half carry the transgene.

To perform the optimization procedure over the parameter *β*, we assumed that the shape of the pollen distribution rapidly reached a steady state. This assumption is reasonable if the orchard is sufficiently small that the honeybees could rapidly spread out in a fairly uniform distribution throughout [42]. Initially, we assume that there is no pollen available, i.e. . The other parameter values are those obtained from the literature: *α* = 1, s^{−1} and . Since a comparison with a classical second-order model would be illustrative, we also consider the case where *α* = 2 and also compute the optimum value of *β* in that case. Model results for both the anomalous and classical diffusion cases are shown in figure 3. For comparison purposes, we also show in figure 3 the results for the model of Tyson *et al.* [42], in which bee movement was modelled using two subpopulations.

When *α* = 1, the best fit is obtained with . When *α* = 2, the best fit is obtained with . While both models provide a good estimate of the proportion of GM pollen close to the source (figure 3*a*), they differ more clearly further away (figure 3*b*). Indeed, beyond a distance of about 40 m from the GM orchard, *Φ*^{GM} quickly vanishes with the classical model (*α* = 2). This is a clear underestimation of pollen dispersal as in this particular dataset, the last non-zero measurement of GM pollen occurs at a distance of about 135 m. With the anomalous diffusion model (*α* = 1), it can be seen that the proportion of GM pollen decreases much more slowly and is still around 0.1% at a distance of 100 m from the GM source. When comparing with the results of the Tyson *et al.* [42] model, we see that the anomalous diffusion with *α* = 1 dispersal curve appears very similar to that obtained using classical diffusion with switching (cf. figure 3*a*(top and bottom panels). The tails however, exhibit very different behaviour (cf. figure 3*b*(top and bottom panels), with the solution from the Tyson *et al.* [42] model showing almost linear decay to zero, while the solution from the anomalous diffusion model decays much more gradually and remains above zero at distances where the Tyson *et al.* [42] model solution becomes negligible.

Finally, it is important to note that the point of non-zero dispersal at about 135 m is due to two seeds out of 752 seeds gathered at that distance. Further work is clearly needed to determine more robust estimates of pollen dispersal distance and pollination effectiveness as a function of dispersal distance. Nonetheless, all the diffusion parameters of the model (i.e. *α*, *λ* and *K*_{α}) were determined from published bee movement data and not from this particular dataset. The only parameter that has been determined from the apple pollen outcrossing data is the pollen deposition rate *β*. This approach ensures that our results are not simply a reflection of the idiosyncrasies of one particular study, but are based on several field studies and so have some robustness built in. In particular, the point of non-zero dispersal at about 135 m does not have undue influence on the parameter values and predictions of the model. The fact that the Lévy model is able to accurately represent the tails of the pollen distribution is, therefore, a qualitative validation of the model.

## 5. Estimating isolation distances

Isolation distances are generally prescribed between GM and non-GM crops to prevent, or at least minimize, conflicts due to transport of GM crop material into neighbouring habitats, conventional fields and/or protected areas. The isolation distance is defined as the minimal required distance between two neighbouring fields or orchards to keep outcrossing below a certain threshold. In the EU, isolation distances are meant to prevent the GM content in the neighbouring conventional field from exceeding 0.9%. For organic crops, this threshold is further reduced to 0.1%. When that threshold is exceeded in a non-GM field, food and feed produced from that field can no longer be declared ‘GM-free’.

Estimating the isolation distance that keeps outcrossing below a given threshold is obviously a difficult task. It involves a delicate compromise between outcrossing risk minimization and cost effectiveness. Indeed, any measure exceeding what is necessary to ensure compliance with the legal threshold would put an extra burden on farmers wishing to adopt GM crops. The adoption of GM crops being a very sensitive political and societal issue, the definition of isolation distances by local government is rarely based solely on scientific knowledge. We nonetheless believe that providing quantitative models that provide a clear estimate of the risks at stake could help clear up the debate.

In this section, we use the model defined by equations (3.1)–(3.2) with the parameter values presented in the previous section to evaluate the isolation distance between a GM and a non-GM orchard in order to keep GM input below a certain threshold. We consider an idealized test case where both orchards have a length of 50 m and are in perfect flowering synchrony. This is, of course, the worse-case scenario as GM crops can often be ‘temporally isolated’ from conventional ones simply through non-overlapping flowering periods. An illustration of the steady state stationary pollen distributions computed with our model as well as the proportion of GM pollen *Φ*^{GM} are shown in figure 4. The isolation distance between both orchards is denoted *l**.

Starting from a set-up where both orchards are next to each other, we can gradually increase the isolation distance *l** and compute the average proportion of GM pollen in the non-GM orchard . That quantity is again computed for *α* = 1 and *α* = 2 in order to highlight how Lévy and Brownian models evaluate the outcrossing risk. Figure 5 shows the evolution of in terms of *l** for both models. By taking the ratio between the values of *l** computed for *α* = 1 and *α* = 2, we can estimate by how much the spatial offset between GM and non-GM orchards should be increased when long-distance pollen flow is taken into account. For orchards of the same size, this ratio increases from about 5.5 to 7.5 for values of ranging from 1% to 0.01%. Such a large increase of the spatial offset requirements highlights the importance of the pollen distribution tails on the outcrossing risk.

Figure 5 also shows at what distances the ‘GM-free’ thresholds are reached if both orchards are the same size. When *α* = 2, our model predicts isolation distances of about 13 m and 27 m for conventional and organic crops, respectively. These distances are in agreement with the distances set by the most liberal European member states. However, when *α* = 1, the thickness of the pollen distribution tail increases and isolation distances of about 72 m and 164 m become necessary. Isolation distances for different orchard size ratios are provided in table 1. These distances are still within the range of isolation distances in application in Europe but they approximately represent a sixfold increase over the predictions of the Brownian model. Our results suggest that the separation distances of several hundreds of metres proposed by some European countries might be unnecessarily large. However, our results also suggest that separating conventional-crop fields by 40 m, as advocated by Riesgo *et al.* [46], might not be sufficient to keep GM outcrossing below the 0.9% threshold.

## 6. Conclusion

GM risk assessment is a central part of the decision-making process surrounding the regulation of GM products. It aims at minimizing environmental impact while preserving the potential economic gains that GM products could bring to farmers. In order to assess the feasibility of coexistence between GM and non-GM crops, quantitative tools are needed to estimate the degree of cross-pollination. Determining appropriate isolation distances has been difficult, however, as the two main modelling paradigms currently available generally provide very different results: Brownian motion models can seriously underestimate the potential for cross-pollination, while pure Lévy movement models have the opposite problem. In this paper, we have shown that a exponentially truncated Lévy movement model could provide policymakers and scientists with a much superior tool. We have focused on bee-mediated pollination, particularly by honeybees, and developed a mechanistic model that explicitly accounts for long-distance pollination.

A number of studies have indeed shown that honeybees can execute displacements over a scale of several hundred metres when looking for alternative food sources or attempting to locate a displaced hive [17,18,43]. These long displacements may result in long-distance pollination and hence increase the risk of adventitious mixing between GM and non-GM crops. Our model explicitly takes that risk into account by allowing long-distance dispersal of GM pollen by pollinating insects. This long-distance dispersal is achieved by replacing the classical assumption that pollinating insects follow a Brownian RW with the broader assumption that they follow a truncated LF. By doing so, we derived a general fractional-order pollen dispersal model of which classical models based on the Brownian motion assumption are an overly conservative special case.

The difference between Brownian and Lévy pollen dispersal models is most apparent in the pollen distribution tails. Removing the Brownian motion assumption increases the amount of long-distance pollen dispersal and hence increases the thickness of the pollen distribution tails. Since more GM pollen is transported faraway from the source crop, our model suggests that isolation distances may be underestimated by models relying on the Brownian motion assumption. Although the dispersal kernel of the Lévy model has heavy tails, it is important to note that all the moments of the kernel are finite. This means that there is a clear upper bound on the distance over which pollen can be transported. Hence the likelihood of outcrossing as predicted by the Lévy model eventually decreases to zero, but the initial decay rate is much slower.

We note here that our anomalous diffusion model also provides a strong alternative to the correlated composite random walk (CCRW)-type model of Tyson *et al.* [42]. There is interest in considering the CCRW paradigm as a replacement of the Lévy movement paradigm; this paper provides a useful comparison of the two in the context of honeybee-mediated pollen dispersal. The best-fitting CCRW model shows essentially no pollen dispersal past a distance of ≈125 m, consistent with the notion that there is a maximum biological dispersal distance, but the data contain a point of non-zero dispersal at ≈135 m. Since accurately estimating long-range dispersal distances is crucial in matters involving GM agriculture, it is important that models capture non-zero dispersal probabilities. In addition, measurement of rare events such as long-distance dispersal distances is difficult and, hence, the actual likelihood of a dispersal event between 135 m and 200 m may be larger than the data suggest. The truncated anomalous diffusion model may provide a more accurate prediction of dispersal distances in this case.

Since models remain approximations of reality, it is noteworthy to point out some potentially important biological details that we have not taken into account. One of them is the influence of the resource distribution on the bees’ foraging patterns. The distance that bees travel depends in a nonlinear fashion on the amount of pollen available. When pollen is plentiful, bees will tend to stay close to the hive, but will travel further when pollen is scarce. However, the Lévy process does not explicitly depend on the spatial distribution of resources. Benhamou [47] and Plank & James [48] have shown that truncated power-law distributions could be generated by other behavioural processes, such as area-restricted search that depends on the memory of locations where resources have been found in the past. Here we did not study the underlying generative processes driving bee foraging patterns. A valuable extension of this study would be to establish how bees respond to resource distributions by following the same approach as the one used by Hills *et al.* [49] to study human foraging.

Other assumptions include the absence of a diurnal cycle for bee activity, and the absence of any homing behaviour that would lead to a higher bee concentration in the vicinity of the hive. The effectiveness of pollination may also decrease with distance from the source, a factor that would reduce the effective dispersal distance. In this study, we assume that all flowers are equally attracting so that bees do not favour a particular part of the modelled region. During their flights, bees obviously lose some pollen, which is therefore not converted into stationary pollen. In addition, bees harvest pollen, storing it in their pollen baskets. We assume that the removal of pollen through loss or harvesting is a spatially uniform effect, and, therefore, does not affect the shape of the distribution of motile pollen available for conversion to stationary pollen. We can thus assume that motile pollen is entirely converted into stationary pollen. Finally, the impact of meteorological or anthropogenic factors on the activity of bees is also neglected. These are factors that would be worth considering in future work.

## Competing interests

We declare we have no competing interests.

## Funding

Funding for this research was provided by NSERC (R.C.T.). E.D. is an honorary Research Associate with the F.R.S.-FNRS.

## Acknowledgments

We gratefully acknowledge the assistance of Okanagan Specialty Fruits who paid for the data presented in figure 3.

## Appendix A. Lévy flights and Lévy walks

It is important to note that Lévy flights (LF) and Lévy walks (LW) are often used interchangeably in the biological literature. Both are random walk (RW) processes with a Lévy probability distribution of step lengths. They are not, however, identical. LF and LW are the heavy-tailed equivalents of the usual uncorrelated and correlated Brownian RW. The former are usually called position-jump processes, while the latter are called velocity-jump processes. In a LF, the particle relocates instantly and hence moves with an infinite velocity while in a LW it moves with a finite velocity. LW are thus slower and seemingly more realistic than LF. However, in the long-time limit, LW become similar to LF as they have more and more time to perform long jumps [50].

For a Brownian process, the density of particles following the uncorrelated RW is the solution of a second-order diffusion equation characterized by an infinite propagation speed. For particles following a correlated RW, the density is the solution of a telegraph equation for which the propagation speed is finite. In the long-time limit, both the uncorrelated and correlated RW processes converge to the same diffusion limit. This correspondence explains why the diffusion equation, and not the more complex telegraph equation, is so widely used to model population redistribution. The same pattern holds for LF and LW. The density of particles following an uncorrelated LF process is the solution of a fractional-order diffusion equation [35]. The density of particles following a correlated LW process is the solution of a more complex fractional-order equation akin to the telegraph equation. The fractional diffusion equation yields the asymptotic behaviour of both LF and LW processes [51].

## Appendix B. Numerical discretization of the tempered fractional diffusion equation

In this section, we present the details of the numerical scheme used to solve equation (2.1). The exponential tempering of the fractional derivative requires a specific numerical treatment inspired by the work of Hanert & Piret [52]. In this study, we have used the finite-element (FE) method to discretize the model equation. With that method, the exact solution is approximated by an expansion in terms of piecewise linear basis functions *ϕ*_{j}(*x*):
In that expression, the unknown nodal values are an approximation of the solution value at a node *j* of a grid (or mesh) that partitions the computational domain into elements. In one dimension, these elements are simply segments of arbitrary size. To handle the exponential tempering, we introduce the following intermediate variables and assume that they have the same discretization as the model solution:
B 1and
B 2where and are the unknown nodal values corresponding to the left and right tempering of the solution.

With the FE method, it is common practice to use the Galerkin method to derive the discrete equations. They are obtained by introducing the discrete solution in equation (2.1) and orthogonalizing the residual with respect to the set of basis functions. The following discrete equations are then obtained:
for *i* = 1, … , *N*. If we replace , and by their expansions, we obtain
B 3where we have introduced the matrices , and . Since the basis functions *ϕ*_{j} are piecewise linear, it is possible to compute their fractional derivative analytically (see Hanert [53] for details).

In order to be able to solve equation (B 3), we still have to express the tempered nodal values and in terms of . To do so, we use the following *L*_{2}-projection which simply amounts to applying a Galerkin formulation to equations (B 1)–(B 2):
and
where we have introduced two additional matrices: and . These matrices allow us to express the vectors of nodal values and in terms of as follows: and . The semi-discrete equations can finally be expressed in matrix form as follows:
B 4The last step is to discretize equation (B 4) in time. This can be done with any standard time integration scheme. In this study, we have used a third-order Adams–Bashforth scheme.

- Received November 6, 2016.
- Accepted January 3, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.