- Email: [email protected]
Simulation modelling of nectar and pollen foraging by honeybees
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
Available online at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/issn/15375110
Research Paper
Simulation modelling of nectar and pollen foraging by honeybees Juan Jose Garcia Adeva a,b,* a b
School of Computer Science and Software Engineering, University of Western Australia, Perth, Australia Cooperative Research Centre for National Plant Biosecurity, Canberra, Australia
article info
This paper describes a simulation model for the foraging activity of honeybees (Apis mel-
Article history:
lifera) as they collect nectar and pollen. Understanding such behaviour and observing the
Received 14 December 2011
resulting foraging spread patterns and distribution are vital from an ecological perspective
Received in revised form
and has useful applications in areas such as agriculture or biosecurity. The simulator is
20 March 2012
based on a spatio-temporal model and implemented following a Web-based architecture.
Accepted 4 May 2012
This simulation model takes into account of detailed floral distributions, weather, hours of
Published online 13 June 2012
daylight, and variations in daily nectar and pollen production by flowers. The produced simulation technology includes a user interface that allows the introduction of the simulation parameters (e.g. spatial resolution and extent, model configuration settings, visual specification of the simulation environment on the map) and receive a detailed results report. The simulation model is evaluated by performing multiple experiments and verifying the results obtained. ª 2012 IAgrE. Published by Elsevier Ltd. All rights reserved.
1.
Introduction
Foraging is the act of honeybees (Apis mellifera) fetching food from sources scattered around their beehive. In this work, we consider foraging to include pollen and nectar from flowers in order to feed the colony, excluding other types of provisions such as water. Nectar provides energy while pollen is the main source of proteins. Foraged pollen and nectar is stored in the beehive, in the form of honey, in order to have enough food for times when flowers do not produce enough nectar or pollen, such as winter (Seeley, 1995). Therefore, being efficient in this foraging activity is key to the colony’s continued survival. Floral patches differ importantly in several factors that affect the profitability of the foraging effort on the patch (Winston, 1991), including the quality and quantity of their
nectar and pollen, the size of the patch, and the distance from the colony. This overall profitability of a floral patch is communicated by individual foragers to other bees in the colony via complex behaviours and other signals (von Frisch, 1967). These other individuals may alter their behaviour accordingly, and their responses may in turn affect the behaviour of the original forager (Seeley, 1985). Thus foraging in honeybees is a dynamic and highly interactive process that is driven by the sharing of information about the location and quality of resources in the surrounding environment. The intensity of this foraging activity and especially its extent over the local flora can have an influence on aspects such as crop pollination and plant biosecurity. There are several reasons for this. On the one hand, the foraging behaviours and hairiness of honeybees mean that they often
* School of Computer Science and Software Engineering, University of Western Australia, Perth, Australia. Tel.: þ61 8 6488 3816; fax: þ61 8 6488 1089. E-mail address: [email protected]. 1537-5110/$ e see front matter ª 2012 IAgrE. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.biosystemseng.2012.05.002
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
inadvertently transfer substances, such as pollen, among plants, hence potentially having important consequences to crop pollination. Honeybees are capable of increasing the yield in 96% of animal-pollinated crops (Potts et al., 2010), and the value of the increased yield and quality achieved by honeybee pollination of crops in the U.S. is estimated as $14.6 billion (Morse & Calderone, 2000). On the other hand, because honeybees can also act as vectors of some plant diseases such as Fire Blight (Erwinia amylovora) and also as vectors of biological control agents (e.g. the microbe Aureobasidium pullulans, which can be used to attack E. amylovora (Moosbeckhofer et al., 2008)), understanding patterns of foraging behaviour will aid in implementing prevention or facilitation measures, respectively, thus contributing to important aspects of disease prevention or eradication in the context of plant biosecurity. This last scenario describing the influence of honeybees in spread of plant diseases is the main motivation for this work. Plant biosecurity is a major concern in a country like Australia where economic viability of plant industries rely to an extent on freedom from pests and diseases. Strategies to test preemptive surveillance and estimate rates of spread of pests or diseases in a spatial environment in real time are not currently available. A ‘what if’ scenario simulation tool is needed to address this shortcoming (Garcı´a Adeva et al., 2010). The aim of the main project, where this current honeybee simulation work resides, is to provide such decision support tools. In consequence, simulating honeybee foraging behaviour will help us to determine the extent and rate at which surrounding flora may be affected by a disease such as Fire Blight. This work consists of a simulation model that is stochastic and discrete in time with 1-trip time steps, which is approximated by difference equations. It is also explanatory (Grant & Swannack, 2008), meaning that its behavioural rules depend on well-established expert knowledge or published empirical data. It can support any number of weather periods, beehives and floral patches over a two dimensional landscape of unrestricted size and varying resolution. Each flower patch can have any shape, size, location, and number of seasons that define the daily nectar and pollen production, therefore extending some previous preliminary work limited to only nectar foraging (Garcı´a Adeva et al., 2011). A weather period is defined by average daily temperature and number of daylight hours. Because the produced technology shares some modelling aspects and resulting features with previous work by the same authors (Garcı´a Adeva, Botha, & Reynolds, 2012), this will be noted where appropriate. From its implementation point of view, the simulator usage is Web based, so that users do not have to install or update software, or worry about enough computational power or platform compatibility. The user interface allows the selection of the simulation parameters (e.g. spatial resolution and size, model configuration settings, visual specification of the simulation environment on the map) and receive a highly detailed report of results. A simulation of several weeks usually takes a few seconds to complete. The structure of this paper is as follows: in Section 2, we discuss some related work in honeybee modelling or Web-
305
based simulation. Section 3 describes the main simulation model and related sub-models. Section 4 offers an evaluation. We finish with some concluding remarks in Section 5.
2.
Related work
There is a number of research efforts in simulating or modelling different aspects of honeybees. This section analyses the most relevant instances of previous honeybee simulation work that we have found. de Vries and Biesmeijer (1998) implemented an agentbased model intended to replicate the foraging experiment that Seeley, Camazine, and Sneyd (1991) had previously performed by using behavioural rules of individual bees. Also based on this same experiment, Camazine and Sneyd (1991) proposed a foraging model based on differential equations for foraging, which considers nectar quality in the decision making process for resource abandonment and recruitment of foraging bees. Cox and Myerscough (2003) proposed a mathematical model of foraging in honeybee colonies. The effort was mainly focused on examining the effect of a small set of foragingrelated parameters (e.g. visit rate, probability of dancing, duration of dancing, probability of abandonment) on the foraging response and consequent nectar intake. This model was intended as an experiment to support the biological hypothesis of ‘genetic variance’ that argues that genetically heterogeneous colonies are better protected against misfortune than homogeneous colonies (Crozier & Page, 1985; Keller & Reeve, 1994). Ghosh and Marshall (2005) proposed a model of learning and collective decision making in honeybees engaged in nectar foraging. However, this model is limited to one of two foraging sites where one site has better nectar than the other in order to understand decision making processes by interacting swarm entities. Bailis, Nagpal, and Werfel (2010) implemented a simulator that models individual foraging bees at a very high level of detail in order to compare the value of sharing food source position information through the waggle dance and private information about previously encountered food sources. The most recent effort consists of a model of honeybee foraging behaviour to estimate almond crop yields as functions of honeybee densities in order to understand the economic behaviour of growers who rely on bees for pollination (Champetier, 2011). What these previous simulators have in common is that they can only run for short periods of real time (e.g. one day (Camazine & Sneyd, 1991) or even a few hours (Bailis et al., 2010)). Most of them are also fairly restrictive in what type of environment (e.g. only one beehive; a specific number of foraging bees or nectar sources) or foraging activity (only either nectar or pollen but not both). None of them support an interactive user interface or GIS integration capabilities. In most cases, the simulation area, its size, and spatial resolution are not configurable. In contrast, we are interested in much longer periods of time (e.g. months) and, most importantly, to help answer the question about how foraging patterns and distributions can affect the local flora, hence potentially
306
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
contributing to the transmission of diseases or parasites, or supporting the distribution of certain biological control agents. From a technical perspective, the simulation models mentioned are locally executed, provide very limited capabilities for user interaction, and do not offer GIS integration or are spatially explicit. It is important to mention that for the last decade or so there has been a strong trend towards remotely-executed simulation, primarily based on the Web (Byrne, Heavey, & Byrne, 2010) and its interactive nature. As such, Web-based simulations have been profusely used in other types of simulation. Instances include the simulation of the glass-pressing process developed by Shi, Zhang, Rao, and Zhang (2008) uses a Web-based client approach that returns a visualisation through the Web browser. The ecological modelling framework EcoNet (Kazanci, 2007) is a noninteractive Web-based simulator which implements a model of population dynamics, although this model did not incorporate a spatial component, instead tracking the change within a single population and producing a plot of population change over time as output. Another example is INDISIM-YEAST, which also chose this approach for modelling yeast populations in a liquid (Ginovart & Can˜adas, 2008); a subject where the spatial component is less relevant due to the small physical scale. The closest effort to ours in terms of technical approach is probably by Xie and Yapa (2006). This is a simulation model of oil spills, which is Web-based using a clientserver architecture, includes a spatial component, and allows interactive, visual editing of the simulation through a Java applet. However this simulator does not incorporate a GIS component.
3.
Simulation model
The main simulation model has four interrelated models that include: i) a spatial model that represents the simulation landscape as a grid of cells; ii) a flora model that represents a patch with a number of flowers and seasons defined by nectar quality and quantity; iii) a weather model that considers time periods of temperature and daily hours of light; iv) a foraging model, which represents the foraging behaviour of a honeybee colony. Figure 1 corresponds to the general conceptual diagram that shows how the main simulation model and its underlying models are connected and what key concepts are shared among them. Our approach does not model the population dynamics of a beehive or the interactions among its individuals, in contrast with other existing work where this individual-based exchange of information leads to a globally intelligent selection of resources (Quijano & Passino, 2010; Tereshko & Loengarov, 2005; Zhang, Ouyang, & Ning, 2010). This is not only for the sake of simplicity, but also because it is not required given the focus and applicability of the simulator. This simulation model is focused on explaining the effects (e.g. total amounts of forage by beehive, visitation rates, foraging distances, etc.) of honeybee foraging in the flora surrounding the beehive on a time scale of days, which can be modelled by using published empirical data as Section 3.4 describes. Honeybees perform a foraging cycle (flying to
Fig. 1 e Conceptual diagram of the simulation model and its related models.
flowers, loading nectar, returning to the beehive, unloading, dancing, etc.) on a time scale of minutes, while recruiting other foragers occurs on a time scale of hours. The consequence from the mathematical point of view is the proportion of bees in each phase of the foraging cycle is essentially at equilibrium on the recruitment time scale (Cox & Myerscough, 2003). Therefore, it is only essential to consider the total number of active foragers at a beehive, based on the total population of its colony. Because traditional mathematical models that provide analytical solutions are only possible for simplified environments, we employ simulation modelling where the abstraction of the particular entity to predict is solved by steps numerically (Grant & Swannack, 2008). The simulation model is based on empirical data and expert knowledge in the field of honeybee ecology and management. Simulation models are often divided between deterministic and stochastic (Odum, 1971). Deterministic simulation models always produce the same result for a given set of input data and model parameters, whereas stochastic ones include one or more random processes, usually as a way of dealing with uncertainty (Law & Kelton, 1991). While the flora model (see Section 3.1), the weather model (see Section 3.2), and the spatial model (see Section 3.3) are deterministic, the foraging model (see Section 3.4) is stochastic. The simulation model specification formalism is DTSS (Discrete Time System Specification) (Nutaro, 2007), where the evolution of the simulation process is approximated by a difference equation. With d representing the duration of the simulation process (e.g. in days), we define Dd as the time-step used to partition d resulting in the first order difference equation yiþ1 ¼ yi þ Dd f ðyi ; xi Þ, which predicts the value of y at time step i þ 1 based on information from its previous time step i. In our particular scenario, Dd corresponds to one foraging trip t, with several of these trips taking place within each simulation day. More specifically, the number of foraging trips ti in a given day di with 0 i < d are calculated by using the
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
average temperature (in C) during the day Ti, the number of daylight hours hi, and the currently available nectar and pollen resources in the surrounding flower patches. Thus, the ¼ m hi , with m ¼ 0.8 maximum number of foraging trips tmax i when 15 < T < 33, m ¼ 0 when T < 8, or m ¼ 0.4 otherwise, while if the the eventual number of foraging trips will be ti tmax i daily nectar is depleted before the end of the day. Figure 2 provides the sequence of tasks that a typical simulation process would consist of. The spatial simulation parameters specified by the user through the Web-based interface, such as geographical bounds and resolution, are used to create the grid representing the simulation area by the spatial model described in Section 3.3. At the beginning of each simulation day di the appropriate average temperature Ti and hours of light hi are selected by applying the model . Also, the described in Section 3.2 and used to calculate tmax i current daily season of each flower patch is chosen using the model described in Section 3.1, and the nectar and pollen production for the day in each cell estimated. For each foraging trip in the day, the foraging model detailed in Section 3.4 is used to describe the honeybee foraging behaviour. At the end of a simulation day, a set of daily data is collected in anticipation for the final results report to be presented to the user through the Web interface. The report is organised by sections, each one corresponding to one simulated day. For each day, the user can inspect: i) how weather has changed in the simulation area; ii) how many foraging trips took place; and iii) technical information that includes: a) the number of cells of the grid representing the simulation area have been used by day; b) the length of time taken for the simulator to process this day; and c) how much RAM the server used. For each foraging trip in a day, the report includes a section with information about: i) the number of different foraging destinations (i.e. cells); ii) the amount of collected nectar; and iii) the amount of collected pollen. Each foraging
307
destination at any time can be examined by the user, hence obtaining a graphical view of visitation rates by honeybees in addition to nectar and pollen production and collection. Moreover, these results can be exported using the wellsupported comma-separated values (CSV) format so that it can be easily imported into popular spreadsheet or data analysis applications.
3.1.
Flora model
A flower patch corresponds to an area containing the following information: total number of flowers in the area and one or more sequential, cyclic seasons. A season indicates for how many days these flowers have the same average production of nectar and pollen per flower. There is an important distinction to make between gross and net nectar. Flowers in a patch may produce a certain amount of (gross) nectar per day. However, this nectar will have a certain ‘quality’ based on its sugar concentration and its subsequent energetic value. While bees will be able to carry a certain load of nectar regardless of its quality, they will prefer to forage flowers that produce higher quality nectar (Hainsworth, 1974). A traditional finite state machine (FSM) is used to represent a sequence of seasons corresponding to the same flower patch. A FSM is a mathematical abstraction useful to model behaviour explained by a finite number of states, transitions between them, and actions. A formal FSM can be defined by a five-tuple M ¼ hI; S; s0 ; O; di, where I is the set of inputs which are used to transition the machine among one of a finite set of states S with s0 ˛S indicating the initial state, O is the set of possible outputs, and d is the transition function d : S I/S. Hence, I O represents the alphabet of M. A state machine will switch to a different state sm ˛S after reading an input in ˛I. Each of the machine states s˛S specifies which state to switch for a given input i˛I. From any state there is only one
Fig. 2 e Overview of tasks performed by the main simulation model.
308
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
outputs O ¼ fh200; 0:7; 20i; h100; 0:9; 30i; h50; 1; 10i; h0; 0; 0ig, and inputs I ¼ f1; 2; .; 66g.
3.2.
Fig. 3 e A finite state machine with five states s0 ; .; s5 , input alphabet I[f0; 1g and outputs O[fl; rg.
transition for any allowed input. Figure 3 shows an example of a FSM known as transducer (or Mealy machine), where outputs depend on both the state and the input, whereas in Moore machines an output is determined only by the state. In this example, when the machine is at state s0 and receives 0 as an input, it will transition to state s1 and produce l as output. However, the state will not change if the received input is 1 while the output value will be r. Considering the set of flower patches in the landscape F ¼ ff0 ; f1 ; .; fm g, each flower patch fi ¼ hli ; Ei i; 0 i m is represented by the total number of flowers in it li and the list of flower seasons Ei ¼ fe0 ; e1 ; .; ek g. Each season ej ; 0 j k is described by a four-tuple hdj ; nj ; qj ; pj i where dj expresses the number of days that this season ej lasts, nj indicates the average amount of nectar production per day per flower (in mg) during this season, qj ˛½0; 1 corresponds to the nectar quality rate value, and Pj indicates the average amount of pollen production per day per flower (in mg) during the season. The FSM representing the seasons of a flower patch has states S ¼ fs0 ; s1 ; .; sk g to represent E, hence jEj ¼ jSj. The input alphabet includes all the consecutive day numbers for E, P thus consisting of I ¼ f1; 2; .; dg, where d ¼ 1 þ ki¼0 di . The output alphabet is given by all the three-tuples extracted from E, producing O ¼ fhn0 ; q0 ; p0 i; hn1 ; q1 ; p1 i; .; hnk ; qk ; pk ig. The transition from one machine state to the next occurs when the input for the current state is the last day of the season represented. A transition from the last state (i.e. flowers season) to the first state is automatically created to represent the cyclic nature of seasons. This transition function has to be triggered on each simulation day for all flower areas. Figure 4 shows the machine representing the four flower patch seasons h10; 200; 0:7; 20i, h5; 100; 0:9; 30i, h20; 50; 1; 10i, and h30; 0; 0; 0i. This means a total natural cycle of 65 days, thus t ¼ 66. This corresponds to four states S ¼ fs0 ; s1 ; s2 ; s4 g, with
Weather model
The simulator uses weather information that is crucial for determining the number of daily foraging trips. This information includes temperature (in C) and number of daylight hours, with both measures applying to the whole simulation area. Thus, given a list of weather periods W ¼ fw0 ; w1 ; .; wk g, an arbitrary weather period wi, with 0 i k, is expressed as a list of three-tuples, each one specifying a number of days along with the temperature and the number of daylight hours during them such as hdi ; Ti ; hi i. The sum of the number of days P for all these pairs d ¼ ki¼0 di corresponds to the total number of days of simulation. Like the Flora Model described in Section 3.1, this model is based on a FSM. In this case, each state si ˛S represents a number of consecutive days di with the same temperature Ti and daylight hours hi. Hence the number of states jSj is determined by the number of temperature periods. The input alphabet includes all the added period days I ¼ f0; 1; 2; .; dg. The output oi ˛O that corresponds to si is the temperature and daylight hours for that period hTi ; hi i. The transition from one state to the next occurs when the input for the current state corresponds to the last day of the period it represents.
3.3.
Spatial model
This section describes the spatial model, which provides a two dimensional abstraction of the simulation area and related entities, including the flower patches, beehives, and honeybee foraging activity. This spatial model is based on the one described in Garcı´a Adeva et al. (2012), which is here extended by adding data to cells in order to track particular aspects of the honeybee foraging activity. In this model, the simulation area consists of a landscape whose bounds are initially defined by the user through the simulator’s user interface. The extent will usually vary from a few hectares to several hundred thousand km2. The landscape can contain areas (i.e. polygons) representing floral patches and locations indicating beehives. Figure 5 provides a graphical representation of this landscape expressed as a rectangular planar region L of width wL and height hL . The boundaries of L are given by two global coordinates ðcbl ; ctr Þ ¼ ððlatb ; lonl Þ; ðlatt ; lonr ÞÞ corresponding to the bottom-left and top-right corners of the extent. To calculate the distance d between two coordinates c1 ¼ ðlat1 ; lon1 Þ and c2 ¼ ðlat2 ; lon2 Þ in L, we define the following equation based on the Spherical Law of Cosines: rðc1 ; c2 Þ ¼ arccos ðsin lat1 sin lat2 þ cos lat1 cos lat2 cos ðlon1 lon2 ÞÞ<
Fig. 4 e The FSM representing the 65-day cycle of a flower patch with four seasons h10; 200; 0:7; 20i, h5; 100; 0:9; 30i, h20; 50; 1; 10i, and h30; 0; 0; 0i.
where < ¼ 6371 is a constant defining the Earth’s radius in km. Discretisation of L is performed by means of a rectangular sparse grid G with its maximum number of horizontal and vertical cells being wG hG . A grid cell is defined by g ¼ ðx; yÞ˛G, where x; y˛Z; 0 x < hG ; and 0 y < hG , with only those non-empty cells being created dynamically when
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
309
Fig. 5 e An extent of landscape L discretised by grid G.
needed. Therefore, wG and hG also indicate the simulation resolution, with its maximum value limited by such a cell size where any given cell contains at least 1 flower. The translation from a coordinate c ¼ ðlat; lonÞ˛L to its cell position g ¼ ðx; yÞ˛G is calculated with
km2, we apply Meister’s formula that calculates the area of a non-self-intersecting polygon using Euclidean geometry: xi ¼ < cos lati cos loni ; yi ¼ < cos lati sin loni ;
x ¼ Prh ðcÞwG wL R; y ¼ Prv ðcÞhG hL R where rv ðcÞ and rh ðcÞ indicate vertical and horizontal distances respectively from c to cbl in L: rh ðc1 Þ ¼ rðcml ; ðlatm ; lon1 ÞÞ; rv ðc1 Þ ¼ rðcbm ; ðlat1 ; lonm ÞÞ: This discretisation process can almost be considered a map cylindrical projection, where the sphere is projected on the walls of a cylinder by mapping circles of longitude (i.e. meridians) as equidistant vertical straight lines and circles of latitude as horizontal straight lines. One of the key differences is how in the traditional cylindrical projection, the projection consists of a cylinder wrapped around the globe at the equator, whereas in this case the cylinder is centred around the midpoints of L. The immediate benefit of this approach is a significant reduction in distortion, hence obtaining a more precise spatial representation. All grid cells are square in shape and have the same size sG ¼ l2 , where sG , wG , and hG are calculated based on the horizontal resolution selected by the simulator’s user (typical values will be in the range 50e500). Thus, obtaining hG and sG in km2 is achieved by
hG ¼
Q S
hL wG wL ;
AðV L Þ ¼
wL hL : wG hG
An area representing a flower patch is expressed as a closed polygon defined by the list of its ordered (either clockwise or counterclockwise) vertex coordinates V L ¼ ðv0 ; v1 ; v2 ; .; vk ; v0 Þ with vi ðlati ; loni Þ˛L; 0 i k. To calculate the extent of V L , in
k1 1X xiþ1 ¼ xi yiþ1 xiþ1 yi : yiþ1 2 i¼0
Discretising an area V L to obtain its corresponding V G simply consists of translating each vertex vi ˛L into corresponding cell positions vi ˛G. Another instance of a frequently required operation is to determine whether a particular cell is included within an area, which we achieve by applying the Jordan curve theorem as proposed by Haines (1994). A cell may contain one or more beehives. In addition, a cell keeps track of how much nectar and pollen it contains by taking into account how many flower patches it is contained by and the season that currently applies to each of these flower patches. In most cases, a cell will be within a single flower patch, although the opposite is possible in cases where flower patches overlap. Hence, given an arbitrary cell c within flower patches F ¼ ff0 ; f1 ; .; fm g, where a single flower patch is represented by fi ¼ hli ; Ei i; 0 i m, li indicates the total number of flowers in it and Ei ¼ fe0 ; e1 ; .; ek g the sequence of seasons, the total amount of forage in c is calculated by lc ¼
m1 X AðcÞ li ; A fi i¼0
nc ¼
m1 X
lc ni ;
i¼0
qc ¼
m1 X nc qc þ ni qi i¼0
sG ¼
k1 1X xi 2 i¼0 yi
pc ¼
m1 X
nc þ ni
;
lc pi ;
i¼0
where lc corresponds to the average number of flowers per cell c, nc corresponds to nectar and pc to pollen, while qc indicates the weighted nectar quality rate that applies to nc.
310
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
3.4.
Foraging model
A key feature of honeybee foraging is how this activity extends over an area larger than 100 km2 around the beehive. The foraging behaviour of honeybees living in nature is characterised by daily changes to their routine and distribution across flower patches and regular foraging at food sources located several km from their beehive (Seeley, 1985). From a given population of forager bees at a beehive, this model aims to determine: i) how many foragers will visit particular flowers depending on proximity, nectar quality, nectar quantity, and pollen quantity; and ii) the number of foraging trips for a particular day. This is a stochastic model due to some of its actions being supported by to the use of weighted random rounding. This stochastic rounding operation happens in multiple aspects of the model including the calculation of foraging population including the different types such as nectar foragers, pollen foragers, and nectar and pollen foragers, the distribution of foragers among cells containing flowers, etc. Thus, we define PxS as the stochastic, weighted rounding operation applied to x, which is based on a Bernoulli random variate, such as PxS ¼
PxR QxS
if ε < fxg if ε fxg
(1)
where x is a positive real number x˛R0 , PxR is the floor function defined by x ¼ max fm˛Zm xg, PxR is the ceiling function defined by x ¼ fn˛Zn xg, fxg is the fractional part of x as given by fxg ¼ x PxR, and ε is the random function with ε : R0 /½0; 1. There is a kernel that defines both the foraging distance and the number of individuals finding a flower patch at a certain distance, by describing the spatial distribution of foragers in the honeybee vicinity as a function of distance from the beehive and the size of flower patches. Foraging is similar to the concept of local dispersal frequently found in insect modelling. There are many kernelbased dispersal models found in the literature including Laplace, Gaussian, Ribbens (Clark, Macklin, & Wood, 1998), and Cauchy (Shaw, 1995). They may be considered arbitrary equations that fit particular dispersal data sets well. The use of these kernels requires that previously parameters for the chosen distribution to be estimated by registering the
a
dispersal distances between occurrences in distribution data from one time period to the next, so that the distribution can be fitted to the frequency distribution of distances (Pitt, 2008). Our dispersal kernel, however, uses an empirical probability distribution to determine foraging distances based on published data (Seeley, 1985; Visscher & Seeley, 1982) and depicted by Fig. 6(a). For instance, 7% of foraging bees will do so within ð0:2; 0:4 km, 28% will do it within ð0:4; 1 km, and so forth. This distribution is based on the assumption that the beehive is surrounded by flower patches rich in forage (i.e. nectar and pollen) and accounts for related activities such as travel times between the beehive and the flower patches, recruitment of other foragers by dancing, etc. Generally, most foragers concentrate their activity in a radius of 2 km from their beehive, but they will travel up to 3 times this distance if needed for survival (Beekman & Lew, 2008). We consider that 30% of the total colony population will be actively foraging during the daylight hours, with 58% of them foraging only nectar, 25% only pollen, and the remaining 17% both pollen and nectar (Winston, 1991). Figure 6(b) provides the graphical representation of the probability of suitable flower patches being detected by nearby foragers based on distance and size. This is based on the assumption that a flower patch will be more detectable by foragers the larger or the closer it is, and vice versa. Therefore, in order to express this probability of detection, we use the probability distribution given by Aðf Þ r 1 g; 1 ; Pf ¼ min 1000 10 where Aðf Þ˛½0; 1000 indicates the flower patch size in m2, r ˛½0; 10 is the distance in km between the flower patch and the beehive, and g > 0 corresponds to the scale parameter, where the larger its value the higher the value of Pf. These two distributions are combined in order to determine the weighted probability that each cell within a flower patch is visited by honeybees. The algorithmic approach for performing this process is based on the grid G, where a cell gi ¼ ðxi ; yi Þ˛G in this grid represents a discrete location in the simulation area, contained by flower patches F ¼ ff0 ; f1 ; .; fm g. The foraging bees from a beehive are represented by B ¼ fb0 ; b1 ; .; bn g.
b
Fig. 6 e Probability distributions of the dispersal kernel used by the foraging model.
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
Table 1 e Basic scenarios and their result to verify the simulation model. Scenario No flower patches around the beehive A flower patch close to the beehive with no nectar or pollen A flower patch rich in nectar and pollen at 20 km from the beehive Temperature below 8 C No daylight hours Beehive has only 1 honeybee
Result No forage activity No forage activity
No forage activity
No forage activity No forage activity Foraging activity 1 only in about of the times 3
By traversing all the possible dispersal distance ranges shown by Fig. 6(a) such as rba ; rcb ; rdc ; . along with the minimum possible distance increment Dr in G (i.e. the length of a cell side) for as many times as possible (i.e. maxfn xi þ 1; xi ; m yi þ 1; yi g), we obtain a list of cells C ¼ ðc0 ; c1 ; .; ck Þ belonging to the ring that corresponds to the first foraging distances interval ½rba ; rcb Þ. This ring of cells is created by applying the well-known midpoint algorithm (van
311
Aken, 1984), thus creating circles of cells as long as rba r < rcb . Those cells from C that are not contained by a flower patch in F are discarded. Also discarded are those cells for which applying the stochastic weighted rounding function defined by Eq. (1) to their corresponding probability of being found Pf results in 0, hence producing the effect that they could not be found by foraging bees due to their location. The remaining cells of C are ordered according to their Pf, to then be used to sequentially distribute the bees from B. The number of bees from B assigned to each remaining cell in C is proportional to the relative amount of forage in the cell. The type of forage considered for this calculation depends on the type of forager: nectar foragers will use the remaining net nectar in the cell, pollen foragers will use the remaining pollen in the cell, while nectar and pollen foragers will use both net nectar and pollen in the cell. The forage in those cells has to be depleted accordingly. As default values, we selected an average load of 40 mg of gross nectar and 20 mg of pollen per bee per trip (Beekman & Lew, 2008; Feuerbacher, Fewell, Roberts, Smith, & Harrison, 2003; Winston, 1991), although these values can be changed by the user through the configuration section of the Web interface. With these values, the depletion process consists of simply removing from the remaining reserve of cell forage the average load of nectar, pollen, or both as many times as
a
b
c
d
Fig. 7 e Location of foraging honeybees on day 1 of the simulation. Dots represent beehives and squares represent a groups of foragers, with the darker the colour the higher the population.
312
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
honeybees chose to forage in the cell such with nc ¼ nc min ðnc ; 40jBjÞ, and pc ¼ pc min ðpc ; 20jBjÞ, while qc remains unchanged.
4.
Evaluation
Evaluating simulation systems a functional point of view is important in order to ensure that the right system has been correctly built (Sojda, 2007). However, because simulation
systems handle complex and poorly structured problems, they are difficult to empirically evaluate. Accordingly, a frequently chosen approach to tackle this difficulty when evaluating a simulation model consists of dividing the effort into two separate tasks: verification and validation (Mihram, 1972). Verification consists of showing the simulation models and software implementation are correct, complete, and coherent, while validation consists of comparing simulation results with observations from the real world. Verification should precede validation (Mihram, 1972).
a
b
c
d
e
f
Fig. 8 e Nectar and pollen depletion after four foraging trips on day 1. Dots represent beehives and squares represent a flower patch section, with the darker the colour the higher the amount of nectar or pollen left.
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
The ideal method of validation would consist of comparing results with a preselected gold standard d a dataset with sufficient accuracy and precision available to the research community. For ecological simulation this gold standard would usually be high-quality historical records. However, it is rare to find one in the natural sciences (Sojda, 2007) and indeed it was the case here d we were not able to find a dataset describing a scenario in a way that can be used to propose a validation strategy. In consequence, our evaluation was limited to verifying the simulation model, by means of
313
performing simulation experiments of a selected set of use cases and then comparing results with expectations. We started the verification process by proposing several basic simulation scenarios that consider extreme use cases. The expected result can be guessed by common sense. Most of the use cases summarised in Table 1 correspond to situations where no foraging activity should be expected, because the temperature is too cold, the flower patch within reach but with no forage, the flower patch has plenty of forage available but it is too far, etc. Although some of these
a
b
c
d
e
f
Fig. 9 e Nectar and pollen depletion after four foraging trips on day 11. Dots represent beehives and squares represent a flower patch section, with the darker the colour the higher the amount of nectar or pollen left.
314
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
a
b
c
d
Fig. 10 e Nectar and pollen depletion after four foraging trips on day 16. Dots represent beehives and squares represent a flower patch section, with the darker the colour the higher the amount of nectar or pollen left.
scenarios were obviously unrealistic (e.g. no daylight hours or a single bee in the beehive), they were still important and useful as a test mechanism to verify that the simulator worked as expected. Next, we proposed additional verification based on a realistic scenario.1 The data generated by running simulations on this scenario and used to produce the graphs shown below is fully available through the Web user interface. The scenario in question comprised a simulation area of 166.28 km2 (16.35 km by 10.17 km) represented by a spatial resolution of 100 by 63 cells (i.e. 0.026 km2 or 2.6 Ha per cell). There are two weather periods: i) 15 days at 22 C with 12 daylight hours; and ii) 15 days at 15 C with 10 daylight hours. There are two beehives, one small with a population of 10,000 bees and a larger one with 20,000 bees. Between the two beehives there is a medium flower patch of about 5 km2 with 500 flowers that produce an average of 40 mg of nectar and 10 mg of pollen per flower per day during the first 15 days, to then reduce their production to 5 mg of nectar (quality rate 0.8) and 2 mg of pollen for the remaining 15 days. East to the larger beehive there is a larger flower patch covering an area of 1
Scenario data available at http://beetle.csse.uwa.edu.au/ honeybee/?bioeng.
21 km2 and containing 5000 flowers. During the first 10 days they produce 50 mg (quality rate 0.5) of nectar and 20 mg of pollen on average per day per flower. During the remaining 20 days the production drops to 10 mg of nectar (quality rate 0.2) and 5 of pollen. The purpose was to use this scenario to examine different types of results produced by the simulation in order to verify the following aspects and effects of the foraging activity: i) foraging bee location and activity over the flower patches; ii) nectar and pollen depletion along several days where some key events took place (e.g. change in weather or in flower production capability); and iii) amount of forage collected by day and during the total duration of the simulation. Figure 7 provides a graphical view of the whereabouts of honeybees on day 1 while performing several foraging trips to the two flower patches, which can be interpreted as an indication of foraging activity along the day. During this day, the maximum number of foraging bees from the small beehive was around 3000. In addition, the maximum expected number from the large beehive was 6000. This accounts for a combined 5220 honeybees foraging nectar, 2250 foraging pollen, and 1530 foraging both nectar and pollen. Foraging bees will prefer to start their work in those parts of the flower patches that are close to the beehive, and this can be verified to happen in
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
Fig. 8(a) and (b). However, as they deplete the nectar and pollen reserves available for the day, these foragers will have to move farther from the beehive, consequently making their work harder since they will have to fly longer distances. Figure 8(c) and (d) confirmed that this is the case by the 5th foraging trip, with it becoming even more prominent by looking at Fig. 8(e) and (f), which correspond to the last foraging trip for the day. These figures also show how the amount of daily resources produced by the smaller flower patch was not enough to provide forage for the whole day, so that by the 5th foraging trip the number of honeybees foraging in this flower patch was very small, while by the end of the day there was virtually no foraging activity left on this flower patch. In contrast, the larger flower patch still offered food based on the small groups of honeybees still foraging on it. However, the honeybees did not continue foraging due to the lack of daylight. Figure 8 offers a graphical view of how nectar and pollen were depleted from the flower patches on day 1 of the simulation. It is interesting to compare these results to those on day 11 (Fig. 9), when the larger flower patch had a season change and its nectar production dropped significantly, and also to those on day 16 (Fig. 10), when there was a change in temperature (thus becoming colder) and number of daylight
315
hours thus producing a reduction in honeybee foraging activity, while at the same time the smaller flower patch had a drop in its nectar production due to seasonal change. Figure 11 shows the remaining nectar at the end of several foraging trips on day 21, where the situation with regards to weather and flower patch seasonality was about the same as on day 16 and hence the results quite similar. As a consequence of weather and flower patch seasonality the number of daily foraging trips between day 1 and day 15 was 9e10, whereas from day 16 to day 30 it was only 4. The last part of this verification focused on amount of collected forage. Figure 12 shows the daily and total amount of collected nectar and pollen by both beehives. Figure 12(a) and (b) show evidence of how the change in weather after day 10 produced a significant reduction in both the daily and total foraged nectar and pollen. Something similar, but with a much lower impact, happened after day 15 as a consequence of the weather change at the time. This effect could also be observed in the total amount of nectar during the 30 days of simulation by looking at Fig. 12(c) and (d). It is important to note that the frequent peaks in Fig. 12(a) and (b) were due to the declining amount of forage collected by the bees as the day progressed. Therefore, the amount of both nectar and pollen in the last foraging trips of the day were
a
b
c
d
Fig. 11 e Nectar and pollen depletion after four foraging trips on day 21. Dots represent beehives and squares represent a flower patch section, with the darker the colour the higher the amount of nectar or pollen left.
316
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
a
b
c
d
Fig. 12 e Daily and total amount of nectar and pollen foraged.
usually smaller than during the first foraging trips of the day, due to availability.
5.
Conclusions
The importance of the ecological role of honeybees by contributing to biodiversity and the production of crops and wild plants through pollination is well understood. Their honey production is also an important economic factor. Even though simulation models can improve our understanding of honeybees, hence helping stakeholders such as biosecurity managers and beekeepers, there are few models available and many are severely limited in their scope. We tackle this by offering a simulation tool that supports long time periods, is highly configurable, offers an interactive user interface, and
produces very detailed results. Its simulation model is stochastic and discrete in time and space, and can support any number of beehives and floral patches over a two dimensional landscape of unrestricted size and resolution. The different underlying models encapsulate well-established ecological knowledge. The verification of the simulation model included several experiments that dealt with nectar and pollen foraging under different conditions of weather and flower patch seasonality. Although it would have been useful to complement this evaluation by validating the simulation results against an existing dataset, the lack of such a dataset made this impossible. We note this as possible future work along with an analysis of the engineering merits of the simulator, including an evaluation of system measures such as performance, time response, and scalability.
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
Acknowledgements The author would like to thank Lori Lach from the Centre for Integrative Bee Research at the University of Western Australia and Rob Manning from the Department of Agriculture and Food Western Australia for pointing at relevant bibliography and providing feedback. The author would like to acknowledge the support of the Australian Government’s Cooperative Research Centre Program.
references
van Aken, J. (1984). An efficient ellipse-drawing algorithm. IEEE Computer Graphics and Applications, 4, 24e35. Bailis, P., Nagpal, R., & Werfel, J. (2010). Positional communication and private information in honeybee foraging models. In Proceedings of the 7th international conference on Swarm intelligence (pp. 263e274). Berlin, Heidelberg: Springer-Verlag. Beekman, M., & Lew, J. B. (2008). Foraging in honeybeesewhen does it pay to dance? Behavioral Ecology, 19(2), 255e261. Byrne, J., Heavey, C., & Byrne, P. J. (2010). A review of Web-based simulation and supporting tools. Simulation Modelling Practice and Theory, 18(3), 253e276. Camazine, S., & Sneyd, J. (1991). A model of collective nectar source selection by honey bees: self-organization through simple rules. Journal of Theoretical Biology, 149(4), 547e571. Champetier, A. (2011). The foraging economics of honey bees in almonds. In Annual meeting of the Agricultural and Applied Economics Association, Pittsburgh, PA. Clark, J. S., Macklin, E., & Wood, L. (1998). Stages and spatial scales of recruitment limitation in southern Appalachian forests. Ecological Monographs, 68(2), 213e235. Cox, M., & Myerscough, M. (2003). A flexible model of foraging by a honey bee colony: the effects of individual behaviour on foraging success. Journal of Theoretical Biology, 223, 179e197. Crozier, R. H., & Page, R. E. (1985). On being the right size: male contributions and multiple mating in social Hymenoptera. Behavioral Ecology and Sociobiology, 18(2), 105e115. Feuerbacher, E., Fewell, J., Roberts, S., Smith, E., & Harrison, J. (2003). Effects of load type (pollen or nectar) and load mass on hovering metabolic rate and mechanical power output in the honey bee Apis mellifera. Journal of Experimental Biology, 206(Pt 11), 1855e1865. von Frisch, K. (1967). The dance language and orientation of bees. Cambridge, MA: Belknap Press of Harvard University Press. Garcı´a Adeva, J. J., Botha, J. H., & Reynolds, M. (2012). A simulation modelling approach to forecast establishment and spread of Bactrocera fruit flies. Ecological Modelling, 227(0), 93e108. Garcı´a Adeva, J. J., Lach, L., & Reynolds, M. (2011). Simulation of honeybee nectar foraging for determining effects on local flora. In 19th International congress on modelling and simulation, the Modelling and Simulation Society of Australia and New Zealand conference, (pp. 2507e2513). Perth, Australia. Garcı´a Adeva, J. J., Sousa-Majer, M. J., Botha, J. H., Hanbury, C. D., Hardie, D. C., & Reynolds, M. (2010). Modelling the establishment and spread of emergency plant pests (EPPs) in Australia: Simulate or suffer. In Global biosecurity conference (p. 101). Brisbane, Australia. Ghosh, S., & Marshall, I. (2005). Simple model of collective decision making during nectar source selection by honey bees.
317
In Workshop on memory and learning mechanisms in autonomous robots (ECAL 2005) p. 10. Ginovart, M., & Can˜adas, J. (2008). INDISIM-YEAST: an individualbased simulator on a website for experimenting and investigating diverse dynamics of yeast populations in liquid media. Journal of Industrial Microbiology & Biotechnology, 35(11), 1359e1366. doi:10.1007/s10295-008-0436-4. Grant, W. E., & Swannack, T. M. (2008). Ecological modeling: A common-sense approach to theory and practice. Wiley-Blackwell. Haines, E. (1994). Graphics Gems IV. Chapter Point in Polygon Strategies. Morgan Kaufmann. Hainsworth, F. R. (1974). Food quality and foraging efficiency. Journal of Comparative Physiology A: Neuroethology, Sensory, Neural, and Behavioral Physiology, 88, 425e431. Kazanci, C. (2007). EcoNet: a new software for ecological modeling, simulation and network analysis. Ecological Modelling, 208(1), 3e8. Keller, L., & Reeve, H. K. (1994). Genetic variability, queen number, and polyandry in social Hymenoptera. Evolution, 48(3), 694e704. Law, A. M., & Kelton, W. D. (1991). Simulation modeling and analysis (2nd ed.). New York, NY: McGraw-Hill. Mihram, G. A. (1972). Some practical aspects of the verification and validation of simulation models. Operational Research Quarterly, 23(1). Moosbeckhofer, R., Loncaric, I., Oberlerchner, J., Persen, U., Ertl, C., & Donat, C. (2008). Use of honeybees (Apis mellifera) as vectors for fire blight antagonists in field experiments. In Acta horticulturae 793: XI International workshop on fire blight (pp. 461e465). Morse, R., & Calderone, N. (2000). The value of honey bees as pollinators of U.S. crops in 2000. Bee Culture Magazine, 128, 1e14. Nutaro, J. (2007). Discrete event simulation of continuous systems. In Handbook of dynamic systems modeling. Chapman & Hall/CRC Press. Odum, E. P. (1971). Fundamentals of ecology. Philadelphia: W.B. Saunders Company. Pitt, J. P. W. (2008). Modelling the spread of invasive species across heterogeneous landscapes. PhD thesis, Lincoln University. Potts, S., Biesmeijer, J., Kremen, C., Neumann, P., Schweiger, O., & Kunin, W. (2010). Global pollinator declines: trends, impacts and drivers. Trends in Ecology and Evolution, 25(6), 345e353. Quijano, N., & Passino, K. M. (2010). Honey bee social foraging algorithms for resource allocation: theory and application. Engineering Applications of Artificial Intelligence, 23(6), 845e861. Seeley, T. D. (1985). Honeybee ecology: A study of adaptation in social life. Princeton, NJ: Princeton University Press. Seeley, T. D. (1995). The wisdom of the hive: The social physiology of honey bee colonies. Cambridge, MA: Harvard University Press. Seeley, T. D., Camazine, S., & Sneyd, J. (1991). Collective decision-making in honey bees: how colonies choose among nectar sources. Behavioral Ecology and Sociobiology, 28(4), 277e290. Shaw, M. W. (1995). Simulation of population expansion and spatial pattern when individual dispersal distributions do not decline exponentially with distance. Proceedings of the Royal Society London Series B, 259, 243e248. Shi, S., Zhang, G., Rao, Y., & Zhang, Y. (2008). A Web-based remote simulation system for the glass pressing process. Simulation, 84(1), 27e40. Sojda, R. (2007). Empirical evaluation of decision support systems: needs, definitions, potential methods, and an example pertaining to waterfowl management. Environmental Modelling and Software, 22(2), 269e277.
318
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
Tereshko, V., & Loengarov, A. (2005). Collective decision making in honey-bee foraging dynamics. Computing and Information Systems, 9(3), 1e7. Visscher, P. K., & Seeley, T. D. (1982). Foraging strategy of honeybee colonies in a temperate deciduous forest. Ecology, 63(6), 1790e1801. de Vries, H., & Biesmeijer, J. C. (1998). Modelling collective foraging by means of individual behaviour rules in honeybees. Behavioral Ecology and Sociobiology, 44, 109e124.
Winston, M. L. (1991). The biology of the honey bee. Cambridge, MA: Harvard University Press. Xie, H., & Yapa, P. D. (2006). Developing a Web-based system for large-scale environmental hydraulics problems with application to oil spill modeling. Journal of Computing in Civil Engineering, 20(3), 197e209. Zhang, C., Ouyang, D., & Ning, J. (2010). An artificial bee colony approach for clustering. Expert Systems with Applications, 37(7), 4761e4767.
Available online at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/issn/15375110
Research Paper
Simulation modelling of nectar and pollen foraging by honeybees Juan Jose Garcia Adeva a,b,* a b
School of Computer Science and Software Engineering, University of Western Australia, Perth, Australia Cooperative Research Centre for National Plant Biosecurity, Canberra, Australia
article info
This paper describes a simulation model for the foraging activity of honeybees (Apis mel-
Article history:
lifera) as they collect nectar and pollen. Understanding such behaviour and observing the
Received 14 December 2011
resulting foraging spread patterns and distribution are vital from an ecological perspective
Received in revised form
and has useful applications in areas such as agriculture or biosecurity. The simulator is
20 March 2012
based on a spatio-temporal model and implemented following a Web-based architecture.
Accepted 4 May 2012
This simulation model takes into account of detailed floral distributions, weather, hours of
Published online 13 June 2012
daylight, and variations in daily nectar and pollen production by flowers. The produced simulation technology includes a user interface that allows the introduction of the simulation parameters (e.g. spatial resolution and extent, model configuration settings, visual specification of the simulation environment on the map) and receive a detailed results report. The simulation model is evaluated by performing multiple experiments and verifying the results obtained. ª 2012 IAgrE. Published by Elsevier Ltd. All rights reserved.
1.
Introduction
Foraging is the act of honeybees (Apis mellifera) fetching food from sources scattered around their beehive. In this work, we consider foraging to include pollen and nectar from flowers in order to feed the colony, excluding other types of provisions such as water. Nectar provides energy while pollen is the main source of proteins. Foraged pollen and nectar is stored in the beehive, in the form of honey, in order to have enough food for times when flowers do not produce enough nectar or pollen, such as winter (Seeley, 1995). Therefore, being efficient in this foraging activity is key to the colony’s continued survival. Floral patches differ importantly in several factors that affect the profitability of the foraging effort on the patch (Winston, 1991), including the quality and quantity of their
nectar and pollen, the size of the patch, and the distance from the colony. This overall profitability of a floral patch is communicated by individual foragers to other bees in the colony via complex behaviours and other signals (von Frisch, 1967). These other individuals may alter their behaviour accordingly, and their responses may in turn affect the behaviour of the original forager (Seeley, 1985). Thus foraging in honeybees is a dynamic and highly interactive process that is driven by the sharing of information about the location and quality of resources in the surrounding environment. The intensity of this foraging activity and especially its extent over the local flora can have an influence on aspects such as crop pollination and plant biosecurity. There are several reasons for this. On the one hand, the foraging behaviours and hairiness of honeybees mean that they often
* School of Computer Science and Software Engineering, University of Western Australia, Perth, Australia. Tel.: þ61 8 6488 3816; fax: þ61 8 6488 1089. E-mail address: [email protected]. 1537-5110/$ e see front matter ª 2012 IAgrE. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.biosystemseng.2012.05.002
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
inadvertently transfer substances, such as pollen, among plants, hence potentially having important consequences to crop pollination. Honeybees are capable of increasing the yield in 96% of animal-pollinated crops (Potts et al., 2010), and the value of the increased yield and quality achieved by honeybee pollination of crops in the U.S. is estimated as $14.6 billion (Morse & Calderone, 2000). On the other hand, because honeybees can also act as vectors of some plant diseases such as Fire Blight (Erwinia amylovora) and also as vectors of biological control agents (e.g. the microbe Aureobasidium pullulans, which can be used to attack E. amylovora (Moosbeckhofer et al., 2008)), understanding patterns of foraging behaviour will aid in implementing prevention or facilitation measures, respectively, thus contributing to important aspects of disease prevention or eradication in the context of plant biosecurity. This last scenario describing the influence of honeybees in spread of plant diseases is the main motivation for this work. Plant biosecurity is a major concern in a country like Australia where economic viability of plant industries rely to an extent on freedom from pests and diseases. Strategies to test preemptive surveillance and estimate rates of spread of pests or diseases in a spatial environment in real time are not currently available. A ‘what if’ scenario simulation tool is needed to address this shortcoming (Garcı´a Adeva et al., 2010). The aim of the main project, where this current honeybee simulation work resides, is to provide such decision support tools. In consequence, simulating honeybee foraging behaviour will help us to determine the extent and rate at which surrounding flora may be affected by a disease such as Fire Blight. This work consists of a simulation model that is stochastic and discrete in time with 1-trip time steps, which is approximated by difference equations. It is also explanatory (Grant & Swannack, 2008), meaning that its behavioural rules depend on well-established expert knowledge or published empirical data. It can support any number of weather periods, beehives and floral patches over a two dimensional landscape of unrestricted size and varying resolution. Each flower patch can have any shape, size, location, and number of seasons that define the daily nectar and pollen production, therefore extending some previous preliminary work limited to only nectar foraging (Garcı´a Adeva et al., 2011). A weather period is defined by average daily temperature and number of daylight hours. Because the produced technology shares some modelling aspects and resulting features with previous work by the same authors (Garcı´a Adeva, Botha, & Reynolds, 2012), this will be noted where appropriate. From its implementation point of view, the simulator usage is Web based, so that users do not have to install or update software, or worry about enough computational power or platform compatibility. The user interface allows the selection of the simulation parameters (e.g. spatial resolution and size, model configuration settings, visual specification of the simulation environment on the map) and receive a highly detailed report of results. A simulation of several weeks usually takes a few seconds to complete. The structure of this paper is as follows: in Section 2, we discuss some related work in honeybee modelling or Web-
305
based simulation. Section 3 describes the main simulation model and related sub-models. Section 4 offers an evaluation. We finish with some concluding remarks in Section 5.
2.
Related work
There is a number of research efforts in simulating or modelling different aspects of honeybees. This section analyses the most relevant instances of previous honeybee simulation work that we have found. de Vries and Biesmeijer (1998) implemented an agentbased model intended to replicate the foraging experiment that Seeley, Camazine, and Sneyd (1991) had previously performed by using behavioural rules of individual bees. Also based on this same experiment, Camazine and Sneyd (1991) proposed a foraging model based on differential equations for foraging, which considers nectar quality in the decision making process for resource abandonment and recruitment of foraging bees. Cox and Myerscough (2003) proposed a mathematical model of foraging in honeybee colonies. The effort was mainly focused on examining the effect of a small set of foragingrelated parameters (e.g. visit rate, probability of dancing, duration of dancing, probability of abandonment) on the foraging response and consequent nectar intake. This model was intended as an experiment to support the biological hypothesis of ‘genetic variance’ that argues that genetically heterogeneous colonies are better protected against misfortune than homogeneous colonies (Crozier & Page, 1985; Keller & Reeve, 1994). Ghosh and Marshall (2005) proposed a model of learning and collective decision making in honeybees engaged in nectar foraging. However, this model is limited to one of two foraging sites where one site has better nectar than the other in order to understand decision making processes by interacting swarm entities. Bailis, Nagpal, and Werfel (2010) implemented a simulator that models individual foraging bees at a very high level of detail in order to compare the value of sharing food source position information through the waggle dance and private information about previously encountered food sources. The most recent effort consists of a model of honeybee foraging behaviour to estimate almond crop yields as functions of honeybee densities in order to understand the economic behaviour of growers who rely on bees for pollination (Champetier, 2011). What these previous simulators have in common is that they can only run for short periods of real time (e.g. one day (Camazine & Sneyd, 1991) or even a few hours (Bailis et al., 2010)). Most of them are also fairly restrictive in what type of environment (e.g. only one beehive; a specific number of foraging bees or nectar sources) or foraging activity (only either nectar or pollen but not both). None of them support an interactive user interface or GIS integration capabilities. In most cases, the simulation area, its size, and spatial resolution are not configurable. In contrast, we are interested in much longer periods of time (e.g. months) and, most importantly, to help answer the question about how foraging patterns and distributions can affect the local flora, hence potentially
306
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
contributing to the transmission of diseases or parasites, or supporting the distribution of certain biological control agents. From a technical perspective, the simulation models mentioned are locally executed, provide very limited capabilities for user interaction, and do not offer GIS integration or are spatially explicit. It is important to mention that for the last decade or so there has been a strong trend towards remotely-executed simulation, primarily based on the Web (Byrne, Heavey, & Byrne, 2010) and its interactive nature. As such, Web-based simulations have been profusely used in other types of simulation. Instances include the simulation of the glass-pressing process developed by Shi, Zhang, Rao, and Zhang (2008) uses a Web-based client approach that returns a visualisation through the Web browser. The ecological modelling framework EcoNet (Kazanci, 2007) is a noninteractive Web-based simulator which implements a model of population dynamics, although this model did not incorporate a spatial component, instead tracking the change within a single population and producing a plot of population change over time as output. Another example is INDISIM-YEAST, which also chose this approach for modelling yeast populations in a liquid (Ginovart & Can˜adas, 2008); a subject where the spatial component is less relevant due to the small physical scale. The closest effort to ours in terms of technical approach is probably by Xie and Yapa (2006). This is a simulation model of oil spills, which is Web-based using a clientserver architecture, includes a spatial component, and allows interactive, visual editing of the simulation through a Java applet. However this simulator does not incorporate a GIS component.
3.
Simulation model
The main simulation model has four interrelated models that include: i) a spatial model that represents the simulation landscape as a grid of cells; ii) a flora model that represents a patch with a number of flowers and seasons defined by nectar quality and quantity; iii) a weather model that considers time periods of temperature and daily hours of light; iv) a foraging model, which represents the foraging behaviour of a honeybee colony. Figure 1 corresponds to the general conceptual diagram that shows how the main simulation model and its underlying models are connected and what key concepts are shared among them. Our approach does not model the population dynamics of a beehive or the interactions among its individuals, in contrast with other existing work where this individual-based exchange of information leads to a globally intelligent selection of resources (Quijano & Passino, 2010; Tereshko & Loengarov, 2005; Zhang, Ouyang, & Ning, 2010). This is not only for the sake of simplicity, but also because it is not required given the focus and applicability of the simulator. This simulation model is focused on explaining the effects (e.g. total amounts of forage by beehive, visitation rates, foraging distances, etc.) of honeybee foraging in the flora surrounding the beehive on a time scale of days, which can be modelled by using published empirical data as Section 3.4 describes. Honeybees perform a foraging cycle (flying to
Fig. 1 e Conceptual diagram of the simulation model and its related models.
flowers, loading nectar, returning to the beehive, unloading, dancing, etc.) on a time scale of minutes, while recruiting other foragers occurs on a time scale of hours. The consequence from the mathematical point of view is the proportion of bees in each phase of the foraging cycle is essentially at equilibrium on the recruitment time scale (Cox & Myerscough, 2003). Therefore, it is only essential to consider the total number of active foragers at a beehive, based on the total population of its colony. Because traditional mathematical models that provide analytical solutions are only possible for simplified environments, we employ simulation modelling where the abstraction of the particular entity to predict is solved by steps numerically (Grant & Swannack, 2008). The simulation model is based on empirical data and expert knowledge in the field of honeybee ecology and management. Simulation models are often divided between deterministic and stochastic (Odum, 1971). Deterministic simulation models always produce the same result for a given set of input data and model parameters, whereas stochastic ones include one or more random processes, usually as a way of dealing with uncertainty (Law & Kelton, 1991). While the flora model (see Section 3.1), the weather model (see Section 3.2), and the spatial model (see Section 3.3) are deterministic, the foraging model (see Section 3.4) is stochastic. The simulation model specification formalism is DTSS (Discrete Time System Specification) (Nutaro, 2007), where the evolution of the simulation process is approximated by a difference equation. With d representing the duration of the simulation process (e.g. in days), we define Dd as the time-step used to partition d resulting in the first order difference equation yiþ1 ¼ yi þ Dd f ðyi ; xi Þ, which predicts the value of y at time step i þ 1 based on information from its previous time step i. In our particular scenario, Dd corresponds to one foraging trip t, with several of these trips taking place within each simulation day. More specifically, the number of foraging trips ti in a given day di with 0 i < d are calculated by using the
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
average temperature (in C) during the day Ti, the number of daylight hours hi, and the currently available nectar and pollen resources in the surrounding flower patches. Thus, the ¼ m hi , with m ¼ 0.8 maximum number of foraging trips tmax i when 15 < T < 33, m ¼ 0 when T < 8, or m ¼ 0.4 otherwise, while if the the eventual number of foraging trips will be ti tmax i daily nectar is depleted before the end of the day. Figure 2 provides the sequence of tasks that a typical simulation process would consist of. The spatial simulation parameters specified by the user through the Web-based interface, such as geographical bounds and resolution, are used to create the grid representing the simulation area by the spatial model described in Section 3.3. At the beginning of each simulation day di the appropriate average temperature Ti and hours of light hi are selected by applying the model . Also, the described in Section 3.2 and used to calculate tmax i current daily season of each flower patch is chosen using the model described in Section 3.1, and the nectar and pollen production for the day in each cell estimated. For each foraging trip in the day, the foraging model detailed in Section 3.4 is used to describe the honeybee foraging behaviour. At the end of a simulation day, a set of daily data is collected in anticipation for the final results report to be presented to the user through the Web interface. The report is organised by sections, each one corresponding to one simulated day. For each day, the user can inspect: i) how weather has changed in the simulation area; ii) how many foraging trips took place; and iii) technical information that includes: a) the number of cells of the grid representing the simulation area have been used by day; b) the length of time taken for the simulator to process this day; and c) how much RAM the server used. For each foraging trip in a day, the report includes a section with information about: i) the number of different foraging destinations (i.e. cells); ii) the amount of collected nectar; and iii) the amount of collected pollen. Each foraging
307
destination at any time can be examined by the user, hence obtaining a graphical view of visitation rates by honeybees in addition to nectar and pollen production and collection. Moreover, these results can be exported using the wellsupported comma-separated values (CSV) format so that it can be easily imported into popular spreadsheet or data analysis applications.
3.1.
Flora model
A flower patch corresponds to an area containing the following information: total number of flowers in the area and one or more sequential, cyclic seasons. A season indicates for how many days these flowers have the same average production of nectar and pollen per flower. There is an important distinction to make between gross and net nectar. Flowers in a patch may produce a certain amount of (gross) nectar per day. However, this nectar will have a certain ‘quality’ based on its sugar concentration and its subsequent energetic value. While bees will be able to carry a certain load of nectar regardless of its quality, they will prefer to forage flowers that produce higher quality nectar (Hainsworth, 1974). A traditional finite state machine (FSM) is used to represent a sequence of seasons corresponding to the same flower patch. A FSM is a mathematical abstraction useful to model behaviour explained by a finite number of states, transitions between them, and actions. A formal FSM can be defined by a five-tuple M ¼ hI; S; s0 ; O; di, where I is the set of inputs which are used to transition the machine among one of a finite set of states S with s0 ˛S indicating the initial state, O is the set of possible outputs, and d is the transition function d : S I/S. Hence, I O represents the alphabet of M. A state machine will switch to a different state sm ˛S after reading an input in ˛I. Each of the machine states s˛S specifies which state to switch for a given input i˛I. From any state there is only one
Fig. 2 e Overview of tasks performed by the main simulation model.
308
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
outputs O ¼ fh200; 0:7; 20i; h100; 0:9; 30i; h50; 1; 10i; h0; 0; 0ig, and inputs I ¼ f1; 2; .; 66g.
3.2.
Fig. 3 e A finite state machine with five states s0 ; .; s5 , input alphabet I[f0; 1g and outputs O[fl; rg.
transition for any allowed input. Figure 3 shows an example of a FSM known as transducer (or Mealy machine), where outputs depend on both the state and the input, whereas in Moore machines an output is determined only by the state. In this example, when the machine is at state s0 and receives 0 as an input, it will transition to state s1 and produce l as output. However, the state will not change if the received input is 1 while the output value will be r. Considering the set of flower patches in the landscape F ¼ ff0 ; f1 ; .; fm g, each flower patch fi ¼ hli ; Ei i; 0 i m is represented by the total number of flowers in it li and the list of flower seasons Ei ¼ fe0 ; e1 ; .; ek g. Each season ej ; 0 j k is described by a four-tuple hdj ; nj ; qj ; pj i where dj expresses the number of days that this season ej lasts, nj indicates the average amount of nectar production per day per flower (in mg) during this season, qj ˛½0; 1 corresponds to the nectar quality rate value, and Pj indicates the average amount of pollen production per day per flower (in mg) during the season. The FSM representing the seasons of a flower patch has states S ¼ fs0 ; s1 ; .; sk g to represent E, hence jEj ¼ jSj. The input alphabet includes all the consecutive day numbers for E, P thus consisting of I ¼ f1; 2; .; dg, where d ¼ 1 þ ki¼0 di . The output alphabet is given by all the three-tuples extracted from E, producing O ¼ fhn0 ; q0 ; p0 i; hn1 ; q1 ; p1 i; .; hnk ; qk ; pk ig. The transition from one machine state to the next occurs when the input for the current state is the last day of the season represented. A transition from the last state (i.e. flowers season) to the first state is automatically created to represent the cyclic nature of seasons. This transition function has to be triggered on each simulation day for all flower areas. Figure 4 shows the machine representing the four flower patch seasons h10; 200; 0:7; 20i, h5; 100; 0:9; 30i, h20; 50; 1; 10i, and h30; 0; 0; 0i. This means a total natural cycle of 65 days, thus t ¼ 66. This corresponds to four states S ¼ fs0 ; s1 ; s2 ; s4 g, with
Weather model
The simulator uses weather information that is crucial for determining the number of daily foraging trips. This information includes temperature (in C) and number of daylight hours, with both measures applying to the whole simulation area. Thus, given a list of weather periods W ¼ fw0 ; w1 ; .; wk g, an arbitrary weather period wi, with 0 i k, is expressed as a list of three-tuples, each one specifying a number of days along with the temperature and the number of daylight hours during them such as hdi ; Ti ; hi i. The sum of the number of days P for all these pairs d ¼ ki¼0 di corresponds to the total number of days of simulation. Like the Flora Model described in Section 3.1, this model is based on a FSM. In this case, each state si ˛S represents a number of consecutive days di with the same temperature Ti and daylight hours hi. Hence the number of states jSj is determined by the number of temperature periods. The input alphabet includes all the added period days I ¼ f0; 1; 2; .; dg. The output oi ˛O that corresponds to si is the temperature and daylight hours for that period hTi ; hi i. The transition from one state to the next occurs when the input for the current state corresponds to the last day of the period it represents.
3.3.
Spatial model
This section describes the spatial model, which provides a two dimensional abstraction of the simulation area and related entities, including the flower patches, beehives, and honeybee foraging activity. This spatial model is based on the one described in Garcı´a Adeva et al. (2012), which is here extended by adding data to cells in order to track particular aspects of the honeybee foraging activity. In this model, the simulation area consists of a landscape whose bounds are initially defined by the user through the simulator’s user interface. The extent will usually vary from a few hectares to several hundred thousand km2. The landscape can contain areas (i.e. polygons) representing floral patches and locations indicating beehives. Figure 5 provides a graphical representation of this landscape expressed as a rectangular planar region L of width wL and height hL . The boundaries of L are given by two global coordinates ðcbl ; ctr Þ ¼ ððlatb ; lonl Þ; ðlatt ; lonr ÞÞ corresponding to the bottom-left and top-right corners of the extent. To calculate the distance d between two coordinates c1 ¼ ðlat1 ; lon1 Þ and c2 ¼ ðlat2 ; lon2 Þ in L, we define the following equation based on the Spherical Law of Cosines: rðc1 ; c2 Þ ¼ arccos ðsin lat1 sin lat2 þ cos lat1 cos lat2 cos ðlon1 lon2 ÞÞ<
Fig. 4 e The FSM representing the 65-day cycle of a flower patch with four seasons h10; 200; 0:7; 20i, h5; 100; 0:9; 30i, h20; 50; 1; 10i, and h30; 0; 0; 0i.
where < ¼ 6371 is a constant defining the Earth’s radius in km. Discretisation of L is performed by means of a rectangular sparse grid G with its maximum number of horizontal and vertical cells being wG hG . A grid cell is defined by g ¼ ðx; yÞ˛G, where x; y˛Z; 0 x < hG ; and 0 y < hG , with only those non-empty cells being created dynamically when
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
309
Fig. 5 e An extent of landscape L discretised by grid G.
needed. Therefore, wG and hG also indicate the simulation resolution, with its maximum value limited by such a cell size where any given cell contains at least 1 flower. The translation from a coordinate c ¼ ðlat; lonÞ˛L to its cell position g ¼ ðx; yÞ˛G is calculated with
km2, we apply Meister’s formula that calculates the area of a non-self-intersecting polygon using Euclidean geometry: xi ¼ < cos lati cos loni ; yi ¼ < cos lati sin loni ;
x ¼ Prh ðcÞwG wL R; y ¼ Prv ðcÞhG hL R where rv ðcÞ and rh ðcÞ indicate vertical and horizontal distances respectively from c to cbl in L: rh ðc1 Þ ¼ rðcml ; ðlatm ; lon1 ÞÞ; rv ðc1 Þ ¼ rðcbm ; ðlat1 ; lonm ÞÞ: This discretisation process can almost be considered a map cylindrical projection, where the sphere is projected on the walls of a cylinder by mapping circles of longitude (i.e. meridians) as equidistant vertical straight lines and circles of latitude as horizontal straight lines. One of the key differences is how in the traditional cylindrical projection, the projection consists of a cylinder wrapped around the globe at the equator, whereas in this case the cylinder is centred around the midpoints of L. The immediate benefit of this approach is a significant reduction in distortion, hence obtaining a more precise spatial representation. All grid cells are square in shape and have the same size sG ¼ l2 , where sG , wG , and hG are calculated based on the horizontal resolution selected by the simulator’s user (typical values will be in the range 50e500). Thus, obtaining hG and sG in km2 is achieved by
hG ¼
Q S
hL wG wL ;
AðV L Þ ¼
wL hL : wG hG
An area representing a flower patch is expressed as a closed polygon defined by the list of its ordered (either clockwise or counterclockwise) vertex coordinates V L ¼ ðv0 ; v1 ; v2 ; .; vk ; v0 Þ with vi ðlati ; loni Þ˛L; 0 i k. To calculate the extent of V L , in
k1 1X xiþ1 ¼ xi yiþ1 xiþ1 yi : yiþ1 2 i¼0
Discretising an area V L to obtain its corresponding V G simply consists of translating each vertex vi ˛L into corresponding cell positions vi ˛G. Another instance of a frequently required operation is to determine whether a particular cell is included within an area, which we achieve by applying the Jordan curve theorem as proposed by Haines (1994). A cell may contain one or more beehives. In addition, a cell keeps track of how much nectar and pollen it contains by taking into account how many flower patches it is contained by and the season that currently applies to each of these flower patches. In most cases, a cell will be within a single flower patch, although the opposite is possible in cases where flower patches overlap. Hence, given an arbitrary cell c within flower patches F ¼ ff0 ; f1 ; .; fm g, where a single flower patch is represented by fi ¼ hli ; Ei i; 0 i m, li indicates the total number of flowers in it and Ei ¼ fe0 ; e1 ; .; ek g the sequence of seasons, the total amount of forage in c is calculated by lc ¼
m1 X AðcÞ li ; A fi i¼0
nc ¼
m1 X
lc ni ;
i¼0
qc ¼
m1 X nc qc þ ni qi i¼0
sG ¼
k1 1X xi 2 i¼0 yi
pc ¼
m1 X
nc þ ni
;
lc pi ;
i¼0
where lc corresponds to the average number of flowers per cell c, nc corresponds to nectar and pc to pollen, while qc indicates the weighted nectar quality rate that applies to nc.
310
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
3.4.
Foraging model
A key feature of honeybee foraging is how this activity extends over an area larger than 100 km2 around the beehive. The foraging behaviour of honeybees living in nature is characterised by daily changes to their routine and distribution across flower patches and regular foraging at food sources located several km from their beehive (Seeley, 1985). From a given population of forager bees at a beehive, this model aims to determine: i) how many foragers will visit particular flowers depending on proximity, nectar quality, nectar quantity, and pollen quantity; and ii) the number of foraging trips for a particular day. This is a stochastic model due to some of its actions being supported by to the use of weighted random rounding. This stochastic rounding operation happens in multiple aspects of the model including the calculation of foraging population including the different types such as nectar foragers, pollen foragers, and nectar and pollen foragers, the distribution of foragers among cells containing flowers, etc. Thus, we define PxS as the stochastic, weighted rounding operation applied to x, which is based on a Bernoulli random variate, such as PxS ¼
PxR QxS
if ε < fxg if ε fxg
(1)
where x is a positive real number x˛R0 , PxR is the floor function defined by x ¼ max fm˛Zm xg, PxR is the ceiling function defined by x ¼ fn˛Zn xg, fxg is the fractional part of x as given by fxg ¼ x PxR, and ε is the random function with ε : R0 /½0; 1. There is a kernel that defines both the foraging distance and the number of individuals finding a flower patch at a certain distance, by describing the spatial distribution of foragers in the honeybee vicinity as a function of distance from the beehive and the size of flower patches. Foraging is similar to the concept of local dispersal frequently found in insect modelling. There are many kernelbased dispersal models found in the literature including Laplace, Gaussian, Ribbens (Clark, Macklin, & Wood, 1998), and Cauchy (Shaw, 1995). They may be considered arbitrary equations that fit particular dispersal data sets well. The use of these kernels requires that previously parameters for the chosen distribution to be estimated by registering the
a
dispersal distances between occurrences in distribution data from one time period to the next, so that the distribution can be fitted to the frequency distribution of distances (Pitt, 2008). Our dispersal kernel, however, uses an empirical probability distribution to determine foraging distances based on published data (Seeley, 1985; Visscher & Seeley, 1982) and depicted by Fig. 6(a). For instance, 7% of foraging bees will do so within ð0:2; 0:4 km, 28% will do it within ð0:4; 1 km, and so forth. This distribution is based on the assumption that the beehive is surrounded by flower patches rich in forage (i.e. nectar and pollen) and accounts for related activities such as travel times between the beehive and the flower patches, recruitment of other foragers by dancing, etc. Generally, most foragers concentrate their activity in a radius of 2 km from their beehive, but they will travel up to 3 times this distance if needed for survival (Beekman & Lew, 2008). We consider that 30% of the total colony population will be actively foraging during the daylight hours, with 58% of them foraging only nectar, 25% only pollen, and the remaining 17% both pollen and nectar (Winston, 1991). Figure 6(b) provides the graphical representation of the probability of suitable flower patches being detected by nearby foragers based on distance and size. This is based on the assumption that a flower patch will be more detectable by foragers the larger or the closer it is, and vice versa. Therefore, in order to express this probability of detection, we use the probability distribution given by Aðf Þ r 1 g; 1 ; Pf ¼ min 1000 10 where Aðf Þ˛½0; 1000 indicates the flower patch size in m2, r ˛½0; 10 is the distance in km between the flower patch and the beehive, and g > 0 corresponds to the scale parameter, where the larger its value the higher the value of Pf. These two distributions are combined in order to determine the weighted probability that each cell within a flower patch is visited by honeybees. The algorithmic approach for performing this process is based on the grid G, where a cell gi ¼ ðxi ; yi Þ˛G in this grid represents a discrete location in the simulation area, contained by flower patches F ¼ ff0 ; f1 ; .; fm g. The foraging bees from a beehive are represented by B ¼ fb0 ; b1 ; .; bn g.
b
Fig. 6 e Probability distributions of the dispersal kernel used by the foraging model.
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
Table 1 e Basic scenarios and their result to verify the simulation model. Scenario No flower patches around the beehive A flower patch close to the beehive with no nectar or pollen A flower patch rich in nectar and pollen at 20 km from the beehive Temperature below 8 C No daylight hours Beehive has only 1 honeybee
Result No forage activity No forage activity
No forage activity
No forage activity No forage activity Foraging activity 1 only in about of the times 3
By traversing all the possible dispersal distance ranges shown by Fig. 6(a) such as rba ; rcb ; rdc ; . along with the minimum possible distance increment Dr in G (i.e. the length of a cell side) for as many times as possible (i.e. maxfn xi þ 1; xi ; m yi þ 1; yi g), we obtain a list of cells C ¼ ðc0 ; c1 ; .; ck Þ belonging to the ring that corresponds to the first foraging distances interval ½rba ; rcb Þ. This ring of cells is created by applying the well-known midpoint algorithm (van
311
Aken, 1984), thus creating circles of cells as long as rba r < rcb . Those cells from C that are not contained by a flower patch in F are discarded. Also discarded are those cells for which applying the stochastic weighted rounding function defined by Eq. (1) to their corresponding probability of being found Pf results in 0, hence producing the effect that they could not be found by foraging bees due to their location. The remaining cells of C are ordered according to their Pf, to then be used to sequentially distribute the bees from B. The number of bees from B assigned to each remaining cell in C is proportional to the relative amount of forage in the cell. The type of forage considered for this calculation depends on the type of forager: nectar foragers will use the remaining net nectar in the cell, pollen foragers will use the remaining pollen in the cell, while nectar and pollen foragers will use both net nectar and pollen in the cell. The forage in those cells has to be depleted accordingly. As default values, we selected an average load of 40 mg of gross nectar and 20 mg of pollen per bee per trip (Beekman & Lew, 2008; Feuerbacher, Fewell, Roberts, Smith, & Harrison, 2003; Winston, 1991), although these values can be changed by the user through the configuration section of the Web interface. With these values, the depletion process consists of simply removing from the remaining reserve of cell forage the average load of nectar, pollen, or both as many times as
a
b
c
d
Fig. 7 e Location of foraging honeybees on day 1 of the simulation. Dots represent beehives and squares represent a groups of foragers, with the darker the colour the higher the population.
312
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
honeybees chose to forage in the cell such with nc ¼ nc min ðnc ; 40jBjÞ, and pc ¼ pc min ðpc ; 20jBjÞ, while qc remains unchanged.
4.
Evaluation
Evaluating simulation systems a functional point of view is important in order to ensure that the right system has been correctly built (Sojda, 2007). However, because simulation
systems handle complex and poorly structured problems, they are difficult to empirically evaluate. Accordingly, a frequently chosen approach to tackle this difficulty when evaluating a simulation model consists of dividing the effort into two separate tasks: verification and validation (Mihram, 1972). Verification consists of showing the simulation models and software implementation are correct, complete, and coherent, while validation consists of comparing simulation results with observations from the real world. Verification should precede validation (Mihram, 1972).
a
b
c
d
e
f
Fig. 8 e Nectar and pollen depletion after four foraging trips on day 1. Dots represent beehives and squares represent a flower patch section, with the darker the colour the higher the amount of nectar or pollen left.
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
The ideal method of validation would consist of comparing results with a preselected gold standard d a dataset with sufficient accuracy and precision available to the research community. For ecological simulation this gold standard would usually be high-quality historical records. However, it is rare to find one in the natural sciences (Sojda, 2007) and indeed it was the case here d we were not able to find a dataset describing a scenario in a way that can be used to propose a validation strategy. In consequence, our evaluation was limited to verifying the simulation model, by means of
313
performing simulation experiments of a selected set of use cases and then comparing results with expectations. We started the verification process by proposing several basic simulation scenarios that consider extreme use cases. The expected result can be guessed by common sense. Most of the use cases summarised in Table 1 correspond to situations where no foraging activity should be expected, because the temperature is too cold, the flower patch within reach but with no forage, the flower patch has plenty of forage available but it is too far, etc. Although some of these
a
b
c
d
e
f
Fig. 9 e Nectar and pollen depletion after four foraging trips on day 11. Dots represent beehives and squares represent a flower patch section, with the darker the colour the higher the amount of nectar or pollen left.
314
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
a
b
c
d
Fig. 10 e Nectar and pollen depletion after four foraging trips on day 16. Dots represent beehives and squares represent a flower patch section, with the darker the colour the higher the amount of nectar or pollen left.
scenarios were obviously unrealistic (e.g. no daylight hours or a single bee in the beehive), they were still important and useful as a test mechanism to verify that the simulator worked as expected. Next, we proposed additional verification based on a realistic scenario.1 The data generated by running simulations on this scenario and used to produce the graphs shown below is fully available through the Web user interface. The scenario in question comprised a simulation area of 166.28 km2 (16.35 km by 10.17 km) represented by a spatial resolution of 100 by 63 cells (i.e. 0.026 km2 or 2.6 Ha per cell). There are two weather periods: i) 15 days at 22 C with 12 daylight hours; and ii) 15 days at 15 C with 10 daylight hours. There are two beehives, one small with a population of 10,000 bees and a larger one with 20,000 bees. Between the two beehives there is a medium flower patch of about 5 km2 with 500 flowers that produce an average of 40 mg of nectar and 10 mg of pollen per flower per day during the first 15 days, to then reduce their production to 5 mg of nectar (quality rate 0.8) and 2 mg of pollen for the remaining 15 days. East to the larger beehive there is a larger flower patch covering an area of 1
Scenario data available at http://beetle.csse.uwa.edu.au/ honeybee/?bioeng.
21 km2 and containing 5000 flowers. During the first 10 days they produce 50 mg (quality rate 0.5) of nectar and 20 mg of pollen on average per day per flower. During the remaining 20 days the production drops to 10 mg of nectar (quality rate 0.2) and 5 of pollen. The purpose was to use this scenario to examine different types of results produced by the simulation in order to verify the following aspects and effects of the foraging activity: i) foraging bee location and activity over the flower patches; ii) nectar and pollen depletion along several days where some key events took place (e.g. change in weather or in flower production capability); and iii) amount of forage collected by day and during the total duration of the simulation. Figure 7 provides a graphical view of the whereabouts of honeybees on day 1 while performing several foraging trips to the two flower patches, which can be interpreted as an indication of foraging activity along the day. During this day, the maximum number of foraging bees from the small beehive was around 3000. In addition, the maximum expected number from the large beehive was 6000. This accounts for a combined 5220 honeybees foraging nectar, 2250 foraging pollen, and 1530 foraging both nectar and pollen. Foraging bees will prefer to start their work in those parts of the flower patches that are close to the beehive, and this can be verified to happen in
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
Fig. 8(a) and (b). However, as they deplete the nectar and pollen reserves available for the day, these foragers will have to move farther from the beehive, consequently making their work harder since they will have to fly longer distances. Figure 8(c) and (d) confirmed that this is the case by the 5th foraging trip, with it becoming even more prominent by looking at Fig. 8(e) and (f), which correspond to the last foraging trip for the day. These figures also show how the amount of daily resources produced by the smaller flower patch was not enough to provide forage for the whole day, so that by the 5th foraging trip the number of honeybees foraging in this flower patch was very small, while by the end of the day there was virtually no foraging activity left on this flower patch. In contrast, the larger flower patch still offered food based on the small groups of honeybees still foraging on it. However, the honeybees did not continue foraging due to the lack of daylight. Figure 8 offers a graphical view of how nectar and pollen were depleted from the flower patches on day 1 of the simulation. It is interesting to compare these results to those on day 11 (Fig. 9), when the larger flower patch had a season change and its nectar production dropped significantly, and also to those on day 16 (Fig. 10), when there was a change in temperature (thus becoming colder) and number of daylight
315
hours thus producing a reduction in honeybee foraging activity, while at the same time the smaller flower patch had a drop in its nectar production due to seasonal change. Figure 11 shows the remaining nectar at the end of several foraging trips on day 21, where the situation with regards to weather and flower patch seasonality was about the same as on day 16 and hence the results quite similar. As a consequence of weather and flower patch seasonality the number of daily foraging trips between day 1 and day 15 was 9e10, whereas from day 16 to day 30 it was only 4. The last part of this verification focused on amount of collected forage. Figure 12 shows the daily and total amount of collected nectar and pollen by both beehives. Figure 12(a) and (b) show evidence of how the change in weather after day 10 produced a significant reduction in both the daily and total foraged nectar and pollen. Something similar, but with a much lower impact, happened after day 15 as a consequence of the weather change at the time. This effect could also be observed in the total amount of nectar during the 30 days of simulation by looking at Fig. 12(c) and (d). It is important to note that the frequent peaks in Fig. 12(a) and (b) were due to the declining amount of forage collected by the bees as the day progressed. Therefore, the amount of both nectar and pollen in the last foraging trips of the day were
a
b
c
d
Fig. 11 e Nectar and pollen depletion after four foraging trips on day 21. Dots represent beehives and squares represent a flower patch section, with the darker the colour the higher the amount of nectar or pollen left.
316
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
a
b
c
d
Fig. 12 e Daily and total amount of nectar and pollen foraged.
usually smaller than during the first foraging trips of the day, due to availability.
5.
Conclusions
The importance of the ecological role of honeybees by contributing to biodiversity and the production of crops and wild plants through pollination is well understood. Their honey production is also an important economic factor. Even though simulation models can improve our understanding of honeybees, hence helping stakeholders such as biosecurity managers and beekeepers, there are few models available and many are severely limited in their scope. We tackle this by offering a simulation tool that supports long time periods, is highly configurable, offers an interactive user interface, and
produces very detailed results. Its simulation model is stochastic and discrete in time and space, and can support any number of beehives and floral patches over a two dimensional landscape of unrestricted size and resolution. The different underlying models encapsulate well-established ecological knowledge. The verification of the simulation model included several experiments that dealt with nectar and pollen foraging under different conditions of weather and flower patch seasonality. Although it would have been useful to complement this evaluation by validating the simulation results against an existing dataset, the lack of such a dataset made this impossible. We note this as possible future work along with an analysis of the engineering merits of the simulator, including an evaluation of system measures such as performance, time response, and scalability.
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
Acknowledgements The author would like to thank Lori Lach from the Centre for Integrative Bee Research at the University of Western Australia and Rob Manning from the Department of Agriculture and Food Western Australia for pointing at relevant bibliography and providing feedback. The author would like to acknowledge the support of the Australian Government’s Cooperative Research Centre Program.
references
van Aken, J. (1984). An efficient ellipse-drawing algorithm. IEEE Computer Graphics and Applications, 4, 24e35. Bailis, P., Nagpal, R., & Werfel, J. (2010). Positional communication and private information in honeybee foraging models. In Proceedings of the 7th international conference on Swarm intelligence (pp. 263e274). Berlin, Heidelberg: Springer-Verlag. Beekman, M., & Lew, J. B. (2008). Foraging in honeybeesewhen does it pay to dance? Behavioral Ecology, 19(2), 255e261. Byrne, J., Heavey, C., & Byrne, P. J. (2010). A review of Web-based simulation and supporting tools. Simulation Modelling Practice and Theory, 18(3), 253e276. Camazine, S., & Sneyd, J. (1991). A model of collective nectar source selection by honey bees: self-organization through simple rules. Journal of Theoretical Biology, 149(4), 547e571. Champetier, A. (2011). The foraging economics of honey bees in almonds. In Annual meeting of the Agricultural and Applied Economics Association, Pittsburgh, PA. Clark, J. S., Macklin, E., & Wood, L. (1998). Stages and spatial scales of recruitment limitation in southern Appalachian forests. Ecological Monographs, 68(2), 213e235. Cox, M., & Myerscough, M. (2003). A flexible model of foraging by a honey bee colony: the effects of individual behaviour on foraging success. Journal of Theoretical Biology, 223, 179e197. Crozier, R. H., & Page, R. E. (1985). On being the right size: male contributions and multiple mating in social Hymenoptera. Behavioral Ecology and Sociobiology, 18(2), 105e115. Feuerbacher, E., Fewell, J., Roberts, S., Smith, E., & Harrison, J. (2003). Effects of load type (pollen or nectar) and load mass on hovering metabolic rate and mechanical power output in the honey bee Apis mellifera. Journal of Experimental Biology, 206(Pt 11), 1855e1865. von Frisch, K. (1967). The dance language and orientation of bees. Cambridge, MA: Belknap Press of Harvard University Press. Garcı´a Adeva, J. J., Botha, J. H., & Reynolds, M. (2012). A simulation modelling approach to forecast establishment and spread of Bactrocera fruit flies. Ecological Modelling, 227(0), 93e108. Garcı´a Adeva, J. J., Lach, L., & Reynolds, M. (2011). Simulation of honeybee nectar foraging for determining effects on local flora. In 19th International congress on modelling and simulation, the Modelling and Simulation Society of Australia and New Zealand conference, (pp. 2507e2513). Perth, Australia. Garcı´a Adeva, J. J., Sousa-Majer, M. J., Botha, J. H., Hanbury, C. D., Hardie, D. C., & Reynolds, M. (2010). Modelling the establishment and spread of emergency plant pests (EPPs) in Australia: Simulate or suffer. In Global biosecurity conference (p. 101). Brisbane, Australia. Ghosh, S., & Marshall, I. (2005). Simple model of collective decision making during nectar source selection by honey bees.
317
In Workshop on memory and learning mechanisms in autonomous robots (ECAL 2005) p. 10. Ginovart, M., & Can˜adas, J. (2008). INDISIM-YEAST: an individualbased simulator on a website for experimenting and investigating diverse dynamics of yeast populations in liquid media. Journal of Industrial Microbiology & Biotechnology, 35(11), 1359e1366. doi:10.1007/s10295-008-0436-4. Grant, W. E., & Swannack, T. M. (2008). Ecological modeling: A common-sense approach to theory and practice. Wiley-Blackwell. Haines, E. (1994). Graphics Gems IV. Chapter Point in Polygon Strategies. Morgan Kaufmann. Hainsworth, F. R. (1974). Food quality and foraging efficiency. Journal of Comparative Physiology A: Neuroethology, Sensory, Neural, and Behavioral Physiology, 88, 425e431. Kazanci, C. (2007). EcoNet: a new software for ecological modeling, simulation and network analysis. Ecological Modelling, 208(1), 3e8. Keller, L., & Reeve, H. K. (1994). Genetic variability, queen number, and polyandry in social Hymenoptera. Evolution, 48(3), 694e704. Law, A. M., & Kelton, W. D. (1991). Simulation modeling and analysis (2nd ed.). New York, NY: McGraw-Hill. Mihram, G. A. (1972). Some practical aspects of the verification and validation of simulation models. Operational Research Quarterly, 23(1). Moosbeckhofer, R., Loncaric, I., Oberlerchner, J., Persen, U., Ertl, C., & Donat, C. (2008). Use of honeybees (Apis mellifera) as vectors for fire blight antagonists in field experiments. In Acta horticulturae 793: XI International workshop on fire blight (pp. 461e465). Morse, R., & Calderone, N. (2000). The value of honey bees as pollinators of U.S. crops in 2000. Bee Culture Magazine, 128, 1e14. Nutaro, J. (2007). Discrete event simulation of continuous systems. In Handbook of dynamic systems modeling. Chapman & Hall/CRC Press. Odum, E. P. (1971). Fundamentals of ecology. Philadelphia: W.B. Saunders Company. Pitt, J. P. W. (2008). Modelling the spread of invasive species across heterogeneous landscapes. PhD thesis, Lincoln University. Potts, S., Biesmeijer, J., Kremen, C., Neumann, P., Schweiger, O., & Kunin, W. (2010). Global pollinator declines: trends, impacts and drivers. Trends in Ecology and Evolution, 25(6), 345e353. Quijano, N., & Passino, K. M. (2010). Honey bee social foraging algorithms for resource allocation: theory and application. Engineering Applications of Artificial Intelligence, 23(6), 845e861. Seeley, T. D. (1985). Honeybee ecology: A study of adaptation in social life. Princeton, NJ: Princeton University Press. Seeley, T. D. (1995). The wisdom of the hive: The social physiology of honey bee colonies. Cambridge, MA: Harvard University Press. Seeley, T. D., Camazine, S., & Sneyd, J. (1991). Collective decision-making in honey bees: how colonies choose among nectar sources. Behavioral Ecology and Sociobiology, 28(4), 277e290. Shaw, M. W. (1995). Simulation of population expansion and spatial pattern when individual dispersal distributions do not decline exponentially with distance. Proceedings of the Royal Society London Series B, 259, 243e248. Shi, S., Zhang, G., Rao, Y., & Zhang, Y. (2008). A Web-based remote simulation system for the glass pressing process. Simulation, 84(1), 27e40. Sojda, R. (2007). Empirical evaluation of decision support systems: needs, definitions, potential methods, and an example pertaining to waterfowl management. Environmental Modelling and Software, 22(2), 269e277.
318
b i o s y s t e m s e n g i n e e r i n g 1 1 2 ( 2 0 1 2 ) 3 0 4 e3 1 8
Tereshko, V., & Loengarov, A. (2005). Collective decision making in honey-bee foraging dynamics. Computing and Information Systems, 9(3), 1e7. Visscher, P. K., & Seeley, T. D. (1982). Foraging strategy of honeybee colonies in a temperate deciduous forest. Ecology, 63(6), 1790e1801. de Vries, H., & Biesmeijer, J. C. (1998). Modelling collective foraging by means of individual behaviour rules in honeybees. Behavioral Ecology and Sociobiology, 44, 109e124.
Winston, M. L. (1991). The biology of the honey bee. Cambridge, MA: Harvard University Press. Xie, H., & Yapa, P. D. (2006). Developing a Web-based system for large-scale environmental hydraulics problems with application to oil spill modeling. Journal of Computing in Civil Engineering, 20(3), 197e209. Zhang, C., Ouyang, D., & Ning, J. (2010). An artificial bee colony approach for clustering. Expert Systems with Applications, 37(7), 4761e4767.