Constructing a Stochastic Model of Bumblebee Flights from Experimental Data

The movement of organisms is subject to a multitude of influences of widely varying character: from the bio-mechanics of the individual, over the interaction with the complex environment many animals live in, to evolutionary pressure and energy constraints. As the number of factors is large, it is very hard to build comprehensive movement models. Even when movement patterns in simple environments are analysed, the organisms can display very complex behaviours. While for largely undirected motion or long observation times the dynamics can sometimes be described by isotropic random walks, usually the directional persistence due to a preference to move forward has to be accounted for, e.g., by a correlated random walk. In this paper we generalise these descriptions to a model in terms of stochastic differential equations of Langevin type, which we use to analyse experimental search flight data of foraging bumblebees. Using parameter estimates we discuss the differences and similarities to correlated random walks. From simulations we generate artificial bumblebee trajectories which we use as a validation by comparing the generated ones to the experimental data.


Foraging Animals
The characteristics of the movement of animals play a key role in a variety of ecologically relevant processes, from foraging and group behaviour of animals [1] to dispersal [2,3] and territoriality [4]. Studying the behaviour of animals, simple random walk models have been proven effective in describing irregular paths [5]. While the first studies on random paths of organisms focused on uncorrelated step sequences [6], in many cases of studies of animal behaviour the directional persistence of the animals suggested a modelling in terms of correlated random walks (CRWs) [7,8]. In many complex environments an intermittent behaviour of animals is observed. In these cases an animal switches either randomly or in reaction to its environment between different movement patterns. The mechanisms which generate, and the factors which influence this switching behaviour have been shown to be important in understanding and modelling complicated animal paths [9][10][11][12]. While there is a source of switching between free flight and food inspections in the experiment we analyse [13], here we concentrate on the former as detailed below. With no clear indication of additional intermittency, we will focus on non-intermittent models in the following.

CRW/Reorientation Model
The planar horizontal movement of an animal is often approximated by a sequence of steps: an angle a(t) describes the current direction of movement in a fixed coordinate frame, while the step length l(t) determines the distance travelled during a time step. The direction a(t), often determined by a specific front direction of the animal, changes each time step by a random turning angle b(t). The description of the dynamics in a co-moving frame, i.e., via the turning angle, turned out to be most useful for analysis of persistent animal movement [7,8]. In many cases b(t) is drawn independently and identically distributed (i.i.d.) from a wrapped normal distribution or a von Mises distribution [14,15] for each time step, giving rise to a persistence in direction depending on how strongly the distribution is concentrated around 0. Usually the step length is taken to be either constant or it is drawn i.i.d. from some distribution. The step length can either be the result of a constant speed and a variable time step or (as in our case below) of a constant time step Dt and a variable speed s(t).
This class of models can generate a variety of different dynamics. Two special cases with a uniform distribution for b and a fixed time step Dt are the standard Gaussian random walk for step lengths l(t)~jz(t)j where z is normally distributed and Lévy flights for power-law tails in the step lengths distributions (l(t)*l {m for 1vmƒ3 and lwl 0 ). Related to Lévy Flights, but using a time step proportional to the step length, are Lévy Walks, which have been of interest as candidates for optimal search behaviour of foraging animals. They have been studied analytically [16], by simulations [10,17], and many experimental data sets have been statistically analysed to determine whether Lévy Walks are suitable to describe the movement of animals (see, e.g., [18][19][20][21][22]).
As Lévy-type models show anomalous diffusive behaviour, in contrast to models with a finite variance of the step length distribution and a fixed time step Dt, only the latter are included in the definition of correlated random walks which are also called reorientation models in the context of animal movement. Apart from pathological cases, CRWs are diffusive in the long time limit according to the central limit theorem.
The estimation of the tortuosity of a trajectory is intimately connected to the distributions of the turning angle and speed [8,14,23]. The relevance of the turning angle distribution for foraging efficiencies when searching in random environments has been analysed, e.g., in [24].

Generalisation of the Model
In the following we will present a generalisation of the CRW above, which we then use to analyse bumblebee flight data. Given movement data with a constant time step Dt, the step length is determined by the speed s(t)~jv(t)j of the animal. As we will be looking at a flying insect in a data recording using a small time step, we may expect to have a deterministic persistence due to the animals momentum. Additionally, the above CRW model assumes that s and b are drawn i.i.d. which is sensible if Dt is large enough. However, for small time steps it cannot be excluded that the decision of the animal to turn left or right takes longer than the time step, which can correlate the turning angles b(t) over a number of time steps. To allow for these possibilities we therefore model the changes in speed and turning angle via two coupled generalized Langevin equations, where we distinguish between the deterministic parts h and g and stochastic terms y andj j s (whose speed dependency will be discussed in the Results section). We assume that the noise processes are stationary with auto-correlation functions which may be non-trivial, and we make no further assumptions for the shape of their stationary distributions. While Eqs. (1,2) represent a timecontinuous description, the turning angle b still yields the change of a according to our fixed time resolution Dt. That is, b(t) relates to a time-continuous angular velocity c of a via b(t)~ð t t{Dt c(t)dt. The animals' position r(t)~(x(t),y(t)) is then given by dx=dt~s cos (a(t)), dy=dt~s sin (a(t)) and da=dt~c(t). Not having experimental access to c, the numerical analysis is done with time-discrete data where the measured turning angle is given by where v(t)~(r(tzDt){r(t))=Dt at times t~nDt, n[N.

Application to Experimental Data
Analysing measured movement data of animals in their natural habitat is intricate due to a variety of factors which may influence the animal's behaviour, ranging from heterogeneous food source distributions [25][26][27] and predation threats [13,28] to individual differences in behaviour within a population [2,3]. Here we analyse data obtained from a small scale laboratory experiment in which single bumblebees forage in an artificial flight arena [29]. The set-up is shown in Fig. 1 together with part of a typical trajectory of a bumblebee on its search for food. Each bumblebee can forage on an artificial flower carpet which is positioned on one of the walls of the arena. In this paper we are not interested in the behaviour resulting from the interaction with the flowers which has been studied in detail in [13]. Instead we only examine the search flights away from the flower carpet. (See section Materials and Methods for details.) We use our generalised stochastic model (Eqs. (1,2)) to describe these flights and to examine in which ways the behaviour deviates from a simple CRW model. Here we will focus on the horizontal movements. By neglecting the slower vertical movements, which are of more interest when analysing the starting and landing behaviour near flowers, we thus restrict ourselves to a two-dimensional model.

Estimation of Drift Terms
Given the experimental data, we start determining the unknown parameters in our model by first estimating the deterministic parts h(b,s) and g(b,s) of the Langevin equation. This is done by numerical estimation [30][31][32][33] of the components of the drift vector field (drift coefficients) where and v:w is the time average over the time series * X conditioned on * X (t)&X , where * X is assumed to be stationary (for a detailed discussion see [32]). The estimation of the drift terms is based on a Markov approximation: only those parts of the  Figure 2 shows the drift vector field, with normalised lengths of the vectors for better visibility. The nearly horizontal vectors show, that the drift quickly pushes the turning angle b towards 0, while the dynamics in the speed s is much slower. As the cross-dependencies of h(b,s) on s and of g(b,s) on b are weak, we can neglect them in our model. Since vector fields are hard to interpret, we will look at the projections in the following.
Examining the drift h(b) of the turning angle in Fig. 3 reveals that the drift term seems linear in b -indeed we find numerically that its slope {k matches exactly to a decay of the turning angle to 0 in a single observation time step Dt by k&1=Dt, disregarding the noise term. This means that by integrating Eq. (1) over a time Dt and approximating the drift h(b) for small Dt by With j s (t) :~Ð While this reduction of the turning angle dynamics from db=dt to b bears similarity to a simple reorientation model, the turning angles are still correlated and speed-dependent, as we will see below.
The speed drift g(s) displayed in Fig. 4 shows that the deterministic part of the speed-Langevin equation alone would have a stable fixed point around s 0~0 :27m=s. Comparing the slopes above and below s 0 reveals that for svs 0 the force towards s 0 is stronger than for sws 0 . This is biologically plausible if one interprets s 0 as a preferred speed: if the bumblebee is slower it   accelerates, but if it is faster it does not rush to decelerate as it would give up the energy spent to reach a high velocity. For very high velocities (over 0.55 m/s) the slope of g(s) increases again. This might be caused by the limited space available to the bumblebee in the flight arena. For our model we approximated g(s) by a piecewise linear function: where d 1 wd 2 w0. As the very high velocities are rare, it made no difference in our model whether we used Eq. (7) or a piecewise linear function with three pieces.

Velocity-Dependent Angle-Noise and Noise Auto-Correlations
What we did not specify before was that the turning angle distribution may depend on the speed of the bumblebees. Given that the force a bumblebee can use to change directions is finite, the largest turning angles have to be smaller when flying with high speeds (see Fig. 5). This is consistent with the absence of simultaneously having high speed and large turning angle in the data -as is evident, e.g., from the data gaps in Fig. 2. However, animals can counteract this geometric dependence by varying the forces used for changing direction with the speed. We approximated the distribution for the turning angles for each speed s by a normal distribution. This approximation works best for low speeds. While there are some deviations for high speeds, it was not possible to reliably fit a better model due to the limited amount of data available. Figure 6 shows how its standard deviation s b depends on the current speed. s b decreases with increasing speed, however it does not decay to 0 as a simple geometric model would predict (see Materials and Methods below).
Instead s b (s) decays roughly exponentially to a constant offset. We therefore model the turning angles as speed-dependent noise with a wrapped normal distribution [14,15]: j s (t)*N (0,s j (s)) with s j (s)~c 1 e {c2s zc 3 . This offset could either be an effect of the boundedness of the flight arena, since the bumblebee has to turn more often to avoid walls when flying fast. Or it could be that the bumblebees use stronger forces for turning during fast flights to maintain their manoeuvrability. It would be interesting to examine free-flight data to check for the cause. In other models in which the momentum of the animal is not important for the observed directional persistence, this cross-dependence is often neglected [7].
For the two stochastic parts of the Langevin equations, we estimated the normalised auto-correlation functions from the data. The turning angle auto-correlation is approximated by a steep power-law as seen in Fig. 7, which in this case is preferable to the alternative fit by a simple exponential decay. By subtraction of our approximation for the deterministic term g(s) from the observed speed changes ds=dt in Eq. (6) we estimated the distribution and auto-correlation of the noise term y(t)~ds(t)=dt{g(s(t)). In   order not to overestimate the noise term, additive discretization errors of an approximate size of s error~D x=Dt 2 due to the finite resolution Dx~10 {3 m of the cameras have been accounted for, giving the variance s 2 y~s 2 y noisy {s 2 error . The noise term y(t) is well approximated by Gaussian noise with an auto-correlation function acf e{e y (t)~ae {l1t z(1{a)e {l2t (see Fig. 8). While an autocorrelation function of the shape of acf p{p y (t)~b(tz1) {p1 z(1{b)(tz1) {p2 can be exluded, a difference between an exponential and a power-law acf e{p y (t)~ce {l3t z(1{c)(tz1) {p2 is not significantly worse than acf e{e y . For our model we chose the simple difference of exponentials acf e{e .
As the observed anti-correlation between delays of 0:1 swtw0:3 s happens on a time scale which is too short to be an effect of the boundedness of the experiment or of residual effects of the presence of the foraging wall [13], it is unclear where the anti-correlation comes from. One could speculate that it might be the result of a stabilising mechanism in the bumblebee dynamics.

Validation
Given all the parameters of the full model (see Materials and Methods) estimated by minimizing the mean squared errors, we used them to generate artificial bumblebee trajectories, as follows: We simulated the dynamics using an Euler-Maruyama scheme with noise terms j s (t),y(t). In rare cases where the Gaussian noise y(t) would lead to a negative speed despite the positive drift g(s) for svs 0 , we enforce a non-negative speed by setting s(t)~0. We correlated the noise terms in advance by modifying their power spectral density in the following way: we take uncorrelated noise of the wanted distribution, multiply its Fourier transform with the Figure 8. Auto-correlation of the non-deterministic speed changes y(t). The auto-correlation function of y(t)~ds(t)=dt{g(s(t)) estimated from the experimental data (dots) with two times the largelag standard error (grey) and three fitted approximations: difference of 2 exponentials (red), difference of 2 power-laws (green), difference of exponential and power-law (blue). The outlier at t~0:02 s is a discretization artifact due to the finite resolution of the data (see [40]). doi:10.1371/journal.pone.0059036.g008   root of the desired power spectral density corresponding to our approximate auto-correlation function and then transform back [34]. To deal with the speed dependence of the turning angle noise j s (t) we first correlate Gaussian noise and afterwards scale with s b (s) at each time step in the integration scheme. While this does not reproduce the auto-correlation of the turning angle exactly, the error made is less than the errors from the estimation of acf b . A sample trajectory of a bumblebee simulated for 200 s using 10 5 time steps is shown in Fig. 9. Using the generated data we checked the validity of the model by comparison to the experimental data of all bumblebees. Figure 10 compares the probability density function pdf(s) of the speed extracted from the simulated data with the experimental data. The auto-correlation functions of the speed and turning angle are shown in Figures 11 and 7. Considering the number of rough approximations we have made for constructing our model, the agreement between simulation results and experimental data is very good.

Summary
We generalised a reorientation model which is often used to describe the correlated random walk of animals by explicitly modelling accelerations via Langevin equations. Analysing movement data from bumblebees, we extracted information on the deterministic and stochastic terms of Eqs. (1,2). Simulations of our model and comparison to the data have shown that the resulting model agrees very well with the experimental data despite the approximations we made for the model. With the estimation of the turning angle drift h(b) we found that while the usual assumption of i.i.d. turning angles is not valid in our case, the lack of a nontrivial drift and the weak auto-correlation of j s are consistent with the usual reorientation model. However, our generalised model exhibits significant differences in the non-trivial deterministic part g(s) of the speed change ds=dt and the speed dependence of the turning angles. In terms of active Brownian particle models [35,36] we described the two-dimensional bumblebee movement by a particle with a non-linear friction term g(s) depending and acting only on the speed, driven by multiplicative coloured noise with different correlations for the angle component and the speed component of the velocity. While this combination of complications might make it difficult to treat the system analytically, progress into this direction has been made [37,38]. We remark that one could ignore the fast decaying auto-correlations of j s and y(t) if one is not interested in the dynamics for short times, thus simplifying the model by using uncorrelated noise terms, since the effect of the noise autocorrelations on the long time dynamics is negligible.
Given that the experiment which yielded our data is rather small and provided the bumblebees with an artificial environment, it would be interesting to apply our new model to free-flying bumblebees to reveal how much the results depend on the specific set-up. This would clarify whether the flight behaviour seen in the laboratory experiment survives as a flight mode for foraging in a patch of flowers in an intermittent model, with an additional flight mode for long flights between flower patches. The analysis of data from other flying insects and birds by using our model could be interesting in order to examine whether the piecewise linear nature of the speed drift and the trivial drift of the turning angle are a common feature. In view of understanding the small-scale biomechanical origin of flight dynamics, our model might serve as a reference point for any more detailed dynamical modelling. That is, we would expect that any more microscopic model should reproduce our dynamics after a suitable coarse graining over relevant degrees of freedom.

Experimental Data
In this experiment 30 bumblebees ( Bombus terrestris) were trained to forage individually in a roughly cubical flight arena with an approximate side length of 0.75 m. Figure 1 shows a diagram of the arena together with data from a typical flight path of a bumblebee. The flight arena included a 4|4 grid of artificial flowers on one of the walls. Each of the 16 flowers consisted of a landing platform, a yellow square floral marker and a replenishing food source where syrup was offered. For the analysis presented in this paper all data in zones (7 cm|9 cm|9 cm) around the flowers has been removed in order to analyse the search behaviour while foraging excluding the interaction with food sources. The 3D flight trajectories of the bumblebees were tracked by two cameras with a temporal resolution of Dt~0:02 s. Each bumblebee was approximated as a point mass with a spatial resolution of 0:1 cm: its position was estimated by the arithmetic mean of all image pixels corresponding to the bumblebee via background subtraction. In total &49000 data points were used for the analysis. For individual bumblebees an average of 51 search trajectories between flower zones have been sampled and analysed. The thorax widths of the bumblebees have a mean of 5:6 mm and a standard deviation of 0:4 mm.
For calculating auto-correlations small gaps in the time series have been interpolated linearly. As the number of gaps was small the correlations for short times were not affected, however, the interpolation increased the usable data for long time delays. Trajectories were split at larger gaps, e.g., when entering a flower zone, to exclude correlations induced by flower visits.
For a discussion of the influence of the boundedness of the flight arena and for the analysis of the foraging dynamics under varying environmental conditions see [13]. More details on the experimental setup can be found in [29,39].

Estimated Model Parameters
The full set of parameters estimated from the data set which was used for the simulation is given here. For the deterministic drift of the speed the change of slope is at s 0~0 :275 m=s while the slopes are d 1~0 :16 and d 2~0 :06. The parameters for the standard deviation s j (s) of the angle noise are c 1~1 26 0 , c 2~1 2 s=m, c 3~1 2:5 0 and its auto-correlation is given by acf b (t)~(tz1) {1:5476 . The non-deterministic changes y(t) of the speed are assumed to be normally distributed with standard deviation s y~3 :52 m=s 2 and auto-correlated according to acf e{e y (t) where a~1:44, l 1~2 5:5 and l 2~1 0:7.

Speed Dependence of Turning Angles
A simple model showing a dependence of the turning angles on the speed (see Fig. 5) is given in the following. Assume that the velocity of an animal changes at each time step Dt by an acceleration vector which is given by a binormal i.i.d. random vector with variance s 2 in both directions. The turning angle b between v t and v tzDt then depends on the quotient g t :~s t =( ffiffi ffi 2 p s) between the former speed s t~j v t j and the noise strength s. By changing to the comoving frame of the animal and integrating out s tzDt the distribution r(b) of the turning angle is given by: z e {g 2 sin 2 (b) 2 ffiffiffi p p g cos (b)(1zerf(g cos (b))) for {pƒbƒp. With vanishing speed s(t)~g(t)~0 the first term gives a uniform distribution as expected, and for g(t)?? the distribution sharply peaks at b~0 with its variance s b approaching 0, similar to the behaviour of the simpler von Mises distribution. As the experimental bumblebee data does not show a decay to s b~0 but to a finite value (see Fig. 6), this simple model does not hold: therefore the accelerations have to be modelled as speed-dependent.