{'fraction_non_alphanumeric': 0.035359899196723894, 'fraction_numerical': 0.016734918884863757, 'mean_word_length': 4.122428398547801, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A marked renewal process model for the size of a honey bee colony

Many areas of agriculture rely on honey bees to provide pollination services and any decline in honey bee numbers can impact on global food security. In order to understand the dynamics of honey bee colonies we present a discrete time marked renewal process model for the size of a colony. We demonstrate that under mild conditions this attains a stationary distribution that depends on the distribution of the numbers of eggs per batch, the probability an egg hatches and the distributions of the times between batches and bee lifetime. This allows an analytic examination of the effect of changing these quantities. We then extend this model to cyclic annual effects where for example the numbers of eggs per batch and {the probability an egg hatches} may vary over the year.


Introduction
The world population reached 7.3 billion as of mid-2015, indicating that the world has added approximately one billion people in the last twelve years and is projected to pass 9 billion by 2050 (UN, 2015). It is vital, therefore, that optimal use is made of the world's finite resources for food production (Collier, 2009). Whilst honey bees are not the only pollinators (Rader et al., 2016), it is estimated that 70% of important food crops are bee-pollinated (Datta et al., 2013) and so the status of pollinators, and bees in particular, is of a key concern in global food security (Lonsdorf et al., 2009;Blaauw and Isaacs, 2014). Many large-scale agricultural businesses now employ peripatetic bee services in order to ensure adequate pollination of crops (Gordon et al., 2014;Bishop et al., 2016).
There is evidence that pollinator populations are declining and this is a cause for concern for policy-makers (Hafi et al., 2012;Breeze et al., 2012;Vanbergen et al., 2014;Chauzat et al., 2014;Potts et al., 2009). This observation motivates our development of a probabilistic model for the colony size. The model allows us to examine how the size depends on the distribution of the numbers of eggs per batch, the chance of hatching and the distributions of the times between batches and bee lifetime. In turn this gives insight into how a disease for example that affects one or more of these quantities affects the colony size.
The dynamics of a bee colony are complex. After her mating flight, the queen remains in the hive laying eggs, fed and cared for by worker bees. The number of eggs laid are influenced by nectar flow into the colony, which in turn, is influenced by the number of bees out foraging and seasonal effects on availability. During winter, in the absence of forage, there is little or no egglaying and the colony relies on the honey and pollen stores accumulated during the summer. In this situation the lifetime of a bee can be significantly longer so that the colony can survive winter (Mattila et al., 2001). Eggs hatch into larvae which pupate after a few days of feeding (Collins, 2004) and emerge as juvenile adult bees about three weeks after the egg is laid and living a further 30-50 days. Juvenile adult worker bees mature as they work within the hive until they become foragers after about three weeks after emergence. Colony sizes can range from 20,000 to 100,000 (Goodman and Kaczynski, 2015;Owen, 2015). Pressures on a colony include short and long term weather effects, for example, the impact of the season on the hatching rate through changes in relative humidity was observed in (Al-Ghamdi et al., 2014). Other effects include diseases (Arundel, 2011;Gordon et al., 2014;Bull et al., 2012;Fürst et al., 2014), pests (Chandler et al., 2000;Cook et al., 2007;Ryabov et al., 2014), predators (Capri et al., 2013), pesticide, fungicide and herbicide use (Botias et al., 2015;Dively et al., 2015;Godfray et al., 2014;Pettis et al., 2013). When these pressures reduce the number of mature foragers, juvenile adults become precocious foragers, that are less efficient and shorter-lived than mature foragers (Schulz et al., 1998). This can lead to collapse of the colony (Perry et al., 2015). We believe that a necessary first step in understanding the dynamics of bee colonies is to understand their behaviour in a stable environment.
The literature contains both probabilistic models examined by simulations and deterministic mathematical models that can be used to examine the effects of environmental changes and disease on colony dynamics (Becher et al., 2013;Khoury et al., 2013;Lever et al., 2014). These models are themselves complex and, whilst the effects are evident in the results, they often give little understanding of the mechanism. For example, the model description for the simulation program BEEHAVE is comprehensive (Becher et al., 2014) but it is complex and mathematically intractable. Our contribution is a parsimonious and analytically tractable probabilistic model for the daily colony size which captures the key features and behaviour of the natural system under stable conditions. This is the first step in the development of more complex models but is of interest as we demonstrate the existence of a stationary distribution for the colony size. In particular, the mean of the stationary distribution is simply expressed in terms of the mean of the batch size, the probability that eggs hatch into bees and the means of bee lifetime and time between batches. Our model can also be viewed as a type of batch immigration-death model: bees arrive in batches of eggs laid by the queen and subsequently die.
The birth and death process with immigration is well-known in the random process literature (Grimmett and Stirzaker, 2001, pp. 276-278;Lawler, 2006, p. 76). Classical immigration-death process models are first order Markovian.
That is, given the population size at time t the population size at time t + 1 is independent of the past of the process. In our case the daily death rate of the bees depends on the age structure of the colony and is not first order Markovian in this sense. In theory if there is an upper bound on bee lifetimes one could obtain a Markov property by extending the state space but we consider a more direct approach based on marked renewal processes.
Our approach requires modelling the distributions of four quantities. The first is the time in days between the days on which the queen lays batches of eggs. Often this will be one day but our results allow this to have a distribution over the positive integers, with daily layings being a special case. The second is the number of eggs in a batch. The third is the probability an egg successfully hatches into a bee. The fourth is the lifetime of a bee. This initial model represents a perfect bee world without seasonal effects where there are always sufficient resources so that the batch sizes on different days have the same distribution and the probability that an egg hatches into a bee remains the same. We then extend this to a cyclic model where the mean batch size and the probability of an egg hatching vary over an annual cycle which allows us to consider seasonal effects.
Our initial hypothesis was that for the given distributions of the time between batches, the number of eggs in a batch, the indicator whether an egg successfully hatches and the lifetime of a bee, with these quantities being mutually independent, there is a stationary distribution for the size of the colony.
The colony size is modelled as a marked renewal process with the number of batches being a renewal process and the number of bees hatching in each batch being the mark. Then coupling techniques (Lindvall, 1992, pp. 21-34) can be used to study its stationary distribution. In Section 2 we describe the model and in Section 3 we give our theoretical results. The main result, Theorem 1, gives an explicit representation of the stationary distribution of the colony size as the number of days becomes large. To our knowledge such a concise construction for a marked renewal process is not noted elsewhere in the literature. Moreover, when the distribution of the time between the batches is nonrandom, which is a typically case for a bee colony, the representation also enables us to obtain the characteristic function of the stationary distribution of the colony without relying on the renewal technique. See Corollary 2. In Theorem 3 we give results on the stationary distribution for the cyclic model. The practical implications of our results are that the colony size will not endlessly increase but will become stable. In effect it gives a baseline for the colony size. Moreover, our results give an explanation for the recovery in the colony size after swarming. If we suppose that the effect of environmental factors or disease is to either reduce the batch size and the hatching probability, increase mortality through reducing bee lifetime by shifting the lifetime distribution to the left, increase the time between batches or a combination of these effects then our results allow us to quantify their effect on the stationary distribution. In Section 4 we examine this for some parametric models. In Section 5 we give some simulations to verify the theoretical results and explore extinction probabilities as well as the time taken to attain the stationary distribution. The proofs of the results in Section 3 are given in the appendices.

Model and Notation
Let τ i be the time between batch i − 1 and batch i, ζ i be the number of eggs in batch i, I ij ∈ {0, 1} be the indicator that the jth egg in batch i hatches and η ij be the lifetime of bee j in batch i, i ∈ N. Then S i = i j=1 τ i is the time that the ith batch of eggs is laid. The number of batches laid by time t is Clearly, the number of bees that hatch in batch i is ξ i = ζi j=1 I ij . Without loss of generality, shifting {S i } by a fixed constant if necessary, we assume that the hatch time from an egg to a bee is zero so that the number of bees in the colony ξi j=1 1 {Si+ηij ≥t} . In this expository work we suppose that {τ i , i ∈ N} are independent and identically distributed (iid) as are {I ij , i, j ∈ N} and {ζ i , i ∈ N}. We further suppose that {η ij : i, j ∈ N} are iid and that the sequences Under our assumptions the times {S i , i ∈ N} at which the batches are laid form a renewal process. We let F τ (x) and F η (x), x ≥ 0, denote the distribution functions of τ 1 and η 11 respectively and we suppose that E(τ 1 ), Var(τ 1 ), E(η 11 ), Var(η 11 ) are all finite. Note that our zero time is arbitrary and in general will not be the time at which a batch of eggs was laid. That is, the time after we start observing until the first batch is laid is τ 1 = τ 1 − u for some u representing the (random) time after the last batch of eggs was laid that our observations begun.
We assume that τ 1 takes positive integer values with distribution f j = P(τ 1 = j), j ≥ 1, where ∞ j=1 f j = 1. If P(τ 1 = 0) > 0, then we can redefine the ζ i 's so that the new renewal process has inter-renewal times taking positive integer values. Let X n be the time till the next renewal from time n, then {X n } is a Markov chain with transition probabilities For simplicity, we also assume that the lifetimes η ij 's only take non-negative integer values. We use d = to stand for two random variables being equal in distribution, d → for the convergence in distribution and ∼ denotes "is distributed as". We let r = P(I 11 = 1) be the probability an egg hatches and set q k =

Results
To state the main result, we need a metric to quantify the speed of convergence of {M n } to its stationary distribution. We define the total variation distance between two integer valued random variables Y 1 , Y 2 as where Z is the space of all integers.
Moreover, EM = rE(ζ 1 )(Eη 11 + 1)/Eτ 1 , where H(n) := EN n satisfies the renewal equation Theorem 1 not only tells us that the marked renewal process {M n } stabilises as n → ∞ but also gives an explicit representation of the stationary distribution.
It enables us to simulate the stationary distribution, i.e., the distribution of M ∞ , without having to simulate the process {M n } for a long time.
It can be observed from Theorem 1 that when τ i 's are not random (i.e. F τ is degenerate) and the lifetimes η ij are bounded by a constant K then M n reaches stationarity for n ≥ K + 1. For a bee colony it is reasonable to take K to be 80 days. An implication of this is that after swarming (where the size of the colony is typically halved) in the stable case that we consider here, the colony size will quickly recover to the stationary size. Thus, at least in the stable case, no other mechanism is required to explain the recovery of the colony size after swarming.
Let Pn(λ) denote the Poisson distribution with mean λ. When τ i 's are not random, we have the following corollary.
Corollary 2. If F τ is degenerate so that we can assume S i , i ≥ 1, are constants, then the stationary distribution of {M n } has the characteristic function where F η is the distribution function of η 11 and ψ ζ (s) := Es ζ1 . In particular, if If the minimum size of a viable colony is ν then using Corollary 2 to approximate P(M n < ν) gives an approximate lower bound on the probability of colony extinction at time n. This is a lower bound as the colony may have become extinct before n (i.e. P(min k≤n M k < ν) ≥ P(M n < ν)). The computation of the extinction probabilities requires derivation of the distribution of min k≤n M k which is beyond the present work and for now we address this through simulations in Section 5. See Figure 4.
To allow the distribution F ζ,i of ζ i , the hatch probability r i of eggs in batch i and the distribution F η,i of η ij , j ≥ 1, to all depend on time i, we propose the periodic model {M lD+i , i = 1, . . . , D}, l ∈ N, to mimic the behaviour over a cycle of D days, where D = 365 days represents a year. That is we assume then (M lD+n ) l∈N reaches stationarity for lD + n > K, and the characteristic function of the stationary distribution is This theorem establishes that the year to year behaviour of the colony size on a given day is stationary.

Applications
We first use our analytic results to examine the effects of changing the distributions of the hatch probability, the batch size, the bee lifetime and the times between batches on the stationary distribution of the colony size. To examine the behaviour of the colony size in more detail we model the lifetimes using a discretization of the skew normal distribution (Azzalini, 2013) truncated at zero as this distribution allows left skewed lifetime distributions.
Left skewness is important as mortality will be low until bees start foraging and this may not occur until they are aged 40 days or so. The skew normal is a three parameter distribution with a location parameter Ξ, a scale parameter ω and a slant parameter α (Azzalini, 2013). Its density is where φ and Φ respectively denotes the probability density function and the cumulative distribution function of the standard normal law. Note that when ω = 0 the distribution of η is not random but constant equal to Ξ. We denote this density SN (Ξ, ω 2 , α). The mean of the skew normal distribution is µ = Ξ + ωδ 2/π where δ = α/ √ 1 + α 2 . If α < 0 then this distribution is skewed to the left. We discretise this distribution by rounding up to the nearest integer, with negative values rounded up to one.
In Figure 1 as an illustration we plot the SN (63, ω 2 , −6) density for varying values of ω. This indicates that the shape of the distribution can become quite skewed. This skewness gives one way of representing an increase in bee mortality, with an increasingly heavy left tail when more bees die at a younger age. Varying the location and slant parameters further change the shape and location of this distribution, which gives the model considerable flexibility. Therefore we use this model to illustrate our results.
For convenience, we assume that we start observations on the first day of the year so that the subindex i represents the ith batch as well as the ith day (replace i with i mod 365 if i > 365) of the year.
In Figure 2

Simulations
To verify the theoretical results we conduct a series of simulations. We initially suppose ζ i ∼ Pn(µ ζ ) and as the probability of an egg hatching is r we have ξ i ∼ Pn(µ ξ ) with µ ξ = rµ ζ . In Figure 3 we report ten simulations of our model with η ∼ SN (63, 30 2 , −6), and E(ξ i ) = r i E(ζ i ) = 1500 × seas(i) + 127  during the first winter. We can see that one of the simulations resulted in the number of bees being under the nominal threshold of 5000 and in that case the colony would have become extinct. To examine this further if the laying times are constants and the distribution of the numbers of eggs laid, the hatching probability and the distribution of the bee lifetimes are specified, we can apply Corollary 2 to determine lower bounds on the extinction probabilities.
In general, we can always use simulations to estimate extinction probabilities, although for many parameter values this will be small. In Figure 4 we plot the extinction probability (supposing that extinction occur when the number of bees in a colony is less than 5000) according to the mean number of eggs laid per day. Here, we simulated the probability of extinction over 30 days as a function of the mean number of eggs laid. This number corresponds to midwinter in our model. Thus our model may be used to determine the effects of changes in the underlying parameters on the extinction probabilities. Note that it only requires a small change in the mean number of eggs laid per day for the extinction probability to increase from zero to one.
In Figure 5 we plot the number of bees in the colony during a period slightly longer than one year using the model described above but incorporating swarming. That is we assumed that when the colony size was more than 80000, half of the colony leave the colony. Clearly, in the right conditions the colony size recovers and as noted above no other mechanism is required to explain the recovery of the colony size.
A further quantity of interest is the time until the colony size attains the stationary distribution. In our case, as the maximum lifetime of a bee is bounded by 80, the colony size attains the stationary distribution after no more than 81 days as long as the number of eggs laid is sufficiently large. This can be observed in Figure 5. Indeed, 80 days after the colony splits, the number of bees in the hive is as if swarming never happened. From an evolutionary perspective by splitting into two, colonies that swarm have an increased probability that one the colonies will survive.

Discussion
Our model is an encouraging first step in modelling the dynamics of a bee colony. We have developed a probabilistic model for the dynamics of a bee colony in an idealised bee world, which may not of course be attained. This degree of simplicity was necessary to obtain analytic results, but nevertheless, our results explain why a typical colony size in summer is around 90,000 bees. Our interpretation is that the stationary distribution represents a state attainable under ideal conditions and thus in many ways is an upper bound. We extended the stationarity results to a more naturalistic seasonal model and showed the marginal distribution of the colony size on a given stage of the season was stationary. Our model also allowed us to give a lower bound for extinction probabilities and demonstrate how the colony size recovers after swarming, as seen in the natural system. The simple expression for the mean colony size is easily computed in terms of the hatching probability and the mean of the numbers of eggs laid, the expected bee lifetime and the expected time between the laying of batches of eggs. A determination of how a factor such as environmental perturbations, disease pressure or weather impacts the colony size (for example through increased mortality or reduced nectar flow and hence the number of eggs laid per day) then allows an examination of its effect on colony size.
There are several ways in which our work can be extended. An extension of our results would first be to consider the joint behaviour of the colony size over several time points which would potentially allow the derivation of exact extinction probabilities in addition to the lower bound. That is to consider the colony size {M lD+i , i = 1, . . . , D} over a year as "functional data" (Ramsay and Silverman, 1997) and examine the behaviour of this function. Corollary 2 examines the marginal distributions of this quantity. One could then consider the effects of long term trends on the colony size and survival. In particular conditional probabilities of survival given a particular colony size could in theory be derived. That is, given the current colony size of say 7000 bees at a given time of year under what conditions would the probability of surviving over winter be 90% for example? Such calculations could indicate whether supplementary feeding was required to increase the number of eggs laid per day to help ensure survival. By relating the effects of environmental changes and disease to these quantities, their effect on colony size may also be examined. Whilst our current results are for stable and cyclic effects we anticipate extending these to even more complex and realistic situations in the future. For example, the availability of nectar in the hive may depend on the numbers of forager bees so that the average batch size could be taken to be dependent on the number of forager bees. In their deterministic models, Khoury et al. (2013) give the steady state population as a function of bee death rate and food collection rate. Bee death rate relates to our bee lifetime distribution. However, we can consider more subtle effects. We do not explicitly include food collection rate in the current model but we have argued this could be included in the models for the mean batch size. The availability of nectar could also depend on the age related efficiency of the forager bees and the average batch size could be taken to be dependent on the ages of the forager bees. These more complex feedback models are technically challenging and beyond the scope of the current work.
Regarding more technical extensions, recall that a lattice distribution is the distribution of a random variable taking possible values of the form {a + bk, a, b are constants, b = 0, k ∈ Z}. The distribution F τ of the renewal times {τ i } in Theorem 1 is a special case of a lattice distribution with b = 1. One can easily adapt the proof of Theorem 1 to obtain the stationary distribution when F τ either follows a general lattice distribution or a non-lattice distribution. We also note that one can also study {M n } as a discrete-time queueing system with batch arrivals (see Takagi (1993)).
We conclude that any change in the parameters of the lifetime distribution or the numbers of eggs laid can affect the population size. This allows a variety of ways in which a disease for example can affect the lifetime distribution and hence colony size. Whilst their properties can be difficult to derive, we believe that probabilistic models for the size of a bee colony are mathematically tractable and provide useful and readily interpretable information on the dependence of the colony size on the times between batches of eggs being laid, the sizes of the batches, the hatching probability and bee lifetime. We further believe First of all, since ξ 1 = ζ1 j=1 I 1j is a mixed binomial, we have and If F is non-degenerate, we write T n = τ 2 + · · · + τ n , then the condition d = 1 implies d T V (T n , T n +1) = O n −1/2 → 0 as n → ∞ (Lindvall, 1992, pp. 41-43).
This ensures as n → ∞. If F τ is degenerate, then S n = S n , giving d T V (S n , S n ) = 0. Hence, we obtain In fact, when F τ is degenerate, we have M n d =M n , so (A.5) is obvious. For the case that F τ is non-degenerate, we have σ τ := Var(τ 1 ) ∈ (0, ∞). Take a n = 0.75n , b n = 0.5n/Eτ 1 , where x is the biggest integer less than or equal to x. By (Barbour et al., 1992, pp. 253-254), enlarging the probability space if necessary, we can construct iid {τ i } having the same distribution as that of τ 1 and independent of {ξ i , η ij : i, j ∈ N} such that S k := k j=1 τ j : k ∈ N satisfies where the second equality is due to (A.4). Set S j := Therefore, it suffices to show that P (M n = M n ) = O(n −1/2 ). To this end, we observe that For the event {S bn = S bn }, (A.6) ensures that P(S bn = S bn ) = O(n −1/2 ). For the event {S bn > a n }, we have P(S bn > a n ) = P 4στ , the second equality is due to (Chen et al., 2011, Theorem 11.2) and the last equality follows from (Chen et al., 2011, (2.82)). Therefore, it remains to show that the final event of (A.7) has probability of order no more than O(n −1/2 ). This can be worked out as follows. On {S bn = S bn ≤ a n , M n = M n }, the two marked renewal processes coincide from time a n , so M n = M n can only happen when at least one of the bees before a n are still alive at time n. That is, the bee must have lived through the interval [a n , n] with lifetime η ij ≥ n − a n , giving {S bn = S bn ≤ a n , M n = M n } ⊂ ∪ bn−1 i=1 ∪ ξn j=1 {η ij ≥ n − a n }.
This yields P(S bn = S bn ≤ a n , M n = M n ) where the last inequality is due to (A.1). This completes the proof of (A.5).