The “ Lévy or Diffusion ” Controversy : How Important Is the Movement Pattern in the Context of Trapping ?

Many empirical and theoretical studies indicate that Brownian motion and diffusion models as its mean field counterpart provide appropriate modeling techniques for individual insect movement. However, this traditional approach has been challenged, and conflicting evidence suggests that an alternative movement pattern such as Lévy walks can provide a better description. Lévy walks differ from Brownian motion since they allow for a higher frequency of large steps, resulting in a faster movement. Identification of the ‘correct’ movement model that would consistently provide the best fit for movement data is challenging and has become a highly controversial issue. In this paper, we show that this controversy may be superficial rather than real if the issue is considered in the context of trapping or, more generally, survival probabilities. In particular, we show that almost identical trap counts are reproduced for inherently different movement models (such as the Brownian motion and the Lévy walk) under certain conditions of equivalence. This apparently suggests that the whole ‘Levy or diffusion’ debate is rather senseless unless it is placed into a specific ecological context, e.g., pest monitoring programs.


Introduction
Pests form a significant threat to agricultural ecosystems worldwide, and therefore effective and reliable monitoring is required to ease the decision making process for intervention.In agro-ecosystems, monitoring is an essential component of integrated pest management programmes (IPM) [1,2], where a control action is implemented if necessary.If the population abundance exceeds a certain predefined threshold level, that is provided the effort and expense is justified, then intervention becomes imminent.Usually the control action takes the form of pesticide application which has many negative implications such as environmental damage in the form of air, soil and water pollution.Such human induced pressures on the environment, often contribute towards bio-diversity loss and affect the functioning of ecosystems [3][4][5].Other major drawbacks which are not necessarily related, include; cancer related diseases for those handling such chemicals [6,7], increased consumer costs [8], poor efficiency in reaching targeted pests [9], pest resistance to regular use [10] and lethal effects on natural enemies [11], possibly leading to a resurgence of the pest population or a secondary pest to emerge.Therefore, in order to avoid unnecessary pesticide application or risk of triggering pest outbreaks, accurate evaluation of population abundance is key [12].Traps are usually installed in the field under controlled experimental conditions as a means to estimate population abundance [13][14][15].They are then exposed for the duration of study, insects are caught, traps are normally emptied at regular intervals and the total number which fall into the trap are counted and form trap counts.It is precisely these counts that are converted to the pest population density at trap locations, and then used to estimate the total pest population size [16,17].
A major ecological challenge is to develop relevant theoretical and mathematical models which can explain patterns and observations obtained from field data [18][19][20].This is primarily due to the fact that inherent complexity found in the behaviour of animals can be difficult to incorporate.However, insects and some invertebrates are easier to model since they have been thought of as non-cognitive.In the case of either single or multiple traps in the same field, individual insect movement can be modelled successfully using a random walk framework [12].The earliest attempts were based on Brownian motion, which provided a framework to characterize patterns of movement with broad applications to conservation [21,22], biological invasions [23][24][25] and, in particular, insect pest monitoring [12,[26][27][28].Theoretical arguments supported by empirical observations suggest that individuals with limited sensory capabilities tend to follow a Brownian movement pattern, more so at large temporal and spatial scales [29][30][31][32][33][34].The corresponding mean field counterpart describes the spatial-temporal population dynamics which is governed by the diffusion equation [12,27,28,35].Recently, both Brownian motion (BM) and the diffusion model have been often criticized and deemed to be an oversimplified description [36].Other revised models have attempted to account for possible intermittent stop-start movement [37] or inherent intensive/extensive behavioural changes [38].Simultaneously, there are also other studies with growing empirical evidence, that support an alternative description which postulates that animal movement can exhibit Lévy walking behaviour [39,40].Lévy walks 1 (LW) differentiate from Brownian movement, since they allow for arbitrarily long step lengths, that is; the probability of executing a larger step is much higher -which results in a faster movement pattern altogether [26,27,42,43].
The Lévy or diffusion controversy has arisen from ongoing debates which provide pro-and-contra arguments for each description [43].Some cases that provide promising evidence for Lévy type movement [44] have been later classified as Brownian [42], and then have been reclassified as Lévy afterwards [45].The confusion arises partly due to different studies providing mixed and often conflicting messages.For example, the movement pattern can switch from Brownian to Lévy in a context specific scenario where resources are scarce [46].In another study, Lévy type characteristics can emerge as a consequence of the fundamental observation that individuals of the same species are non-identical [47].It is also possible that the underlying movement pattern can be misidentified, since variation in individual walking behavior of diffusive insects can create the impression of a Lévy flight [48].Also, composite correlated random walks can produce similar movement patterns as Lévy walks, therefore current methods fail to reliably differentiate between these two models (although recent attempts have been made to address this issue [49]).On the other hand, diffusive properties can appear for a population of Lévy walking individuals due to strong interactions e.g if movement is stopped when individuals encounter each other [50].Even more recently, the diffusion and Brownian framework has been revisited and shown to be in excellent agreement with field data [28].In either case, it is still unclear which type of movement is adopted by insects, what conditions can alter the pattern and which mathematical framework is most efficient -hence the controversy persists.
The idea of using a modified diffusion model as an equivalent framework for Lévy walking individuals was first introduced by Ahmed [35] and later further developed [26,27].In particular, it was demonstrated that in the context of pest monitoring, trap count patterns could be reproduced when comparing a type of Lévy walk to time dependent diffusion.Further discussions have highlighted that the ecological basis behind incorporating time dependent diffusion is not clearly understood, and how this is linked to the type of mechanisms [26,51,52].Our study has been instigated by these fruitful discussions leading to an extension to the previous study by Ahmed and 1 To avoid confusion, it must be noted that what is called a Lévy walk in ecological literature often corresponds to a Lévy flight in physics [41].
Petrovskii [27].Our aim is two-fold; that is to propose a diffusion model which consists of parameters that are of ecological significance and can be interpreted with reference to the diffusive properties.Also, to investigate if trap counts can be effectively reproduced for a broader class of Lévy walks using diffusion.If so, we question the relative importance of identifying the underlying movement pattern?From a cost perspective, it may make sense to concentrate more on the geometry and design of the experiment rather than the particular movement model.

Brownian Motion
Individual based models provide a complementary cost-effective methodology to field experiments, and can be used to simulate movement and analyse trap counts [12,18,53].The idea is to replicate these experiments through a virtual environment where supplementary and even alternative information is sought [54].In the specific case of low-density populations, the magnitude of stochastic fluctuations can be quite large, and an individual based modelling framework can be particularly useful to describe the movement dynamics i.e if dangerous pests are present.Since our interest is primarily based on understanding the underlying mechanisms which govern movement, it is sufficient to focus on a 1D conceptual scenario.Despite the fact that this case is hardly realistic in terms of modelling movement in a real field setting, it does provide a theoretical background for the more realistic 2D case [55].Also, any unnecessary additional complexity which would arise due to the effects of trap and field geometries are then avoided.
The basic idea in 2D is to model walking or crawling insect movement along a continuous curvilinear path which can be mapped to a broken line with position recorded at discrete times [18].In mathematical terms, the 1D analogue over unbounded space for a population of N individuals, records the position X (n) i of the n th individual at time t, which is discretized as t i = {t 0 = 0, t 1 , t 2 ..., t S = T}, i = 0, 1, .., S, n = 1, 2, .., N, so that Here, the script (n) is included to differentiate between movement tracks of different individuals of the same insect.We assume that positions are recorded at regular intervals with constant time increment where the first observation is recorded at time t = 0 with total number of steps S and total time t S = T.Each individual moves from position X and average velocity defined over each step.The n th individual placed at some position X (n) i can move either to the right or to the left with resulting position X (n) i+1 .Assuming that each subsequent step is completely random 2 and generated according to a predefined probability distribution φ(ξ), then each further position can be determined by The randomness of animal movement is obviously an idealization which, however, is well justified under certain conditions, e.g.see [18] for a detailed discussion of this issue.
where ξ is a random variable.Trap counts can then be obtained if the individuals which fall within a well pre-defined region are removed from the system at regular intervals and counted [12,35].Generally, animal movement is anisotropic due to the mere fact that animals have a front and rear end [56], resulting in a correlated random walk [30].For simpler movement modes, such as that of insects, we can assume that the walk is uncorrelated and the direction of movement is completely independent of previous directions.Each position is then solely dependent on the previous location, and therefore essentially Markovian [57,58].Under this assumption, the resulting movement is isotropic and there is no preferential direction i.e no advection or drift of any kind -which can arise in the presence of an attractant.In relation to the step distribution, φ(ξ) is a symmetric probability density function (pdf) with zero mean, that is; φ(ξ) = φ(−ξ).In the case of Brownian motion, the corresponding pdf is normal which reads with scale parameter σ, which can possibly be dependent on time, (the subscript n refers to normal).Alternatively, we write ξ ∼ N(0, σ 2 ) which denotes ξ is randomly drawn from a normal distribution with mean zero and variance σ 2 .In realistic ecological applications, many insects are released into the field instead of a single individual.Using the 1D scenario as a baseline case, the corresponding analogy is to initially distribute N individuals along a finite spatial interval x ∈ (0, L).Common release methods used in trap count studies are either of two types; that is uniform or a point source.
In the uniform case, we prescribe initial position as , where U(a, b) denotes the uniform distribution over the interval a < x < b.For a point source we have X (n) 0 = x for all n = 1, 2, ..., N at t = 0 with x = x as the centralised location of the release point.The resulting movement pattern is then completely determined by the type of initial condition and step distribution φ(ξ).Note that, the population distribution essentially becomes uniform for larger times, and therefore identical trap counts are obtained irrespective of the type of initial distribution.This effect is realised due to the inherent random movement of individuals in the field.For more detailed information on how the shape of the trap count profile is affected by the type of initial condition, the reader is redirected to [35].Generally, in most field applications, a uniform initial population distribution can be reasonably assumed.Even in the case when the true release distribution is characterized by multiple point source releases, the effect on trap count variation is somewhat minimal [28].Henceforth, all simulations in this paper adopt the uniform distribution as an initial condition.
To describe the individual based model fully, boundary conditions are enforced and defined as follows; an impermeable stop-go 'sticky' type boundary is installed at the external boundary at x = L [59], such that at any instant in time if the individual position exceeds this boundary; that is if X (n) i > L then it remains at location x = L.The next position in the process is determined purely by (1), and the individual continues to interplay with the dynamics of the system provided 0 < X (n) i+1 < L at the next step, otherwise it either remains at the external boundary, if The trap boundary at x = 0 introduces a perturbation to the movement, and is incorporated in the following way; if X (n) i < 0 at any instant in time then the individual is removed from the system and the total trap count increases by one, functioning as an absorbing boundary.Cumulative trap counts J t are obtained at each step over discrete times t = i∆t, resulting in a stochastic trap count trajectory.It may be important to mention that, the choice of time step ∆t has some significance, since it is known that time scale invariance is lost in the presence of an absorbing and/or stop-go boundary, possibly leading to noticeable differences in trap count recordings [59].Therefore, ∆t is chosen small enough so that trap counts are in line with counts obtained using alternative methods, such as the mean field analytic or numerical solutions (see §3.1 and Appendix A for details).Also, at the same time, ∆t is chosen to be sufficiently large, so that the assumption that subsequent steps are uncorrelated is feasible [12,31].

Condition of equivalence
In many ecological applications (such as, for instance, conservation and monitoring), it is important to know the probability that a given animal will remain inside a certain domain or area.Since the movement is described by the probability distribution function (pdf) φ(ξ).One can expect that the probability of remaining inside the domain depends on the properties of φ, in particular on its rate of decay at large distances.
We first consider the case where φ(ξ) is fat-tailed, φ(ξ) ∼ |ξ| −(α+1) with 0 < α < 2, which is the characteristic exponent for Lévy walks (see §4.1 for more details later).We focus on the special case α = 2 as it is thought, basing both on observations of movement patterns [60,61] and on some evolutionary argument [50] to be ecologically most relevant.In this case, the pdf for the Lévy walk is described by the Cauchy distribution: where γ is a scale parameter and the subscript c refers to Cauchy.In case at t = 0 the animal is at x = 0, function (3) gives the probability density of its position after one step.It is straightforward to see that the pdf of the animal position after i steps, i.e. at time t i = i∆t, is given by the same distribution but with a re-scaled value of the parameter γ, that is From equation ( 4) we may define a symmetric interval of interest x i (γ, ), such that the integral in this interval is always equal to a certain quantity , and determine the limits of this interval as a function of the parameters of the process and the arbitrary probability : Our goal is to obtain an alternative stochastic process, composed by the sum of random variables from the Gaussian family, that may be comparable to this one in the sense that it replicates the same probability over the interval of interest (e.g.see [28]).To do that, we will consider the sum Ȳ = Y 1 + Y 2 + • • • + Y i of normally distributed random variables, defined by the pdf where the variance ∆ 2 i is just the sum of the variances from each random variable Y k , As the relation between the two process is obtained by equating the integrals on the same domain [−x i , x i ], we may compute the probability over the Gaussian process using equation ( 7), and, introducing the inverse error function as Φ (•) ≡ erf −1 (•) (see footnote 3 ), we relate the two processes by the following relation: In order to determine the behavior of the increments σ 2 k (γ, ), i.e., the variance of each additional Gaussian variable required to be comparable to the original Cauchy process, we only need to write equation ( 8) as which immediately leads to the expression Now we recall that i is the time of the movement (measured in steps).From Eq. ( 12), we therefore conclude that the probability to confine the insect performing the Lévy walk over the spatial domain −L < x < L coincides exactly with the probability of the same event in the case where the insect performs the Brownian motion, provided the variance of the Brownian motion (i.e., essentially, the diffusion coefficient) increases linearly with time, σ 2 t ∼ 2t.

Time Dependent Diffusion
Insect movement is inherently more complex in nature, due to a contribution from both external and internal factors.Typical external factors include environmental effects or stimuli, which can be quite challenging to incorporate from a modelling perspective.Since our interest lies in the actual mechanisms at play, we assume homogeneity in the sense that external factors are absent.In terms of the underlying movement mechanisms, examples of internal factors include; individual variation, composite and/or intermittent movement or even time-density dependent diffusive behaviour [62,63].The main challenge is then to develop a coherent model which can include these different processes and accurately describe the population dynamics.Obviously, the issue becomes more difficult if a combination of these features are present.In the context of insect pest monitoring, diffusion models have shown to provide a good theoretical framework and the means for trap count interpretation [12,27,35].In particular, time dependent diffusion provides an adequate description for more complicated behaviour, at least where standard diffusion fails [26,27,35,64].The notion of insect movement with time dependent diffusivity is not new, and has been observed in field studies [19].
The 1D diffusion equation for the population density u(x, t) with time dependent diffusion coefficient D = D(t) over the semi-infinite domain 0 < x < ∞, with initial density u(x, t = 0) = u 0 (x) and zero density condition u(x = 0, t) = 0 at the trap boundary, reads By introducing a change of time variable Peer-reviewed version available at Mathematics 2018, 6, 77; doi:10.3390/math6050077 the system of equations ( 13) transforms to The general solution [26,27,65] is then given by where is the fundamental solution of the diffusion equation in (15), which reduces to in the specific case of constant diffusivity.The diffusive flux through the boundary at x = 0 corresponds to trap counts j(τ), which can be determined by Fick's law, that is j with cumulative trap counts J(t) (total flux) given by where erfc(z is the complimentary error function.Therefore, the total number of trap counts J(t) for the system (13) in normal time t is given by, In the case of a uniform distribution u 0 (x) = U 0 , (19) reduces to with in the special case with constant diffusivity D, corresponding to standard diffusion.

Equivalence of Trap Counts: Brownian Motion vs Diffusion in a Semi-bounded Space
For the diffusion equation ( 15), the mean location and mean squared displacement (MSD) are useful statistics that characterize the movement, where F(x, τ) is defined by (17).It is well known that diffusion is the macroscopic description of Brownian motion [66], where the MSD is equal to the variance of the step distribution φ(ξ), From this we obtain a link between the scale parameter in (2) and diffusion coefficient, For a discrete time model, one can expect that this remains valid, at least approximately, for a small but finite value of ∆t, that is In the case of standard diffusion, the MSD grows linearly with time and is related to the scale parameter by which is known as the hallmark of Brownian motion [12,18,67].More generally, for anomalous diffusion, the MSD grows according to some power law relationship where H is the Hurst exponent.Here, H = 1 2 corresponds to standard diffusion (26), 1  2 < H < 1 corresponds to super-diffusion and H = 1 corresponds to ballistic or wavelike motion.A full comprehensive summary of movement properties with reference to H can be found in [58].
To demonstrate equivalence between Brownian motion and an anomalous diffusion model, consider as a baseline case where D 0 is the initial diffusivity and D 1 controls the effect of time dependency for larger time.This structure is chosen as an example, so that the scale parameter ( 24) is in accordance with (27), i.e The analytical solution for the model with an initial uniform distribution can be derived from (20), which reads and approximates to the flux for constant diffusion (21) for small time.25) and (28), that is )∆t with D 0 , D 1 , H given above.Each individual executes a total of S = 5000 steps with constant time step increment ∆t = 0.001 and total time T = S∆t = 5.Individuals initially uniformly distributed X (n) 0 ∼ U(0, L = 5).Trap installed at position x = 0 and simulations are conducted with external boundary condition described in §(2.1).Trap counts are replicated and averaged over ten realizations to reduce the effect of stochasticity.Numerical solution: (Green dashed) represent the mean field numerical solution using the method of explicit finite differences.See Appendix A for further details on the numerical scheme.Plot (b) Absolute relative error: A(t) plotted at times t = hk, h = 0.1, k = 0, 1, 2, ..., 50, with average Ā = 0.306 (red), 0.157 (blue), 0.167 (black), to illustrate the magnitude of the discrepancy between the analytic solution and trap counts.(For interpretation of the references to colour in this figure legend, and all subsequent figures, the reader is referred to the web version of the article.) In Figure ( 1) Plot (a), we find that there is almost identical agreement between trap counts obtained from the Brownian individual based model and mean field diffusion model, as expected from theory.More specifically, it is shown here that the diffusive flux can be used to reproduce trap count patterns for standard diffusion H = 1 2 , super diffusion H = 3 4 and ballistic movement H = 1, to a high level of accuracy.Intuitively, we expect that this holds for diffusion coefficients which have a more complicated time-dependency.The diffusion coefficient consists of three unknown parameters, namely; D 0 , D 1 and H.In terms of usage, if initial diffusivity can be measured through experiments, then other parameters can be estimated using the tools outlined in [27], i.e by approximating the flux rate in the limit t → 0 and relating it to the expected number of individuals trapped after one time step.Note that, since an analytical solution cannot be obtained for the diffusion model ( 15) over a finite domain, the numerical solution is also shown for consistency.See Appendix A for details on the explicit finite difference scheme used, and how the flux is computed at the trap boundary.
We introduce the average absolute error Ā as a means to quantify the discrepancy between the models.Although advanced statistical tools exist to measure differences between stochastic and deterministic processes, for our purposes this simple statistical metric will suffice, and will later prove to be effective.The absolute error (relative to the total population) evaluated at discrete times t = hk is defined by, with time increment h and total time T = hK.The errors are then averaged over (K + 1) differences, Table 1.Tabulated values of the average absolute relative error Ā as defined by (32), to compare the fit between the anomalous diffusion model (30) and trap counts obtained from Brownian motion (see Figure 1 ).
Standard diffusion H = Plot (b) illustrates the discrepancy using the absolute error, which lies within 0.6 − 0.7% of all total trap counts.Theoretically, the Brownian and corresponding diffusion model are equivalent, and the errors can be dismissed as somewhat negligible, partly due to accumulation of small computational errors.We expect this error to tend to zero, with magnitude of stochastic fluctuations decreasing as N − 1 2 in the limit N → ∞, [68].Also, longer time simulations (not shown) demonstrate that the discrepancy increases as the effect of external boundary encounters is realised.Therefore, we require that the finite domain is large enough, that is L 2 D is much larger than the typical characteristic trapping time.The quantified errors in Table 1 are probably not useful on their own, however, for now they function as benchmark values which indicate an extremely good fit -indicating equivalence.As a general rule of thumb, we classify the level of fit as; equivalent 0 < Ā ≤ 0.5, good fit 0.5 < Ā ≤ 1, moderate fit 1 < Ā ≤ 1.5 and poor fit Ā > 1.5.These will be useful later as a point of reference, when comparisons are made between various diffusion models, in an attempt to reproduce Lévy trap count data.

Stable Laws
In the case of Brownian motion, the step distribution is Normal (2), and the end tails decay exponentially fast (thin tail).A large number of studies have shown that animal movement can follow a more complicated movement pattern where the step distribution decays much slower according to some type of inverse power law (heavy or fat tail), known as Lévy walks [32,44,69].As a result, individuals have a greater chance of executing 'rare' large steps, and therefore the properties of the random walk are altered.The biological consequence is such that the overall movement pattern is faster in comparison to what is typically observed in Brownian motion.Lévy walks can be characterized by Lévy α-stable distributions, simply known as stable laws.A distribution is said to be stable if the sum (or, more generally, a linear combination with positive weights) of two independent random variables has the same distribution up to a scaling factor and shift, [70,71].The mechanisms behind the resulting movement are governed by the step distribution, which is completely described using four parameters, namely, a tail index α ∈ (0, 2], skewness parameter β ∈ [−1, 1], scale parameter γ ∈ (0, ∞) and a location parameter δ ∈ R. The asymptotic behaviour of the end tails is, where α ∈ (0, 2] determines the rate at which the tails of the distribution taper off.For α ≤ 0 the distribution cannot be normalized, and therefore the pdf cannot be defined.For α ≥ 2 the end tails decay sufficiently fast at large |ξ|, ensuring that all moments exist and the central limit theorem (CLT) applies, that is the probability density of the walker after S steps converges to the Normal distribution as S → ∞.The generalized central limit theorem (gCLT) states that the sum of identically distributed random variables with distributions having inverse power law tails converges to one of the stable laws, of which the Normal distribution is a special case.For the range 0 < α < 2 the gCLT applies and the condition on the second moment is relaxed, that is; second moments diverge and the tails are asymptotically equivalent to a Pareto law.Since their first introduction, usage of stable laws has been overlooked and somewhat neglected, mainly due to difficulties arising from an infinite variance.However, there are now well developed and readily available algorithms that can be exploited for simulation runs [72,73], and thus stable laws are increasingly being considered, particularly in movement ecology [40].
Stable laws can be parametrized in Z different but equivalent ways, currently there exists at least eleven different variations which has led to much confusion [74,75].Each type has an advantage over the others and the parameter Z is often chosen based on purpose for use, i.e simulation based studies, data fitting, or study of algebraic/analytic properties.Since our focus is primarily based on obtaining trap counts from simulations, we choose the Z = 0 parametrization and henceforth use this setting.We adopt the notation introduced by [70]; that is the random variable ξ is drawn from the stable distribution S(α, β, γ, δ; Z ).
Since explicit pdfs are not available for all values of α ∈ (0, 2] the distribution is often described in terms of the characteristic function -the inverse Fourier transform of the pdf, For an isotropic random walk, the pdf of the step distribution is symmetrical and therefore both the skewness and location parameters are fixed with β = δ = 0.The resulting distribution is known as an α-stable symmetric Lévy distribution, which is then completely characterized solely by the index α and scale parameter γ.For brevity, we adopt the notation where it is understood that all parameters are zero except α and γ.The resulting characteristic function in (34) reads which is a useful way to mathematically describe all (symmetric) stable distributions, since a closed form or analytical expression does not exist for all indices α, with the exception of the Normal α = 2 and Cauchy α = 1 cases; where, the subscripts n, h, c, l have been included to distinguish between the different distributions.
For the normal distribution, σ is the standard deviation shown in the step distribution (2), with γ n defined in this way due to the choice of parametrization.This relation can be easily derived through the characteristic function (35).Note that, some authors use Lévy distribution to refer to stable laws, however, more commonly it refers to α = 1 2 , β = 1 which is a skewed distribution, defined for ξ ≥ 0 [70].Since we consider symmetric distributions i.e. β = 0, we will refer to (39) as the 'Symmetric'-Lévy distribution.In some cases, the pdf can be expressed analytically, even if it cannot be written in closed form e.g the Holtsmark distribution (37) can be written using hyper-geometric functions (if symmetric) or the Whittaker function (if skewed).Also, the Symmetric-Lévy distribution can be expressed in terms of special functions, such as Fresnel integrals [74].However, these representations are bulky and not useful in the context of this study.

Equivalence of Trap Counts: Cauchy Walk vs. Diffusion
Standard diffusion provides an oversimplified description [76], partly due to an underlying assumption that all individuals are identical in terms of their movement capabilities.In reality, this is not entirely true, and it is known that individuals in the field do not possess an equal ability for self-motion.Even in the case of a population of identical insect specie type, distinct traits can significantly affect movement abilities -such as body mass, length of wings, or more generally shape and size [77].To overcome this assumption, modified diffusion models have been successfully introduced.For example, [47] take into account that the diffusion coefficient can vary according to some type of diffusivity distribution function, rather than being constant.As a result, it is shown that by introducing the concept of a statistically structured population, the fat tails which are inherent in Lévy walks, can appear due to the fact that individuals of the same species are non-identical.Therefore, the mechanism of fat tails formation in a real population is always present, even if it can sometimes be induced by a mixture of other processes.This is not the only approach which attempts to explain the phenomena of fat tails appearing.We propose an alternative model (see later Peer-reviewed version available at Mathematics 2018, 6, 77; doi:10.3390/math6050077§4.3),where diffusivity varies continuously with time.On an individual level, the interpretation is such that distinct diffusive rates are adopted, and therefore the model takes into account individual variation.However, at the population level, when rates are aggregated, the movement is governed by time dependent diffusion.
In §3.1 it was demonstrated that equivalent trap counts can be obtained for Brownian movement using an anomalous diffusion model (30).Here, we test whether this same model can explain trap counts from non-Brownian movement, with particular interest in the level of discrepancy.Figure ( 3) Plot (a) compares trap counts J t from the Cauchy walk 4 against the diffusive flux J(t) for exponents; standard diffusion H = 1 2 , super-diffusion H = 3 4 and ballistic diffusion H = 1.A non-linear curve fitting tool is used to determine the best-fit parameters D 0 , D 1 by fitting the analytic solution (30) against trap count recordings in the least squares sense, (see Table B.1 for complete list of trap count recordings).Plot (b) illustrates the discrepancy between the trap counts and diffusion model.The relative error is shown (instead of the absolute relative error used previously in Figure 1) to distinguish between the time intervals when trap counts are either over or under estimated, here, a positive relative error corresponds to the diffusive flux forming an overestimation, and vice versa.Table (2) quantifies the fit and compares whether the diffusion model can reproduce trap counts as effectively as what was previously seen in the Brownian case (see Figure 1).4 This specific type of random walk is of significant interest in foraging theory since an inverse square power-law distribution of flight lengths provides an optimal strategy to detect target sites provided that the sites are sparse and can be revisited [60].Also, see §2.2.(32), to compare the fit between the anomalous diffusion model (30) and trap counts obtained from Brownian motion (see Figure 1) and the Cauchy walk (see Figure 3).
Standard diffusion H = From Figure 3 Plot (b) it is clear that standard diffusion fails to predict trap counts adequately, with a maximum relative error of about 5%.This is expected since standard diffusion corresponds to the random walk model with a normal step distribution whose end tails decay exponentially, unlike the Cauchy distribution.On comparison, Table 2 shows that both Super diffusion and ballistic diffusion provide a good/moderate fitting, respectively.Petrovskii et al. [26] demonstrated that in the case of some simple diffusion models, such as linear dependency on time D(t) = at + b and its variation D(t) = at + bt 1 3 (a, b constant), both of which constitute ballistic motion, provide a reasonable fit to trap count data.Figure (3) also confirms a reasonable fit for the ballistic case where the relative error lies within approx.3%, however, still an evident discrepancy is noticed.Trap counts are much better reproduced by super-diffusion, in this case H = 3  4 with relative error within approx.2%.This is somewhat expected, since it is well known that generally, super diffusion has long been acquainted with Lévy walks [41,78].Despite this promising accuracy, note that the diffusion flux tends to alternate; in the sense that trap counts are overestimated for small time, underestimated for intermediate time and overestimated again on a larger time scale -for both the super diffusive and ballistic case.This phenomena is also realised for other Hurst exponents in the super-diffusive and ballistic regime (simulations not shown here), even when a variety of scale parameters are considered.Although, the accuracy of the matching is somewhat ecologically acceptable in either case, a diffusion coefficient which yields trap counts to a higher level of precision is sought.
In a recent study, Ahmed and Petrovskii [27] demonstrated that the time-dependency of the diffusion coefficient could be inherently more complex than that proposed by anomalous diffusion.Subsequently, a time dependent diffusion model was developed to provide an alternative framework for a Cauchy walk with step distribution (38).In particular, passive trap counts were reproduced effectively and the study indicated that in the case of a Cauchy walk, the problem of trap count interpretation can be addressed with a high precision based on the diffusion equation.However, some drawbacks with this model 5 include, firstly; the complicated structure of the diffusion coefficient leads to practicality issues, since it is not expressed in closed form.Secondly; the model consists of multiple unknown parameters with little room for interpretation i.e the ecological significance of parameters or how they relate to the movement pattern is unclear.Finally; the study is constrained to Cauchy walks, and it is not known whether the diffusion model is effective in predicting trap counts for a broader range of tail indices α.With this background, we attempt to address the following: Can trap counts obtained from a system of genuine Lévy walkers be accurately reproduced using the diffusion equation, in particular, with a greater accuracy than what is observed for anomalous diffusion in Figure (3)?If so, what is the structure of the diffusion coefficient and how can the behaviour of the resulting diffusion profiles be explained from an ecological viewpoint in relation to any parameters?

Proposed Diffusion Coefficient
Observations of trap count patterns (such as those typically observed in Figure 3) suggest that the coefficient proposed should consist of some type of growth function G(t), which should behave 5 See equation (4.6) in Ahmed and Petrovskii [27] for a detailed description of the model previously proposed.as a controlling mechanism for diffusivity on a short and/or intermediate time scale.In addition to this, a suitable decay function should also be introduced to induce a dampening effect to ensure that trap counts are not overestimated for larger time, typically observed when the movement process grows faster than standard diffusion.Intuitively, we propose the following structure Exponential decay (40) with G(0) = 0, i.e growth is zero at t = 0 so that initial diffusivity is defined as D(0) = D 0 , with obvious meaning.Here, the growth function is subject to exponential decay causing the diffusivity to be damped over larger times with damping coefficient ν.Subsequently, the diffusivity returns to an initial state D 0 in the large time limit as t → ∞, provided the growth function is not faster than exponential growth, that is lim t→∞ G(t)e −νt = 0.The corresponding diffusive flux for an initial uniform population U 0 across a semi-infinite domain x > 0 with zero density condition at x = 0 (described in §3) can be derived using (20), which reads A number of possible candidates for the growth function exist in the literature, but are often applied to model population dynamics.Examples of such include; logistic, Gompertz, von Bertalanffy, generalized or hyper-logistic growth [79].The simplest of these is the logistic type, and an example of application is the Rosenzweig and MacArthur [80] model for predator-prey interactions with logistically growing prey.We propose that the growth function G(t) in ( 40) grows logistically, as a means to model diffusivity, rather than the typical use of modelling populations, so that with corresponding diffusion coefficient The movement dynamics are then completely governed by set parameters {D 0 , D 1 , k, ν}, all positive.This growth function is a parabolic function of time, which increases from G = 0 at time t = 0, until a maximum G max = 1 4 D 1 k is attained at time t = k 2 , and then decreases until the growth diminishes, i.e G = 0 at time t = k.The corresponding diffusion coefficient is still valid for negative growth over the interval k < t < t * provided D(t) > 0, where t * is a solution of ln D 1 (t * − k) − νt * = ln D 0 k.For fixed D 1 , the value of k controls the maximum growth capacity and also determines the instant in time when growth alternates from positive to negative.In the special case, for sufficiently small time with large k, the term t k in ( 43) is negligible and the growth function is approximately linear G(t) ≈ D 1 t, and in the limiting case as k → ∞, with corresponding diffusion coefficient, which now depends on three parameters {D 0 , D 1 , ν}.The corresponding diffusive flux can be derived for the logistic model ( 42) from (41), which results in Model 1: and simplifies to Model 2: in the reduced linear case (44).Henceforth, we will refer to this as Model 1 & 2, respectively.
We presume that this increase in diffusion rate can induce faster movement, which can be comparable (at some level) to the pattern inherent in Lévy walks.Following this, the diffusivity begins to decrease until it reaches a minimum level D min at time t = k 2 + 1 ν + k 2 1 + 4 (kν) 2 , and then asymptotically approaches the initial state D 0 in the large time limit.In the special case of linear growth function, insect diffusivity reaches D max = D 0 + D 1 eν at time t = 1 ν with no subsequent local minimum, and same asymptotic behaviour as t → ∞. Figure 4 Plot (c) illustrates the flux for each corresponding diffusion coefficient in Plot (b), it is precisely these diffusion Models 1 & 2 (46 -47) which will be tested in §4.4 against Lévy trap count data.

Reproducing Lévy Trap counts using Diffusion
In this section we test Models 1 (46) & 2 (47) against Lévy trap count data (see Figure 5  laws, defined in (37 -39).Also, the assumptions that the walk is uncorrelated and unbiased in a homogeneous environment still apply.In a system of N individuals executing a Lévy walk, the position of the n th individual at the (i + 1) th step can be described by The trap count is expected to grow faster with time, compared to what is usually recorded for Brownian movement.This is due to the frequency of long jumps increasing, and therefore the contribution from remote parts of the population to the trap count also increases.For our purposes, we simulate trap counts for the tail indices α = 3 2 , 1, 1 2 , referring to the Holtsmark (37), Cauchy (38) and Symmetric-Lévy (39) distributions, previously introduced in §4.1.The movement dynamics are completely governed by the scale parameters γ h , γ c , γ l .Although, comparing random walks prior to simulation runs can reveal information on parameter selection [81] (also see §2.2), for our purposes it suffices to arbitrarily select three distinct scale parameters for each case i.e. γ h = 0.01, 0.02, 0.04, γ c = 0.0005, 0.002, 0.003 and γ l = 1 × 10 −6 , 4 × 10 −6 , 2 × 10 −5 .Here, parameters are chosen so that trap count data is obtained with reasonable level of variation (see Table B.1).The diffusive flux curves J(t) given by Model 1 ( 46) & Model 2 (47) are then fitted (in the least squares sense) against these trap counts using a non-linear curve fitting tool and best-fit parameters are estimated, listed in Table 3. (n) i < 0 at any instant in time then the individual is removed from the system and the accumulated trap count increases by one.Consequently, the number of individuals in the population decrease as time flows, and an increasing stochastic trap count trajectory is formed.Trap counts: Bold dots depict cumulative trap counts J t recorded for the cases (a) Holtsmark, (b) Cauchy, (c) Symmetric-Lévy at times t = 0, 0.1, 0.2, ..., 3..Different scale parameters are considered for each respective case.Also, trap counts are averaged over five realizations to reduce the effect of stochasticity.(For full list of recordings, see Table B.1). Diffusive flux: Curves J(t) shown for Model 2 in all three cases (red, blue, black curves).Model 1 shown only for the case corresponding to Holtsmark S(α = 1.5, γ h = 0.04), (magenta curve).All best fit parameters are listed in Table 3.  Figure 5 illustrates the fitting between the diffusive flux J(t) and Lévy trap count data J t , for the Holtsmark α = 3 2 , Cauchy α = 1 and symmetric-Lévy α = 1 2 distributions, respectively.Model 2 is shown and found to form an almost identical fit, with the exception of the Holtsmark case with γ h = 0.04.In this special case, we find that Model 2 eventually overestimates trap counts, which is more apparent for larger γ h .Here, we find that Model 1 forms a better fit and this can also be realised on inspection of best fit parameters in Table 3.The parameter k is significant (compare boxed value with other k), and therefore the term t k behaves as some type of correction term which slows the diffusivity rate.The corresponding growth function is of logistic type, with diffusion coefficient In all other cases, k is relatively large, and therefore on a short time scale the term t k in this diffusion coefficient is negligible.As a result, the growth function is then approximately linear and Model 1 reduces to Model 2 with diffusion coefficient D(t) = D 0 + D 1 te −νt .Figure 6 shows a plot of the diffusion coefficients, corresponding to each case in Figure 5.The diffusive profiles tend to follow a particular pattern.Typically, the diffusivity begins at some initial value D 0 and begins to increase until it peaks at D max , and then subsequently decays to re-approach the initial value for larger time (where the latter is not typically essential as the interest is on short time dynamics).In the special case, seen in Plot (a), the rate of decay for Model 2 (black curve) is slower than required, resulting in larger diffusivity and explains the overestimation previously observed in Figure 5 Plot (a) for the Holtsmark case with γ h = 0.04.On comparing the diffusion coefficients for both models, we find that trap counts are better estimated using a profile with larger initial diffusivity, smaller peak and faster decay (see dashed curve).Table 4. Tabulated values of the average absolute relative error Ā as defined by (32), to compare the fit between models 1 (46) and 2 (47) and trap counts.Ā is also included for the anomalous diffusion model, see Figure 3 and Table 2. Figure 7 illustrates that a high level of accuracy is maintained, where the error lies roughly within 1% using; (i) Model 1 for the Holtsmark case with γ h = 0.04 (magenta circles in Plot (a)), and (ii) Model 2 for all other cases.On comparing the absolute relative error in Table 4, we find that the proposed diffusion models significantly improve trap count prediction, than what is obtained from super diffusion.Moreover, the numerical values in Table 4 are indicative of equivalence (compare boxed values to others), since values of the same order of magnitude (0 < Ā ≤ 0.5) were obtained when comparing standard diffusion to Brownian motion -which are theoretically equivalent movement processes (see Table 1).Evidently, these proposed diffusion models can be used effectively to reproduce trap counts for a system of Lévy walking individuals to a remarkable level of accuracy -yielding almost identical counts.

Discussion
The biological basis behind time dependent diffusion, and how this is linked to the type of mechanisms involved has not been clearly understood [51].Here, the parameters introduced have ecological significance with obvious meaning, in the sense that they are not all arbitrary and are related to the underlying movement dynamics.One possible explanation of the type of diffusive pattern seen in Figure 6 can be arrived at through the concept of differential energetics [82][83][84].If an individual executes a large step, then this will incur higher energy costs than a small step, so energy expenditure for different step lengths is a given.In relation to basic individual traits, such as body mass, a heavier body requires a larger force and a larger energy expense to change direction or execute a larger step.Therefore, it may be expected that the frequency of moving long distances is lower for heavier individuals.Also, for those individuals that do prefer to take large steps, this could possibly induce larger rest pauses contributing towards intermittent behaviour [85], from which decrease in the rate of diffusivity can be explained.On the other hand, the mechanisms behind the movement behaviour could involve a degree of spatial synchronisation of individual variation.In the case of movement with a physiological origin, this could result from a diffusivity distribution as indicated by Petrovskii and Morozov [47] or alternatively, even time dependent diffusion.Also, with reference to behavioural effects, a sudden event may trigger change in behaviour in one individual that results in swarming behaviour, leading to a change at the population level [86].It must also be realised that the build up of insect population and following changes in diffusivity can also be a result of transient environmental factors, e.g temperature [87], as sudden temperature changes can excite movement or even lead to erratic behaviour.If one or more of these factors can bring about diffusive movement varying with time, then as a consequence, the pattern produced by an insect population performing Brownian motion may be indistinguishable from the pattern produced by insects performing Lévy walks.Altogether, it seems that identifying a genuine Lévy walk may be more challenging than previously thought.In light of this, it is not surprising that the Lévy or diffusion controversy has been persistent, with strong evidences arguing for either side -in this study, we have demonstrated that both sides are somewhat equivalent, at least in the context of trapping.
On a final note, we would like to mention some limitations and suggest possible further research directions.Although, this study is purely theoretical and attempts to answer some important issues in movement ecology, it is limited to the 1D spatial scale.In a more practical scenario, it would be interesting to see a similar study in a more realistic 2D domain, which would be more relevant to field studies with trapping of walking/crawling insects [12,28].The problem would then have enhanced complexity due to the introduction of a single trap of different possible shapes and sizes.The movement pattern would then be altered, as the effects of trap and field boundaries are realised by enforcing a confined or restricted space [59].In particular, interest would lie in how the time dependent diffusion model could be further refined in order to produce trap counts to a high level of accuracy within these different geometries.The diffusion model can then be checked and tested for good agreement against insect trap counts in the real field, given that there is evidence of Lévy type movement beforehand.Developing a corresponding model in 2D would be the next obvious step to take.
Attempting to model more realistic scenarios would further increase the level of model complexity.For example, in the real field, very rarely a single trap is installed, rather, a multiple trapping system is implemented [88].In terms of modelling, the geometry of the domain becomes intrinsically complicated, and the effects of perturbations from each trap may become difficult to unravel.In addition, removing the assumption that the random walk is unbiased would result in directed movement.This can be related to baited traps, which are widely used in practice to increase the frequency of captures [89,90].Application of a certain agent can induce behavioural responses and attract insects to the trap, such as light, odours or pheromone.Since installation of baited traps considerably alter the insect behavior, they are much more difficult to model [91]  corresponding theory is largely absent.In this case, the convenient mathematical framework would consist of advection-diffusion equations, where possibly, a time dependent diffusivity may be used.

Concluding Remarks
In this paper, we show that the diffusion coefficients (43) and the reduced version (45) can be incorporated into a time-dependent diffusion model, which can then be used to predict, almost identically, trap counts from a system of individuals who undergo a Lévy walk.Moreover, we show that this can be achieved for a broad range of Lévy tail indices.Also, we find that these proposed diffusion models are much more accurate and effective when compared to super diffusion.Alongside the development of the models, we explore the biological basis for time-dependent diffusion in more detail, and interpret parameters in relation to diffusive patterns.We argue that, if these inherently different movement models yield almost identical trap counts, then how important is the movement pattern in the context of trapping?This study suggests that the movement pattern is not that important after all, rather emphasis should be on the ecological context, at least in integrated pest monitoring programmes.

Figure 4
Figure 4 Plot (a) illustrates the logistic growth function as a parabolic profile for different values of k, with linear growth as k → ∞. Figure 4 Plot (b) shows how the diffusion coefficient behaves for particular parameter values.In ecological terms, the mechanistic process is such that insect diffusivity increases from the initial value D 0 until maximum diffusivity D max is attained at time t = k 2 + 1 ν −

Figure 5 .
Figure 5. Simulation details: In accordance with the simulation setting in §2.1, N = 1000 individuals are initially uniformly distributed along a 1D spatial domain 0 < x < L = 5.After one time step ∆t = 0.001, each individual executes a single step, with subsequent position defined by the recurrence relation (48).A total number of s = 3000 steps are executed, with total time of exposure T = st = 3. Prior to the simulation run, an impermeable external boundary is installed at x = 5 ensuring that no individual can escape or enter the system at this end (no-migration/immigration), by forcing the condition; if X (n) i > 5 at any instant in time then X (n) i = 5.The point trap at x = 0 functions

Figure 6 .
Figure 6.Diffusion coefficients given by (45) for Model 2 (red, blue black curves) with listed best-fit parameters in Table3.Diffusion coefficient given by(43) for Model 1 shown only in Plot (a) (dashed curve).

Figure 7 .
Figure 7. Absolute relative error between trap counts and diffusive flux for the cases (a) Holtsmark, (b) Cauchy and (c) Symmetric-Lévy.Each colour corresponds to those cases with scale parameters shown in Figure 5.

Table 2 .
Tabulated values of the average absolute relative error Ā as defined by

Table 3 .
(46)-fit parameters using a non-linear curve fitting tool (in the least squares sense) by fitting Model 1(46)and Model 2 (47) against cumulative trap counts (see TableB.1 for complete list of recordings).The diffusion coefficients in Figure6are those plotted with highlighted parameters in the Table below.