tiny[table]font=tiny

Modelling  daily weight variation in honey bee hives

Karina Arias-Calluari1,2 , Theotime Colin2, Tanya Latty2, Mary Myerscough1,
Eduardo G. Altmann1

1 School of Mathematics and Statistics, The University of Sydney, Sydney, New
South Wales, Australia

2 School of Life and Environmental Sciences, The University of Sydney, Sydney,
New South Wales, Australia

These authors contributed equally to this work.

## Abstract

A quantitative understanding of the dynamics of bee colonies is important to
support global efforts to improve bee health and enhance pollination services.
Traditional approaches focus either on theoretical models or data-centred
statistical analyses. Here we argue that the combination of these two
approaches is essential to obtain interpretable information on the state of
bee colonies and show how this can be achieved in the case of time series of
intra-day weight variation. We model how the foraging and food processing
activities of bees affect global hive weight through a set of ordinary
differential equations and show how to estimate  reliable ranges for the ten
parameters of this model from measurements  on a single day. Our analysis of
10 hives at different times shows that crucial indicators of the health of
honey bee colonies are estimated robustly and fall in ranges compatible with
previously reported results. The indicators include the amount of food
collected (foraging success) and the number of active foragers, which may be
used to develop early warning indicators of colony failure.

## Author summary

Honey bees are under threat and dying at an alarming rate due to pesticides,
parasites, and other stressors. Obtaining information about the health of bee
colonies is essential to understand how this happens and to identify measures
that can prevent this from happening. Herein, we built a mathematical model of
the daily dynamics of hives that allows such information to be extracted
without detailed and expensive measurements. Based only on measurements of how
the weight of the hive changes during the day, our model can be used to
estimate how many bees are collecting food, how successful they are, and how
much time they spend outside the hive. Due to its simplicity, the model
presented here can be applied to a wide range of hive scale systems and help
beekeepers track how healthy and productive their bees are.

##  1 Introduction

There is growing concern that crop pollination by honey bees might no longer
be a sustainable operation in the near future [1]. Single and combined stress
from pesticides, parasites and diseases have led to an increase in bee colony
mortality worldwide[2, 3, 4, 5, 6, 7]. Despite beekeeper’s efforts to control
pests and diseases, annual colony losses of about 30% are frequently reported
in countries across Europe [8] and North America[9, 10, 11]. Automated methods
to track the status of hives can play a key role in helping beekeepers prevent
the death of colonies[12, 13, 14, 15, 16].

The social organization of a bee colony is complex and relies on positive and
negative feedback loops to maintain homeostasis. Feedback loops can help
buffer the effect of stressors or cause snowballing effects which in some
cases precipitate colony death. A well-known feedback loop that has been
hypothesized to accelerate the failure of the colony is the one that regulates
the ontogeny of foraging[17]. In a bee colony, workers typically spend the
first few days of their adult life performing various maintenance, defence,
and care tasks before transitioning into foragers[18]. The presence of
sufficiently large number of foragers inhibits the transition of younger bees
into foragers[19, 17]. Under ideal conditions, this ensures enough bees are
performing in-hive tasks, and that the colony is able to process all the food
brought back by the foragers and to raise the next generation of workers. If
an environmental stress reduces the number of foragers, this same mechanism
ensures the rapid replacement of the foraging force. With fewer foragers, the
transition of younger bees into foragers is increased and the foragers’
population is rapidly replenished. Under sustained stress however, it has been
hypothesized that this buffering mechanism can lead to generations of workers
starting to forage at a younger age, which reduces their life expectancy and
the foragers’ population increasingly fast until colony death[20]. Thus,
information about the number of forager bees and their activities can provide
important guidance about the health of the hive.

Fig 1: (a) The weight of a bee hive shows systematic variations throughout a
day. The blue dots show the weight of a hive measured every minute for four
consecutive days. (b) Each daily hive weight’s curve present consistent
patterns at specific time ranges of morning, afternoon, evening and night
time.

Directly measuring the number of forager bees is difficult, but the use of new
technologies provides opportunities to estimate them from more accessible data
[21, 12]. Recently developed affordable and precise balances can easily be
fitted to hives in the field to study the long-term weight changes[22] and
shorter within-day weight variations of hives[12, 23]. Within-day hive weight
variations are clearly visible and exhibit similar patterns that are
consistent across multiple day as shown in Fig. 1. In the morning, rapid
weight loss occurs because foragers leave the hive [24]. For this reason,
previous methods have focused on measuring the slope and low point of the
morning departure segment of daily hive weight’s curves[24], and have used
this measure as likely representing forager activity[25]. However, these
phenomenological methods provide limited information about the daily internal
dynamics of the hive. Other methods perform piece-wise linear regression
analysis where up to five breaking points are suggested after an abrupt change
of slope [24, 26]. In particular, these methods do not provide quantitative
estimations of the relevant bee activities such as the number of forager bees
or the amount of food brought to the hive because the different bee activities
are not explicitly modeled and because the simultaneous departures and returns
of bees prevent trivial calculation of the population of forager bees.

In this paper, we develop a mechanistic mathematical model  of the daily
dynamics of hives which, in combination with hive weight time series,
provides quantitative estimations of key indicators of the state of the hive
such as the total number of active forager bees and the average mass of food
they bring to the hive in each trip. We start in Sec. 2.1 with our model
presentation, its qualitative description, and its main properties. In Sec.
2.2 we present how the parameters can be inferred from data. In Sec. 2.3 we
apply our methodology to time series of ten different hives measured across a
year and show how in most cases the results lead to biologically meaningful
estimations. Finally, in Sec. 3 we summarize our conclusions. A detailed
description of our data and computational methods is presented in the
Materials and Method section. The data and codes used in our analysis are
available in a repository [27].

##  2 Results

###  2.1 Model

Our aim is to obtain a mechanistic model that contains the essential processes
necessary to describe how the weight $W$ of a bee hive changes throughout time
$t$ over a daily cycle.

#### Qualitative description.

Our model will incorporate the following well-known processes and honey-bee
behaviours that are also schematically illustrated in  Fig 2:

  1. 1.

In the morning, the weight of the hive rapidly falls as foragers leave the
hive[24]. This initial weight drop usually lasts for a few hours, indicating
that the weight loss is only partly counterbalanced by forager bees returning
to the hive. A successful forager will return with pollen, nectar, water or
propolis. The weight of a forager bee is expected to be higher when they
return than when they leave, if the resources they seek are available in their
environment.

  2. 2.

Around midday, departures and arrivals reach roughly an equilibrium. Hives in
environments with food resources start gaining weight, possibly indicating
that the weight of departing foragers is less than the weight of effective
foragers returning with food to the colony[24]. The weight increases until all
foragers have returned at dusk.

  3. 3.

In the late afternoon, the number of departures reduce until they stop
completely around dusk. Once all the foragers have returned, the hive  stops
gaining weight.

  4. 4.

At night, the weight of a hive declines slowly untill dawn due to
respiration, food consumption, and evaporation. The metabolic activity of a
bee colony remains high at night because bees actively maintain the
temperature of the brood at about 35 degrees, as well as humidity and $CO_{2}$
concentration to optimize development. Bees actively warm the colony by
shivering, or cool down and ventilate the hive by evapotranspiration and
fanning. These activities require a great deal of movement and, thus, causes
honey bees to consume significant amounts of honey at night. Sugars in the
honey are broken down to produce energy and exhaled as $CO_{2}$. A colony also
evaporates a large amount of water, either to facilitate biological processes,
concentrate sugars in food stores, or cool down the hive. Moreover, it is
likely that a small amount of the weight loss at night results from cleaning
activities, with bees likely getting rid of pests, dead bees and chewed brood
cell caps throughout the night. The cycle resumes the next day again.

Fig 2:  Schematic representation of the daily cycle of bee activities that
impact in modeling colony weight, where time flows in a clockwise direction
(blue line). The main time-dependent elements we model are the number of
forager bees $N(t)$ and the food $F(t)$ (in green) inside the hive. The
foragers are represented as active bees that depart the hive (black arrow) and
return to it (red arrow) with food (in orange). The ratio of forager bees
inside the hive at a particular time is represented by $x(t)$, which is
constant before and after the dynamics of forager bees because the $N_{max}$
foragers are all inside the hive. The static weight $W_{0}$ is composed of the
hive structure and bees that do not leave the hive.

The net outcome of the activities described above depends on the hive, the
environment, and the  weather conditions. In favourable environments and at
its peak population in summer, a single bee colony can harvest a few kilograms
of food in a single day. Most of this weight comes in the form of nectar and
honeydew, which generally contains more than 50% water and diluted sugars[28].
Hives in environments with little or no resources may  lose more weight than
they gain in a day as a consequence of honey consumption inducing respiration
and evaporation (i.e., foragers expend more resources than they are able to
bring back to the colony ).

#### Mathematical model.

We consider the time dependent weight $W\in\mathbb{R}^{+}$ (measured in $g$)
to depend on the number of forager bees $N\in\mathbb{Z},0\leq N\leq N_{max}$
(measured in bees) and on the food $F\in\mathbb{R}^{+}$ (measured in g) as

$W(t)=F(t)+wN(t)+W_{0},$ (1)

where $W_{0}$ is the constant weight of the hive and $w$ is the (average)
weight of a bee. Next we discuss the variation in $N(t)$ and $F(t)$.

For mathematical convenience, we model the variation in number of foragers
bees $N(t)$ by focusing on the fraction $x\in\mathbb{R}$ (nondimensional) of
all $N_{max}$ foragers bees that are inside the hive:

$x\equiv\frac{N}{N_{max}}\Rightarrow N(t)=N_{max}x(t).$ (2)

We write the rate of change of $x$ as

$\frac{dx}{dt}=a(t)(1-x)-d(t)x,$ (3)

where $a(t)$ , measured in $1/s$, is the arrival rate of the fraction of
forager bees outside the hive and $d(t)$ , measured in $1/s$, is the departure
rate of the fraction of forager bees inside the hive. The inverse of these
rates can be thought of as the typical period a bee spends outside and inside
the hive, respectively. These terms multiply the fraction of bees outside
($1-x$) and inside ($x$), respectively. The rate of bees arriving is

$A(t)=a(t)(N_{max}-N(t))=a(t)N_{max}(1-x(t))$

and the rate of bees departing is

$D(t)=d(t)N(t)=d(t)N_{max}x(t),$

where $A$ and $D$ are  measured in $bees/s$.

The food inside the hive $F(t)$ increases due to the mass $m$ (measured in
$g$) of food brought in by arriving bees and decreases due to evaporation,
respiration, and food consumption, which we model as a constant loss $\ell$
(measured in $g/s$)111 More sophisticated models could consider a time
dependent loss (e.g., due to temperature or weather variations) and a term
proportionally to $N$ (or $x$), indicating a rate of consumption proportional
to the number of (forager) bee.  so that

$\frac{dF}{dt}=-\ell+mA(t)=N_{max}(ma(t)(1-x))-\ell.$ (4)

Notice that $F(t)$ is linear in $x$ and can be computed directly from $x(t)$
and $a(t)$ by integrating the above equation. Substituting $F(t)$ and $x(t)$
in Eq. (1) we obtain the time-dependent weight $W(t)$ of the hive, our
measurable quantity. This implies that the equations can be solved
sequentially, starting from $x(t)$. The final step of our model is to specify
the departure $d(t)$ and arrival $a(t)$ rates defined in Eq. (3) and then
solve for $x(t)$. We will consider a very simple choice for $a(t)$ and $d(t)$
– piecewise constant in $t$ – which allows for an explicit solutions for
$x(t)$ and, thus, $W(t)$. More sophisticated models could use more elaborated
functions $a(t)$ and $d(t)$ or even delay equations.

#### Day and night cycles.

Without loss of generality, in our model, we consider:

  * $t=0$

to be the time in the morning (around sunrise) at which the bees start exiting
the hives.

  * $t=1$

to be the time in the late afternoon (around sunset) when bees stop leaving
the hive.

Therefore, in all cases below we assume $d(t)=0$ for $t<0$ and $t>1$. The
matching of the model time $t$ to the clock time requires the estimation from
data of two clock times: $t_{0}$ (corresponding to $t=0$) and $t_{1}$
(corresponding to $t=1$).

#### Night dynamics.

Let us start with the (easiest) night case $t>1$, when bees do not depart
$d(t)=0$ and arrive at a constant rate $a(t)=a$. Assuming that at the
beginning of this regime at $t=1$ there are $x_{1}$ bees inside, the
particular solution of Eq. (3) is

$x(t)=1-(1-x_{1})e^{-a(t-1)},$ (5)

which corresponds to an exponential decay towards all the bees being inside
the hive ($x=1$) with a rate $a$. Therefore, we can consider $1/a$
proportional to the duration for which bees remain outside the hive.

#### Day dynamics.

Let us consider the simple case where $d(t)=d$ and $a(t)=a$, both constant in
time. The initial condition in this case is $x=1$ at $t=0$ because we assume
all forager bees spent the night in the hive (i.e., $1/a\ll$ night period and
all previous-day foragers return). In this case, Eq. (3) leads to

$x(t)=\frac{a}{d+a}+\frac{d}{d+a}e^{-(d+a)t}.$ (6)

#### Explicit solutions.

We are now able to combine the results above to write explicit solutions for
the fraction of foragers inside the hive $x(t)$, the weight of food inside the
hive $F(t),$ and the total weight of the hive $W(t)$. We consider the arrival
rates during the day $a_{1}$ and at dusk $a_{2}$ to be different:

$a(t)=\left\\{\begin{array}[]{ll}a_{1}&\text{ for }t\leq 1,\\\ a_{2}&\text{
for }t>1,\end{array}\right.$ (7)

where we expect to retrieve $a_{2}>a_{1}$ (bees spend less time out when it is
getting dark).

Introducing these arrival rates and a constant departure rate $d$ in Eqs. (5)
and (6) we find an explicit solution divided in three time intervals for
$x(t)$:

$x(t)=\left\\{\begin{array}[]{ll}1&\text{ for }t<0,\\\
\frac{a_{1}}{d+a_{1}}+\frac{d}{d+a_{1}}e^{-(d+a_{1})t}&\text{ for }0\leq t\leq
1,\\\ 1-(1-x_{1})e^{-a_{2}(t-1)}&\text{ for }t>1\end{array}\right.$ (8)

where $x_{1}=\frac{a_{1}}{d+a_{1}}+\frac{d}{d+a_{1}}e^{-(d+a_{1})}$ is the
value of $x(t=1)$ and ensures the continuity of $x(t)$.

We can now turn to the dynamics of food inside the hive $F(t)$. From Eq. (4)
we obtain

$F(t)=C+(\tilde{m}-\ell)t-\tilde{m}\int x(t)dt,$

where $\tilde{m}=mN_{max}a$ and $C$ is an integration constant to be fixed
using the initial condition $F(0)=F_{0}$ and $F(1)=F_{1}$. The solution for
$F(t)$ is then solved by integrating Eq. (8), which can be done analytically
in each of the three branches leading us to

$F(t)=\left\\{\begin{array}[]{ll}F_{0}-\ell t&\text{ for }t<0,\\\
F_{0}-\frac{\tilde{m}d}{(d+a_{1})^{2}}+(\tilde{m}(1-\frac{a_{1}}{d+a_{1}})-\ell)t+\tilde{m}\frac{d}{(d+a_{1})^{2}}e^{-(d+a_{1})t}&\text{
for }0\leq t\leq 1,\\\
F_{1}+\ell+\tilde{m}\frac{1-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}x_{1}}}{a_{2}}-\ell
t-\tilde{m}\frac{1-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}x_{1}}}{a_{2}}e^{-a_{2}(t-1)}&\text{
for }t>1,\end{array}\right.$ (9)

where $F_{0}\equiv F(t=0)$ – for this model an arbitrary constant that is
incorporated into the static weight of the hive $W_{0}$ in Eq. (1),
consequently, without loss of generality, we take $F_{0}=0$ – and $F_{1}\equiv
F(t=1)$ is computed setting $t=1$ in the second line of Eq. (9) so that

$F_{1}=F_{0}-\frac{\tilde{m}d}{(d+a_{1})^{2}}+\left(\tilde{m}\left(1-\frac{a_{1}}{a_{1}+d}\right)-\ell\right)+\frac{\tilde{m}d}{(d+a_{1})^{2}}e^{-(d+a_{1})}.$

Finally, the weight of the hive $W(t)$ is obtained from Eq. (1) as the sum of
Eq. (8) (times $wN_{max}$), Eq. (9), and the constant weight of the hive
$W_{0}$, leading to an explicit piece-wise continuous function $W(t)$ that
will be compared to the measured data. A summary of the parameters used in the
model is given in Table 1.

Parameter | Description | Prior
---|---|---
|  | Range
$w:$ |  the average weight of a single bee, introduced in Eq. (1). |  0.113 g/bee [29, 30]
$W_{0}:$ | the static weight of the hive that includes the |  $[5000\,,\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}max(W)}]$g
|  structure, food stores, non-foraging bees, etc., see Eq. (1) |
$N_{max}:$ | the total number of foragers active on that day, |  $(0\,,\,80000]$ bee
|  introduced in Eq. (2) |
$\ell:$ |  the continuous loss, introduced in Eq. (4) |  $(0\,,\,\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}max(W)}/h]$ g/h
$m:$ |  the average mass of food brought by foragers, introduced in Eq. (4) |  $(0\,,\,0.78w]$ g/bee [31, 32, 33]
$a_{1},a_{2}:$ |  the rates of arrivals, see Eq. (7), |  $[0.10\,,\,4.80]$ 1/h [34, 35, 36]
| inversely proportional to the time bees spend out- and |
$d:$ |  the rates of departures, see Eq. (7), |  $[0.81\,,\,3.01]$ 1/h [34, 35, 36]
| inversely proportional to the time bees spend out- and |
| in-side the hive, respectively. |
$t_{0},t_{1}:$ |  the clock time corresponding to $t=0$ and $t=1$, respectively |  $[6.50\,\,9.00]$ h
|  |  $[14.5\,\,18.5]$ h
Table 1: The parameters $\theta$ of our model. For simplicity, the average
weight of bees $w$ is fixed based on previously reported results. The
remaining parameters are inferred from data within a prior range of admissible
parameters (right column).  The boundaries of the range of flat priors were
chosen based on the cited references or on values smaller/larger than
reasonably possible values (e.g., the maximum weight max(W) is chosen to be
the maximum total weight of the hive at all times). In all cases reported
below, the inferred parameters remained away from the boundaries, indicating
that our choice of boundaries has no influence on our results.

#### Effective description of the model solution W(t).

The mathematical model described above is divided in three regimes –
$t<0,0\leq t\leq 1,t>1$ – with similar characteristics as shown in Fig. 3. The
three regimes contain segments of linear dependence of $W(t)$ – a decay in the
first and last regime, and typically an increase in the intermediate regime –
with a transition period (shaded region) between the regimes:

  * $t<0$

The first regime corresponds to all bees inside the hive $(x=1)$ and a linear
decay of $W$ on time $t$ with rate $\ell$– evaporation and respiration rate
defined in Eq. (4)– and with an intercept $A\equiv W(t=0)$.

  * $0\leq t\leq 1$

The second regime starts with an exponential relaxation of $x(t)$ from $0$ to
the fixed point value $x^{*}\equiv\frac{a}{d+a}$, as described in Eq. (6),
which leads to a sharp decay of $W(t)$ until an inflection time $t_{c}$.
Shortly after $t_{c}$, $W(t)$ shows a linear dependency with rate $\alpha$ and
constant $B$.

  * $t>1$

The third regime shows an exponential relaxation of $x(t)$ from $x^{*}$ to
$x=1$ as described in Eq. (5), which leads to a sharp growth in $W(t)$ as all
the bees return and before the linear decay of the first regime starts again
completing the cycle.

The critical regime describing bee activities is the second regime that takes
place during the day $0\leq t\leq 1$. It can be effectively described by the
parameters defined in Table 1, which can be rewritten in terms of the robust
parameters of the model as:

Fig 3: Graphical representation of the solution of our model. (a) The weight
as a function of time, $W(t)$, in which the three linear regimes and the two
transition regimes (shaded areas) are visible. (b-d) The three key elements
that compose $W(t)$: (b) the number of foragers inside the hive $x(t)$ that
converges to the constant value $x^{*}$; (c) the amount of food $F(t)$ inside
the hive; and (d) the constant weight of the hive $W_{0}$. In all four plots,
five curves are shown for different model parameters
$\theta=\\{N_{max},W_{0},a_{1},a_{2},d,m,\ell=18.74\,\,g/h,w=0.113\,\,g/bee\\}$
that lead to the same effective parameters $A,B,t_{c},\alpha$ shown in panel
(a). The values of the remaining parameters in numerical order for $N_{max}$
are
$W_{0}[g]=26230,\,\,26173,\,\,26117,\,\,26060,\,\,26004;\,\,a_{1}[1/h]=0.882,\,\,0.914,\,\,0.941,\,\,0.964,\,\,,0.984;\,\,a_{2}[1/h]=2.329,\,\,2.328,\,\,2.328,\,\,2.327,\,\,2.327;\,\,d[1/h]=0.469,\,\,0.424,\,\,0.387,\,\,0.356,\,\,0.330;and\,\,m[g/bee]=0.038,\,\,0.037,\,\,0.036,\,\,0.035,\,\,0.034$.
The two grey spans highlight the time intervals of exponential decay and
growth in $x(t)$, respectively.

  * A:

The total weight of the hive at $W({t=0})$ defined as: $A=N_{max}w+W_{0}$

  * $\alpha$:

The rate of weight gain due to foraging, from Eq. (9):
$\alpha=\tilde{m}(1-\frac{a_{1}}{d+a_{1}})-\ell$

  * $B:$

The intersection point of weight gain at $t=0$:
$B=-(\alpha+\ell)\left(\frac{w}{ma_{1}}+\frac{1}{d+a_{1}}\right)+A$

  * $t_{c}:$

The time of the minimum weight:
$t_{c}={\frac{-1}{({d+a_{1}})}\ln\left(\frac{\alpha}{dN_{max}w+\alpha+\ell}\right)}$

Fig. 3 shows an important mathematical property of our model: different
choices of parameters lead to nearly identical curves $W(t)$ making our model
effectively underdetermined by $W(t)$. Mathematically, this underdetermination
happens because the _four_ effective parameters introduced above
($A,B,\alpha,t_{c}$), which provide a simplified quantitative description of
the solution $W(t)$, depend on _five_ model parameters
($N_{max},W_{0},m,a_{1},d$, note that $\tilde{m}\equiv mN_{max}$ and the mass
of a bee $w$ is taken as fixed). This means that there are different choices
of the parameters $\theta$ of our model that lead to the same four effective
parameters describing the second regime $W(t)$ (even if the curves $W(t)$ are
mathematically distinct, their difference is indistinguishable in Fig. 3 and
much smaller than the experimental precision of the data). Biologically, this
underdetermination appears due to our definition of which bees in the hive
count as forager bees. Our model parameters are fixed and can be thought as an
average over a spectrum of bees that ranges from very active to those that do
only very few forager excursions in a day. Depending on whether we include the
least active bees as foragers or not, we obtain different model parameters but
the same effective parameters defined above. For instance, the same initial
weight $A\equiv W(t=0)=N_{max}w+W_{0}$ can be obtained by including the less
active bees as foragers – thus increasing $N_{max}$ – or as part of the non-
foragers bees that spend most of the time in the hive – thus contributing to
the constant $W_{0}$. Similarly, counting the least active bees as foragers
reduces the average mass $m$ of food brought to the hive, increases the
average time they spend in the hive (i.e., decreasing $d$), decreases the time
spent outside the hive (i.e., increase $a_{1},a_{2}$), and leads to a larger
fraction of forager bees inside the hive ($x^{*}$ increases). As shown in Fig.
3, solutions of $W(t)$ with very different $N_{max},W_{0}$ in panels (b) and
(d) – which imply also different $d,a_{1},a_{2},\textrm{and}\,\,m$ – lead to
essentially the same curves for the food $F(t)$ in panel (c) and weight $W(t)$
in panel (a).

###  2.2 Connecting model and data

Here we show how to combine our model with data to extract information about
the state of hives. We describe the theoretical and computational methodology
used to infer the model parameters $\theta$ from time series of hive weights
$W^{\prime}(t)$.

#### Theory.

We assume that the measured weight $W^{\prime}$ deviates from the actual
weight $W$ as

$W^{\prime}=W+\varepsilon,$

where $\varepsilon$ is an independent Gaussian distributed random variable
with zero mean and standard deviation $\sigma$, i.e.,
$\varepsilon\sim\mathcal{N}(0,\sigma^{2})$, that accounts for measurement
errors and random effects not accounted in our model.

If $\theta$ are the model parameters generating $W\equiv W_{\theta}$, at an
arbitrary time $t$, as described in the model of the previous section, the
probability likelihood of an observation $W^{\prime}$ is a Gaussian
distribution centred at $W_{\theta}$ given by

$P(W^{\prime}|\theta)=\mathcal{N}(W_{\theta};\sigma^{2}).$

Assuming that observations are performed at times $q_{1},q_{2},\ldots q_{n}$
with $q_{i}=q_{0}+n\delta q$ and are (conditionally) independent from each
other, the log-likelihood of the observations is given by:

$\log
P(W^{\prime}|\theta)=\sum_{i=1}^{n}\log\mathcal{N}(W_{\theta}(q_{i});\sigma)=\sum_{i=1}^{n}-\frac{(W^{\prime}(q_{i})-W_{\theta}(q_{i}))^{2}}{2\sigma^{2}}-\log(\sqrt{2\pi}\sigma)$
(10)

We can include further information on the parameters $\theta$ by considering a
prior distribution $P(\theta)$. Here we consider flat priors for all
parameters, i.e., $P(\theta)=cst.$ for $\theta\in[\theta_{1},\theta_{2}]$ as
specified in Tab. 1 above. Effectively, this restricts the range of admissible
parameters without changing the probability of parameters within the range.
The estimation of the parameters $\theta$ can be done maximizing the posterior
distribution $P(\theta|D)$ which in our case is equivalent to the maximization
of the log-likelihood $(\ref{eq.L})$ or to the minimisation of the squares of
$W^{\prime}-W$ as

$L=\sum_{i=1}^{n}(W^{\prime}(\tau_{i})-W(\tau_{i}))^{2}$ (11)

in the range of admissible parameters set by the priors.

#### Computation.

The problem of obtaining the ten parameters introduced in our model – see
Table 1 – has been thus reduced to the optimization problem of finding the
minimun of $L$ in a ten-dimensional space. The computational method we use for
this optimization is centred around the Levenberg–Marquardt algorithm as
implemented in the library SciPy (See Materials and Methods Sec-4.1). However,
a robust estimation of the parameters cannot be obtained directly through the
naïve application of this method both because of the high-dimensionality of
the problem and the effective underdetermination of our model reported at the
end of last section (which leads to essentially flat regions of $L$ in the
parameter space with multiple local minima). Here we use properties of bee
hives and our model to simplify this problem and obtain a robust estimation of
parameters.

We start by considering  four parameters:

  * $w$

For simplicity, we fix the average weight of bees to $w=0.113\,g/bee$[29, 30].

  * $t_{0},t_{1}$

We infer the times $t_{0}$ and $t_{1}$ in a grid of plausible values
(separated by five minutes). During this minimization we do not distinguish
between model parameters that lead to the same effective parameters
$A,B,\alpha,\,\,\textrm{and}\,\,t_{c}$ because $t_{0}$ and $t_{1}$ are not
involved in them. The estimation of $t_{0}$ and $t_{1}$ allow us to re-scale
the data to the model time scale, facilitating the numerical estimation of the
remaining parameters.

  * $\ell$

The continuous loss can be calculated directly from the negative linear
decrease when $t<0$.

These parameters, and the effective parameters ($A,B,\alpha,t_{c}$), are not
strongly affected by the underdetermination problem discussed above and are
thus part of our robust estimation. The remaining six parameters –
$N_{max},W_{0},m,a_{1},a_{2},d$ –  are directly connected to the
underdetermination problem and need to be estimated with care. We start by
obtaining a robust estimation of $A\,\,,B$,$\,\,\alpha$, and $t_{c}$ by direct
minimization of $L$. The mapping of the effective parameters to the model
parameters is done in two different ways:

  * (i)

We  seek a range of plausible model parameters determined as the values of
$W_{0}$ and $N_{max}$ that keep the four effective parameters constant and
that maintain the parameters $\ell,m,a_{1},a_{2}$ and $d$ within their prior
range (set in Table. 1). In practice, this is done by defining different
linear combinations between $N_{max}$ and $W_{0}$ that are introduced in the
system of six equations that connect the parameters. For each possible
$N_{max}$, we obtain a particular $a_{1}$, $d$, and $m$ by solving the systems
of equations given by $A\,\,,B$ and $t_{c}$. The second rate of arrival,
$a_{2}$, is obtained by the conservation of departed and returned bees. This
method leads to the estimation of an interval of possible parameters.

  * (ii)

Alternatively, by pre-defining one parameter, we calculate the $N_{max}$ value
directly because the system of equations presents a unique solution. We focus
on the departure rate $d$ because  of the consistency among different reported
studies on the time interval between trips that a forager bee spends inside
the hive [37] and because $d$ can be directly connected to the simplifying
modeling decision of dividing bees into foragers and those that do not leave
the hive (i.e., we consider as foragers the most active bees so that
collectively their average departure rate is $d$, the other bees are
considered to effectively stay inside the hive, and their weight is counted as
part of $W_{0}$). In practice, we choose $d$ so that the average time spent by
forager bees inside the hive between trips is  $\tau_{d}=0.816$ h 222the half-
life $\tau$ is derived from the arrival or departure rate, and defined in Eq.
(12) below. [36].

Even if the first approach does not lead to a point estimation of the
parameters, in practice most solutions can be restricted to relatively narrow
intervals by considering physically-based constraints. Some of these
constraints are the average mass of food carried by a single bee $m$, the
interval of time between trips $\tau$, and the number of trips needed to
satisfy food requirements. Details on the computational methods appear in
Materials and Methods Sec.-4.2 and the codes are available in the repository
[27].

###  2.3 Results from data analysis

In this section we report and discuss the results obtained by applying the
model defined in Sec. 2.1 to data using the inference methods explained in
Sec. 2.2.

#### Data.

We  used 24h time series of weights $W^{\prime}(t)$ sampled every $\delta t=1$
minute for ten different hives on different days. The data was collected
between 27th of November 2017 and 2nd of March 2019, from hives located at
Macquarie Park, NSW, Australia (33° 46′ 06.6″ S 151° 06′ 43.8″ E, see [34]).
Here we present data from the 7th of April in Table 2 and 5 additional dates
in Figure 5 (b) and (c). In all cases, the daily maximum temperature was above
16∘C (which is sufficient for honeybees to forage [38]), and no rain was
recorded (ensuring that the recorded weights are not affected by
precipitation). This dataset also contains an independent estimation of the
total number of bees in each of the hives, obtained from measures of the total
weight of the hive with and without bees (performed two days after the day we
use to infer the parameters of our model) [34]. As a pre-processing step to
the analysis of weight time series, we removed outliers that correspond to
exogenous influences as described in Materials and Methods Sec.-4.1.

Fig 4: Comparison between the model and data. Data corresponds to measurements
on 2018-4-7 for ‘Hive 6’. By fixing $\tau_{d}=49$ min, the estimated
parameters are $N_{max}=5830$ bees, $W_{0}=24068\,$g, $m=0.03$ g/bee,
$a_{1}=0.99$ $h^{-1}$, $a_{2}=2.1$ $h^{-1}$ , $d=0.85$ $h^{-1}$, $l=20.66$
g/h, $t_{0}=8.40\,$h and $t_{1}=17.93\,$h.

#### Representative case.

Fig 4 shows the comparison between our model with inferred parameters and data
for one specific hive. We see how the combination of the three piece-wise
continuous curves in our model (for $t<0$, $0\leq t\leq 1$, and $t>1$,
respectively) provides a reasonable account of the characteristic pattern of
within-day hive weight variations. The main advantage of our methodology is
the mechanistic generative model associated to this curve which ensures the
parameters describing the curve have a relevant meaning. For instance, we
estimate that on this day there were $N_{max}=5,830$ active forager bees, that
they brought on average $0.03$ g of food on each trip, and that trips were on
average  $\tau_{a_{1}}=0.70\,\textrm{h}$ 2 long.

Fig 5: Inference of activity times $t_{0}$ and $t_{1}$. (a) Contour plot of
$L$ obtained performing multiple minimizations at fixed $t_{0}$ and $t_{1}$
(‘Hive 6’ at 2018-4-7). The minimum $L$ is at $t_{0}=8.40\,\mathrm{h}$ and
$t_{1}=17.93\,\mathrm{h}$. (b)  Estimated $t_{0}$ and $t_{1}$ fluctuate around
similar values for ten different hives located in close proximity when
analysing data sets during the same day. E.g. at 2018-4-7, the different
activity times oscillate around $t_{0}=8.4\,\mathrm{h}$ and
$t_{1}=17.85\,\mathrm{h}$, respectively (dashed lines). (c) Time evolution for
the optimum $t_{0}$ and $t_{1}$ for ‘Hive 6’ only. Transitioning from Autumn
to Winter,  we observe a linear increase for $t_{0}$ and a linear decrease for
$t_{1}$ (linear regression).

#### Activity time.

As a further test of the robustness of our inference method, and the
meaningfulness of its outcome, we look in detail at the inferred times in
which foragers start ($t_{0}$) and stop ($t_{1}$) departing the hive on
different days and different hives. The results summarized in Fig. 5 show
that: (a) a single well-defined optimal parameter value exists, confirming the
first and crucial step of our inference procedure; (b) similar values are
obtained across multiple hives in the same day and similar location, in
agreement with the expectation that the initial times are similar for hives at
similar locations; and (c) the estimated values for the same hive varies
throughout the year, with the period of activity becoming shorter as the days
become shorter (in the southern-hemisphere winter). These results confirm the
success of our approach which is to treat these times as additional parameters
to be inferred from the data simultaneously with the other parameters.

Hives
---
$\theta$ |  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 |  10
Robust estimation
$t_{0}\,$(h) |  8.90 |  8.30 |  8.70 |  8.20 |  8.10 |  8.40 |  8.40 |  8.50 |  8.40 |  8.10
$t_{1}\,$(h) |  17.77 |  17.36 |  18.42 |  18.26 |  17.36 |  17.93 |  17.93 |  17.77 |  17.77 |  17.93
$l\,\,$(g/h) |  9.762 |  40.805 |  5.600 |  16.628 |  33.183 |  20.663 |  15.016 |  3.974 |  7.687 |  4.399
$A$(g) |  17,503 |  23,384 |  17,878 |  26,809 |  23,412 |  24,731 |  21,244 |  18,057 |  19,485 |  19,364
$B$(g) |  17,409 |  22,.849 |  17,808 |  26,562 |  22,905 |  24,369 |  20,940 |  17,925 |  19,342 |  19,262
$\alpha\,$(g/h) |  16.594 |  95.081 |  4.595 |  38.529 |  72.577 |  59.560 |  24.780 |  4.721 |  23.754 |  4.476
$t_{c}\,$(h) |  9.521 |  9.598 |  9.852 |  9.633 |  9.565 |  9.706 |  10.525 |  13.485 |  9.953 |  11.423
Interval estimation
$N_{max}\,$ (bees) |  [1950, 5400] |  [5050, 8200] |  [940,  2500] |  [2310, 3450] |  [7180, 7650] |  [4150, 6030] |  $\emptyset$ |  $\emptyset$ |  $\emptyset$ |  [840,  940]
$W_{0}$(g) |  [16892, 17282] |  [22451, 22794] |  [17596, 17772] |  [26413, 26537] |  [22542, 22594] |  [24045, 24254] |  $\emptyset$ |  $\emptyset$ |  $\emptyset$ |  [19257, 19267]
$a_{1}$(1/h) |  [3.59, 4.79] |  [0.61, 0.92] |  [1.44, 2.63] |  [0.54, 0.79] |  [0.81, 0.85] |  [0.76, 1.01] |  $\emptyset$ |  $\emptyset$ |  $\emptyset$ |  [0.20, 0.25]
$\tau_{a1}$(h) |  [0.14, 0.19] |  [0.75, 1.13] |  [0.26, 0.48] |  [0.88, 1.28] |  [0.82, 0.86] |  [0.69, 0.92] |  $\emptyset$ |  $\emptyset$ |  $\emptyset$ |  [2.75, 3.40]
$a_{2}$(1/h) |  [1.34, 1.34] |  [1.01, 1.02] |  [3.20, 3.21] |  [2.43, 2.47] |  [1.07, 1.07] |  [2.10, 2.11] |  $\emptyset$ |  $\emptyset$ |  $\emptyset$ |  [2.36, 2.35]
$\tau_{a2}$(h) |  [0.52, 0.52] |  [0.68, 0.69] |  [0.22, 0.22] |  [0.28, 0.29] |  [0.65, 0.65] |  [0.33, 0.33] |  $\emptyset$ |  $\emptyset$ |  $\emptyset$ |  [0.29, 0.30]
$d$(1/h) |  [0.82, 2.38] |  [0.82, 1.41] |  [0.82, 2.28] |  [0.82, 1.26] |  [0.82, 0.87] |  [0.82, 1.24] |  $\emptyset$ |  $\emptyset$ |  $\emptyset$ |  [0.82, 0.94]
$\tau_{d}$(h) |  [0.29, 0.85] |  [0.49, 0.85] |  [0.30, 0.85] |  [0.55, 0.85] |  [0.79, 0.85] |  [0.56, 0.85] |  $\emptyset$ |  $\emptyset$ |  $\emptyset$ |  [0.74, 0.85]
$m\left(\frac{\times 10^{-3}g}{bee}\right)$ |  [7.00, 10.00] |  [37.40, 60.00] |  [7.00, 12.00] |  [38.50, 59.80] |  [32.40, 34.10] |  [28.80, 39.60] |  $\emptyset$ |  $\emptyset$ |  $\emptyset$ |  [39.00, 45.00]
Estimation assuming  $\tau_{d}=0.816$ h
$N_{max}\,$ (bees) |  5220$\pm$0 |  7930$\pm$0 |  2410$\pm$0 |  3330$\pm$0 |  7400$\pm$0 |  5830$\pm$0 |  3540$\pm$0 |  970$\pm$0 |  1260$\pm$0 |  910$\pm$0
$W_{0}$(g) |  16912$\pm$3 |  22481$\pm$7 |  17606$\pm$5 |  26426$\pm$3 |  22570$\pm$4 |  24068$\pm$5 |  20840$\pm$1 |  17946$\pm$4 |  19339$\pm$2 |  19260$\pm$1
$a_{1}$(1/h) |  4.77$\pm$0.6 |  0.90$\pm$0.0 |  2.60$\pm$0.2 |  0.77$\pm$0.0 |  0.83$\pm$0.0 |  0.99$\pm$0.1 |  0.50$\pm$0.0 |  0.11$\pm$0.0 |  0.44$\pm$0.0 |  0.24$\pm$0.0
$a_{2}$(1/h) |  1.34$\pm$0.1 |  1.02$\pm$0.1 |  3.21$\pm$0.3 |  2.47$\pm$0.2 |  1.07$\pm$0.1 |  2.11$\pm$0.1 |  2.00$\pm$0.3 |  1.56$\pm$0.1 |  1.84$\pm$0.1 |  2.35$\pm$0.2
$d$(1/h) |  0.85$\pm$0.1 |  0.85$\pm$0.0 |  0.85$\pm$0.1 |  0.85$\pm$0.1 |  0.85$\pm$0.0 |  0.85$\pm$0.1 |  0.85$\pm$0.0 |  0.85$\pm$0.0 |  0.85$\pm$0.1 |  0.85$\pm$0.0
$m\left(\frac{\times 10^{-3}g}{bee}\right)$ |  10.0$\pm$1 |  38.0$\pm$2 |  10.0$\pm$1 |  40.0$\pm$3 |  33.0$\pm$1 |  30.0$\pm$3 |  33.0$\pm$1 |  53.0$\pm$2 |  77.0$\pm$3 |  41.0$\pm$2
Table 2: Inferred parameters for 10 different hives on the same day
(2018-4-7). Values $\emptyset$ in the interval estimation indicates that there
was no interval of parameters satisfying all imposed constraints for this
hive.The values of $\tau$ next to the parameters $a_{1},a_{2},$ and $d$
correspond to the time spent outside (respectively, inside) the hive and were
computed using Eq. (12). The uncertainties in the values of the estimation
assuming  $\tau_{d}=0.816$ h minutes, reported under the curly brackets, were
computed using bootstrapping, see Materials and Methods Sec.-4.3.

#### Systematic analysis.

Motivated by the success of our methodology in the cases above, we applied it
systematically to 10 different hives on the same day (2018-4-7). The results
obtained from our two inference methods are reported in Table 2. The
comparison of the inferred parameters of these hives with expectations about
the behaviour of bees and hives provide further evidence that our mathematical
model is capturing meaningful information from the data. In particular:

  * •

The time forager bees start ($t_{0}$) and stop ($t_{1}$) departing the hives
are aligned with sunlight.

  * •

The time of minimum weight was inferred to be around $t_{c}=9:45$ a.m. shortly
after $t_{0}$ (for all hives except ‘Hive 8’, which presents a later minimum).
This result is compatible with the minimum daily weights reported in
experimental data [39], which reveal minimum weights around $9:00$ a.m.
Recently, results from segmented linear regression showed minimum losses can
take place even 4 hours after sunrise [40].

  * •

The estimation of the number of forager bees $N_{max}$ is around a few
thousand bees, in line with alternative methods of estimation.  In particular,
Figure 6(a) shows that the total number of bees per hive correlates positively
with our estimation of $N_{max}$ ($R^{2}=0.596$, p=0.036)333 In all our
analysis of correlation between variables, the reported p-value p is the
probability of obtaining a positive coefficient of determination $R^{2}$ equal
or larger than the reported $R^{2}$ under the null hypothesis that the
variables are independent from each other.. The fraction of forager bees
(assuming $\tau_{d}=0.816h$) is estimated to be around $18\%$, comparable to
previous estimations (around $25\%-30\%$ in Refs. [41, 42]). Additionally, we
obtained a positive correlation between $N_{max}$ and  rate of weight gained
during the day $\alpha$, ($R^{2}$=0.778, p=0.0056) consistent with the
previously reported link between the size of the population and the amplitude
of weight variations during the active foraging period [40].

  * •

The average amount of food brought back to the hive by a bee was inferred to
be typically in the range  $m\in[10,40]\times 10^{-3}$ g/bee, which is
comparable to reported crop loads carried by active foraging bees when
returning to the hive with loads of pollen [31], water [32], and different
sweet solutions [33]. The two exceptions (Hives 8 and 9) show values up to
$m=77\times 10^{-3}$ g/bee and possible correspond to sugar concentrations
with total dissolved solids (TDS) of 60% [33].

  * •

In 90% of the cases, $a_{1}<a_{2}$, which is in line with the expectation that
bees return faster to the hive at dusk (after they stop leaving).

  * •

The typical time spent by forager bees inside and outside the hive is computed
through an analogy with the half-life, $\tau$, which represents the time that
half of the foragers’ colony will need to go outside the hive or vice-versa.
The half-life, $\tau$ depends on the departure or arrival rate
$r\in\\{d,a_{1},a_{2}\\}$ in hours as

$\tau=(t_{1}-t_{0})\frac{\ln 2}{r}.$ (12)

From the calculated values of $d,a_{1,2},$ and $t_{1}-t_{0}$, we determine
that the time in the hive between trips, $\tau_{d}$, ranges from 0.29h to
0.85h and the typical time of a foraging trip,$\tau_{a}$, is in between 0.14h
to 3.40h. These estimations align with the results in Refs. [37, 36, 42]. The
only hive that shows values up to 6.25 h for trip duration is Hive 8, this
hive shows behaviour outside the ones predicted by our model, as evidenced by
$\emptyset$ in the interval estimation. However, some reports points that bees
are able to be outside the hive 6.25 h at the age of 36 days [35].

Fig 6:  Comparing our model to independent measurements. (a) Our estimation of
the number of active foragers $N_{max}$ correlates positively with the total
number of bees ($R^{2}=0.596,p=0.036)$. The dashed straight line corresponds
to 18.2% of all bees being foragers, constant across all hives. The total
number of bees was estimated through an independent measurement performed on
2018-04-9, two days after the data used in our inferrence [34]. (b)The plot
represents the difference between the number of bees arriving and bees
departing a honey-bee hive during a 5 minutes interval.  The model prediction
(blue line) was computed integrating $A(t)-D(t)$ over each time interval. The
real data was measured using bee-tracking methods.

#### Arrivals and departures.

Finally, we test the predictions of our model against independent experimental
measures that tracked the number of bees that arrive and depart from a hive at
a time interval [37]. Figure 6(b) compares these measurements to the
assumption in our model and reveals the limitations of our simplifying choice
of constant and piece-wise linear arrivals and departures in Eq. (7).  While
our model captures the main characteristics of the data, unsurprisingly, bee
activities do not start or end abruptly and the measured data shows smoother
variations than the prediction of our model. A similar behaviour can be
observed for $W(t)$ – such as in Fig. 3 – where the measured data do not show
the abrupt changes seen in the model curve at $t=0$ and $t=1$ (discontinuous
derivative). This limitation can be partially overcome within our modeling
framework – Eqs. (1) to (4) – by changing the assumptions of constant $a(t)$
and $d(t)$ – e.g., in Eq. (7). For instance, the next simplest alternative is
to consider that the arrival remains constant $a(t)=a$ but that the departures
varies continuously from $0$ at time $t=0$ to a maximum value $d_{max}$ at
time $t=0.5$ and back to $0$ at $t=1$. The simplest polynomial with these
characteristics is

$d(t)=4d_{max}t(1-t).$ (13)

In this case, Eq. (3) is a linear first-order equation and the solution is:

$x(t)=\frac{C}{e^{r(t)}}+a\int e^{r(t)dt}dt,$ (14)

where $r(t)=at+2d_{max}t^{2}-\frac{4}{3}d_{max}t^{3}$ and $C$ is the arbitrary
constant to be fixed by the initial condition $x(0)=1$. This alternative model
has the same number of parameters ($d_{max}$ instead of $d$), but requires
additional computational resources to numerically find the solution for $x(t)$
through the integration of Eq. (14). This example illustrate the flexibility
of our model to consider more realistic settings and also the trade-off
between model complexity and the need of computational approaches.

##  3 Discussion

We introduced a simple mathematical model that describes the hive processes
affecting within-day weight variations. Simplistic assumptions of the
departure and arrival of bees allow us to obtain closed form solutions for the
model that result in a continuous, piece-wise differentiable, curve $W(t)$
with $10$ biologically-relevant parameters. We inferred the parameters of
$W(t)$ for $10$ different hives and for different days across the year.
Overall, the results reveal that the parameters we infer are compatible with
previous knowledge of bee behavior.

The main advantage of our approach when compared to previous works is that it
combines mechanistic mathematical models with data analysis to obtain
interpretable quantitative information about the condition of the hive.
Additional advantages of our model approach include: (i) indications from the
inference of when results cannot be trusted (either because they are not
capturing meaningful information or because optimal results lie outside of the
allowed error range); (ii) that our framework allows for the proposal of more
sophisticated models, e.g., by considering additional terms contributing to
the food dynamics in Eq. (4) or more detailed modes of arrival $a(t)$ and
departure $d(t)$ such as in Eq. (13). An important methodological expansion of
our work would be to consider and compare different type of functional forms
in our models. This will require considering cases in which solutions are not
obtained in closed forms and the implementation of a Bayesian-model comparison
that can compare models of different complexities, in line with modern
simulation-based and Bayesian inference that are increasingly applied to
dynamical systems of biological interest [43, 44, 45].

The success of our model in extracting meaningful information about bee hives
opens the possibility of using it for a systematic analysis of the health of
hives across different hives and across time. The following quantities
provided by our model are of particular interest:

  * •

The number of active foragers in a bee colony $N$ is a crucial parameter of
the capacity of a colony at sustaining itself. Environmental stressors such as
pesticides can cause a rapid depletion of the population of active foragers
[46]. Estimations of $N$ provide a rapid assessment of the daily foragers
population and allow quick intervention and stress reduction to prevent colony
collapse [20].

  * •

The time spent by foragers outside the hive $\tau_{a}$ can be connected to the
efficiency of the colony because long times outside the hive lead to a
reduction of the life expectancy of foragers [37].

  * •

The amount of food collected by bees in each trip $m$ allows for a better
assessment of foraging difficulty and thus of environments that can best
support beekeeping operations [22].

  * •

The rate of weight grow of a hive $\alpha$ gives a reliable indicator of
whether foragers are supplying their colony properly (leading to optimal
colony development) or if instead the colony is endangered (suggesting the
existence of stressors or environmental limitations).

Tracking the daily evolution of these parameters opens the perspective of
developing weight-based early-warning indicators of bee colony failure that
could lead to new strategies for risk control.

##  4 Materials and Methods

###  4.1 Computational methods

Here we describe the computational methods used to estimate ranges of the
parameters $\theta$ of our model presented in Table 1 from the within day time
series of weights $W^{\prime}(t)$. In Sec.(2.2) of this manuscript this issue
has been formulated as the problem of finding the parameters $\theta$ that
minimize the sum of squares $L$ in Eq. (11). The ten parameters proposed in
this model were reduced to nine parameters after considering the weight of a
bee $w$ as $0.113$ g/bee.

The minimisation was performed in a nine-dimensional space through a least-
square method composed of three different approximation levels. First, the
clock times when bees started departing $t_{0}$, and when they stopped leaving
$t_{1}$. Next, a robust estimation of the characteristic features from the
time series of the weight fluctuations was performed. Finally, the remaining
seven parameters were obtained from the connection with the robust parameters.

For the estimation of clock times for $t_{0}$ and $t_{1}$ we created a
squared-grid for the times $t_{0}\in[$6:30am,9:00am$]$ and
$t_{1}\in[$2:30pm,6:30pm$]$ spaced every $5$ minutes. These times combinations
provides 1440 different scaling options. Then, we applied to each ordered
pair, ($t_{0},t_{1}$), the Levenberg–Marquardt algorithm for square
minimization. We used the python library SciPy (through
“scipy.optimize.curve_fit” Pythons’ function) with the inputs described in
Tab. 3 to apply the algorithm.

Parameter | Description
---|---
$y_{data}:$ |  the measured weight $W^{\prime}$ per minute over a full day
$x_{data}:$ |  the time per minute scaled for the corresponding measured weight $W^{\prime}$
$f:$ |  the model function, introduced in Eq. (1)
$P_{0}:$ |  the initial guess for the $\theta$ parameters that are not bee-facts
| $W_{0}=15000\,g,\,\,N_{max}=15000\,\,\mathrm{bees},\,\,l=17g/h$
|  $m=0.01\,g/bee,a_{1}=0.21\,h^{-1}$, $a_{2}=0.41h^{-1}$, $\,\,d=1h^{-1}$
$Bounds$: | lower and upper bounds on parameters
| $W_{0}=[5000\quad max(W^{\prime}))\,g,\>\>\>N_{max}=(0\quad
80000]\,\,\mathrm{bees},\qquad l=(0\quad max(W^{\prime}))g/s,$
|  $\quad\quad\quad m=(0\quad 0.78w]\,g,\>\>\>\>\>\>\qquad\qquad
a_{1}=[0.10\quad 4.80]\,h^{-1}$,
|  $\qquad\quad a_{2}=[0.10\quad 4.80]h^{-1}$, $d=[0.81\quad 3.01]h^{-1}$
$maxfev:$ | the maximum number of evaluations allowed,
| for this function, 5000 iterations were used
$t_{0},t_{1}:$ |  the clock time corresponding to $t_{0}$ and $t_{1}$, were obtained previously
Table 3: The input required for the “scipy.optimize.curve_fit” Pythons’
function.

We choose the combination of $t_{0}$ and $t_{1}$ that minimizes the error as
the inferred scale time parameters. Then, after re-scaling our data set on the
time axis, we proceeded with the robust estimation of $A$, $B$, $t_{c}$, and
$\alpha$ from the daily-hive data. We infered the remaining seven model
parameters from a system of six equations, which consist of; (i) four
equations presented as the “effective description of the model”. (ii) One
equation from the conservation of the number of departed and returned bees.
(iii) The last equation is the rate of evaporation $\ell$, which can be
calculated directly from the time series after the re-scaling process when
$t<0$. Then, we performed a simultaneous calculation with six equations and
seven variables; consequently, the system is undetermined with no unique
solution due to an extra degree of freedom. We proposed that our dependent
system of equations must be solved for any possible $N_{max}$ value. Note that
the total weight of the hive at $W(t=0)$ can be decomposed in two terms only:
the weight given by the total number of forager bees $N_{max}$, and the
initial constant weight $W_{0}$. Consequently, our solution will depend on
picking a value for $N_{max}$ and finding the corresponding $W_{0}$ value.
Thus, an infinite number of solutions will be obtained and then reduced to a
confidence interval as can be observed in Figure 7 in the Supporting
information. Alternatively, by pre-defining one value, such as departure rate
$d$, we can reduce the variables of our equation’s system and obtain a unique
solution for the parameters $\hat{\theta}$. The solution has extended to 10
independent and different hives during the day 2019-4-7 and is presented in
Sec 2.3.

###  4.2 Pre-processing of data

Before any computation, we removed outliers that do not reflect actual change
in hive weight and that would otherwise affect our analysis. These outliers
appeared due to equipment failure, animal or human interactions with the
hives, or rainfall. In practice, we removed all measurement points that
corresponded to a difference of 20g or more within two consecutive minutes, a
conservative threshold which ensure that strong oscillations were excluded.
From the ten different hives analysed, the hive that presented the highest
number of outliers was hive ‘10’, with 89 outliers, so at most only 6.1 % of
the data were excluded from this hive  (Figure 8).

###  4.3 Bootstrapping

Bootstrapping is used for assessing the effects of data uncertainty. We
applied a bootstrapping method to estimate the uncertainty around the
estimated parameters $\hat{\theta}$. For each daily data of 10 different
hives, $200$ re-samplings were generated, and the minimisation procedure
described above was applied to them (for the fixed value of $t_{0},t_{1}$
obtained in the full dataset). Suppose the inference fails or the number of
iterations exceeds $5,000$ evaluation, and the optimal parameters are still
not found. In that case, the re-sample is discarded (for the 10 analysed hives
listed in Tab. 2 only 5% of the 200 re-samplings were discarded). The inferred
parameters of the accepted re-samples are used to estimate the range of
plausible parameters. The uncertainty reported in the main manuscript
corresponds to the standard deviation over the parameters estimated in the re-
sampled data. For $10$ hives reported in Tab. 2 the uncertainty of the 7
inferred parameters corresponds to an area of 95.4 % under the curve of the
probability density function. Fig. 9 in the Supporting information section
displays an example of the probability density to calculate the range of
values for the seven parameters of ‘Hive 6’ on 2018-4-7.

## Supporting Information

#### Figure 7:

Multiple solutions are obtained from the Hive 6 on $2018-4-7$, and then
reduced to a confidence interval.

#### Figure 8:

Example of outlier detection for the Hive 10 on $2018-4-7$

#### Figure 9:

Probability Density Function for the parameters displayed at Table 1 after the
bootstrapping process of the Hive 6 on $2018-4-7$.

#### Data reporting:

Data is available at the Supplementary Information of Ref. [34] and in our
repository [27].

#### Codes:

The implementation of the proposed model and data analysis methods presented
in this paper are available in the GitHub repository (Click Here), and in the
Zenodo repository (Click Here) [27].

Fig 7: Results obtained from the ‘Hive 6’ the day $2018-4-7$ (a) The time
spent inside, $\tau_{d}$, and outside the hive during the day, $\tau_{a_{1}}$
is compared with reported intervals of time inside the hive taking by the
forager bees between trips, and an average of time outside the hive per
flying. Then, a valid interval is selected considering previous reports. (b)
The maximum reported food load collected per trip (crop load) is $0.077$ g.
The shaded area represents a region where an acceptable crop load is obtained.
(c) Representation of two characteristic parameters, the rate of gained weight
and the intersection with $t=0$, as $N_{max}$ decreases $B$ and $\alpha$
become unstable. (d) Error calculated as the sum of the squared difference
between the model and the real data. The shaded area represents a region where
an acceptable error occurs based on bootstrapping simulation.  Fig 8: Example
of outlier detection for the ‘Hive 10’ on 2018 − 4 − 7. The outliers detected
are in red color and were removed before the data analysis.

Fig 9:  PDF for the fitting parameters: $N_{max}$[bees], $W_{0}$[kg], $m$[g],
$a_{1}$[$h^{-1}$], $a_{2}$[$h^{-1}$], $d$[$h^{-1}$], $l$[g/h] with
$t_{0}=8.4$h and $t_{1}=17.93$h after bootstrapping process of ‘Hive 6’ during
the 2018-4-7.

## Author Contribution

Conceptualization: KAC, TC, MM, TL, EGA, Experiments: TC, Data Analysis: KAC,
TC, Mathematical Model: KAC, TC, EGA, Coding: KAC, TC, Writing – Original
Draft Preparation: KAC, TC, EGA, and Writing – Review & Editing: KAC, TC, MM,
TL, EGA

## Acknowledgement

We acknowledge the Australian Research Council grant DP190101994. T.C.
acknowledges funding from the Lord Mayors Charitable Foundation and Eldon and
Anne Foote Trust. K.A.C. thanks The Sydney Informatics Hub at The University
of Sydney for providing access to HPC-Artemis for data processing. We thank
David Riley for helpful comments.

## References

  * 1. Goulson D, Nicholls E, Botías C, Rotheray EL.  Combined stress from parasites, pesticides and lack of flowers drives bee declines.  Science. 2015;347(6229):1255957.
  * 2. Goulson D, Nicholls E, Botías C, Rotheray EL.  Bee declines driven by combined stress from parasites, pesticides, and lack of flowers.  Science. 2015;347(6229).
  * 3. Hladik ML, Main AR, Goulson D. Environmental risks and challenges associated with neonicotinoid insecticides; 2018.
  * 4. Siviter H, Bailes EJ, Martin CD, Oliver TR, Koricheva J, Leadbeater E, et al.  Agrochemicals interact synergistically to increase bee mortality.  Nature. 2021;596(7872):389–392.
  * 5. Le Conte Y, Ellis M, Ritter W.  Varroa mites and honey bee health: can Varroa explain part of the colony losses?  Apidologie. 2010;41(3):353–363.
  * 6. Powney GD, Carvell C, Edwards M, Morris RK, Roy HE, Woodcock BA, et al.  Widespread losses of pollinating insects in Britain.  Nature communications. 2019;10(1):1–6.
  * 7. Woodcock BA, Bullock J, Shore R, Heard M, Pereira M, Redhead J, et al.  Country-specific effects of neonicotinoid pesticides on honey bees and wild bees.  Science. 2017;356(6345):1393–1395.
  * 8. Corbet SA, Williams IH, Osborne JL.  Bees and the pollination of crops and wild flowers in the European Community.  Bee world. 1991;72(2):47–59.
  * 9. Lee KV, Steinhauer N, Rennich K, Wilson ME, Tarpy DR, Caron DM, et al.  A national survey of managed honey bee 2013–2014 annual colony losses in the USA.  Apidologie. 2015;46(3):292–305.
  * 10. Currie RW, Pernal SF, Guzmán-Novoa E.  Honey bee colony losses in Canada.  Journal of Apicultural Research. 2010;49(1):104–106.
  * 11. Jacques A, Laurent M, Consortium E, Ribière-Chabert M, Saussac M, Bougeard S, et al.  A pan-European epidemiological study reveals honey bee colony survival depends on beekeeper education and disease control.  PLoS one. 2017;12(3):e0172591.
  * 12. Meikle W, Holst N.  Application of continuous monitoring of honeybee colonies.  Apidologie. 2015;46(1):10–22.
  * 13. Francis RM, Nielsen SL, Kryger P.  Varroa-virus interaction in collapsing honey bee colonies.  PloS one. 2013;8(3):e57540.
  * 14. Morawetz L, Köglberger H, Griesbacher A, Derakhshifar I, Crailsheim K, Brodschneider R, et al.  Health status of honey bee colonies (Apis mellifera) and disease-related risk factors for colony losses in Austria.  PloS one. 2019;14(7):e0219293.
  * 15. Horn J, Becher MA, Johst K, Kennedy PJ, Osborne JL, Radchuk V, et al.  Honey bee colony performance affected by crop diversity and farmland structure: a modeling framework.  Ecological Applications. 2021;31(1):e02216.
  * 16. Charreton M, Decourtye A, Henry M, Rodet G, Sandoz JC, Charnet P, et al.  A locomotor deficit induced by sublethal doses of pyrethroid and neonicotinoid insecticides in the honeybee Apis mellifera.  PloS one. 2015;10(12):e0144879.
  * 17. Leoncini I, Crauser D, Robinson G, Le Conte Y.  Worker-worker inhibition of honey bee behavioural development independent of queen and brood.  Insectes Sociaux. 2004;51(4):392–394.
  * 18. Seeley TD, Kolmes SA.  Age polyethism for hive duties in honey bees—illusion or reality?  Ethology. 1991;87(3-4):284–297.
  * 19. Robinson GE.  Regulation of honey bee age polyethism by juvenile hormone.  Behavioral ecology and sociobiology. 1987;20(5):329–338.
  * 20. Perry CJ, Søvik E, Myerscough MR, Barron AB.  Rapid behavioral maturation accelerates failure of stressed honey bee colonies.  Proceedings of the National Academy of Sciences. 2015;112(11):3427–3432.
  * 21. Henry M, Beguin M, Requier F, Rollin O, Odoux JF, Aupinel P, et al.  A common pesticide decreases foraging success and survival in honey bees.  Science. 2012;336(6079):348–350.
  * 22. Quinlan GM, Sponsler D, Gaines-Day HR, McMinn-Sauder HB, Otto CR, Smart AH, et al.  Grassy–herbaceous land moderates regional climate effects on honey bee colonies in the Northcentral US.  Environmental Research Letters. 2022;.
  * 23. Meikle W, Weiss M, Stilwell A.  Monitoring colony phenology using within-day variability in continuous weight and temperature of honey bee hives.  Apidologie. 2016;47(1):1–14.
  * 24. Meikle WG, Holst N, Colin T, Weiss M, Carroll MJ, McFrederick QS, et al.  Using within-day hive weight changes to measure environmental effects on honey bee colonies.  PloS one. 2018;13(5):e0197589.
  * 25. Colin T, Meikle WG, Paten AM, Barron AB.  Long-term dynamics of honey bee colonies following exposure to chemical stress.  Science of the Total Environment. 2019;677:660–670.
  * 26. Meikle WG, Rector BG, Mercadier G, Holst N.  Within-day variation in continuous hive weight data as a measure of honey bee colony activity.  Apidologie. 2008;39(6):694–707.
  * 27. Arias-Calluari K, Altmann EG. Codes for modelling weight of bee hive; 2022.  Available from: https://doi.org/10.5281/zenodo.7272480.
  * 28. Chalcoff VR, Aizen MA, Galetto L.  Nectar concentration and composition of 26 species from the temperate forest of South America.  Annals of botany. 2006;97(3):413–421.
  * 29. Thompson H.  Extrapolation of acute toxicity across bee species.  Integrated environmental assessment and management. 2016;12(4):622–626.
  * 30. Thompson HM, Pamminger T.  Are honeybees suitable surrogates for use in pesticide risk assessment for non-Apis bees?  Pest management science. 2019;75(10):2549–2557.
  * 31. García-García MC, Ortiz PL, Dapena MJD.  Variations in the weights of pollen loads collected by Apis mellifera L.  Grana. 2004;43(3):183–192.
  * 32. Visscher PK, Crailsheim K, Sherman G.  How do honey bees (Apis mellifera) fuel their water foraging flights?  Journal of Insect Physiology. 1996;42(11-12):1089–1094.
  * 33. Afik O, Dag A, Shafir S.  The effect of avocado (Persea americana) nectar composition on its attractiveness to honey bees (Apis mellifera).  Apidologie. 2006;37(3):317–325.
  * 34. Colin T, Forster CC, Westacott J, Wu X, Meikle WG, Barron AB.  Effects of late miticide treatments on foraging and colony productivity of European honey bees (Apis mellifera).  Apidologie. 2021;52(2):474–492.
  * 35. He X, Wang W, Qin Q, Zeng Z, Zhang S, Barron AB.  Assessment of flight activity and homing ability in Asian and European honey bee species, Apis cerana and Apis mellifera, measured with radio frequency tags.  Apidologie. 2013;44(1):38–51.
  * 36. Rodney S, Purdy J.  Dietary requirements of individual nectar foragers, and colony-level pollen and nectar consumption: a review to support pesticide exposure assessment for honey bees.  Apidologie. 2020;51(2):163–179.
  * 37. Colin T, Warren RJ, Quarrell SR, Allen GR, Barron AB.  Evaluating the foraging performance of individual honey bees in different environments with automated field RFID systems.  Ecosphere. 2022;13(5):e4088.
  * 38. Clarke D, Robert D.  Predictive modelling of honey bee foraging activity using local weather conditions.  Apidologie. 2018;49(3):386–396.
  * 39. Hambleton JI.  The effect of weather upon the change in weight of a colony of bees during the honey flow.  1339\. US Department of Agriculture; 1925.
  * 40. Holst N, Meikle WG.  Breakfast canyon discovered in honeybee hive weight curves.  Insects. 2018;9(4):176.
  * 41. Van Der Steen J.  The foraging honey bee.  BBKA News-The British Bee Journal. 2015;2015(February):43–46.
  * 42. Rodney S, Kramer VJ.  Probabilistic assessment of nectar requirements for nectar-foraging honey bees.  Apidologie. 2020;51(2):180–200.
  * 43. Cranmer K, Brehmer J, Louppe G.  The frontier of simulation-based inference.  PNAS. 2020;doi:10.1073/pnas.1912789117.
  * 44. Roda WC.  Bayesian inference for dynamical systems.  Infectious Disease Modelling. 2020;5:221–232.  doi:10.1016/j.idm.2019.12.007.
  * 45. Barnes CP, Silk D, Sheng X, Stumpf MP.  Bayesian design of synthetic biological systems.  Proceedings of the National Academy of Sciences. 2011;108(37):15190–15195.
  * 46. Colin T, Meikle WG, Wu X, Barron AB.  Traces of a neonicotinoid induce precocious foraging and reduce foraging performance in honey bees.  Environmental science & technology. 2019;53(14):8252–8261.


HP2 | 2.2249 |  R20ee$e$Value was not explicitly included in the text but appeared in the RV phase-folded plot.
HT1 | $2.2244\pm 0.0010$ | C20
HT2 | $2.2225\pm 0.0044$ | C20
G1 + G2 | $2.2246\pm 0.0219$ |  C20ff$f$A narrowly constrained Gaussian prior of $2.225\pm 0.020$ d was used to detect a 3$\sigma$ upper limit of $K=23$ m s-1.
HARPS + HIRES | $2.22567^{+0.00026}_{-0.00011}$ | T18
HARPS | 2.22704 |  Rei2013gg$g$Value based off of the most prominent peak in the periodogram.
HARPS | $2.2267\pm 0.0001$ | Bon2013
Tomography
ESPaDOnS + NARVAL | $2.2399\pm 0.0006$ | Mor2008
252525 $b$$b$footnotetext: See Sect. 2.2 for the analysis setup.

Bon2013: Bonfils et al. (2013); Eng2009: Engle et al. (2009); Hun2012: Hunt-
Walker et al. (2012); Mor2008: Morin et al. (2008); Rei2013: Reiners et al.
(2013).

Table 5: CARMENES VIS RV data of AD Leo.

BJD (TDB**$*$Barycentric dynamical time.) | RV (m s-1) |  $\sigma_{\textnormal{RV}}$ (m s-1)
---|---|---
2458199.31070 | $17.26$ | $3.18$
2458199.43446 | $-0.85$ | $2.68$
2458199.53410 | $-7.02$ | $2.70$
2458200.30763 | $-13.30$ | $1.61$
2458200.33569 | $-12.72$ | $1.91$
2458200.52737 | $-0.81$ | $2.85$
2458205.33771 | $18.22$ | $1.62$
2458205.43139 | $23.51$ | $1.59$
2458205.57006 | $24.40$ | $2.42$
2458209.35609 | $-6.66$ | $1.37$
2458209.47549 | $3.29$ | $1.75$
2458210.31450 | $21.81$ | $7.25$
2458211.33225 | $-15.82$ | $4.55$
2458211.33586 | $-12.85$ | $3.16$
2458212.36090 | $32.09$ | $2.45$
2458212.56685 | $25.08$ | $2.78$
2458213.41097 | $-23.28$ | $2.34$
2458213.51416 | $-22.32$ | $2.31$
2458215.31521 | $-22.16$ | $2.03$
2458215.44196 | $-31.14$ | $2.08$
2458215.54733 | $-25.31$ | $2.17$
2458216.36673 | $5.91$ | $3.27$
2458217.35751 | $-7.13$ | $3.32$
2458218.35496 | $-8.17$ | $4.46$
2458225.34239 | $9.43$ | $1.61$
2458225.57485 | $24.92$ | $2.99$
2458893.39518 | $21.10$ | $2.58$
2458893.39739 | $18.67$ | $2.96$
2458893.56436 | $18.29$ | $2.49$
2458893.56616 | $16.85$ | $2.35$
2458893.69412 | $14.18$ | $2.18$
2458893.69600 | $14.00$ | $1.91$
2458894.39429 | $-14.04$ | $2.12$
2458894.39632 | $-9.74$ | $2.55$
2458894.40312 | $-12.89$ | $2.22$
2458894.40435 | $-12.00$ | $2.12$
2458894.55514 | $-15.42$ | $1.90$
2458894.55639 | $-17.37$ | $1.86$
2458894.70318 | $-16.27$ | $2.50$
2458894.70436 | $-15.28$ | $2.15$
2458895.38983 | $11.95$ | $2.90$
2458895.39120 | $7.18$ | $2.46$
2458895.54067 | $14.16$ | $2.47$
2458895.54207 | $16.96$ | $2.62$
2458895.71912 | $21.55$ | $3.84$
2458895.72044 | $16.96$ | $3.15$
262626

Table 6: CARMENES NIR RV data of AD Leo.

BJD (TDB**$*$Barycentric dynamical time.) | RV (m s-1) |  $\sigma_{\textnormal{RV}}$ (m s-1)
---|---|---
2458199.31074 | $-9.74$ | $9.77$
2458199.43432 | $-1.48$ | $9.80$
2458199.53395 | $-21.53$ | $7.61$
2458200.30763 | $-16.58$ | $5.70$
2458200.33556 | $-24.46$ | $6.47$
2458200.52732 | $5.43$ | $24.79$
2458205.33748 | $13.42$ | $6.57$
2458205.43117 | $27.81$ | $5.84$
2458205.57007 | $9.75$ | $5.85$
2458209.35584 | $-5.05$ | $5.92$
2458209.47579 | $-11.00$ | $7.33$
2458210.31440 | $-5.64$ | $23.89$
2458211.33294 | $-24.69$ | $13.57$
2458211.33596 | $-21.72$ | $9.17$
2458212.36111 | $26.48$ | $5.11$
2458212.56689 | $7.25$ | $7.70$
2458213.41095 | $-14.93$ | $6.10$
2458213.51413 | $-16.77$ | $5.91$
2458215.31521 | $-7.97$ | $4.84$
2458215.44212 | $-6.86$ | $5.62$
2458215.54747 | $-11.64$ | $7.33$
2458216.36658 | $-7.24$ | $9.65$
2458217.35744 | $6.59$ | $23.42$
2458218.35494 | $-9.44$ | $11.83$
2458225.34234 | $10.32$ | $7.20$
2458225.57469 | $-1.85$ | $10.24$
2458893.39513 | $17.18$ | $5.81$
2458893.39739 | $31.12$ | $12.16$
2458893.56431 | $8.59$ | $7.67$
2458893.56616 | $6.39$ | $10.36$
2458893.69395 | $-2.56$ | $11.05$
2458893.69590 | $15.99$ | $13.10$
2458894.39425 | $-1.52$ | $7.95$
2458894.39614 | $-3.07$ | $12.02$
2458894.40311 | $15.02$ | $8.90$
2458894.40428 | $0.83$ | $6.93$
2458894.55510 | $-2.68$ | $7.32$
2458894.55633 | $-13.38$ | $7.02$
2458894.70320 | $-27.15$ | $11.60$
2458894.70442 | $-15.44$ | $39.11$
2458895.38980 | $18.25$ | $6.22$
2458895.39117 | $2.50$ | $8.46$
2458895.54064 | $10.29$ | $7.53$
2458895.54203 | $-1.03$ | $10.48$
2458895.71911 | $13.87$ | $14.45$
2458895.72021 | $5.28$ | $18.57$
272727

Table 7: HARPS RV and accompanying data for AD Leo used in this paper, first
processed using serval (Zechmeister et al. 2018) and then corrected for
nightly zero points (Trifonov et al. 2020).

BJD (TDB**$*$Barycentric dynamical time.) | RV (m s-1) |  $\sigma_{\textnormal{RV}}$ (m s-1)
---|---|---
2452986.859 | $28.07$ | $1.89$
2453511.548 | $-19.61$ | $1.01$
2453520.521 | $-15.62$ | $1.00$
2453543.482 | $13.69$ | $1.26$
2453544.453 | $-9.50$ | $0.92$
2453550.460 | $39.00$ | $2.32$
2453728.866 | $39.05$ | $1.03$
2453758.755 | $-21.83$ | $0.74$
2453760.755 | $-12.09$ | $0.76$
2453761.781 | $18.34$ | $0.67$
2453783.726 | $-7.40$ | $0.70$
2453785.727 | $-18.55$ | $0.69$
2453809.661 | $-2.50$ | $0.56$
2453810.677 | $10.93$ | $0.68$
2453811.676 | $8.04$ | $0.66$
2453812.664 | $-9.53$ | $0.96$
2453813.659 | $20.68$ | $0.66$
2453814.655 | $-20.46$ | $0.85$
2453815.571 | $26.14$ | $0.72$
2453815.622 | $25.40$ | $0.66$
2453815.736 | $24.43$ | $0.63$
2453816.544 | $-14.43$ | $0.74$
2453816.656 | $-19.40$ | $0.65$
2453816.722 | $-21.99$ | $0.68$
2453817.551 | $24.31$ | $0.63$
2453817.676 | $26.22$ | $0.67$
2453829.615 | $-1.07$ | $0.62$
2453830.540 | $-7.15$ | $0.79$
2453831.670 | $8.61$ | $0.57$
2453832.650 | $-11.77$ | $0.62$
2453833.628 | $23.40$ | $0.57$
2453834.611 | $-27.06$ | $0.63$
2453835.656 | $27.56$ | $0.63$
2453836.617 | $-19.54$ | $0.59$
2453861.594 | $-1.22$ | $0.69$
2453863.569 | $-22.08$ | $0.69$
2453864.535 | $29.03$ | $0.66$
2453867.542 | $-13.93$ | $0.69$
2453868.518 | $22.20$ | $0.70$
2453871.563 | $19.01$ | $0.62$
2456656.850 | $-36.79$ | $0.87$
2456656.861 | $-35.19$ | $0.86$
2456657.853 | $39.91$ | $0.88$
2456658.865 | $-45.31$ | $0.88$
2456658.876 | $-45.60$ | $0.86$
2456659.859 | $43.23$ | $1.17$
2456797.513 | $-29.24$ | $0.57$
282828

Table 8: HIRES RV data for AD Leo used in this paper produced by Tal-Or et al.
(2019).

BJD (TDB**$*$Barycentric dynamical time.) | RV (m s-1) |  $\sigma_{\textnormal{RV}}$ (m s-1)
---|---|---
2452064.85709 | $-4.47$ | $2.24$
2452334.01660 | $28.75$ | $1.89$
2452602.13013 | $-17.43$ | $1.94$
2452652.10228 | $8.91$ | $1.97$
2452653.01007 | $-5.16$ | $1.98$
2452711.93243 | $-4.57$ | $1.80$
2452712.95569 | $-5.29$ | $2.08$
2452804.84095 | $-14.89$ | $2.40$
2452828.76276 | $5.11$ | $2.12$
2453044.88295 | $17.33$ | $2.10$
2453339.12059 | $13.73$ | $2.02$
2453340.13691 | $-24.87$ | $2.03$
2453398.91738 | $16.05$ | $1.94$
2453398.92487 | $12.70$ | $1.82$
2453724.00531 | $22.43$ | $2.11$
2453724.01159 | $20.59$ | $1.88$
2453748.00269 | $-21.62$ | $1.73$
2453749.89640 | $-21.39$ | $1.72$
2453749.90307 | $-24.10$ | $1.83$
2453750.91326 | $15.30$ | $1.55$
2453751.93398 | $-12.65$ | $3.19$
2453753.00706 | $17.96$ | $1.75$
2453753.01294 | $28.22$ | $2.16$
2453753.94933 | $-12.15$ | $1.71$
2453753.95561 | $-9.32$ | $1.73$
2453776.07390 | $-0.55$ | $1.69$
2453776.08012 | $-0.77$ | $1.59$
2453777.06086 | $-22.01$ | $1.65$
2453777.06716 | $-17.44$ | $1.73$
2453777.88577 | $8.69$ | $1.56$
2453778.14548 | $0.71$ | $1.79$
2453778.82453 | $-26.82$ | $1.71$
2453837.86778 | $19.00$ | $1.99$
2453837.87490 | $17.74$ | $1.83$
2453841.86395 | $10.58$ | $1.91$
2453841.87031 | $8.42$ | $1.79$
2454130.10052 | $-36.01$ | $1.90$
2454130.10671 | $-37.32$ | $1.90$
2454130.99591 | $5.96$ | $2.17$
2454490.99221 | $-31.71$ | $1.93$
2455197.95670 | $26.65$ | $2.73$
2455905.11046 | $-31.07$ | $2.00$
2456641.07001 | $-58.13$ | $1.94$
292929

##  Appendix B Aliasing in HIRES data

The periodogram of the HIRES data set as a whole exhibits the strongest peak
at 2.07 d. instead of at 2.23 d (see Fig. 8). This can be understood as the
window function has a peak at 1 month, more specifically, 29.53 d; 2.07 d and
2.23 d are thus an alias pair with respect to this period. As in Dawson &
Fabrycky (2010) and Stock et al. (2020), a dealiasing approach using the
AliasFinder303030https://github.com/JonasKemmer/AliasFinder (Stock & Kemmer
2020) was performed. However, probably due to the fact that the underlying
signal may not be a simple sinusoid, but rather a quasi-periodic signal with
varying amplitudes and phaseshifts, the results were not fully conclusive. To
further investigate, we performed a stacked Bayesian GLS (s-BGLS) periodogram
(see Fig. 8) where the stacking enables us to determine the coherence of a
signal with increasing number of observations (Mortier et al. 2015; Mortier &
Collier Cameron 2017). As Mortier & Collier Cameron (2017), we normalized all
s-BGLS periodograms to their respective minimum values, with the minimum
probability set to 1. We found that until about 35 data points, which
corresponds to the start of the sparse data, the 2.23 d signal seems to be the
most prominent. Afterwards, the alias at 2.07 d gains significance. This
behavior is in accordance to what we expect from an unstable, incoherent
signal: when one signal loses significance, another one gains it (as also
demonstrated Mortier & Collier Cameron 2017, with the Sun as an example). We
conclude that the HIRES data are fully consistent with an RV modulation with a
2.23 d period.

Figure 8: S-BGLS periodograms of the HIRES data of AD Leo, zoomed into two
frequency ranges (left, 2.23 d, and right, 2.07 d) for the HIRES data. The
nominal value $\log{P}$ does not carry significance, rather the relative
values of $\log{P}$ are of importance.

##  Appendix C Stellar activity

Figure 9: Correlation plots with the RVs (left) and GLS periodograms (right)
of the various stellar activity indicators from the CARMENES VIS spectroscopic
data for AD Leo. For the correlation plots, the circles and triangles
represent the first and second subset of CARMENES VIS data, respectively. Data
points are color-coded with the rotation phase. The Pearson-$r$ correlation
coefficient combining both subsets is shown within each panel (R). For the
periodograms, only the first season of the CARMENES VIS data was considered
for plotting. The orange vertical solid and dashed lines represent the
rotation period at $P=2.23$ d and its daily alias at 1.81 d. The horizontal
dotted, dot-dashed, and dashed blue lines represent the 10 %, 1 %, and 0.1 %
FAP levels.

Figure 10: Same as in Fig. 9, but for the CARMENES NIR channel.

Figure 11: Same as in Fig. 9, but for the HARPS instrument.


Next, we have the similar result as Lemma 17

$\left\|\omega_{0}^{t}-\omega_{i}^{t}\right\|\leq\left\|\omega_{0}^{t-1}-\omega_{0}^{t}\right\|+\frac{L_{\Phi}}{\mu_{1}}\left\|\omega_{i}^{t-1}-\omega_{i}^{t}\right\|+\frac{L_{12}}{\mu_{1}}\left\|\psi_{i}^{t}-\psi^{\star}\left(\omega_{i}^{t}\right)\right\|+\frac{L_{12}}{\mu_{1}}\left\|\psi_{i}^{t-1}-\psi^{\star}(\omega_{i}^{t-1})\right\|+\frac{1}{\mu_{1}}\left\|e_{\omega,i}^{t}\right\|+\frac{1}{\mu_{1}}\left\|e_{\omega,i}^{t-1}\right\|.$
(189)

Then we bound $\left\|\nabla\Psi\left(\omega_{0}^{t}\right)\right\|^{2}$ with
the following result similar to Lemma 18

$\begin{split}\left\|\nabla\Phi\left(\omega_{0}^{t}\right)\right\|\leq&L_{\Phi}\left\|\omega^{t}_{0}-\omega^{t-1}_{0}\right\|+\frac{\mu_{1}^{2}+L_{\Phi}^{2}}{N\mu_{1}}\sum_{i=1}^{N}\left\|\omega^{t-1}_{i}-\omega^{t}_{i}\right\|+\frac{L_{12}(L_{\Phi}+\mu_{1})}{N\mu_{1}}\sum_{i=1}^{N}\left\|\psi_{i}^{t-1}-\psi^{\star}(\omega_{i}^{t-1})\right\|\\\
&+\frac{L_{12}L_{\Phi}}{N\mu_{1}}\sum_{i=1}^{N}\left\|\psi_{i}^{t}-\psi^{\star}\left(\omega_{i}^{t}\right)\right\|+\frac{\mu_{1}+L_{\Phi}}{N\mu_{1}}\sum_{i=1}^{N}\left\|e_{\omega,i}^{t}\right\|+\frac{L_{\Phi}}{N\mu_{1}}\sum_{i=1}^{N}\left\|e_{\omega,i}^{t-1}\right\|.\end{split}$
(190)

By applying the Cauchy-Schwarz inequality and summing up $t$ from $1$ to $T$,
we have a similar result as Eqn. (145), with the particular form as

$\begin{split}\sum_{t=1}^{T}\left\|\nabla\Phi\left(\omega_{0}^{t}\right)\right\|^{2}\leq&4L_{\Phi}^{2}\sum_{t}\left\|\omega^{t}_{0}-\omega^{t-1}_{0}\right\|^{2}+\frac{4(\mu_{1}^{2}+L_{\Phi}^{2})^{2}}{N\mu_{1}^{2}}\sum_{i=1}^{N}\sum_{t=1}^{T}\left\|\omega^{t-1}_{i}-\omega^{t}_{i}\right\|^{2}\\\
&+\frac{8L_{12}^{2}\left((L_{\Phi}+\mu_{1})^{2}+L_{\Phi}^{2}\right)}{N\mu_{1}^{2}}\sum_{i=1}^{N}\sum_{t=1}^{T}\left\|\psi_{i}^{t}-\psi^{\star}\left(\omega_{i}^{t}\right)\right\|^{2}+8\frac{(\mu_{1}+2L_{\Phi})^{2}}{\mu_{1}^{2}}T\epsilon.\\\
\end{split}$ (191)

Then by substituting Eqn. (184) into the previous equation, we have the
following result

$\begin{split}\sum_{t=1}^{T}\left\|\nabla\Phi\left(\omega_{0}^{t}\right)\right\|^{2}\leq&\frac{4L^{2}_{\Phi}\mu^{2}_{1}N+8C_{1}L^{2}_{12}\left((L_{\Phi}+\mu_{1})^{2}+L_{\Phi}^{2}\right)}{N\mu_{1}^{2}}\sum_{t=1}^{T}\left\|\omega^{t}_{0}-\omega^{t-1}_{0}\right\|^{2}\\\
&+\frac{4(\mu_{1}^{2}+L_{\Phi}^{2})^{2}}{N\mu_{1}^{2}}\sum_{i=0}^{N}\sum_{t=1}^{T}\left\|\omega^{t-1}_{i}-\omega^{t}_{i}\right\|^{2}\\\
&+\frac{8C_{2}L^{2}_{12}\left((L_{\Phi}+\mu_{1})^{2}+L_{\Phi}^{2}\right)}{N\mu_{1}^{2}}\sum_{i=1}^{N}\left\|\psi_{i}^{0}-\psi^{\star}(\omega_{i}^{0})\right\|^{2}\\\
&+\frac{8C_{3}L^{2}_{12}\left((L_{\Phi}+\mu_{1})^{2}+L_{\Phi}^{2}\right)}{N\mu_{1}^{2}}\sum_{i=1}^{N}\sum_{j\neq
i}\left\|\omega^{0}_{i}-\omega^{0}_{j}\right\|^{2}\\\
&+\left(\frac{8(\mu_{1}+2L_{\Phi})^{2}}{\mu_{1}^{2}}+\frac{8C_{4}L^{2}_{12}\left((L_{\Phi}+\mu_{1})^{2}+L_{\Phi}^{2}\right)}{N\mu_{1}^{2}}\right)T\epsilon.\end{split}$
(192)

Finally, we reach the convergence result, which is the counterpart of Theorem
2 by following similar manipulations in our final proof of Theorem
2(substituting the above Eqn into Eqn (188)), which is given by

$\Phi(\omega_{0}^{0})-\Psi(\omega_{0}^{T})\leq-
E_{1}\sum_{t=1}^{T}\left\|\nabla\Phi\left(\omega_{0}^{t}\right)\right\|^{2}+E_{2}\sum_{i=1}^{N}\left\|\psi_{i}^{0}-\psi^{\star}(\omega_{i}^{0})\right\|^{2}+E_{3}\sum_{i=1}^{N}\sum_{j\neq
i}\left\|\omega^{0}_{i}-\omega^{0}_{j}\right\|^{2}+E_{4}T\epsilon.$ (193)

where

$\displaystyle E_{1}$
$\displaystyle=\min\left\\{\frac{D_{3}}{D_{1}},\frac{D_{4}}{D_{2}}\right\\},$
(194) $\displaystyle E_{2}$
$\displaystyle=\frac{8E_{1}C_{2}\left((L_{\Phi}+\mu_{1})^{2}+L_{\Phi}^{2}\right)L_{12}^{2}}{N\mu_{1}^{2}}+\frac{C_{2}(5\mu_{1}+32L_{\Phi})L_{12}^{2}}{NL_{\Phi}\mu_{1}},$
(195) $\displaystyle E_{3}$
$\displaystyle=\frac{8E_{1}C_{3}\left((L_{\Phi}+\mu_{1})^{2}+L_{\Phi}^{2}\right)L_{12}^{2}}{N\mu_{1}^{2}}+\frac{C_{3}(5\mu_{1}+32L_{\Phi})L_{12}^{2}}{NL_{\Phi}\mu_{1}},$
(196) $\displaystyle E_{4}$
$\displaystyle=\frac{8NE_{1}(\mu_{1}+2L_{\Phi})^{2}+8E_{1}C_{4}L^{2}_{12}\left((L_{\Phi}+\mu_{1})^{2}+L_{\Phi}^{2}\right)}{N\mu_{1}^{2}}+C_{4}(2+\frac{8}{\mu_{1}}),$
(197) $\displaystyle D_{1}$
$\displaystyle=\frac{4L^{2}_{\Phi}\mu^{2}_{1}N+4C_{1}L^{2}_{12}(L_{\Phi}+\mu_{1})^{2}}{N\mu_{1}^{2}},$
(198) $\displaystyle D_{2}$
$\displaystyle=\frac{8C_{4}L^{2}_{12}\left((L_{\Phi}+\mu_{1})^{2}+L_{\Phi}^{2}\right)}{N\mu_{1}^{2}},$
(199) $\displaystyle D_{3}$
$\displaystyle=\frac{\mu_{1}^{2}-2L_{\Phi}\mu_{1}-4L_{\Phi}^{2}}{2\mu_{1}N}-\frac{C_{1}^{2}(5\mu_{1}+32L_{\Phi})L_{12}^{2}}{NL_{\Phi}\mu_{1}},$
(200) $\displaystyle D_{4}$ $\displaystyle=\frac{\mu_{1}}{4}.$ (201)

Similarly, since $\Phi(\omega_{0}^{0})-\Psi(\omega_{0}^{T})$ is lower-bounded
by a constant $C$. Rearranging terms on the RHS of (193), we have

$\sum_{t=1}^{T}\left\|\nabla\Phi\left(\omega_{0}^{t}\right)\right\|^{2}\leq\frac{E_{2}\sum_{i=1}^{N}\left\|\psi_{i}^{0}-\psi^{\star}(\omega_{i}^{0})\right\|^{2}+E_{3}\sum_{i=1}^{N}\sum_{j\neq
i}\left\|\omega^{0}_{i}-\omega^{0}_{j}\right\|^{2}-C}{E_{1}}+\frac{E_{4}T\epsilon}{E_{1}}.$

Dividing both sides by $T$ and taking $\limsup_{T\to\infty}$, we obtain

$\displaystyle\limsup_{T\to\infty}\frac{\sum_{t=1}^{T}\left\|\nabla\Phi\left(\omega_{0}^{t}\right)\right\|^{2}}{T}$
$\displaystyle\leq\limsup_{T\to\infty}\frac{E_{2}\sum_{i=1}^{N}\left\|\psi_{i}^{0}-\psi^{\star}(\omega_{i}^{0})\right\|^{2}+E_{3}\sum_{i=1}^{N}\sum_{j\neq
i}\left\|\omega^{0}_{i}-\omega^{0}_{j}\right\|^{2}-C}{TE_{1}}+\frac{E_{4}\epsilon}{E_{1}}$
$\displaystyle=\frac{E_{4}\epsilon}{E_{1}},$

which implies that
$\sum_{t=1}^{T}\left\|\nabla\Phi\left(\omega_{0}^{t}\right)\right\|^{2}=O(T\epsilon)$,
and for sufficiently large $t$,
$\left\|\nabla\Phi\left(\omega_{0}^{t}\right)\right\|^{2}=O(\epsilon)$,
completing the proof.


##  A Further examples