Impact of high wind speed on blooming plants-honeybees-honey production model

Shireen Jawad; Zainab Hayder Abid AL-Aali

doi:10.12688/f1000research.172134.1

Home Browse Impact of high wind speed on blooming plants-honeybees-honey production...

ALL Metrics

-

Views

-

Downloads

Get PDF

Get XML

Export

▬

✚

Research Article

Impact of high wind speed on blooming plants-honeybees-honey production model 

[version 1; peer review: awaiting peer review]

Shireen Jawad ¹, Zainab Hayder Abid AL-Aali¹

PUBLISHED 26 Dec 2025

Author details Author details

¹ Mathematics, University of Baghdad, Baghdad, Baghdad Governorate, Iraq

Shireen Jawad
Roles: Methodology, Supervision, Writing – Review & Editing

Zainab Hayder Abid AL-Aali
Roles: Formal Analysis, Investigation, Writing – Original Draft Preparation

OPEN PEER REVIEW

REVIEWER STATUS AWAITING PEER REVIEW

This article is included in the Fallujah Multidisciplinary Science and Innovation gateway.

Abstract

Background

Local ecosystems and global agriculture are contingent upon the mutualistic relationship between pollinators and floral plants. In symbiosis, pollinators increase agricultural production by improving plant cross-pollination, genetic variety, crop quality, and yield. The potential impact on plant reproduction is particularly alarming due to the decline of pollinating insects. Habitat loss, diseases, climate change, pesticides, and predation have all contributed to the decline of pollinator species. High-speed wind is a significant factor that impacts the mutualistic relationship between plants and pollinators.

Methods

Studying the dynamics of interactions between blooming plants and honeybee populations is crucial for addressing honeybee decline and ensuring sustainable ecosystems. This work employs mathematical modeling to analyze the dynamics of a blooming plant, honeybee population, and honey production symbiosis, with a special emphasis on the effect of high-speed wind flow.

Results

The stability of various ecological equilibria has been investigated using dynamical system theory. Bifurcation phenomena, such as transcritical and Hopf bifurcations, have been discovered using bifurcation theory. Furthermore, the numerical results show that high wind flow can cause the extinction of the honeybee population and honey production.

Conclusions

Due to the rapid depletion of flowering plants and the high rate of wind speed, the populations of honeybees and blossoming plants are at risk of becoming unsustainable. However, the combination of reduced wind flow and increased symbiotic strengths can bolster the stability and sustainability of blooming plant-honeybee-honey production ecosystems. These findings inform conservation policies targeted toward protecting honeybees and increasing biodiversity.

Keywords

Wind speed, blooming plants-honeybees model, mutualistic relationship, Beddington-DeAngelis functional response, dynamical systems, stability analysis, bifurcation.

Corresponding author: Shireen Jawad

Competing interests: No competing interests were disclosed.

Grant information: The author(s) declared that no grants were involved in supporting this work.

Copyright: © 2025 Jawad S and AL-Aali ZHA. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite: Jawad S and AL-Aali ZHA. Impact of high wind speed on blooming plants-honeybees-honey production model  [version 1; peer review: awaiting peer review]. F1000Research 2025, 14:1459 (https://doi.org/10.12688/f1000research.172134.1) First published: 26 Dec 2025, 14:1459 (https://doi.org/10.12688/f1000research.172134.1) Latest published: 26 Dec 2025, 14:1459 (https://doi.org/10.12688/f1000research.172134.1)

Introduction

The mutualistic relationship between blooming plants and honeybee populations is a crucial ecological interaction for the sustainability of both local ecosystems and global agricultural systems.¹ Animal pollinators, including honeybees, provide pollination services through symbiosis, which is essential for the successful reproduction of approximately 300,000 plant species worldwide. In symbiosis, pollinators are vital for promoting plant cross-pollination, genetic diversity, and crop quality and yield, substantially contributing to agricultural productivity.² Consequently, this is a critical research area in conservation biology and ecology. A significant number of studies have been undertaken regarding the symbiotic relationships between plants and pollinator systems.³ Hadani⁴ suggested that the symbiosis can be characterized by the Beddington-DeAngelis response function, which incorporates competition for resource exploitation among pollinators and the obligatory relationship between the plant and the pollinator. This relationship has been observed to maintain a steady state, provided that the initial population level is sufficiently substantial.⁵ Biswas et al. examined a plant–pollinator model to investigate the impact of predation on pollinator species. They conclude that the pollinator is at risk of extinction if the predation rate cannot be controlled. Moreover, this hypothesis has prompted the advancement of extensive studies examining the impacts of nectar theft and ants on the plant-pollinator system.⁶ The reduction of biodiversity is a widespread issue, although the decrease of pollinating insects is especially alarming due to its possible effects on plant reproduction.⁷ A recent report on the global reduction of honeybee and bumblebee populations has highlighted pollination’s ecological and economic significance.⁸ The reduction of pollinator species can be ascribed to various ecological and environmental reasons, including habitat loss, illnesses, climate change, pesticides, and predation.^9–12 Invasive predators can significantly affect pollinators by diminishing their quantity, altering plant reproductive success, and undermining the plant-pollinator relationship.¹³

One factor that negatively affects the mutualism in plant-pollinator systems is high-speed wind.^5,14–17 High winds disrupt flower fragrance messages, reducing honeybee attraction. Strong winds quickly distribute flowers’ aromatic chemicals in different directions, reducing their “scent signal”. As aroma dispersion reduces, honeybees’ ability to discover food sources decreases, resulting in fewer trips to flowers.¹⁸ In addition, flying in severe winds demands more energy for balance and advancement. Honeybees may stay in the hive or fly less to preserve energy. Pollination opportunities decrease as fewer flowers are visited. Strong winds can break flowers’ petals or open their parts abnormally, decreasing honeybees’ access to their reproductive organs. This temporal shift may diminish pollination when flowers are most pollinating. High winds can cause honeybees to crash, fall, or be blown off course, destroying the colony. Lost workers or persistent stress can harm the colony and its capacity to provide enough workers for visits.¹⁹

In this study, we discuss the effect of high wind speed on a three-dimensional blooming plant–honeybee–honey production mathematical model ( $phn$ model) that takes into account saturated mutualism between blooming plants and honeybees by the Beddington–DeAngelis response function. The primary focus is to investigate the dynamics of the plant–pollinator system while considering the impacts of wind speed on the mutualism between blooming plants and honeybees. The study aims to understand the interactions between mutualism and wind speed, and how these interactions affect the overall ecological balance and sustainability of the blooming plant–honeybee–honey production system.

Methods

Assumptions of the model

In this section, a $phn$ model is formulated to describe the interaction among blooming plants $p (t)$ , honey honeybees $h (t)$ and the production of honey $n (t)$ at time $t$ . Then the blooming plant, the honey honeybees, and the production of honey system can be depicted by

(1)

\begin{array}{l} \frac{dp}{dt} = r_{1} p (1 - \frac{p}{k_{1}}) + \frac{α_{1} hp}{1 + w + (ah + bp)} - γ_{1} p, \\ \frac{dh}{dt} = r_{2} h (1 - \frac{h}{k_{2}}) + \frac{α_{2} hp}{1 + w + (ah + bp)} - γ_{2} h, \\ \frac{dn}{dt} = \frac{α_{3} h}{1 + w + ch} - βnh - γ_{3} n . \end{array}

1. The blooming plants are assumed to grow in the absence of honeybees at the intrinsic growth rate $r_{1}$ , depletion rate $γ_{1}$ and carrying capacity $k_{1}$ . Since blooming plants provide nutrients for honeybees, honeybees offer pollination services to blooming plants; hence, their relationship is mutualistic. Beddington-DeAngelis functional response $(\frac{α_{1} hp}{1 + w + (ah + bp)})$ can be used to express the mutualistic relationship, where $α_{1}$ denotes the positive effect of honeybees (a kind of pollinator) on plants, $a$ refers to the undepleted equilibrium rate for the blooming plant–honeybee interaction, which incorporates travel and unloading durations at a central location along with individual-level blooming plant–honeybee interactions, and $b$ indicates the intensity of competition among honeybees for floral resources.
2. High winds break plants’ stems and branches, causing flowers to collapse and damage their petals. This damage reduces the flowers’ attractiveness to pollinating insects. Also, wind speeds reduce the honeybees’ ability to search in a windy environment. Let $ϑ (w) = \frac{1}{1 + w}$ be the efficiency of wind, which satisfied the following
- • $ϑ (0) = 1$ means in the absence of wind, the mutualistic interaction between blooming plants and honeybees remains as before, i.e., $(\frac{α_{1} hp}{1 + (ah + bp)})$ .
- • $ϑ (w) > 1$ means the efficiency of mutualistic interaction decreases in a high, windy environment.
3. The honeybee population are expected to grow from external resources at the intrinsic growth rate $r_{2}$ , depletion rate $γ_{2}$ and carrying capacity $k_{2}$ . $(\frac{α_{2} hp}{1 + w + (ah + bp)})$ represents the nutrients that blooming plants provide to honeybees, where $α_{2}$ stands for the corresponding value of honeybee nutrients from blooming plants.
4. The term $(\frac{α_{3} h}{1 + w + ch})$ represents the honey production in the colonies where $α_{3}$ is the rate of production that is contingent upon the quantity of honeybees, and $c$ is the half-saturation rate. Wind speed negatively affects the amount of honey produced since it diminishes the honeybees’ capacity to search for nutrition.
5. During winter or drought, honeybees rely on their stored honey to survive. This represents a natural consumption process but diminishes the quantity available for harvest. Therefore, $β$ signifies the rate at which honeybees consume honey to survive.
6. The amount of honey produced decreases due to many factors, such as absorbing moisture from the atmosphere. In humid environments, honey absorbs moisture, which reduces its concentration and makes it susceptible to fermentation and spoilage. High temperatures cause the honey content to evaporate. Therefore, $γ_{3}$ denotes the rate of natural causes of honey loss.

Further, the schematic sketch of $phn$ system is illustrated in Figure 1.

Figure 1. Flowchart of the $phn$ model.

Model analysis

Before analyzing our model, it is pertinent to invoke the following lemmas, the demonstration of which is available in Refs. 20-22.

Lemma 1.

If $H, ℊ > 0$ and $\dot{H} ⩾ (⩽) H (H - ℊ H^{α})$ , where $a$ is a positive constant, when $t ⩾ 0$ and $H (0) > 0$ , then

H (t) ⩾ (⩽) {(\frac{H}{ℊ})}^{1 / α} {[1 + (\frac{H H^{- a} (0)}{ℊ} - 1) e^{- H αt}]}^{- 1 / α}

Lemma 2.

(Comparison lemma) Suppose that $H, ℊ > 0,$ with $Ρ (0) > 0$ . Then for $\frac{dΡ}{dt} \leq Ρ (t) [H - ℊ Ρ (t)]$ , then $lim_{t \to \infty} sup P (t) \leq \frac{H}{ℊ}$ , and if $\frac{d Ρ}{dt} \geq Ρ (t) [H - ℊ Ρ (t)]$ , then $lim_{t \to \infty} inf P (t) \geq \frac{H}{ℊ}$ .

Uniqueness.

Since the right side of the $phn$ model $C^{1} (R_{+}^{3})$ , so they satisfy the Lipschitz condition. Therefore, the solution to $phn$ model that starts in $R_{+}^{3}$ exists and is unique.

Positivity

Theorem 1.

The solution $(p (t), h (t), n (t))$ of $phn$ system with the initial condition $(p_{0}, h_{0}, n_{0})$ is positive.

Proof.

According to the blooming plants and honeybees’ equations of $phn$ system, we get

p (t) = p_{0} e^{\int_{0}^{t} (r_{1} (1 - \frac{p (τ)}{k_{1}}) + \frac{α_{1} h (τ)}{1 + w + (ah (τ) + bp (τ))} - γ_{1}) dτ} \geq 0

h (t) = h_{0} e^{\int_{0}^{t} (r_{2} (1 - \frac{h (τ)}{k_{2}}) + \frac{α_{2} p (τ)}{1 + w + (ah (τ) + bp (τ))} - γ_{2}) dτ} \geq 0

From the equation for honey production, we attain

\frac{dn}{dt} \geq - n (β h_{0} e^{\int_{0}^{t} (r_{2} (1 - \frac{h (τ)}{k_{2}}) + \frac{α_{2} p (τ)}{1 + w + (ah (τ) + bp (τ))} - γ_{2}) dτ} + γ_{3}), which implies that n (t) \geq n_{0} e^{\int_{0}^{t} (β h_{0} e^{\int_{0}^{t} (r_{2} (1 - \frac{h (τ)}{k_{2}}) + \frac{α_{2} p (τ)}{1 + w + (ah (τ) + bp (τ))} - γ_{2}) dτ} + γ_{3}) dτ} \geq 0

Therefore, the solution $(p (t), h (t), n (t))$ will remain positive for all $t \geq 0$ .

Boundedness

Theorem 2.

All solutions of $phn$ system are uniformly bounded.

Proof.

From the blooming plants equation of the $phn$ system, we obtain

\frac{dp}{dt} = r_{1} p (1 - \frac{p}{k_{1}}) + \frac{α_{1} hp}{1 + w + (ah + bp)} - γ_{1} p \leq p [(r_{1} - γ_{1}) - (\frac{r_{1} p}{k_{1}}) + α_{1} h]

.

After using the honeybees’ carrying capacity, we get

\frac{dp}{dt} \leq p [(r_{1} + α_{1} k_{2} - γ_{1}) - (\frac{r_{1} p}{k_{1}})]

Using Lemma 1, we deduced that.

p (t) \leq {[\frac{(r_{1} + α_{1} k_{2} - γ_{1}) k_{1}}{r_{1}}]}^{\frac{1}{2}} {(1 + [\frac{(r_{1} + α_{1} k_{2} - γ_{1}) k_{1} p^{- 2} (0)}{r_{1}} - 1] e^{- 2 (r_{1} + α_{1} k_{2} - γ_{1}) t})}^{- \frac{1}{2}}

It is obtained after taking $t \to \infty$

p (t) \leq {[\frac{(r_{1} + α_{1} k_{2} - γ_{1}) k_{1}}{r_{1}}]}^{\frac{1}{2}} = p_{m}

Similarly, for the honey honeybee equation of the $phn$ system, we find

h (t) \leq {[\frac{k_{1} (r_{2} + α_{2} p_{m} - γ_{2})}{r_{2}}]}^{\frac{1}{2}} = h_{m}

Finally, from the amount of honey produced equation, we attain

\frac{dn}{dt} = \frac{α_{3} h}{1 + w + ch} - βnh - γ_{3} n \leq α_{3} h_{m} - γ_{3} n

Then, by Lemma 2, we get

0 \leq lim_{t \to \infty} sup [n (t)] \leq \frac{α_{3} h_{m}}{γ_{3}} = n_{m}

Therefore, any solutions of $phn$ system will be attracted to $Σ = {(h, p, n) \in R_{+}^{3} : p (t) \leq p_{m}, h (t) \leq h_{m}, n (t) \leq n_{m}}$ .

Persistence

The $phn$ system is said to be persistent if all its components survive in future times.

Theorem 3.

All $phn$ system components are persistent if

(2)

r_{1} > γ_{1},

(3)

r_{2} > γ_{2} .

Proof:

First, to prove that $p (t)$ If it is persistent, we have to show that $lim_{t \to \infty} inf p (t) > 0$ i.e. $p (t)$ will not decay to zero.

From the blooming plants equation, we get

\frac{dp}{dt} = r_{1} p (1 - \frac{p}{k_{1}}) + \frac{α_{1} hp}{1 + w + (ah + bp)} - γ_{1} p \geq p [(r_{1} - γ_{1}) - (\frac{r_{1}}{k_{1}}) p] .

By Lemma 2, we have $lim_{t \to \infty} inf p (t) \geq L_{1}$ , where $L_{1} = \frac{k_{1} (r_{1} - γ_{1})}{r_{1}}$ provided $r_{1} > γ_{1}$ . Thus, for a small $E_{1} > 0$ , $\exists$ a positive number $t_{1} > 0$ such that

p (t) \geq L_{1} - E_{1} \forall t > t_{1}

Applying the same strategy for the honeybee equation, we get

h (t) \geq L_{2} - E_{2} \forall t > t_{2},

where

L_{2} = \frac{k_{2} (r_{2} - γ_{2})}{r_{2}}

provided

r_{2} > γ_{2}

.

From the honey production equation, we attain

\frac{dn}{dt} = \frac{α_{3} h}{1 + w + ch} - βnh - γ_{3} n \geq \frac{α_{3} (L_{2} - E_{2})}{1 + w + c (L_{2} - E_{2})} - n (β h_{m} - γ_{3})

Since

• $g (h) = \frac{α_{3} h}{1 + w + ch}$ is an increasing function, which means if $h (t) \geq L_{2} - E_{2}$ , then $g (h) \geq g (L_{2} - E_{2})$
• $h (t) \leq h_{m}$

Then, by Lemma 2, $lim_{t \to \infty} inf n (t) \geq \frac{α_{3} (L_{2} - E_{2})}{[1 + w + c (L_{2} - E_{2})] (β h_{m} - γ_{3})}$ , where $L_{3} = \frac{α_{3} (L_{2} - E_{2})}{[1 + w + c (L_{2} - E_{2})] (β h_{m} - γ_{3})}$

Thus, for a small $E_{3} > 0$ , $\exists$ a positive number $t_{3} > 0$ such that

n (t) \geq L_{3} - E_{3} \forall t > t_{3}

Therefore, $phn$ system is uniformly persistent.

Remark 1:

Conditions 2 and 3 indicate that all $phn$ system components will survive in the future times if the intrinsic growth rates of the blooming plants and honeybees exceed their depletion rates.

Equilibrium analysis

The possible equilibria of $phn$ system are

1) The extinction point $S_{1} = (0, 0, 0)$ .
2) The blooming plants point $S_{2} = (p_{1}, 0, 0)$ , where $p_{1} = \frac{k_{1} (r_{1} - γ_{1})}{r_{1}}$ . For $p_{1}$ to be positive, condition (2) must be satisfied.
3) The honeybee point $S_{3} = (0, h_{2}, 0)$ , where $h_{2} = \frac{k_{2} (r_{2} - γ_{2})}{r_{2}}$ . Clearly $h_{2} > 0$ if condition (3) is satisfied.
4) The blooming plants free point $S_{4} = (0, h_{3}, n_{3})$ , where $h_{3} = \frac{k_{2} (r_{2} - γ_{2})}{r_{2}}$ and $n_{3} = \frac{α_{3} h_{3}}{(1 + w + c h_{3}) (β h_{3} + γ_{3})} > 0$ . Clearly $h_{3} > 0$ if the condition (3) is satisfied.
5) The coexistence point $S_{5} = (p_{4}, h_{4}, n_{4})$ , here $n_{4} = \frac{α_{3} h_{4}}{(1 + w + c h_{4}) (β h_{4} + γ_{3})}$ , $h_{4} = \frac{r_{1} b p_{4}^{2} + z_{1} p_{4} + z_{2}}{z_{3} - a r_{1} p_{4}}$ where,
$z_{1} = r_{1} (1 + w - k_{1} b) + γ_{1} k_{1} b,$
$z_{2} = k_{1} (γ_{1} (1 + w) - (r_{1} + w)),$
$z_{3} = k_{1} (α_{1} + a (r_{1} - γ_{1} - p)),$

and $p_{4}$ is the root of the following equation

(4)

A_{0} p^{3} + A_{1} P^{2} + A_{2} P + A_{3} = 0,

where,

$A_{0} = r_{1} [r_{1} a (r_{2} b + r_{2} bw + α_{2} a k_{2}) - r_{2} b (b z_{3} + a z_{1})],$

$A_{1} = (r_{2} - γ_{2}) [r_{1} ({ak}_{2} (r_{1} + wa r_{1} - b z_{3} - a z_{1})] + r_{2} (a z_{1} + aw z_{1} - bw - bw z_{3} - ab z_{2}) - r_{2} b z_{1} z_{3} - r_{2} a z_{1}^{2}$ $- a r_{1} α_{2} z_{3} k_{2} - a α_{2} r_{1} z_{3} k_{2},$

$A_{2} = (r_{2} - γ_{2}) [z_{3} k_{2} (z_{3} b - r_{1} - 2 r_{1} aw + a z_{1} - r_{1} a) - k_{2} r_{1} a^{2} z_{2}] + r_{2} (z_{2} (a (r_{1} - 2 z_{1} + w r_{1}) - b z_{3})$ $+ z_{3} (k_{2} α_{2} - r_{1} - z_{1} w)),$

$A_{3} = (r_{2} - γ_{2}) [z_{3} (z_{3} - k_{2} (w z_{3} - a z_{2}))] - r_{2} z_{2} (z_{3} + w z_{3} + a z_{2}) .$

By Descartes’s rule of signs, Equation (4) has a unique positive root $p_{4}$ , if one of the following conditions holds:

\begin{array}{l} A_{0} > 0, and A_{2}, A_{3} < 0, \\ A_{0}, A_{1} > 0, and A_{3} < 0, \\ A_{0} < 0, and A_{2}, A_{3} > 0, \\ A_{0}, A_{1} < 0, and A_{2}, A_{3} > 0 . \end{array}

Further, $h_{4} > 0$ if one of the following conditions holds:

\begin{matrix} r_{1} b p_{4}^{2} + z_{1} p_{4} + z_{2} > 0 and z_{3} > a r_{1} p_{4} \\ r_{1} b p_{4}^{2} + z_{1} p_{4} + z_{2} < 0 and z_{3} < a r_{1} p_{4} \end{matrix}

Local stability of equilibrium points

To investigate the local stability, one needs to determine the Jacobian matrix at any point. Thus, the Jacobian matrix at any point $(p, h, n)$ is

J (p, h, n) = [\begin{matrix} a_{11} & \frac{α_{1} p (1 + w + bp)}{{(1 + w + (ah + bp))}^{2}} & 0 \\ \frac{α_{2} h (1 + w + ah)}{{(1 + w + (ah + bp))}^{2}} & a_{22} & 0 \\ 0 & \frac{α_{3} (1 + w)}{{(1 + w + ch)}^{2}} - βn & - βh - γ_{3} \end{matrix}],

where

a_{11} = \frac{- 2 r_{1} p}{k_{1}} + \frac{α_{1} h (1 + w + ah)}{{(1 + w + (ah + bp))}^{2}} + (r_{1} - γ_{1})

, and

a_{22} = \frac{- 2 r_{2} h}{k_{2}} + \frac{α_{2} p (1 + w + bp)}{{(1 + w + (ah + bp))}^{2}} + (r_{2} - γ_{2})

.

Theorem 4.

The extinction point $S_{1} = (0, 0, 0)$ is locally asymptotically stable if

(5)

γ_{1} > r_{1},

(6)

γ_{2} > r_{2} .

Proof:

The Jacobian matrix at $S_{1} = (0, 0, 0)$ is

(7)

J (S_{1}) = [\begin{matrix} r_{1} - γ_{1} & 0 & 0 \\ 0 & r_{2} - γ_{2} & 0 \\ 0 & \frac{\propto_{3}}{1 + w} & - γ_{3} \end{matrix}],

Then, the eigenvalues of $J (S_{1})$ are $λ_{11} = r_{1} - γ_{1}$ , $λ_{12} = r_{2} - γ_{2}$ and $λ_{13} = - γ_{3} < 0$ .

Therefore, $S_{1}$ is asymptotically stable under conditions 5 and 6.

Remark 2:

The biological interpretation of conditions 5 and 6 indicates $phn$ system reaches asymptotically the extinction point when the depletion rates of the blooming plants and honeybees exceed their intrinsic growth rates.

Theorem 5.

The blooming plants point $S_{2} = (p_{1}, 0, 0)$ is locally asymptotically stable if

(8)

r_{2} + \frac{α_{2} p_{1}}{(1 + w + b p_{1})} < γ_{2},

Proof:

$J (S_{2}) = J (p_{1}, 0, 0)$ is given by

(9)

J (S_{2}) = [\begin{matrix} - (r_{1} - γ_{1}) & \frac{α_{1} p_{1}}{(1 + w + b p_{1})} & 0 \\ 0 & \frac{α_{2} p_{1}}{(1 + w + b p_{1})} + (r_{2} - γ_{2}) & 0 \\ 0 & \frac{α_{3}}{1 + w} & - γ_{3} \end{matrix}] .

So, the eigenvalues of $S_{2}$ are

$λ_{21} = - (r_{1} - γ_{1}) < 0$ , under the existence condition of the blooming plants point.

$λ_{22} = \frac{α_{2} p_{1}}{(1 + w + b p_{1})} + (r_{2} - γ_{2}),$ and $λ_{23} = - γ_{3} < 0$ .

Thus, $S_{2}$ is asymptotically stable if condition 8 is satisfied $.$

Theorem 5.

The honeybee point $S_{3} = (0, h_{2}, 0)$ is locally asymptotically stable if

(10)

r_{1} + \frac{α_{1} h_{2}}{(1 + w + a h_{2})} < γ_{1},

(11)

(r_{2} - γ_{2}) < \frac{2 r_{2} h_{2}}{k_{2}}

Proof:

$J (S_{3}) = J (0, h_{2}, 0)$ is given by

(12)

J (S_{3}) = [\begin{matrix} \frac{α_{1} h_{2}}{(1 + w + a h_{2})} + (r_{1} - γ_{1}) & 0 & 0 \\ \frac{α_{2} h_{2}}{1 + w + a h_{2}} & (r_{2} - γ_{2}) - \frac{2 r_{2} h_{2}}{k_{2}} & 0 \\ 0 & \frac{α_{3} (1 + w)}{{(1 + w + c h_{2})}^{2}} & - β h_{2} - γ_{3} \end{matrix}] .

So, the eigenvalues of $J (S_{3})$ are $λ_{31} = \frac{α_{1} h_{2}}{(1 + w + a h_{2})} + (r_{1} - γ_{1}), λ_{32} = (r_{2} - γ_{2}) - \frac{2 r_{2} h_{2}}{k_{2}}$ , and $λ_{33} = - β h_{2} - γ_{3} < 0$ . Thus, $S_{3}$ is asymptotically stable if the conditions 10 and 11 are satisfied.

Theorem 6.

The blooming plants free point $S_{4} = (0, h_{3}, n_{3})$ is locally asymptotically stable if

(13)

r_{1} + \frac{α_{1} h_{3}}{(1 + w + a h_{3})} < γ_{1},

(14)

(r_{2} - γ_{2}) < \frac{r_{2} h_{3}}{k_{2}}

Proof:

$J (S_{4}) = J (0, h_{3}, n_{3})$ is given by

(15)

J (S_{4}) = [\begin{matrix} \frac{α_{1} h_{3}}{(1 + w + a h_{3})} + (r_{1} - γ_{1}) & 0 & 0 \\ \frac{α_{2} h_{3}}{(1 + w + a h_{3})} & (r_{2} - γ_{2}) - \frac{2 r_{2} h_{3}}{k_{2}} & 0 \\ 0 & \frac{α_{3} (1 + w)}{(1 + w + c h_{3})} - β n_{3} & - β h_{3} - γ_{3} \end{matrix}] .

The eigenvalues of $J (S_{4})$ are $λ_{41} = \frac{α_{1} h_{3}}{(1 + w + a h_{3})} + (r_{1} - γ_{1})$ , $λ_{42} = (r_{2} - γ_{2}) - \frac{2 r_{2} h_{3}}{k_{2}}$ , and $λ_{43} = - β h_{3} - γ_{3} < 0$ . So, $S_{4}$ is asymptotically stable if conditions 13 and 14 are satisfied.

Theorem 7.

The coexistence point $S_{5} = (p_{4}, h_{4}, n_{4})$ is locally asymptotically stable if

(16)

a_{ii} < 0, i = 1, 2,

(17)

a_{11} a_{22} > [\frac{α_{1} α_{2} p_{4} h_{4} (1 + w + b p_{4}) (1 + w + a h_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{4}}] .

Proof:

$J (S_{5}) = J (p_{4}, h_{4}, n_{4})$ is given by

(18)

J (S_{5}) = [\begin{matrix} a_{11} & \frac{α_{1} p_{4} (1 + w + b p_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}} & 0 \\ \frac{α_{2} h_{4} (1 + w + a h_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}} & a_{22} & 0 \\ 0 & \frac{α_{3} (1 + w)}{{(1 + w + c h_{4})}^{2}} - β n_{4} & - β h_{4} - γ_{3} \end{matrix}],

where

a_{11} = (r_{1} - γ_{1}) - \frac{2 r_{1} p_{4}}{k_{1}} + \frac{α_{1} h_{4} (1 + w + a h_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}}

, and

a_{22} = (r_{2} - γ_{2}) - \frac{2 r_{2} h_{4}}{k_{2}} + \frac{α_{2} p_{4} (1 + w + b p_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}}

.

Then, the characteristic equation of $J (S_{5})$ is given by:

(19)

(- β h_{4} - γ_{3} - λ) [λ^{2} - Trλ + Det] = 0,

where,

\begin{array}{l} λ_{51} = - β h_{4} - γ_{3}, \\ Tr = a_{11} + a_{22} = (r_{1} - γ_{1}) - \frac{2 r_{1} p_{4}}{k_{1}} + \frac{α_{1} h_{4} (1 + w + a h_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}} - \frac{2 r_{2} h_{4}}{k_{2}} + \frac{α_{2} p_{4} (1 + w + b p_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}} + (r_{2} - γ_{2}), \\ Det = a_{11} a_{22} - [\frac{α_{1} α_{2} p_{4} h_{4} (1 + w + b p_{4}) (1 + w + a h_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{4}}] = ((r_{1} - γ_{1}) - \frac{2 r_{1} p_{4}}{k_{1}} + \frac{α_{1} h_{4} (1 + w + a h_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}}) ((r_{2} - γ_{2}) - \frac{2 r_{2} h_{4}}{k_{2}} + \frac{α_{2} p_{4} (1 + w + b p_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}}) - [\frac{α_{1} α_{2} p_{4} h_{4} (1 + w + b p_{4}) (1 + w + a h_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{4}}] . \end{array}

Thus, $S_{5}$ exhibits local stability if conditions 16 and 17 are fulfilled.

Global stability

In this section, the Lyapunov method is used to illustrate the global stability of the previous points, as shown in the following theorems.

Theorem 8.

The extinction point $S_{1} = (0, 0, 0)$ is a global asymptotic stability (GAS) if the following conditions are met.

(20)

r_{1} + \frac{h_{m} (α_{1} + α_{2})}{1 + w} < γ_{1}

(21)

r_{2} + \frac{α_{3}}{1 + w} < γ_{2}

Proof:

Let $E_{1} (p, h, n) = p + h + n$ , where $E_{1} : R_{+}^{3} \to R$ , which satisfies $E_{1} (0, 0, 0) = 0$ and $E_{1} (p, h, n) > 0$ for all $(p, h, n) \in R_{+}^{3}$ with $(p, h, n) \neq (0, 0, 0)$ . Then

\frac{d E_{1}}{dt} = \frac{dp}{dt} + \frac{dh}{dt} + \frac{dn}{dt} = (r_{1} p (1 - \frac{p}{k_{1}}) + \frac{α_{1} hp}{1 + w + (ah + bp)} - γ_{1} p) + (r_{2} h (1 - \frac{h}{k_{2}}) + \frac{α_{2} hp}{1 + w + (ah + bp)} - γ_{2} h) + (\frac{α_{3} h}{1 + w + ch)} - βnh - γ_{3} n) = r_{1} p - r_{1} \frac{p^{2}}{k_{1}} + \frac{hp (α_{1} + α_{2})}{1 + w + (ah + bp)} - γ_{1} p + r_{2} h - r_{2} \frac{h^{2}}{k_{2}} - γ_{2} h + \frac{α_{3} h}{1 + w + ch)} - βnh - γ_{3} n

Then, by using the upper bound of the honeybees’ population, we get

\frac{d E_{1}}{dt} \leq p (r_{1} - γ_{1} + \frac{h_{m} (α_{1} + α_{2})}{1 + w}) + h (r_{2} - γ_{2} + \frac{α_{3}}{1 + w}) - r_{1} \frac{p^{2}}{k_{1}} - r_{2} \frac{h^{2}}{k_{2}} - βnh - γ_{3} n .

The first two terms are negative definite if conditions 20 and 21 are satisfied. Hence, $d E_{1} / dt$ is a negative definite. Therefore, the extinction point $S_{1}$ is GAS.

Theorem 9.

The blooming plants point $S_{2} = (p_{1}, 0, 0)$ is GAS if

(22)

r_{2} + \frac{α_{2} p_{m} + α_{3}}{1 + w} < γ_{2},

(23)

{(\frac{α_{1}}{D})}^{2} \leq \frac{r_{1} r_{2}}{k_{1} k_{2}},

where,

D = 1 + w + (ah + bp)

.

Proof:

Let $E_{2} = (p - p_{1} - p_{1} ln (\frac{p}{p_{1}})) + h + n$ , where $E_{2} : R_{+}^{3} \to R$ , which satisfies $E_{2} (p_{1}, 0, 0) = 0$ and $E_{2} (p, h, n) > 0$ for all $(p, h, n) \in R_{+}^{3}$ with $(p, h, n) \neq (p_{1}, 0, 0)$ , then

\frac{d E_{2}}{dt} = (\frac{p - p_{1}}{p}) \frac{dp}{dt} + \frac{dh}{dt} + \frac{dn}{dt} = \frac{- r_{1}}{k_{1}} {(p - p_{1})}^{2} + \frac{α_{1} h (p - p_{1})}{D} + h r_{2} (1 - \frac{h}{k_{2}}) + \frac{α_{2} ph}{D} - γ_{2} h + \frac{α_{3} h}{1 + w + ch} - βnh - γ_{3} n

Then, by using the upper bound of the blooming plants population, we get

\frac{d E_{2}}{dt} \leq - (\frac{r_{1}}{k_{1}} {(p - p_{1})}^{2} - \frac{α_{1}}{D} h (p - p_{1}) + \frac{r_{2}}{k_{2}} h^{2}) + h (r_{2} + \frac{α_{2} p_{m} + α_{3}}{1 + w} - γ_{2}) - βnh - γ_{3} n .

Thus,

\frac{d E_{2}}{dt} \leq - {(\sqrt{\frac{r_{1}}{k_{1}}} (p - p_{1}) + \sqrt{\frac{r_{2}}{k_{2}}} h)}^{2} + h (r_{2} + \frac{α_{2} p_{m} + α_{3}}{1 + w} - γ_{2}) - βnh - γ_{3} n .

The first two terms are negative definite if conditions 22 and 23 are satisfied. Hence, $d E_{2} / dt$ is a negative definite. Therefore, the blooming plants point $S_{2} = (p_{1}, 0, 0)$ is GAS.

Theorem 10.

The honeybee point $S_{3} = (0, h_{2}, 0)$ , is GAS if

(24)

r_{1} + \frac{α_{1} h_{m}}{1 + w} < γ_{1},

(25)

{(\frac{α_{2}}{D})}^{2} \leq \frac{r_{1} r_{2}}{k_{1} k_{2}},

(26)

\frac{α_{3}}{β (1 + w)} < n,

where,

D = 1 + w + (ah + bp)

.

Proof:

Let $E_{3} = p + (h - h_{2} - h_{2} ln \frac{h}{h_{2}}) + n$ , where $E_{3} : R_{+}^{3} \to R$ , which satisfies $E_{3} (0, h_{2}, 0) = 0$ and $E_{3} (p, h, n) > 0$ for all $(p, h, n) \in R_{+}^{3}$ with $(p, h, n) \neq (0, h_{2}, 0)$ , then

\frac{d E_{3}}{dt} = \frac{dp}{dt} + (\frac{h - h_{2}}{h}) \frac{dh}{dt} + \frac{dn}{dt} = (r_{1} p (1 - \frac{p}{k_{1}}) + \frac{α_{1} hp}{1 + w + (ah + bp)} - γ_{1} p) + (h - h_{2}) ((r_{2} (1 - \frac{h}{k_{2}}) + \frac{α_{2} p}{1 + w + (ah + bp)} - γ_{2})) + (\frac{α_{3} h}{1 + w + ch} - βnh - γ_{3} n) \leq r_{1} p - \frac{r_{1} p^{2}}{k_{1}} + \frac{α_{1} hp}{1 + w} - γ_{1} p - \frac{r_{2}}{k_{2}} {(h - h_{2})}^{2} + \frac{α_{2} (h - h_{2}) p}{D} + \frac{α_{3} h}{1 + w} - βnh - γ_{3} n

Then, by using the upper bound of the honeybees’ population, we get

\frac{d E_{3}}{dt} \leq - (\frac{r_{1}}{k_{1}} p^{2} - \frac{α_{2} (h - h_{2}) p}{D} + (\frac{r_{2}}{k_{2}}) {(h - h_{2})}^{2}) + p (r_{1} + \frac{α_{1} h_{m}}{1 + w} - γ_{1}) + h (\frac{α_{3}}{1 + w} - βn) - γ_{3} n

Thus,

\frac{d E_{3}}{dt} \leq - {(\sqrt{\frac{r_{1}}{k_{1}} p} + \sqrt{\frac{r_{2}}{k_{2}}} (h - h_{2}))}^{2} + p (r_{1} + \frac{α_{1} h_{m}}{1 + w} - γ_{1}) + h (\frac{α_{3}}{1 + w} - βn) - γ_{3} n .

The first three terms are negative definite if conditions 24-26 are satisfied. Hence, $d E_{3} / dt$ is a negative definite. Therefore, the blooming plants point $S_{3} = (0, h_{2}, 0)$ is GAS.

The blooming plants free point $S_{4} = (0, h_{3}, n_{3})$

Theorem 11.

The blooming plants free point $S_{4} = (0, h_{3}, n_{3})$ is GAS if

(27)

{(\frac{α_{2}}{D})}^{2} \leq \frac{r_{1} r_{2}}{2 k_{1} k_{2}},

(28)

{(\frac{α_{3} (1 + w)}{D_{1} D_{2}} - βn)}^{2} \leq \frac{r_{2} (β h_{3} + γ_{3})}{2 k_{2}},

where,

D = 1 + w + (ah + bp), D_{1} = (1 + w + ch)

and

D_{2} = (1 + w + c h_{3})

.

Proof:

Let $E_{4} = p + (h - h_{3} - h_{3} ln \frac{h}{h_{3}}) + {(\frac{n - n_{3}}{2})}^{2}$ , where $E_{4} : R_{+}^{3} \to R$ , which satisfies $E_{4} (0, h_{3}, n_{3}) = 0$ and $E_{4} (p, h, n) > 0$ for all $(p, h, n) \in R_{+}^{3}$ with $(p, h, n) \neq (0, h_{3}, n_{3})$ , then

\frac{d E_{4}}{dt} = \frac{dp}{dt} + (\frac{h - h_{3}}{h}) \frac{dh}{dt} + (n - n_{3}) \frac{dn}{dt} = (r_{1} p (1 - \frac{p}{k_{1}}) + \frac{α_{1} hp}{1 + w + (ah + bp)} - γ_{1} p) + (h - h_{3}) (r_{2} (1 - \frac{h}{k_{2}}) + \frac{α_{2} p}{1 + w + (ah + bp)} - γ_{2}) + (n - n_{3}) (\frac{α_{3} h}{1 + w + ch} - βnh - γ_{3} n) = r_{1} p - \frac{r_{1} p^{2}}{k_{1}} + \frac{α_{1} hp}{1 + w + (ah + bp)} - γ_{1} p + (h - h_{3}) (- \frac{r_{2}}{k_{2}} (h - h_{3}) + \frac{α_{2} p}{1 + w + (ah + bp)}) + (n - n_{3}) (\frac{α_{3} h}{1 + w + ch} - \frac{α_{3} h_{3}}{1 + w + c h_{3}} - βnh - β n_{3} h_{3} - γ_{3} (n - n_{3}))

Then, by using the upper bound of the honeybees’ population, we get

\frac{d E_{4}}{dt} \leq - (\frac{r_{1}}{k_{1}} p^{2} - \frac{α_{2} (h - h_{3}) p}{1 + w + (ah + bp)} + (\frac{r_{2}}{2 k_{2}}) {(h - h_{3})}^{2}) + p (r_{1} + \frac{α_{1} h_{m}}{1 + w} - γ_{1}) - [(\frac{r_{2}}{2 k_{2}}) {(h - h_{3})}^{2} - (\frac{α_{3} (1 + w)}{(1 + w + ch) (1 + w + c h_{3})} - βn) (h - h_{3}) (n - n_{3}) + (β h_{3} + γ_{3}) {(n - n_{3})}^{2}]

Thus,

\frac{d E_{4}}{dt} \leq - {(\sqrt{\frac{r_{1}}{k_{1}} p} + \sqrt{\frac{r_{2}}{2 k_{2}}} (h - h_{3}))}^{2} + p (r_{1} + \frac{α_{1} h_{m}}{1 + w} - γ_{1}) - {(\sqrt{\frac{r_{2}}{2 k_{2}}} (h - h_{3}) + \sqrt{(β h_{3} + γ_{3})} (n - n_{3}))}^{2}

The first and the third terms are negative definite if conditions 27 and 28 are satisfied, while the second term is negative under condition 24. Hence, $d E_{4} / dt$ is a negative definite. Therefore, the blooming plants free point $S_{4} = (0, h_{3}, n_{3})$ is GAS.

Theorem 12.

The coexistence point $S_{5} = (p_{4}, h_{4}, n_{4})$ is GAS if

(29)

{(\frac{(α_{1} + α_{2}) (1 + w) + α_{1} b p_{4} + α_{2} a h_{4}}{N_{1} N_{2}})}^{2} \leq \frac{1}{2} (\frac{r_{1}}{k_{1}} + \frac{α_{1} b h_{4}}{N_{1} N_{2}}) (\frac{r_{2}}{k_{2}} + \frac{α_{2} a p_{4}}{N_{1} N_{2}}),

(30)

{(\frac{α_{3} (1 + w)}{D_{1} D_{2}} - βn)}^{2} \leq \frac{1}{2} (\frac{r_{2}}{k_{2}} + \frac{α_{2} a p_{4}}{N_{1} N_{2}}) (β h_{4} + γ_{3}),

where,

N_{1} = 1 + w + (ah + bp)

,

N_{2} = 1 + w + (a h_{4} + b p_{4}), D_{1} = (1 + w + ch)

and

D_{2} = (1 + w + c h_{3})

.

Proof: Let $E_{5} = (p - p_{4} - p_{4} ln \frac{p}{p_{4}}) + (h - h_{4} - h_{4} ln \frac{h}{h_{4}}) + {(\frac{n - n_{4}}{2})}^{2}$ , where $E_{5} : R_{+}^{3} \to R$ , which satisfies $E_{5} (p_{4}, h_{4}, n_{4}) = 0$ and $E_{5} (p, h, n) > 0$ for all $(p, h, n) \in R_{+}^{3}$ with $(p, h, n) \neq (p_{4}, h_{4}, n_{4})$ , then

\frac{d E_{5}}{dt} = (\frac{p - p_{4}}{p}) \frac{dp}{dt} + (\frac{h - h_{4}}{h}) \frac{dh}{dt} + (n - n_{4}) \frac{dn}{dt} = (p - p_{4}) (r_{1} (1 - \frac{p}{k_{1}}) + \frac{α_{1} h}{1 + w + (ah + bp)} - γ_{1}) + (h - h_{4}) (r_{2} (1 - \frac{h}{k_{2}}) + \frac{α_{2} p}{1 + w + (ah + bp)} - γ_{2}) + (n - n_{4}) (\frac{α_{3} h}{1 + w + ch} - βnh - γ_{3} n) \frac{d E_{5}}{dt} = - [(\frac{r_{1}}{k_{1}} + \frac{α_{1} b h_{4}}{N_{1} N_{2}}) {(p - p_{4})}^{2} - (\frac{(α_{1} + α_{2}) (1 + w) + α_{1} b p_{4} + α_{2} a h_{4}}{N_{1} N_{2}}) (p - p_{4}) (h - h_{4}) + \frac{1}{2} (\frac{r_{2}}{k_{2}} + \frac{α_{2} a p_{4}}{N_{1} N_{2}}) {(h - h_{4})}^{2}] - [\frac{1}{2} (\frac{r_{2}}{k_{2}} + \frac{α_{2} a p_{4}}{N_{1} N_{2}}) {(h - h_{4})}^{2} - (\frac{α_{3} (1 + w)}{D_{1} D_{2}} - βn) (n - n_{4}) (h - h_{4}) + (β h_{4} + γ_{3}) {(n - n_{4})}^{2}]

Thus,

\frac{d E_{5}}{dt} \leq - {(\sqrt{(\frac{r_{1}}{k_{1}} + \frac{α_{1} b h_{4}}{N_{1} N_{2}})} (p - p_{4}) + \sqrt{\frac{1}{2} (\frac{r_{2}}{k_{2}} + \frac{α_{2} a p_{4}}{N_{1} N_{2}})} (h - h_{4}))}^{2} - {(\sqrt{\frac{1}{2} (\frac{r_{2}}{k_{2}} + \frac{α_{2} a p_{4}}{N_{1} N_{2}})} (h - h_{4}) + \sqrt{(β h_{4} + γ_{3})} (n - n_{4}))}^{2}

Hence, $d E_{5} / dt$ is a negative definite under conditions 29 and 30. Therefore, the coexistence point $S_{5} = (p_{4}, h_{4}, n_{4})$ is GAS.

Bifurcation

This section explores the probability of occurrence of transcritical (TB) and Hopf bifurcation (HB) around the non-hyperbolic equilibrium points. For more details, see Refs. 23-26.

Theorem 13.

For $r_{2}^{*} = γ_{2}$ , the $phn$ model faces TB at the extinction point $S_{1} = (0, 0, 0)$ .

Proof:

According to $J (S_{1})$ given by Eq. (7), the $phn$ system at $S_{1}$ has a zero eigenvalue $λ_{12} = 0$ , at $r_{2}^{*} = γ_{2}$ , and $J (S_{1})$ at $r_{2}^{*} = γ_{2}$ becomes

J^{*} (S_{1}) = [\begin{matrix} r_{1} - γ & 0 & 0 \\ 0 & 0 & 0 \\ 0 & \frac{α_{3}}{1 + w} & - γ_{3} \end{matrix}]

Now, suppose that $ϑ^{[1]} = {(ϑ_{1}^{[1]}, ϑ_{2}^{[1]}, ϑ_{3}^{[1]})}^{T}$ , and ${(T^{[1]})}^{T} = {(t_{1}^{[1]}, t_{2}^{[1]}, t_{3}^{[1]})}^{T}$ be eigenvectors to $λ_{22} = 0$ of $J^{*} (S_{1})$ , and $J^{* T} (S_{1})$ , respectively. The calculation gives $ϑ^{[1]} = (0, 1, \frac{α_{3}}{γ_{3} (1 + w)}$ ), and ${(T^{[1]})}^{T} = (0, 1, 0)$ by solving $(J^{*} (S_{1}) - λ_{12} I) ϑ^{[1]}$ , and $(J^{* T} (S_{1}) - λ_{12} I) T^{[1]}$ for $ϑ^{[1]}$ and $T^{[1]}$ .

Further,

F_{r_{2}} (S, r_{2}) = (0, h (1 - \frac{h}{k_{2}}), 0) \Rightarrow F_{r_{2}} (S_{1}, r_{2}^{*}) = (0, 0, 0)

{(T^{[1]})}^{T} F_{r_{2}} (S_{1}, r_{2}^{*}) = (0, 1, 0) (0, 0, 0) = 0

{(T^{[1]})}^{T} D F_{r_{2}} (S_{1}, r_{2}^{*}) ϑ^{[1]} = (0, 1, 0) [\begin{matrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{matrix}] = 1 \neq 0

{(T^{[1]})}^{T} D^{2} F_{r_{2}} (S_{1}, r_{2}^{*}) (ϑ^{[1]}, ϑ^{[1]}) = (0, 1, 0) [\begin{matrix} 0 \\ \frac{- 2 r_{2}^{*}}{k_{2}} \\ 0 \end{matrix}] = \frac{- 2 r_{2}^{*}}{k_{2}} \neq 0

Therefore, there is a TB around $S_{1}$ with the parameter $r_{2}^{*} = γ_{2}$ .

Theorem 14.

For $r_{1}^{*} = γ_{1}$ , the $phn$ model faces TB at the extinction point $S_{2} = (p_{1}, 0, 0)$ if

(30)

p_{1} = k_{1},

(31)

{(T^{[2]})}^{T} D^{2} F_{r_{1}} (S_{2}, r_{1}^{*}) (ϑ^{[2]}, ϑ^{[2]}) \neq 0

Proof:

According to $J (S_{2})$ given by Eq. (9), the $phn$ system at $S_{2}$ has a zero eigenvalue $λ_{22} = 0$ , at $r_{1}^{*} = γ_{1}$ , and $J (S_{2})$ at $r_{1}^{*} = γ_{1}$ becomes

J^{*} (S_{2}) = [\begin{matrix} 0 & 0 & 0 \\ 0 & r_{2} - γ_{2} & 0 \\ 0 & \frac{α_{3}}{1 + w} & - γ_{3} \end{matrix}],

Now, suppose that $ϑ^{[2]} = {(ϑ_{1}^{[2]}, ϑ_{2}^{[2]}, ϑ_{3}^{[2]})}^{T}$ and ${(T^{[1]})}^{T} = {(t_{1}^{[2]}, t_{2}^{[2]}, t_{3}^{[2]})}^{T}$ be eigenvectors with respect to $λ_{22} = 0$ of $J^{*} (S_{2})$ and ${J^{*}}^{T} (S_{2})$ respectively. Solving $(J^{*} (S_{2}) - λ_{22} I) ϑ^{[2]} = 0$ , and $(J^{* T} (S_{2}) - λ_{22} I) T^{[2]} = 0$ for $ϑ^{[2]}$ , and $T^{[2]}$ gives $ϑ^{[2]} = {(1, 1, \frac{α_{3}}{γ_{3 (1 + w)}})}^{T}$ and ${(T^{[2]})}^{T} = {(1, \frac{α_{3}}{(1 + w) (γ_{2} - r_{2})}, 1)}^{T}$ , where $(γ_{2} - r_{2}) \neq 0$ .

Further,

F_{r_{1}} (S, r_{1}) = (p (1 - \frac{p}{k_{1}}), 0, 0) ⟹ F_{r_{1}} (S_{2}, r_{1}^{*}) = (p_{1} (1 - \frac{p_{1}}{k_{1}}), 0, 0)

{(T^{[2]})}^{T} F_{r_{1}} (S_{2}, r_{1}^{*}) = (1, \frac{\propto_{3}}{(1 + w) (γ_{2} - r_{2})}, 1) (p_{1} - \frac{p_{1}^{2}}{k_{1}}, 0, 0) = p_{1} (1 - \frac{p_{1}}{k_{1}}) = 0 under condition 30 .

{(T^{[2]})}^{T} {DF}_{r_{1}} (S_{2}, r_{1}^{*}) (ϑ^{[2]}) = (1, \frac{α_{3}}{(1 + w) (γ_{2} - r_{2})}, 1) [\begin{matrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}] = 1 \neq 0

${(T^{[2]})}^{T} D^{2} F_{r_{1}} (S_{2}, r_{1}^{*}) (ϑ^{[2]}, ϑ^{[2]}) = (1, \frac{α_{3}}{(1 + w) (γ_{2} - r_{2})}, 1) [\begin{matrix} X_{11}^{[2]} + X_{12}^{[2]} \\ X_{21}^{[2]} + X_{22}^{[2]} \\ X_{32}^{[2]} + X_{33}^{[2]} \frac{\propto_{3}}{γ_{3} (1 + w)} \end{matrix}]$

$= (X_{11}^{[2]} + X_{12}^{[2]}) + \frac{α_{3}}{(1 + w) (γ_{2} - r_{2})} (X_{21}^{[2]} + X_{22}^{[2]}) + (X_{32}^{[2]} + X_{33}^{[2]} \frac{α_{3}}{γ_{3} (1 + w)}) \neq 0$

under condition 31, here, $X_{11}^{[2]} = \frac{- 2 r_{1}}{k_{1}} + [\frac{α_{1} [(1 + w)]}{{(1 + w + b p_{1})}^{2}}]$ , $X_{12}^{[2]} = [\frac{α_{1} (1 + w)}{{(1 + w + b p_{1})}^{2}}] + [\frac{- 2 a α_{1} p_{1}}{{(1 + w + b p_{1})}^{2}}]$ , $X_{21}^{[2]} = [\frac{α_{2} (1 + w)}{{(1 + w + b p_{1})}^{2}}]$ , $X_{22}^{[2]} = [\frac{α_{2} [(1 + w)]}{{(1 + w + b p_{1})}^{2}}] + [\frac{- 2 r_{2}}{k_{2}} - \frac{2 a α_{2} p_{1}}{{(1 + w + b p_{1})}^{2}}]$ , $X_{32}^{[2]} = \frac{- 2 c α_{3}}{{(1 + w)}^{2}} - \frac{β α_{3}}{γ_{3} (1 + w)}, and X_{33}^{[2]} = \frac{- β α_{3}}{γ_{3} (1 + w)} .$

Therefore, there is a TB around $S_{2}$ with the parameter $r_{1}^{*} = γ_{1}$ .

Theorem 15.

For ${r_{2}}^{* *} = \frac{γ_{2} k_{2}}{k_{2} - 2 h_{2}}$ , the $phn$ model faces TB at the honeybee point $S_{3} = (0, h_{2}, 0)$ if

(32)

h_{2} = k_{2},

Proof:

According to $J (S_{3})$ given by Eq. (11), the $phn$ system at $S_{3}$ has a zero eigenvalue $λ_{32} = 0$ , at ${r_{2}}^{* *}$ , and $J (S_{3})$ at ${r_{2}}^{* *}$ becomes

J^{*} (S_{3}) = [\begin{matrix} \frac{α_{1} h_{2}}{1 + w + a h_{2}} + (r_{1} - γ_{1}) & 0 & 0 \\ \frac{α h_{2}}{1 + w + a h_{2}} & 0 & 0 \\ 0 & \frac{(1 + w) α_{3}}{{(1 + w + c h_{2})}^{2}} & - (β h_{2} + γ_{3}) \end{matrix}]

.

Now, Suppose that $ϑ^{[3]} = (ϑ_{1}^{[3]}, ϑ_{2}^{[3]}, ϑ_{3}^{[3]})$ and ${(T^{[3]})}^{T} = {(t_{1}^{[3]}, t_{2}^{[3]}, t_{3}^{[3]})}^{T}$ be eigenvectors to $λ_{32} = 0$ of $J^{*} (S_{3})$ and ${J^{*}}^{T} (S_{3})$ respectively. Solving $(J^{*} (S_{3}) - λ_{32} I) ϑ^{[3]} = 0$ , and $(J^{* T} (S_{3}) - λ_{32} I) T^{[3]} = 0$ for $ϑ^{[3]}$ and $T^{[3]}$ gives $V^{[3]} = (0, 1, \frac{α_{3} (1 + w)}{(β h_{2} + γ_{3}) {(1 + w + c h_{2})}^{2}})$ and ${(T^{[3]})}^{T} = (\frac{- α_{2} h_{2}}{(α_{1} h_{2} + (r_{1} - γ_{1}) (1 + w + a h_{2}))}, 1, 0)$ , where $[α_{1} h_{2} + [(r_{1} - γ_{1}) (1 + w + a h_{2})] \neq 0$ .

Further

\begin{matrix} F_{r_{2}} (S_{3}, r_{2}) = (0, h - \frac{h}{k_{2}}, 0) \Rightarrow F_{r_{2}} (S_{3}, r_{2}^{* *}) = (0, h_{2} (1 - \frac{h_{2}}{k_{2}}), 0) \\ {(T^{[3]})}^{T} F_{r_{2}} (S_{3}, r_{2}^{* *}) = (\frac{- α_{2} h_{2}}{(α_{1} h_{2} + (r_{1} - γ_{1}) (1 + w + a h_{2}))}, 1, 0) (0, h_{2} (1 - \frac{h_{2}}{k_{2}}), 0) \\ \begin{matrix} {(T^{[3]})}^{T} F_{r_{2}} (S_{3}, r_{2}^{* *}) = h_{2} (1 - \frac{h_{2}}{k_{2}}) = 0 under condition 32 . \\ \begin{matrix} {(T^{[3]})}^{T} {DF}_{r_{2}} (S_{3}, r_{2}^{* *}) (ϑ^{[3]}) = (\frac{- α_{2} h_{2}}{(α_{1} h_{2} + (r_{1} - γ_{1}) (1 + w + a h_{2}))}, 1, 0) [\begin{matrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{matrix}] = 1 \neq 0 \\ {(T^{[3]})}^{T} D^{2} F_{r_{2}} (S_{3}, r_{2}^{* *}) (ϑ^{[3]}, ϑ^{[3]}) = (\frac{- α_{2} h_{2}}{(α_{1} h_{2} + (r_{1} - γ_{1}) (1 + w + a h_{2}))}, 1, 0) [\begin{matrix} 0 \\ \frac{- 2 r_{2}^{* *}}{k_{2}} \\ X_{32}^{[3]} - β \frac{α_{3} (1 + w)}{(β h_{2} + γ_{2}) {(1 + w + c h_{2})}^{2}} \end{matrix}] = \frac{- 2 r_{2}^{* *}}{k_{2}} \neq 0, \end{matrix} \end{matrix} \end{matrix}

where

X_{32}^{[3]} = \frac{- 2 c α_{3} (1 + w)}{{(1 + w + c h_{2})}^{3}} - \frac{β α_{3} (1 + w)}{(β h_{2} + γ_{2}) {(1 + w + c h_{2})}^{2}}

Therefore, there is a TB around $S_{3}$ with the parameter $r_{2}^{* *}$ .

Theorem 16.

For ${r_{2}}^{* * *} = \frac{γ_{2} k_{2}}{k_{2} - 2 h_{3}}$ , the $phn$ model faces TB at the honeybee point $S_{4} = (0, h_{3}, n_{3})$ if

(33)

h_{3} = k_{2},

Proof:

According to $J (S_{4})$ given by Eq. (15), the $phn$ system at $S_{4}$ has a zero eigenvalue $λ_{42} = 0$ , at ${r_{2}}^{* * *}$ , and $J (S_{4})$ at ${r_{2}}^{* * *}$ becomes

J^{*} (S_{4}) = [\begin{matrix} a_{11} & 0 & 0 \\ a_{21} & 0 & 0 \\ 0 & a_{32} & - β h_{3} - γ_{3} \end{matrix}],

where,

a_{11} = \frac{α_{1} h_{3}}{(1 + w + a h_{3})} + (r_{1} - γ_{1})

,

a_{21} = \frac{α_{2} h_{3}}{1 + w + a h_{3}}

and

a_{32} = \frac{α_{3} (1 + w)}{(1 + w + c h_{3})} - β n_{3}

.

Now, Suppose that $ϑ^{[4]} = (ϑ_{1}^{[4]}, ϑ_{2}^{[4]}, ϑ_{3}^{[4]})$ , and ${(T^{[4]})}^{T} = (t_{1}^{[4]}, t_{2}^{[4]}, t_{3}^{[4]})$ be an eigenvector of to $λ_{42} = 0$ of $J^{*} (S_{4})$ and $J^{* T} (S_{4})$ , respectively, which gives $ϑ^{[4]} = (0, 1, \frac{a_{32}}{β h_{3} + γ_{3}})$ and ${(T^{[4]})}^{T} = (\frac{- a_{21}}{a_{11}}, 1, 0)$ .

Further,

\begin{matrix} F_{r_{2}} (S, r_{2}) = (0, h (1 - \frac{h}{k_{2}}), 0) \Rightarrow F_{r_{2}} (S_{4}, r_{2}^{* * *}) = (0, h_{2} - \frac{h_{2}}{k_{2}}, 0) \\ {(T^{[4]})}^{T} F_{r_{2}} (S_{3}, r_{2}^{*}) = (\frac{- a_{21}}{a_{11}}, 1, 0) (0, h_{3}^{2} - \frac{h_{3}^{2}}{k_{2}}, 0) = h_{3} (1 - \frac{h_{3}}{k_{2}}) = 0 under condition 33 . \\ \begin{matrix} {(T^{[4]})}^{T} {DF}_{r_{2}} (S_{4}, r_{2}^{*}) (ϑ^{[4]}) = (\frac{- a_{21}}{a_{11}}, 1, 0) [\begin{matrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{matrix}] = 1 \neq 0 \\ {(T^{[4]})}^{T} D^{2} {F_{r}}_{2} (S_{4}, r_{2}^{*}) (ϑ^{[4]}, ϑ^{[4]}) = (\frac{- a_{21}}{a_{11}}, 1, 0) [\begin{matrix} 0 \\ \frac{- 2 r_{2}^{*}}{k_{2}} \\ \frac{- 2 c α_{3} (1 + w)}{{(1 + w + c h_{3})}^{3}} + β {(\frac{a_{32}}{β h_{3} + γ_{3}})}^{2} \end{matrix}] = \frac{- 2 r_{2}^{*}}{k_{2}} \neq 0, \end{matrix} \end{matrix}

Therefore, there is a TB around $S_{4}$ with the parameter $r_{2}^{* * *}$ .

Theorem 17.

The $phn$ undergoes a Hopf bifurcation at the coexistence point $S_{5}$ to the bifurcation parameter $w^{*}$ if

(34)

{[Det] |}_{w = w^{*}} > 0,

(35)

b α_{1} p_{4} h_{4} + a α_{2} p_{4} h_{4} \neq α_{1} h_{4} (1 + w^{*} + a h_{4}) + α_{2} p_{4} (1 + w^{*} + b p_{4}),

where the formula of

w^{*}

is given in the proof.

Proof.

The Jacobian matrix at $S_{5}$ with $w^{*}$ is given by

J (S_{5}, w^{*}) = [\begin{matrix} a_{11} & \frac{α_{1} p_{4} (1 + w^{*} + b p_{4})}{{(1 + w^{*} + (a h_{4} + b p_{4}))}^{2}} & 0 \\ \frac{α_{2} h_{4} (1 + w^{*} + a h_{4})}{{(1 + w^{*} + (a h_{4} + b p_{4}))}^{2}} & a_{22} & 0 \\ 0 & \frac{α_{3} (1 + w^{*})}{{(1 + w^{*} + c h_{4})}^{2}} - β n_{4} & - β h_{4} - γ_{3} \end{matrix}],

where $a_{11} = (r_{1} - γ_{1}) - \frac{2 r_{1} p_{4}}{k_{1}} + \frac{α_{1} h_{4} (1 + w^{*} + a h_{4})}{{(1 + w^{*} + (a h_{4} + b p_{4}))}^{2}}$ , and $a_{22} = (r_{2} - γ_{2}) - \frac{2 r_{2} h_{4}}{k_{2}} + \frac{α_{2} p_{4} (1 + w^{*} + b p_{4})}{{(1 + w^{*} + (a h_{4} + b p_{4}))}^{2}}$ .

The Hopf bifurcation occurred if the following conditions are satisfied.

1. ${[Tr] |}_{w = w^{*}} = 0$ ,
2. ${[Det] |}_{w = w^{*}} > 0$ ,
3. $\frac{\partial}{∂w} {[Re (λ_{1, 2})]}_{w = w^{*}} \neq 0$ (Transversality condition),

where $Tr$ and $Det$ are defined in the characteristic equation given by (19). Now, we set $Tr = 0$ to find $w^{*}$ , which gives

Tr = 0 ⟹ (r_{1} - γ_{1}) - \frac{2 r_{1} p_{4}}{k_{1}} + \frac{α_{1} h_{4} (1 + w + a h_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}} - \frac{2 r_{2} h_{4}}{k_{2}} + \frac{α_{2} p_{4} (1 + w + b p_{4})}{{(1 + w + (a h_{4} + b p_{4}))}^{2}} + (r_{2} - γ_{2}) = 0 .

Let $b_{1} = (r_{1} - γ_{1}) + (r_{2} - γ_{2}) - \frac{2 r_{1} p_{4}}{k_{1}} - \frac{2 r_{2} h_{4}}{k_{2}}$ , $b_{2} = (a h_{4} + b p_{4})$ , $b_{3} = α_{1} a h_{4}^{2} + α_{2} b p_{4}^{2}$ , $b_{4} = (α_{1} h_{4} + α_{2} p_{4})$ , and $v = (1 + w)$ . Then

Tr = 0 ⟹ b_{1} + \frac{b_{4} v + b_{3}}{{(v + b_{2})}^{2}} = 0

Solving the above equation for $v$ , gives

{(v + b_{2})}^{2} b_{1} + b_{4} v + b_{3} = 0

Then,

\begin{matrix} {(v + b_{2})}^{2} b_{1} + b_{4} v + b_{3} = 0 \\ b_{1} v^{2} + (2 b_{1} b_{2} + b_{4}) v + (b_{1} b_{2}^{2} + b_{3}) = 0 . \end{matrix}

The discriminant of the above equation is

Δ = {(2 b_{1} b_{2} + b_{4})}^{2} - 4 b_{1} (b_{1} b_{2}^{2} + b_{3}) = b_{4}^{2} + 4 b_{1} (b_{2} b_{4} - b_{3}) .

Thus,

v = \frac{- (2 b_{1} b_{2} + b_{4}) \pm \sqrt{Δ}}{2 b_{1}}, w^{*} = v - 1, b_{1} \neq 0

That means condition (1) is satisfied at $w^{*}$ , and the characteristic equation given by (19) can be rewritten as

(- β h_{4} - γ_{3} - λ) [λ^{2} + Det] = 0

Solving the above equation yields $λ_{1} = - β h_{4} - γ_{3}$ , $λ_{2, 3} = \pm i \sqrt{Det}$ . Clearly $λ_{2}$ and $λ_{3}$ are complex conjugates if condition 34 is satisfied. In addition, the general roots of Eq. (19) in the neighborhood of $w^{*}$ as $λ_{2, 3} = \frac{Tr \pm i \sqrt{{Tr}^{2} - 4 Det}}{2}$ , then

\frac{\partial}{∂w} {[Re (λ_{2, 3})]}_{w = w^{*}} = \frac{b α_{1} p_{4} h_{4} + a α_{2} p_{4} h_{4} - α_{1} h_{4} (1 + w^{*} + a h_{4}) - α_{2} p_{4} (1 + w^{*} + b p_{4})}{{(1 + w^{*} + (a h_{4} + b p_{4}))}^{3}} \neq 0 under condition 35 .

So, $phn$ system undergoes a Hopf bifurcation at $S_{5}$ with the bifurcation parameter $w^{*}$ .

Numerical simulations and discussion

A numerical confirmation is carried out to complete the analytical results for $phn$ system using MATLAB. The simulations are conducted by using the following set of parameters.

(36)

r_{1} = 0.56, k_{1} = 10, α_{1} = 0.04, w = 3, a = 0.02, b = 0.03 γ_{1} = 0.04, r_{2} = 0.18, k_{2} = 6, α_{2} = 0.08 = 0.4, γ_{2} = 0.03, α_{3} = 0.01, β = 0.00157, γ_{3} = 0.0018, c = 0.01 .

For the parameters listed above in Eq. (36), the nullclines of $phn$ system are indicated in Figure 2. The figure depicts the coexistence point $S_{5} = (p_{4}, h_{4}, n_{4}) = (0.21, 0.11, 0.23)$ , while Figure 3 illustrates the global behavior of $S_{5}$ . That means the parameters listed in Eq. (36) corroborate the findings of Theorem 12, which shows that all other points act as saddle points, except for $S_{5}$ . These findings also confirm the uniform persistence for $phn$ system, which confirms the output of Theorem 3.

Figure 2. The nullclines of $phn$ system with the parameters listed in Eq. (36).

Figure 3. The global stability of $S_{5} = (0.21, 0.11, 0.23)$ .

An important parameter to investigate is the effect of the wind level $(w)$ on the interaction of blooming plants $p (t)$ , honey honeybees $h (t),$ and honey production $n (t)$ . We address two aspects of $w$ : first, the extent to which it influences the species densities in the inner equilibrium, and second, how it can alter $phn$ system’s stability. For a low level of wind flow, i.e., $0 < w < 0.017,$ the solution of the $phn$ system approaches a chaotic attractor, see Figure 4(a). While in the interval $0.017 < w \leq 0.74$ , the solution of the $phn$ system converges to a limit cycle, see Figure 4(b) and Figure 5. Consequently, for the region $0.74 < w < 3.4$ , the solution stabilized at the coexistence point $S_{5}$ , see Figure 2, which was plotted when $w = 3$ . Finally, for a high level of wind speed, i.e., for $w > 3.4$ , $phn$ system loses two of its components, and the solution in this case settles down to the blooming plants’ point $S_{2} = (0.16, 0, 0)$ , see Figure 4(c). This indicates that when the wind speed is light, the wind helps pollinate or carries the seeds of flowering plants from one place to another, aiding their reproduction. Therefore, the presence of flowering plants contributes to the system’s persistence. In contrast, for a high wind flow, we see that the populations of honeybees and the production of honey face extinction.

Figure 4. The effect of various wind speeds $w$ .

Figure 5. Hopf bifurcation at $w = 0.74$ .

Now, we have the intrinsic growth rate of the honeybee population ( $r_{2})$ , which is an essential quantity to discuss since it can potentially influence the population’s densities of blooming plants, honeybees, and honey production. It is clear from Figure 6 that the solution of $phn$ system settles down to the coexistence point $S_{5}$ for $r_{2} > 0.03$ , while it stabilized at the extinction point $S_{1} = (0, 0, 0)$ for $r_{2} \leq 0.03$ . This result confirms the occurrence of a transcritical bifurcation at $r_{2} = r_{2}^{TB} = 0.03$ , which confirms the output of Theorem 13, see Figure 7. Further, Figure 8 indicates the global behavior of the extinction point $S_{1}$ . This result confirms the global stability condition of $S_{1}$ which has honeybee stated in Theorem 8. Further, this result shows that $r_{2}$ is a critical parameter that impacts the continuity of the whole system’s coexistence.

Figure 6. The effect of varying $r_{2}$ .

Figure 7. Transcritical bifurcation at $r_{2} = 0.03$ .

Figure 8. The global stability of $S_{1}$ .

The effect of varying the intrinsic growth rate $r_{1}$ of the population of blooming plants is investigated in Figure 9. The figure indicates that the solution of $phn$ system converges to the blooming plants point $S_{2} = (0.28, 0, 0)$ for $r_{1} \leq 0.04$ , which means the system faces an occurrence of transcritical bifurcations at $r_{1} = r_{1}^{TB}$ = 0.04, which confirms the output of Theorem 14. So, $phn$ system loses two of its components for $r_{1} \leq 0.04$ . While for $r_{1} > 0.04$ , $phn$ system converges to the coexistence point $S_{5}$ . Further, the global stability of $S_{2}$ is illustrated in Figure 10. We can conclude that the intrinsic growth rate $r_{1}$ of the blooming plants population is a critical parameter affecting the honeybees’ persistence and honey production.

Figure 9. The effect of varying $r_{1}$ .

Figure 10. The global stability of $S_{2}$ .

The mutualistic rates between the honeybee population and the blossoming plants, $α_{1}$ and $α_{2}$ , are examined in Figures 11 and 12. The density of blossoming plants, the honeybee population, and honey production are all improved as a result of the increase in mutualistic rates.

Figure 11. The effect of varying $α_{1}$ .

Figure 12. The effect of varying $α_{2}$ .

A rise in the depletion rate of the population of blossoming plants population $γ_{1}$ (when $γ_{1} = 0.57)$ leads to losses of the honeybee population and honey in $phn$ system, and the solution stabilized at the honeybee point $S_{3} = (0, 0.52, 0)$ from various initial conditions. This result confirms the global stability theorem of $S_{3}$ which was stated in Theorem 10. As a result, the continued presence of a blossoming plant population significantly impacts the persistence of honey production since the flowering plants are the primary source of sustenance for pollinators. See Figure 13.

Figure 13. The global stability of $S_{3}$ when $γ_{1} = 0.57$ .

To establish the effect of $β$ on the dynamics of $phn$ , Figure 14 has been drawn with three values of the rate at which honeybees consume honey to survive, i.e., $β$ . The figure shows the solutions settling down to the coexistence point $S_{5}$ . The increase in $β$ substantially results in a decline in honey output.

Figure 14. The effect of varying $β$ .

In order to investigate the sensitivity of the coexistence point $S_{5} = (p_{4}, h_{4}, n_{4})$ of $phn$ system, we implement partial rank correlation coefficients (PRCC). The parameters $r_{1}, k_{1}, α_{1}, w, a, b, γ_{1}, r_{2}, k_{2}, α_{2}, γ_{2}, β, γ_{3}, c$ and $α_{3}$ serve as input parameters, whereas $p_{4}, h_{4},$ and $n_{4}$ the output variables. We subsequently generate Figure 15 by utilizing the parameter set in Eq. (36). Blossoming plants and the honeybee population demonstrate heightened sensitivity to honeybees’ carrying capacity $k_{2}$ and the corresponding value of honeybee nutrients from blooming plants $α_{2}$ . While the honey production determines heightened sensitivity to $α_{3}$ . On the other hand, the wind flow, i.e., $w,$ significantly influences $p_{4}$ , $h_{4}$ and $n_{4}$ . The wind flow significantly reduces the blooming plants, honeybee population, and honey production. It can be inferred that the wind flow is a critical parameter that influences the coexistence of $p_{4}, h_{4},$ and $n_{4}$ , see Figure 15.

Figure 15. The sensitivity of the data given in Eq. 36 relative to the $S_{5} = (p_{4}, h_{4}, n_{4})$ .

Conclusion

Wind flow profoundly affects blossoming plants, honeybee populations, and honey production dynamics, and impacts ecosystem stability. An ODE mathematical model has been studied to understand these dynamics. The solution of $phn$ system has been established to possess the fundamental attributes, such as positivity and persistence, boundedness, local and global stability, and bifurcation. Numerical results indicated that the honeybee species may be extinct due to increased wind velocities within specific parameter ranges. Furthermore, the coexistence equilibrium becomes unstable as a result of a Hopf bifurcation when a low wind flow induces periodic oscillations. Further, the simulations indicated that the threshold values for the transcritical bifurcation have been precisely determined at a decreased honeybee and blooming plant growth rate. However, the system may reach a point where both the blooming plant and the honeybee populations are no longer viable due to an elevated mortality rate of flowering plants. On the other hand, the increase in mutualistic rates between the honeybee population and the blooming plants has a regenerative effect, supporting the sustainability of the honeybee–honey production system.

Data availability

This study relies on numerical data generated from the proposed mathematical model. This data includes initial coefficients, initial conditions, and numerical simulation outputs (tables and figures). All of this data is not derived from field measurements but was generated programmatically for the purposes of theoretical analysis and numerical simulation of the system under study. This aligns with the general trend of sharing research data to enhance replication and reuse within the scientific community. All data underlying the results are available as part of the article, and no additional source data are required

References

1. Dean AM: A simple model of mutualism. Am. Nat. 1983; 121(3): 409–417. Publisher Full Text
2. Ollerton J, Winfree R, Tarrant S: How many flowering plants are pollinated by animals? Oikos. 2011; 120(3): 321–326. Publisher Full Text
3. Lee Y-D, Yokoi T, Nakazawa T: A pollinator crisis can decrease plant abundance despite pollinators being herbivores at the larval stage. Sci. Rep. 2024; 14(1): 18523. PubMed Abstract | Publisher Full Text | Free Full Text
4. Fishman MA, Hadany L: Plant–pollinator population dynamics. Theor. Popul. Biol. 2010; 78(4): 270–277. Publisher Full Text
5. Jawad S, Thirthar AA, Nisar KS: The impact of climate change on flowering plants-bees-Vespa orientalis model. Results Control Optim. 2025; 20: 100583. Publisher Full Text
6. Biswas A, Medda R, Pal S: Dynamics of predatory effect on saturated plant–pollinator mutualistic relationship. Chaos An Interdiscip. J. Nonlinear Sci. 2025; 35(2). PubMed Abstract | Publisher Full Text
7. Hakeem E, Jawad S, Ali AH, et al.: How mathematical models might predict desertification from global warming and dust pollutants. MethodsX. 2025; 14: 103259. PubMed Abstract | Publisher Full Text | Free Full Text
8. Lever JJ, van Nes EH , Scheffer M, et al.: The sudden collapse of pollinator communities. Ecol. Lett. 2014; 17(3): 350–359. PubMed Abstract | Publisher Full Text
9. Al Nuaimi M, Jawad S: Modelling and stability analysis of the competitional ecological model with harvesting. Commun. Math. Biol. Neurosci. 2022; 2022: Article-ID.
10. Shalan RN, Shireen R, Lafta AH: Discrete an SIS model with immigrants and treatment. J. Interdiscip. Math. 2020; 24: 1201–1206. Publisher Full Text
11. Ali A, Jawad S: Stability analysis of the depletion of dissolved oxygen for the Phytoplankton-Zooplankton model in an aquatic environment. Iraqi J. Sci. 2024; 2736–2748. Publisher Full Text
12. Nezar F, Jawad S, Winter M, et al.: Stability analysis of excessive carbon dioxide gas emission model through following reforestation policy in low-density forest biomass. Baghdad Sci. J. 2025; 22(4): 1335–1353.
13. Cresswell JE: A demographic approach to evaluating the impact of stressors on bumble bee colonies. Ecol. Entomol. 2017; 42(2): 221–229. Publisher Full Text
14. Thirthar AA, Panja P, Abdeljawad T: Impact of alarm signals and mutualistic interactions in a food chain model of oxpeckers, zebras, and lions. Partial Differ. Equations Appl. Math. 2025; 14: 101189. Publisher Full Text
15. Jawad S, Roy S, Thirthar AA, et al.: Deterministic and stochastic risks assessment of excessive CO2 emission on forest biomass under weak allee effect. J. Appl. Math. Comput. 2025; 71: 9129–9156. Publisher Full Text
16. Thirthar AA, Alaoui AL, Roy S, et al.: Fractional and stochastic dynamics of predator–prey systems: The role of fear and global warming. Eur. Phys. J. B. 2025; 98(7): 1–21. Publisher Full Text
17. Ahmed M, Jawad S, Das D, et al.: Impact of dust storms on plant biomass: Model structure and dynamic study. Alex. Eng. J. 2025; 126: 605–622. Publisher Full Text
18. Balfour NJ, Ratnieks FLW: Wind Alters Plant-Pollinator Community Structure, Bee Foraging Rate & Movements Between Plants. Behav. Ecol. 2025; 36: araf067. PubMed Abstract | Publisher Full Text | Free Full Text
19. Friedman J, Barrett SCH: Wind of change: new insights on the ecology and evolution of pollination and mating in wind-pollinated plants. Ann. Bot. 2009; 103(9): 1515–1527. PubMed Abstract | Publisher Full Text | Free Full Text
20. Chen F: On a nonlinear nonautonomous predator–prey model with diffusion and distributed delay. J. Comput. Appl. Math. 2005; 180(1): 33–49. Publisher Full Text
21. Place CM: Dynamical Systems: Differential Equations, Maps, and Chaotic Behaviour. Routledge; 2017.
22. Aamer Z, Jawad S, Batiha B, et al.: Evaluation of the Dynamics of Psychological Panic Factor, Glucose Risk and Estrogen Effects on Breast Cancer Model. Computation. 2024; 12(8): 160. Publisher Full Text
23. Ahmed M, Jawad S: Bifurcation analysis of the role of good and bad bacteria in the decomposing toxins in the intestine with the impact of antibiotic and probiotics supplement. AIP Conference Proceedings. AIP Publishing; 2024.
24. Thirthara AA, Jawadb S, Shahc K, et al.: How does media coverage affect a COVID-19 pandemic model with direct and indirect transmission?2024.
25. Javaid Y, Jawad S, Ahmed R, et al.: Dynamic complexity of a discretized predator-prey system with Allee effect and herd behaviour. Appl. Math. Sci. Eng. 2024; 32(1): 2420953. Publisher Full Text
26. Roy S, AL-Jaf DS, Thirthar AA, et al.: The impact of human shields in autonomous and non-autonomous prey-predator models with modified Cosner functional response. Math. Comput. Simul. 2025; 242: 54–73. Publisher Full Text

Comments on this article Comments (0)

Version 1

VERSION 1 PUBLISHED 26 Dec 2025

Author details Author details

¹ Mathematics, University of Baghdad, Baghdad, Baghdad Governorate, Iraq

Shireen Jawad
Roles: Methodology, Supervision, Writing – Review & Editing

Zainab Hayder Abid AL-Aali
Roles: Formal Analysis, Investigation, Writing – Original Draft Preparation

Competing interests

No competing interests were disclosed.

Grant information

The author(s) declared that no grants were involved in supporting this work.

Article Versions (1)

version 1

Published: 26 Dec 2025, 14:1459

https://doi.org/10.12688/f1000research.172134.1

Copyright

© 2025 Jawad S and AL-Aali ZHA. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Download

Export To

metrics

	Views	Downloads
F1000Research	-	-
PubMed Central Data from PMC are received and updated monthly.	-	-

Citations

0

SEE MORE DETAILS

CITE

how to cite this article

Jawad S and AL-Aali ZHA. Impact of high wind speed on blooming plants-honeybees-honey production model  [version 1; peer review: awaiting peer review]. F1000Research 2025, 14:1459 (https://doi.org/10.12688/f1000research.172134.1)

NOTE: If applicable, it is important to ensure the information in square brackets after the title is included in all citations of this article.

track

receive updates on this article

Track an article to receive email alerts on any updates to this article.

Open Peer Review

Comments on this article Comments (0)

Version 1

VERSION 1 PUBLISHED 26 Dec 2025

Open Peer Review

Reviewer Status

AWAITING PEER REVIEW

Comments on this article

All Comments(0)

Add a comment

Sign up for content alerts

Browse by related subjects

[1] 1. Dean AM: A simple model of mutualism. Am. Nat. 1983; 121(3): 409–417. Publisher Full Text

[2] 2. Ollerton J, Winfree R, Tarrant S: How many flowering plants are pollinated by animals? Oikos. 2011; 120(3): 321–326. Publisher Full Text

[3] 3. Lee Y-D, Yokoi T, Nakazawa T: A pollinator crisis can decrease plant abundance despite pollinators being herbivores at the larval stage. Sci. Rep. 2024; 14(1): 18523. PubMed Abstract | Publisher Full Text | Free Full Text

[4] 4. Fishman MA, Hadany L: Plant–pollinator population dynamics. Theor. Popul. Biol. 2010; 78(4): 270–277. Publisher Full Text

[5] 5. Jawad S, Thirthar AA, Nisar KS: The impact of climate change on flowering plants-bees-Vespa orientalis model. Results Control Optim. 2025; 20: 100583. Publisher Full Text

[6] 6. Biswas A, Medda R, Pal S: Dynamics of predatory effect on saturated plant–pollinator mutualistic relationship. Chaos An Interdiscip. J. Nonlinear Sci. 2025; 35(2). PubMed Abstract | Publisher Full Text

[7] 7. Hakeem E, Jawad S, Ali AH, et al.: How mathematical models might predict desertification from global warming and dust pollutants. MethodsX. 2025; 14: 103259. PubMed Abstract | Publisher Full Text | Free Full Text

[8] 8. Lever JJ, van Nes EH , Scheffer M, et al.: The sudden collapse of pollinator communities. Ecol. Lett. 2014; 17(3): 350–359. PubMed Abstract | Publisher Full Text

[9] 9. Al Nuaimi M, Jawad S: Modelling and stability analysis of the competitional ecological model with harvesting. Commun. Math. Biol. Neurosci. 2022; 2022: Article-ID.

[10] 10. Shalan RN, Shireen R, Lafta AH: Discrete an SIS model with immigrants and treatment. J. Interdiscip. Math. 2020; 24: 1201–1206. Publisher Full Text

[11] 11. Ali A, Jawad S: Stability analysis of the depletion of dissolved oxygen for the Phytoplankton-Zooplankton model in an aquatic environment. Iraqi J. Sci. 2024; 2736–2748. Publisher Full Text

[12] 12. Nezar F, Jawad S, Winter M, et al.: Stability analysis of excessive carbon dioxide gas emission model through following reforestation policy in low-density forest biomass. Baghdad Sci. J. 2025; 22(4): 1335–1353.

[13] 13. Cresswell JE: A demographic approach to evaluating the impact of stressors on bumble bee colonies. Ecol. Entomol. 2017; 42(2): 221–229. Publisher Full Text

[14] 14. Thirthar AA, Panja P, Abdeljawad T: Impact of alarm signals and mutualistic interactions in a food chain model of oxpeckers, zebras, and lions. Partial Differ. Equations Appl. Math. 2025; 14: 101189. Publisher Full Text

[15] 15. Jawad S, Roy S, Thirthar AA, et al.: Deterministic and stochastic risks assessment of excessive CO2 emission on forest biomass under weak allee effect. J. Appl. Math. Comput. 2025; 71: 9129–9156. Publisher Full Text

[16] 16. Thirthar AA, Alaoui AL, Roy S, et al.: Fractional and stochastic dynamics of predator–prey systems: The role of fear and global warming. Eur. Phys. J. B. 2025; 98(7): 1–21. Publisher Full Text

[17] 17. Ahmed M, Jawad S, Das D, et al.: Impact of dust storms on plant biomass: Model structure and dynamic study. Alex. Eng. J. 2025; 126: 605–622. Publisher Full Text

[18] 18. Balfour NJ, Ratnieks FLW: Wind Alters Plant-Pollinator Community Structure, Bee Foraging Rate & Movements Between Plants. Behav. Ecol. 2025; 36: araf067. PubMed Abstract | Publisher Full Text | Free Full Text

[19] 19. Friedman J, Barrett SCH: Wind of change: new insights on the ecology and evolution of pollination and mating in wind-pollinated plants. Ann. Bot. 2009; 103(9): 1515–1527. PubMed Abstract | Publisher Full Text | Free Full Text

[20] 20. Chen F: On a nonlinear nonautonomous predator–prey model with diffusion and distributed delay. J. Comput. Appl. Math. 2005; 180(1): 33–49. Publisher Full Text

[21] 21. Place CM: Dynamical Systems: Differential Equations, Maps, and Chaotic Behaviour. Routledge; 2017.

[22] 22. Aamer Z, Jawad S, Batiha B, et al.: Evaluation of the Dynamics of Psychological Panic Factor, Glucose Risk and Estrogen Effects on Breast Cancer Model. Computation. 2024; 12(8): 160. Publisher Full Text

[23] 23. Ahmed M, Jawad S: Bifurcation analysis of the role of good and bad bacteria in the decomposing toxins in the intestine with the impact of antibiotic and probiotics supplement. AIP Conference Proceedings. AIP Publishing; 2024.

[24] 24. Thirthara AA, Jawadb S, Shahc K, et al.: How does media coverage affect a COVID-19 pandemic model with direct and indirect transmission?2024.

[25] 25. Javaid Y, Jawad S, Ahmed R, et al.: Dynamic complexity of a discretized predator-prey system with Allee effect and herd behaviour. Appl. Math. Sci. Eng. 2024; 32(1): 2420953. Publisher Full Text

[26] 26. Roy S, AL-Jaf DS, Thirthar AA, et al.: The impact of human shields in autonomous and non-autonomous prey-predator models with modified Cosner functional response. Math. Comput. Simul. 2025; 242: 54–73. Publisher Full Text

Impact of high wind speed on blooming plants-honeybees-honey production model

Abstract

Background

Methods

Results

Conclusions

Keywords

Introduction

Methods

Assumptions of the model

(1)

Figure 1. Flowchart of the phn model.

Model analysis

Lemma 1.

Lemma 2.

Positivity

Theorem 1.

Proof.

Boundedness

Theorem 2.

Proof.

Persistence

Theorem 3.

(2)

(3)

Proof:

Remark 1:

Equilibrium analysis

(4)

Local stability of equilibrium points

Theorem 4.

(5)

(6)

Proof:

(7)

Remark 2:

Theorem 5.

(8)

Proof:

(9)

Theorem 5.

(10)

(11)

Proof:

(12)

Theorem 6.

(13)

(14)

Proof:

(15)

Theorem 7.

(16)

(17)

Proof:

(18)

(19)

Global stability

Theorem 8.

(20)

(21)

Proof:

Theorem 9.

(22)

(23)

Proof:

Theorem 10.

(24)

(25)

(26)

Proof:

Theorem 11.

(27)

(28)

Proof:

Theorem 12.

(29)

(30)

Bifurcation

Theorem 13.

Proof:

Impact of high wind speed on blooming plants-honeybees-honey production model 

Figure 1. Flowchart of the $phn$ model.

Figure 2. The nullclines of $phn$ system with the parameters listed in Eq. (36).

Figure 3. The global stability of $S_{5} = (0.21, 0.11, 0.23)$ .

Figure 4. The effect of various wind speeds $w$ .

Figure 5. Hopf bifurcation at $w = 0.74$ .

Figure 6. The effect of varying $r_{2}$ .

Figure 7. Transcritical bifurcation at $r_{2} = 0.03$ .

Figure 8. The global stability of $S_{1}$ .

Figure 9. The effect of varying $r_{1}$ .

Figure 10. The global stability of $S_{2}$ .

Figure 11. The effect of varying $α_{1}$ .

Figure 12. The effect of varying $α_{2}$ .

Figure 13. The global stability of $S_{3}$ when $γ_{1} = 0.57$ .

Figure 14. The effect of varying $β$ .

Figure 15. The sensitivity of the data given in Eq. 36 relative to the $S_{5} = (p_{4}, h_{4}, n_{4})$ .