Comment on “An approach to the modeling of honey bee colonies” by Romero-Leiton et al. (2022)

The results of a recent paper on the modeling of a honey bee colony have been supplemented and strengthened by showing that the domain of biological validity of a mathematical model with unstable singularity removed is contained in the domain of attraction of positive stable equilibrium. The results go beyond these of the commented paper.

1 Introduction

Due to ecological role of honey bees, it is important to more accurately assess the stability of honey bee colonies. Recently, Romero-Leiton et al. (2022) (Sect. 3) proposed a mathematical model of honey bee colony dynamics in the form of a 3-dimensional nonlinear autonomous system. The asymptotic stability analysis has been local and examined by the linearization.

Observe that their model can be decomposed into the population part:

$$\begin{array}{}\text{(1)}& \dot{B}=\mathit{\beta }f\left(T\right)-bB,\dot{T}=\mathit{\omega }B-aT\end{array}$$

and the equation of honey production M is as follows:

$$\begin{array}{}\text{(2)}& \dot{M}=\mathit{\rho }{\displaystyle \frac{T}{T+u}}-\mathit{\alpha }M-\mathit{\delta }TM,\end{array}$$

which significantly facilitates an analysis. Equation (1) is classified as the age/structure model (Chen et al., 2021).

B and T are populations of immature and adult honey bees, respectively; β stands for immature bee reproduction rate; ω is adult bee maturation rate; $a:=\mathit{\sigma }+{\mathit{\mu }}_T$, where σ is adult bee death rate from a stressful factor; μ_T is adult bee natural death rate; $b:={\mathit{\mu }}_B+\mathit{\omega }$, and μ_B denotes immature bee natural death rate. The eclosion function f has the form of Hille's sigmoidal function of the first order $f\left(y\right)=\frac{y}{y+\mathit{\nu }}$, where ν stands for an average saturate rate (number of adult bees needed for immature bees to reach half of its maximum number). Notice that f is strictly increasing.

M is honey production; ρ is honey production rate; α denotes rate of honey loss due to natural causes; δ stands for honey bee consumption rate in adult state, and u is an average saturate rate.

Positive quadrant ℝ²⁺ (octant ℝ³⁺) is the domain of biological validity of Eq. (1) (Eqs. 1 and 2).

The system Eq. (1) has always an equilibrium at the origin and a desired positive equilibrium (B₀,T₀), and clearly

$$\begin{array}{}\text{(3)}& \mathrm{0}=\mathit{\beta }f\left(T_{\mathrm{0}}\right)-b{B}_{\mathrm{0}},\mathrm{0}=\mathit{\omega }{B}_{\mathrm{0}}-a{T}_{\mathrm{0}}.\end{array}$$

The equilibrium point of Eq. (1) at zero ought to be an unstable saddle point, which is the case if and only if

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\nu }}}>{\displaystyle \frac{ab}{\mathit{\omega }\mathit{\beta }}}:=k,\end{array}$$

where k⁻¹ is called ecological threshold. Then the one-dimensional stable eigenspace of the linearized version of Eq. (1) at (0,0) is the supporting line of ℝ²⁺. By translation of coordinates $x_{\mathrm{1}}:=B-B_{\mathrm{0}}$ and $x_{\mathrm{2}}:=T-T_{\mathrm{0}}$, we reduce Eq. (1) to

$$\begin{array}{}\text{(5)}& {\displaystyle }\begin{array}{l}\begin{array}{l}{\dot{x}}_{\mathrm{1}}=\mathit{\beta }g\left(x_{\mathrm{2}}\right)-b{x}_{\mathrm{1}},{\dot{x}}_{\mathrm{2}}=\mathit{\omega }{x}_{\mathrm{1}}-a{x}_{\mathrm{2}},\\ g\left(y\right):=f(y+T_{\mathrm{0}})-f\left(T_{\mathrm{0}}\right)=\frac{\mathit{\nu }y}{(T_{\mathrm{0}}+\mathit{\nu }+y)(T_{\mathrm{0}}+\mathit{\nu })},\end{array}\end{array}\text{(6)}& {\displaystyle }g\left(\mathrm{0}\right)=\mathrm{0},g(-T_{\mathrm{0}})=-{\displaystyle \frac{T_{\mathrm{0}}}{T_{\mathrm{0}}+\mathit{\nu }}}\stackrel{\left(\text{Eq.}\mathrm{1.3}\right)}{\overbrace{=}}-{\displaystyle \frac{b{B}_{\mathrm{0}}}{\mathit{\beta }}}=-k{T}_{\mathrm{0}}.\end{array}$$

The translation maps ℝ²⁺, the domain of biological validity of the model (Eq. 1), onto $S_{\mathrm{0}}:=[-B_{\mathrm{0}},\mathrm{\infty })\times [-T_{\mathrm{0}},\mathrm{\infty })$.

2 Preliminary analysis

Combining Eq. (5) we establish that x₂ satisfies the Liénard equation (LaSalle and Lefschetz, 1961, p. 59):

$${\ddot{x}}_{\mathrm{2}}=-(a+b){\dot{x}}_{\mathrm{2}}-\mathit{\omega }\mathit{\beta }\left[k{x}_{\mathrm{2}}-g\left(x_{\mathrm{2}}\right)\right].$$

Whence the energetic Lyapunov functional (weakly decreasing along solutions)

$$V={\displaystyle \frac{{\dot{x}}_{\mathrm{2}}^{\mathrm{2}}}{\mathrm{2}\mathit{\omega }\mathit{\beta }}}+\underset{\mathrm{0}}{\overset{x_{\mathrm{2}}}{\int }}\left[ky-g\left(y\right)\right]\mathrm{d}y,\dot{V}=-{\displaystyle \frac{a+b}{\mathit{\omega }\mathit{\beta }}}{\dot{x}}_{\mathrm{2}}^{\mathrm{2}}\le \mathrm{0}$$

is readily available (LaSalle and Lefschetz, 1961, p. 60).

A set is invariant if solutions map this set into itself and strongly invariant if they map this set onto itself.

Lemma 2.1. The sector S₀ is an invariant subset for solutions of Eq. (5).

Proof. With Eq. (6), on a semi–straight line ($x_{\mathrm{1}}=-B_{\mathrm{0}}$, $x_{\mathrm{2}}\ge -T_{\mathrm{0}}$), it holds that

$${\dot{x}}_{\mathrm{1}}=\mathit{\beta }g\left(x_{\mathrm{2}}\right)-b{x}_{\mathrm{1}}\ge \mathit{\beta }g(-T_{\mathrm{0}})+b{B}_{\mathrm{0}}=-{\displaystyle \frac{ab{T}_{\mathrm{0}}}{\mathit{\omega }}}+b{\displaystyle \frac{a{T}_{\mathrm{0}}}{\mathit{\omega }}}=\mathrm{0},$$

and trajectories weakly enter S₀. Similarly on a semi-straight line ($x_{\mathrm{1}}\ge -B_{\mathrm{0}}$, $x_{\mathrm{2}}=-T_{\mathrm{0}}$), one has

$${\dot{x}}_{\mathrm{2}}=\mathit{\omega }{x}_{\mathrm{1}}-a{x}_{\mathrm{2}}\ge -\mathit{\omega }{B}_{\mathrm{0}}+a{T}_{\mathrm{0}}=\mathrm{0},$$

so therein trajectories weakly enter S₀.

3 Main results

Theorem 3.1. Assume that Eq. (4) holds. Every solution of Eq. (5) starting from $S_{\mathrm{0}}\setminus \mathit{\left\{}\right(B_{\mathrm{0}},-T_{\mathrm{0}}\left)\mathit{\right\}}$ tends to zero.

Proof. Thanks to Lemma 2.1 and since V is a Lyapunov functional, for every l≥0, the set $\mathrm{\Omega }:=\mathit{\{}x\in S_{\mathrm{0}}:V(x)\le l\mathit{\}}$ is invariant and compact (closed and bounded). The strongly invariant limit set contained in $\left\{x\in \mathrm{\Omega }:\dot{V}={\dot{x}}_{\mathrm{2}}=\mathit{\omega }{x}_{\mathrm{1}}-a{x}_{\mathrm{2}}=\mathrm{0}\right\}$ consists of two equilibrium points. Indeed ${\dot{x}}_{\mathrm{2}}\equiv \mathrm{0}$ implies ${\ddot{x}}_{\mathrm{2}}\equiv \mathrm{0}$, which with the Liénard equation leads to kx₂≡g(x₂). Now, with Eq. (6), either x₂≡0 or $x_{\mathrm{2}}\equiv -T_{\mathrm{0}}$. In the first case x₁≡0, and in the second one $x_{\mathrm{1}}=-B_{\mathrm{0}}$. The second case is eliminated by excluding the unstable equilibrium point. Modifying the result slightly (LaSalle and Lefschetz, 1961, Theorem VI, p. 58) or applying directly the result of Grabowski (2020, Theorem 1.4), we get the claim.

Importing T=T(t) from Eq. (1), we can regard Eq. (2) as a linear nonautonomous equation.

$$\dot{M}\left(t\right)=-\left[\mathit{\alpha }+\mathit{\delta }T\left(t\right)\right]M+\mathit{\rho }{\displaystyle \frac{T\left(t\right)}{T\left(t\right)+u}},t\ge \mathrm{0}$$

with the solution $M\left(t\right)=M_{\mathrm{1}}\left(t\right)+M_{\mathrm{2}}\left(t\right)$, where

$$\begin{array}{l}{M}_{\mathrm{1}}\left(t\right):=\exp (-{\int }_{\mathrm{0}}^t\left[\mathit{\alpha }+\mathit{\delta }T\left(\mathit{\tau }\right)\right]\mathrm{d}\mathit{\tau })M\left(\mathrm{0}\right),\\ M_{\mathrm{2}}\left(t\right):=\frac{{\int }_{\mathrm{0}}^t\frac{\mathit{\rho }T\left(\mathit{\tau }\right)}{T\left(\mathit{\tau }\right)+u}\exp \left({\int }_{\mathrm{0}}^{\mathit{\tau }}\left[\mathit{\alpha }+\mathit{\delta }T\left(s\right)\right]\mathrm{d}s\right)\mathrm{d}\mathit{\tau }}{\exp \left({\int }_{\mathrm{0}}^t\left[\mathit{\alpha }+\mathit{\delta }T\left(s\right)\right]\mathrm{d}s\right)}.\end{array}$$

M₁(t) exponentially decays to 0. With Theorem 3.1, for $\left(B\right(\mathrm{0}),T(\mathrm{0}\left)\right)\in {\mathbb{R}}^{\mathrm{2}+}\setminus \mathit{\left\{}\mathrm{0}\mathit{\right\}}$, it holds that T(t)⟶T₀ as t→∞. Then, applying L'Hôspital's rule, we get

$$\underset{t\to \mathrm{\infty }}{lim}M\left(t\right)={\displaystyle \frac{\mathit{\rho }{T}_{\mathrm{0}}}{(T_{\mathrm{0}}+u)(\mathit{\alpha }+\mathit{\delta }{T}_{\mathrm{0}})}}:=M_{\mathrm{0}},$$

so M(t)⟶M₀ as t→∞ and M(t)>0, provided that M(0)>0. It is easy to see that M₀ is an equilibrium of Eq. (2). Thus, from Theorem 3.1 we obtain the following corollary.

Theorem 3.2. Assume that Eq. (4) holds. The domain of biological validity ℝ³⁺ of the full model (Eqs. 1 and 2) with unstable equilibrium at 0 removed is the domain of attraction of the positive equilibrium.

4 Conclusions

Theorem 3.2 describes the global behavior of the full system and goes beyond the results of Romero-Leiton et al. (2022, Sect. 3), which were local ones, in a neighborhood of positive equilibrium points. Basic facts are the following: (i) the original model admits a decomposition, and (ii) Eq. (1) is related to the Liénard equation for which (iii) some construction of Lyapunov functionals are known (LaSalle and Lefschetz, 1961, p. 60, p. 115 and Grabowski, 2020, Sect. 2).

Let us indicate ecological aspects of our results:

• i.

  The colony stability is much more robust than it would be following the linearization; i.e., larger deviations of initial conditions from the steady state are allowed (deviations of initial conditions can be interpreted as the Dirac-type impulses additively disturbing the system).

• ii.

  The colony reaches this via the internal stabilizing feedback realized by the eclosion.

• iii.

  The eclosion function may have a different form than Eq. (5), provided that $\frac{g\left(y\right)}{y}<k$ for y>0 and the line ky is not an asymptote for g (this is needed to have level sets of V bounded as required in the proof of Theorem 3.2).