F1000Research F1000Research 2046-1402 F1000 Research Limited London, UK 10.12688/f1000research.172134.2 Research Article Articles Impact of high wind speed on blooming plants-honeybees-honey production model 

[version 2; peer review: 2 approved]

 AL-Aali Zainab Hayder Abid Formal Analysis Investigation Writing – Original Draft Preparation 1 Jawad Shireen Methodology Supervision Writing – Review & Editing https://orcid.org/0000-0002-3090-8357 a 1 Mathematics, University of Baghdad, Baghdad, Baghdad Governorate, Iraq shireen.jawad@sc.uobaghdad.edu.iq

No competing interests were disclosed.

 29 4 2026 2025 14 1459 17 4 2026 Copyright: © 2026 AL-Aali ZHA and Jawad S 2026 This is an open access article distributed under the terms of the Creative Commons Attribution Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 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.

 Wind speed blooming plants-honeybees model mutualistic relationship Beddington-DeAngelis functional response dynamical systems stability analysis bifurcation. The author(s) declared that no grants were involved in supporting this work. Amendments from Version 1

Reviewers' remarks prompted minor modifications. The introduction included more details on references 5 and 7 and on the model assumptions and provided a biological interpretation of w. The numerical simulations section presented more confirmatory results, and the conclusion included a paragraph on sustainable honey production and a biological summary.  Finally, version 2 includes the correct author order.

 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. 1 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. 2 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. 3 Hadani 4 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. In, 5 the researcher discusses that the symbiotic relationship is adversely affected by two factors: global warming and the predation of Vespa orientalis on bees. The researchers elaborated on these issues. They observed that the bees' foraging is impeded and the flowering season is adversely affected by the ongoing rise in global temperatures. Consequently, the symbiotic relationship with extinction is at risk. 6 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. In addition, the researchers demonstrated, using mathematical modelling, that a disease outbreak between a plant and its pollinator is associated with declines in pollinator populations worldwide. They also concluded that a better understanding of how land-use change affects species interactions is critical for conserving pollinator health. 7 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. 8 A recent report on the global reduction of honeybee and bumblebee populations has highlighted pollination’s ecological and economic significance. 9 The reduction of pollinator species can be ascribed to various ecological and environmental reasons, including habitat loss, illnesses, climate change, pesticides, and predation. 10– 13 Invasive predators can significantly affect pollinators by diminishing their quantity, altering plant reproductive success, and undermining the plant-pollinator relationship. 14

One factor that negatively affects the mutualism in plant-pollinator systems is high-speed wind. 5, 15– 18 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. 19 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. 20

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 dp dt = r 1 p ( 1 − p k 1 ) + α 1 hp 1 + w + ( ah + bp ) − γ 1 p , dh dt = r 2 h ( 1 − h k 2 ) + α 2 hp 1 + w + ( ah + bp ) − γ 2 h , dn dt = α 3 h 1 + w + ch − βnh − γ 3 n .

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 ( α 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 and w signifies the wind flow.

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 ) = 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., ( α 1 hp 1 + ( ah + bp ) ) .

ϑ ( w ) > 1 means the efficiency of mutualistic interaction decreases in a high, windy environment.

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 . ( α 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.

The term ( α 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.

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.

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.

Flowchart of the <inline-formula> <mml:math display="inline"> <mml:mi>phn</mml:mi> </mml:math> </inline-formula> model. Model analysis

Before analyzing our model, it is pertinent to invoke the following lemmas, the demonstration of which is available in Refs. 21 – 24

If H , ℊ > 0 and H ̇ ⩾ ( ⩽ ) H ( H − ℊ H α ) , where a is a positive constant, when t ⩾ 0 and H ( 0 ) > 0 , then H ( t ) ⩾ ( ⩽ ) ( H ℊ ) 1 / α [ 1 + ( H H − a ( 0 ) ℊ − 1 ) e − H αt ] − 1 / α

(Comparison lemma) Suppose that H , ℊ > 0 , with Ρ ( 0 ) > 0 . Then for dΡ dt ≤ Ρ ( t ) [ H − ℊ Ρ ( t ) ] , then lim t → ∞ sup P ( t ) ≤ H ℊ , and if d Ρ dt ≥ Ρ ( t ) [ H − ℊ Ρ ( t ) ] , then lim t → ∞ inf P ( t ) ≥ 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

The solution ( p ( t ) , h ( t ) , n ( t ) ) of phn system with the initial condition ( p 0 , h 0 , n 0 ) is positive.

According to the blooming plants and honeybees’ equations of phn system, we get p ( t ) = p 0 e ∫ 0 t ( r 1 ( 1 − p ( τ ) k 1 ) + α 1 h ( τ ) 1 + w + ( ah ( τ ) + bp ( τ ) ) − γ 1 ) dτ ≥ 0 h ( t ) = h 0 e ∫ 0 t ( r 2 ( 1 − h ( τ ) k 2 ) + α 2 p ( τ ) 1 + w + ( ah ( τ ) + bp ( τ ) ) − γ 2 ) dτ ≥ 0

From the equation for honey production, we attain

dn dt ≥ − n ( β h 0 e ∫ 0 t ( r 2 ( 1 − h ( τ ) k 2 ) + α 2 p ( τ ) 1 + w + ( ah ( τ ) + bp ( τ ) ) − γ 2 ) dτ + γ 3 ) , which implies that n ( t ) ≥ n 0 e ∫ 0 t ( β h 0 e ∫ 0 t ( r 2 ( 1 − h ( τ ) k 2 ) + α 2 p ( τ ) 1 + w + ( ah ( τ ) + bp ( τ ) ) − γ 2 ) dτ + γ 3 ) dτ ≥ 0

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

Boundedness

All solutions of phn system are uniformly bounded.

From the blooming plants equation of the phn system, we obtain dp dt = r 1 p ( 1 − p k 1 ) + α 1 hp 1 + w + ( ah + bp ) − γ 1 p ≤ p [ ( r 1 − γ 1 ) − ( r 1 p k 1 ) + α 1 h ] .

After using the honeybees’ carrying capacity, we get dp dt ≤ p [ ( r 1 + α 1 k 2 − γ 1 ) − ( r 1 p k 1 ) ]

Using Lemma 1, we deduced that. p ( t ) ≤ [ ( r 1 + α 1 k 2 − γ 1 ) k 1 r 1 ] 1 2 ( 1 + [ ( r 1 + α 1 k 2 − γ 1 ) k 1 p − 2 ( 0 ) r 1 − 1 ] e − 2 ( r 1 + α 1 k 2 − γ 1 ) t ) − 1 2

It is obtained after taking t → ∞ p ( t ) ≤ [ ( r 1 + α 1 k 2 − γ 1 ) k 1 r 1 ] 1 2 = p m

Similarly, for the honey honeybee equation of the phn system, we find h ( t ) ≤ [ k 1 ( r 2 + α 2 p m − γ 2 ) r 2 ] 1 2 = h m

Finally, from the amount of honey produced equation, we attain dn dt = α 3 h 1 + w + ch − βnh − γ 3 n ≤ α 3 h m − γ 3 n

Then, by Lemma 2, we get 0 ≤ lim t → ∞ sup [ n ( t ) ] ≤ α 3 h m γ 3 = n m

Therefore, any solutions of phn system will be attracted to Σ = { ( h , p , n ) ∈ R + 3 : p ( t ) ≤ p m , h ( t ) ≤ h m , n ( t ) ≤ n m } .

Persistence

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

All phn system components are persistent if r 1 > γ 1 , r 2 > γ 2 .

First, to prove that p ( t ) If it is persistent, we have to show that lim t → ∞ inf p ( t ) > 0 i.e. p ( t ) will not decay to zero.

From the blooming plants equation, we get dp dt = r 1 p ( 1 − p k 1 ) + α 1 hp 1 + w + ( ah + bp ) − γ 1 p ≥ p [ ( r 1 − γ 1 ) − ( r 1 k 1 ) p ] .

By Lemma 2, we have lim t → ∞ inf p ( t ) ≥ L 1 , where L 1 = k 1 ( r 1 − γ 1 ) r 1 provided r 1 > γ 1 . Thus, for a small E 1 > 0 , ∃ a positive number t 1 > 0 such that p ( t ) ≥ L 1 − E 1 ∀ t > t 1

Applying the same strategy for the honeybee equation, we get h ( t ) ≥ L 2 − E 2 ∀ t > t 2 , where L 2 = k 2 ( r 2 − γ 2 ) r 2 provided r 2 > γ 2 .

From the honey production equation, we attain dn dt = α 3 h 1 + w + ch − βnh − γ 3 n ≥ α 3 ( L 2 − E 2 ) 1 + w + c ( L 2 − E 2 ) − n ( β h m − γ 3 )

Since

g ( h ) = α 3 h 1 + w + ch is an increasing function, which means if h ( t ) ≥ L 2 − E 2 , then g ( h ) ≥ g ( L 2 − E 2 )

h ( t ) ≤ h m

Then, by Lemma 2, lim t → ∞ inf n ( t ) ≥ α 3 ( L 2 − E 2 ) [ 1 + w + c ( L 2 − E 2 ) ] ( β h m − γ 3 ) , where L 3 = α 3 ( L 2 − E 2 ) [ 1 + w + c ( L 2 − E 2 ) ] ( β h m − γ 3 )

Thus, for a small E 3 > 0 , ∃ a positive number t 3 > 0 such that n ( t ) ≥ L 3 − E 3 ∀ t > t 3

Therefore, phn system is uniformly persistent.

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

The extinction point S 1 = ( 0 , 0 , 0 ) .

The blooming plants point S 2 = ( p 1 , 0 , 0 ) , where p 1 = k 1 ( r 1 − γ 1 ) r 1 . For p 1 to be positive, condition (2) must be satisfied.

The honeybee point S 3 = ( 0 , h 2 , 0 ) , where h 2 = k 2 ( r 2 − γ 2 ) r 2 . Clearly h 2 > 0 if condition (3) is satisfied.

The blooming plants free point S 4 = ( 0 , h 3 , n 3 ) , where h 3 = k 2 ( r 2 − γ 2 ) r 2 and n 3 = α 3 h 3 ( 1 + w + c h 3 ) ( β h 3 + γ 3 ) > 0 . Clearly h 3 > 0 if the condition (3) is satisfied.

The coexistence point S 5 = ( p 4 , h 4 , n 4 ) , here n 4 = α 3 h 4 ( 1 + w + c h 4 ) ( β h 4 + γ 3 ) , h 4 = 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 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: 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 .

Further, h 4 > 0 if one of the following conditions holds: 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

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 ) = [ a 11 α 1 p ( 1 + w + bp ) ( 1 + w + ( ah + bp ) ) 2 0 α 2 h ( 1 + w + ah ) ( 1 + w + ( ah + bp ) ) 2 a 22 0 0 α 3 ( 1 + w ) ( 1 + w + ch ) 2 − βn − βh − γ 3 ] , where a 11 = − 2 r 1 p k 1 + α 1 h ( 1 + w + ah ) ( 1 + w + ( ah + bp ) ) 2 + ( r 1 − γ 1 ) , and a 22 = − 2 r 2 h k 2 + α 2 p ( 1 + w + bp ) ( 1 + w + ( ah + bp ) ) 2 + ( r 2 − γ 2 ) .

The extinction point S 1 = ( 0 , 0 , 0 ) is locally asymptotically stable if γ 1 > r 1 , γ 2 > r 2 .

The Jacobian matrix at S 1 = ( 0 , 0 , 0 ) is J ( S 1 ) = [ r 1 − γ 1 0 0 0 r 2 − γ 2 0 0 ∝ 3 1 + w − γ 3 ] ,

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.

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.

The blooming plants point S 2 = ( p 1 , 0 , 0 ) is locally asymptotically stable if r 2 + α 2 p 1 ( 1 + w + b p 1 ) < γ 2 ,

J ( S 2 ) = J ( p 1 , 0 , 0 ) is given by J ( S 2 ) = [ − ( r 1 − γ 1 ) α 1 p 1 ( 1 + w + b p 1 ) 0 0 α 2 p 1 ( 1 + w + b p 1 ) + ( r 2 − γ 2 ) 0 0 α 3 1 + w − γ 3 ] .

So, the eigenvalues of S 2 are

λ 21 = − ( r 1 − γ 1 ) < 0 , under the existence condition of the blooming plants point.

λ 22 = α 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 .

The honeybee point S 3 = ( 0 , h 2 , 0 ) is locally asymptotically stable if r 1 + α 1 h 2 ( 1 + w + a h 2 ) < γ 1 , ( r 2 − γ 2 ) < 2 r 2 h 2 k 2

J ( S 3 ) = J ( 0 , h 2 , 0 ) is given by J ( S 3 ) = [ α 1 h 2 ( 1 + w + a h 2 ) + ( r 1 − γ 1 ) 0 0 α 2 h 2 1 + w + a h 2 ( r 2 − γ 2 ) − 2 r 2 h 2 k 2 0 0 α 3 ( 1 + w ) ( 1 + w + c h 2 ) 2 − β h 2 − γ 3 ] .

So, the eigenvalues of J ( S 3 ) are λ 31 = α 1 h 2 ( 1 + w + a h 2 ) + ( r 1 − γ 1 ) , λ 32 = ( r 2 − γ 2 ) − 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.

The blooming plants free point S 4 = ( 0 , h 3 , n 3 ) is locally asymptotically stable if r 1 + α 1 h 3 ( 1 + w + a h 3 ) < γ 1 , ( r 2 − γ 2 ) < r 2 h 3 k 2

J ( S 4 ) = J ( 0 , h 3 , n 3 ) is given by J ( S 4 ) = [ α 1 h 3 ( 1 + w + a h 3 ) + ( r 1 − γ 1 ) 0 0 α 2 h 3 ( 1 + w + a h 3 ) ( r 2 − γ 2 ) − 2 r 2 h 3 k 2 0 0 α 3 ( 1 + w ) ( 1 + w + c h 3 ) − β n 3 − β h 3 − γ 3 ] .

The eigenvalues of J ( S 4 ) are λ 41 = α 1 h 3 ( 1 + w + a h 3 ) + ( r 1 − γ 1 ) , λ 42 = ( r 2 − γ 2 ) − 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.

The coexistence point S 5 = ( p 4 , h 4 , n 4 ) is locally asymptotically stable if a ii < 0 , i = 1 , 2 , a 11 a 22 > [ α 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 ] .

J ( S 5 ) = J ( p 4 , h 4 , n 4 ) is given by J ( S 5 ) = [ a 11 α 1 p 4 ( 1 + w + b p 4 ) ( 1 + w + ( a h 4 + b p 4 ) ) 2 0 α 2 h 4 ( 1 + w + a h 4 ) ( 1 + w + ( a h 4 + b p 4 ) ) 2 a 22 0 0 α 3 ( 1 + w ) ( 1 + w + c h 4 ) 2 − β n 4 − β h 4 − γ 3 ] , where a 11 = ( r 1 − γ 1 ) − 2 r 1 p 4 k 1 + α 1 h 4 ( 1 + w + a h 4 ) ( 1 + w + ( a h 4 + b p 4 ) ) 2 , and a 22 = ( r 2 − γ 2 ) − 2 r 2 h 4 k 2 + α 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: ( − β h 4 − γ 3 − λ ) [ λ 2 − Trλ + Det ] = 0 ,

where, λ 51 = − β h 4 − γ 3 , Tr = a 11 + a 22 = ( r 1 − γ 1 ) − 2 r 1 p 4 k 1 + α 1 h 4 ( 1 + w + a h 4 ) ( 1 + w + ( a h 4 + b p 4 ) ) 2 − 2 r 2 h 4 k 2 + α 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 − [ α 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 ) − 2 r 1 p 4 k 1 + α 1 h 4 ( 1 + w + a h 4 ) ( 1 + w + ( a h 4 + b p 4 ) ) 2 ) ( ( r 2 − γ 2 ) − 2 r 2 h 4 k 2 + α 2 p 4 ( 1 + w + b p 4 ) ( 1 + w + ( a h 4 + b p 4 ) ) 2 ) − [ α 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 ] .

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.

The extinction point S 1 = ( 0 , 0 , 0 ) is a global asymptotic stability (GAS) if the following conditions are met. r 1 + h m ( α 1 + α 2 ) 1 + w < γ 1 r 2 + α 3 1 + w < γ 2

Let E 1 ( p , h , n ) = p + h + n , where E 1 : R + 3 → R , which satisfies E 1 ( 0 , 0 , 0 ) = 0 and E 1 ( p , h , n ) > 0 for all ( p , h , n ) ∈ R + 3 with ( p , h , n ) ≠ ( 0 , 0 , 0 ) . Then d E 1 dt = dp dt + dh dt + dn dt = ( r 1 p ( 1 − p k 1 ) + α 1 hp 1 + w + ( ah + bp ) − γ 1 p ) + ( r 2 h ( 1 − h k 2 ) + α 2 hp 1 + w + ( ah + bp ) − γ 2 h ) + ( α 3 h 1 + w + ch ) − βnh − γ 3 n ) = r 1 p − r 1 p 2 k 1 + hp ( α 1 + α 2 ) 1 + w + ( ah + bp ) − γ 1 p + r 2 h − r 2 h 2 k 2 − γ 2 h + α 3 h 1 + w + ch ) − βnh − γ 3 n

Then, by using the upper bound of the honeybees’ population, we get d E 1 dt ≤ p ( r 1 − γ 1 + h m ( α 1 + α 2 ) 1 + w ) + h ( r 2 − γ 2 + α 3 1 + w ) − r 1 p 2 k 1 − r 2 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.

The blooming plants point S 2 = ( p 1 , 0 , 0 ) is GAS if r 2 + α 2 p m + α 3 1 + w < γ 2 , ( α 1 D ) 2 ≤ r 1 r 2 k 1 k 2 , where, D = 1 + w + ( ah + bp ) .

Let E 2 = ( p − p 1 − p 1 ln ( p p 1 ) ) + h + n , where E 2 : R + 3 → R , which satisfies E 2 ( p 1 , 0 , 0 ) = 0 and E 2 ( p , h , n ) > 0 for all ( p , h , n ) ∈ R + 3 with ( p , h , n ) ≠ ( p 1 , 0 , 0 ) , then d E 2 dt = ( p − p 1 p ) dp dt + dh dt + dn dt = − r 1 k 1 ( p − p 1 ) 2 + α 1 h ( p − p 1 ) D + h r 2 ( 1 − h k 2 ) + α 2 ph D − γ 2 h + α 3 h 1 + w + ch − βnh − γ 3 n

Then, by using the upper bound of the blooming plants population, we get d E 2 dt ≤ − ( r 1 k 1 ( p − p 1 ) 2 − α 1 D h ( p − p 1 ) + r 2 k 2 h 2 ) + h ( r 2 + α 2 p m + α 3 1 + w − γ 2 ) − βnh − γ 3 n .

Thus, d E 2 dt ≤ − ( r 1 k 1 ( p − p 1 ) + r 2 k 2 h ) 2 + h ( r 2 + α 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.

The honeybee point S 3 = ( 0 , h 2 , 0 ) , is GAS if r 1 + α 1 h m 1 + w < γ 1 , ( α 2 D ) 2 ≤ r 1 r 2 k 1 k 2 , α 3 β ( 1 + w ) < n , where, D = 1 + w + ( ah + bp ) .

Let E 3 = p + ( h − h 2 − h 2 ln h h 2 ) + n , where E 3 : R + 3 → R , which satisfies E 3 ( 0 , h 2 , 0 ) = 0 and E 3 ( p , h , n ) > 0 for all ( p , h , n ) ∈ R + 3 with ( p , h , n ) ≠ ( 0 , h 2 , 0 ) , then d E 3 dt = dp dt + ( h − h 2 h ) dh dt + dn dt = ( r 1 p ( 1 − p k 1 ) + α 1 hp 1 + w + ( ah + bp ) − γ 1 p ) + ( h − h 2 ) ( ( r 2 ( 1 − h k 2 ) + α 2 p 1 + w + ( ah + bp ) − γ 2 ) ) + ( α 3 h 1 + w + ch − βnh − γ 3 n ) ≤ r 1 p − r 1 p 2 k 1 + α 1 hp 1 + w − γ 1 p − r 2 k 2 ( h − h 2 ) 2 + α 2 ( h − h 2 ) p D + α 3 h 1 + w − βnh − γ 3 n

Then, by using the upper bound of the honeybees’ population, we get d E 3 dt ≤ − ( r 1 k 1 p 2 − α 2 ( h − h 2 ) p D + ( r 2 k 2 ) ( h − h 2 ) 2 ) + p ( r 1 + α 1 h m 1 + w − γ 1 ) + h ( α 3 1 + w − βn ) − γ 3 n

Thus, d E 3 dt ≤ − ( r 1 k 1 p + r 2 k 2 ( h − h 2 ) ) 2 + p ( r 1 + α 1 h m 1 + w − γ 1 ) + h ( α 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 )

The blooming plants free point S 4 = ( 0 , h 3 , n 3 ) is GAS if ( α 2 D ) 2 ≤ r 1 r 2 2 k 1 k 2 , ( α 3 ( 1 + w ) D 1 D 2 − βn ) 2 ≤ 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 ) .

Let E 4 = p + ( h − h 3 − h 3 ln h h 3 ) + ( n − n 3 2 ) 2 , where E 4 : R + 3 → R , which satisfies E 4 ( 0 , h 3 , n 3 ) = 0 and E 4 ( p , h , n ) > 0 for all ( p , h , n ) ∈ R + 3 with ( p , h , n ) ≠ ( 0 , h 3 , n 3 ) , then d E 4 dt = dp dt + ( h − h 3 h ) dh dt + ( n − n 3 ) dn dt = ( r 1 p ( 1 − p k 1 ) + α 1 hp 1 + w + ( ah + bp ) − γ 1 p ) + ( h − h 3 ) ( r 2 ( 1 − h k 2 ) + α 2 p 1 + w + ( ah + bp ) − γ 2 ) + ( n − n 3 ) ( α 3 h 1 + w + ch − βnh − γ 3 n ) = r 1 p − r 1 p 2 k 1 + α 1 hp 1 + w + ( ah + bp ) − γ 1 p + ( h − h 3 ) ( − r 2 k 2 ( h − h 3 ) + α 2 p 1 + w + ( ah + bp ) ) + ( n − n 3 ) ( α 3 h 1 + w + ch − α 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 d E 4 dt ≤ − ( r 1 k 1 p 2 − α 2 ( h − h 3 ) p 1 + w + ( ah + bp ) + ( r 2 2 k 2 ) ( h − h 3 ) 2 ) + p ( r 1 + α 1 h m 1 + w − γ 1 ) − [ ( r 2 2 k 2 ) ( h − h 3 ) 2 − ( α 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, d E 4 dt ≤ − ( r 1 k 1 p + r 2 2 k 2 ( h − h 3 ) ) 2 + p ( r 1 + α 1 h m 1 + w − γ 1 ) − ( r 2 2 k 2 ( h − h 3 ) + ( β 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.

The coexistence point S 5 = ( p 4 , h 4 , n 4 ) is GAS if ( ( α 1 + α 2 ) ( 1 + w ) + α 1 b p 4 + α 2 a h 4 N 1 N 2 ) 2 ≤ 1 2 ( r 1 k 1 + α 1 b h 4 N 1 N 2 ) ( r 2 k 2 + α 2 a p 4 N 1 N 2 ) , ( α 3 ( 1 + w ) D 1 D 2 − βn ) 2 ≤ 1 2 ( r 2 k 2 + α 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 p p 4 ) + ( h − h 4 − h 4 ln h h 4 ) + ( n − n 4 2 ) 2 , where E 5 : R + 3 → R , which satisfies E 5 ( p 4 , h 4 , n 4 ) = 0 and E 5 ( p , h , n ) > 0 for all ( p , h , n ) ∈ R + 3 with ( p , h , n ) ≠ ( p 4 , h 4 , n 4 ) , then d E 5 dt = ( p − p 4 p ) dp dt + ( h − h 4 h ) dh dt + ( n − n 4 ) dn dt = ( p − p 4 ) ( r 1 ( 1 − p k 1 ) + α 1 h 1 + w + ( ah + bp ) − γ 1 ) + ( h − h 4 ) ( r 2 ( 1 − h k 2 ) + α 2 p 1 + w + ( ah + bp ) − γ 2 ) + ( n − n 4 ) ( α 3 h 1 + w + ch − βnh − γ 3 n ) d E 5 dt = − [ ( r 1 k 1 + α 1 b h 4 N 1 N 2 ) ( p − p 4 ) 2 − ( ( α 1 + α 2 ) ( 1 + w ) + α 1 b p 4 + α 2 a h 4 N 1 N 2 ) ( p − p 4 ) ( h − h 4 ) + 1 2 ( r 2 k 2 + α 2 a p 4 N 1 N 2 ) ( h − h 4 ) 2 ] − [ 1 2 ( r 2 k 2 + α 2 a p 4 N 1 N 2 ) ( h − h 4 ) 2 − ( α 3 ( 1 + w ) D 1 D 2 − βn ) ( n − n 4 ) ( h − h 4 ) + ( β h 4 + γ 3 ) ( n − n 4 ) 2 ]

Thus, d E 5 dt ≤ − ( ( r 1 k 1 + α 1 b h 4 N 1 N 2 ) ( p − p 4 ) + 1 2 ( r 2 k 2 + α 2 a p 4 N 1 N 2 ) ( h − h 4 ) ) 2 − ( 1 2 ( r 2 k 2 + α 2 a p 4 N 1 N 2 ) ( h − h 4 ) + ( β 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. 24 – 27

For r 2 ∗ = γ 2 , the phn model faces TB at the extinction point S 1 = ( 0 , 0 , 0 ) .

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 ) = [ r 1 − γ 0 0 0 0 0 0 α 3 1 + w − γ 3 ]

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 , α 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 − h k 2 ) , 0 ) ⇒ 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 ) [ 0 0 0 0 1 0 0 0 0 ] = 1 ≠ 0 ( T [ 1 ] ) T D 2 F r 2 ( S 1 , r 2 ∗ ) ( ϑ [ 1 ] , ϑ [ 1 ] ) = ( 0 , 1 , 0 ) [ 0 − 2 r 2 ∗ k 2 0 ] = − 2 r 2 ∗ k 2 ≠ 0

Therefore, there is a TB around S 1 with the parameter r 2 ∗ = γ 2 .

For r 1 ∗ = γ 1 , the phn model faces TB at the extinction point S 2 = ( p 1 , 0 , 0 ) if p 1 = k 1 , ( T [ 2 ] ) T D 2 F r 1 ( S 2 , r 1 ∗ ) ( ϑ [ 2 ] , ϑ [ 2 ] ) ≠ 0

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 ) = [ 0 0 0 0 r 2 − γ 2 0 0 α 3 1 + w − γ 3 ] ,

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 , α 3 γ 3 ( 1 + w ) ) T and ( T [ 2 ] ) T = ( 1 , α 3 ( 1 + w ) ( γ 2 − r 2 ) , 1 ) T , where ( γ 2 − r 2 ) ≠ 0 .

Further, F r 1 ( S , r 1 ) = ( p ( 1 − p k 1 ) , 0 , 0 ) ⟹ F r 1 ( S 2 , r 1 ∗ ) = ( p 1 ( 1 − p 1 k 1 ) , 0 , 0 ) ( T [ 2 ] ) T F r 1 ( S 2 , r 1 ∗ ) = ( 1 , ∝ 3 ( 1 + w ) ( γ 2 − r 2 ) , 1 ) ( p 1 − p 1 2 k 1 , 0 , 0 ) = p 1 ( 1 − p 1 k 1 ) = 0 under condition 30 . ( T [ 2 ] ) T DF r 1 ( S 2 , r 1 ∗ ) ( ϑ [ 2 ] ) = ( 1 , α 3 ( 1 + w ) ( γ 2 − r 2 ) , 1 ) [ 1 0 0 0 0 0 0 0 0 ] = 1 ≠ 0

( T [ 2 ] ) T D 2 F r 1 ( S 2 , r 1 ∗ ) ( ϑ [ 2 ] , ϑ [ 2 ] ) = ( 1 , α 3 ( 1 + w ) ( γ 2 − r 2 ) , 1 ) [ X 11 [ 2 ] + X 12 [ 2 ] X 21 [ 2 ] + X 22 [ 2 ] X 32 [ 2 ] + X 33 [ 2 ] ∝ 3 γ 3 ( 1 + w ) ]

= ( X 11 [ 2 ] + X 12 [ 2 ] ) + α 3 ( 1 + w ) ( γ 2 − r 2 ) ( X 21 [ 2 ] + X 22 [ 2 ] ) + ( X 32 [ 2 ] + X 33 [ 2 ] α 3 γ 3 ( 1 + w ) ) ≠ 0

under condition 31, here, X 11 [ 2 ] = − 2 r 1 k 1 + [ α 1 [ ( 1 + w ) ] ( 1 + w + b p 1 ) 2 ] , X 12 [ 2 ] = [ α 1 ( 1 + w ) ( 1 + w + b p 1 ) 2 ] + [ − 2 a α 1 p 1 ( 1 + w + b p 1 ) 2 ] , X 21 [ 2 ] = [ α 2 ( 1 + w ) ( 1 + w + b p 1 ) 2 ] , X 22 [ 2 ] = [ α 2 [ ( 1 + w ) ] ( 1 + w + b p 1 ) 2 ] + [ − 2 r 2 k 2 − 2 a α 2 p 1 ( 1 + w + b p 1 ) 2 ] , X 32 [ 2 ] = − 2 c α 3 ( 1 + w ) 2 − β α 3 γ 3 ( 1 + w ) , and X 33 [ 2 ] = − β α 3 γ 3 ( 1 + w ) .

Therefore, there is a TB around S 2 with the parameter r 1 ∗ = γ 1 .

For r 2 ∗ ∗ = γ 2 k 2 k 2 − 2 h 2 , the phn model faces TB at the honeybee point S 3 = ( 0 , h 2 , 0 ) if h 2 = k 2 ,

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 ) = [ α 1 h 2 1 + w + a h 2 + ( r 1 − γ 1 ) 0 0 α h 2 1 + w + a h 2 0 0 0 ( 1 + w ) α 3 ( 1 + w + c h 2 ) 2 − ( β h 2 + γ 3 ) ] .

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 , α 3 ( 1 + w ) ( β h 2 + γ 3 ) ( 1 + w + c h 2 ) 2 ) and ( T [ 3 ] ) T = ( − α 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 ) ] ≠ 0 .

Further F r 2 ( S 3 , r 2 ) = ( 0 , h − h k 2 , 0 ) ⇒ F r 2 ( S 3 , r 2 ∗ ∗ ) = ( 0 , h 2 ( 1 − h 2 k 2 ) , 0 ) ( T [ 3 ] ) T F r 2 ( S 3 , r 2 ∗ ∗ ) = ( − α 2 h 2 ( α 1 h 2 + ( r 1 − γ 1 ) ( 1 + w + a h 2 ) ) , 1 , 0 ) ( 0 , h 2 ( 1 − h 2 k 2 ) , 0 ) ( T [ 3 ] ) T F r 2 ( S 3 , r 2 ∗ ∗ ) = h 2 ( 1 − h 2 k 2 ) = 0 under condition 32 . ( T [ 3 ] ) T DF r 2 ( S 3 , r 2 ∗ ∗ ) ( ϑ [ 3 ] ) = ( − α 2 h 2 ( α 1 h 2 + ( r 1 − γ 1 ) ( 1 + w + a h 2 ) ) , 1 , 0 ) [ 0 0 0 0 1 0 0 0 0 ] = 1 ≠ 0 ( T [ 3 ] ) T D 2 F r 2 ( S 3 , r 2 ∗ ∗ ) ( ϑ [ 3 ] , ϑ [ 3 ] ) = ( − α 2 h 2 ( α 1 h 2 + ( r 1 − γ 1 ) ( 1 + w + a h 2 ) ) , 1 , 0 ) [ 0 − 2 r 2 ∗ ∗ k 2 X 32 [ 3 ] − β α 3 ( 1 + w ) ( β h 2 + γ 2 ) ( 1 + w + c h 2 ) 2 ] = − 2 r 2 ∗ ∗ k 2 ≠ 0 , where X 32 [ 3 ] = − 2 c α 3 ( 1 + w ) ( 1 + w + c h 2 ) 3 − β α 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 ∗ ∗ .

For r 2 ∗ ∗ ∗ = γ 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 h 3 = k 2 ,

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 ) = [ a 11 0 0 a 21 0 0 0 a 32 − β h 3 − γ 3 ] , where, a 11 = α 1 h 3 ( 1 + w + a h 3 ) + ( r 1 − γ 1 ) , a 21 = α 2 h 3 1 + w + a h 3 and a 32 = α 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 , a 32 β h 3 + γ 3 ) and ( T [ 4 ] ) T = ( − a 21 a 11 , 1 , 0 ) .

Further, F r 2 ( S , r 2 ) = ( 0 , h ( 1 − h k 2 ) , 0 ) ⇒ F r 2 ( S 4 , r 2 ∗ ∗ ∗ ) = ( 0 , h 2 − h 2 k 2 , 0 ) ( T [ 4 ] ) T F r 2 ( S 3 , r 2 ∗ ) = ( − a 21 a 11 , 1 , 0 ) ( 0 , h 3 2 − h 3 2 k 2 , 0 ) = h 3 ( 1 − h 3 k 2 ) = 0 under condition 33 . ( T [ 4 ] ) T DF r 2 ( S 4 , r 2 ∗ ) ( ϑ [ 4 ] ) = ( − a 21 a 11 , 1 , 0 ) [ 0 0 0 0 1 0 0 0 0 ] = 1 ≠ 0 ( T [ 4 ] ) T D 2 F r 2 ( S 4 , r 2 ∗ ) ( ϑ [ 4 ] , ϑ [ 4 ] ) = ( − a 21 a 11 , 1 , 0 ) [ 0 − 2 r 2 ∗ k 2 − 2 c α 3 ( 1 + w ) ( 1 + w + c h 3 ) 3 + β ( a 32 β h 3 + γ 3 ) 2 ] = − 2 r 2 ∗ k 2 ≠ 0 ,

Therefore, there is a TB around S 4 with the parameter r 2 ∗ ∗ ∗ .

The phn undergoes a Hopf bifurcation at the coexistence point S 5 to the bifurcation parameter w ∗ if [ Det ] | w = w ∗ > 0 , 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 ) , where the formula of w ∗ is given in the proof.

The Jacobian matrix at S 5 with w ∗ is given by J ( S 5 , w ∗ ) = [ a 11 α 1 p 4 ( 1 + w ∗ + b p 4 ) ( 1 + w ∗ + ( a h 4 + b p 4 ) ) 2 0 α 2 h 4 ( 1 + w ∗ + a h 4 ) ( 1 + w ∗ + ( a h 4 + b p 4 ) ) 2 a 22 0 0 α 3 ( 1 + w ∗ ) ( 1 + w ∗ + c h 4 ) 2 − β n 4 − β h 4 − γ 3 ] ,

where a 11 = ( r 1 − γ 1 ) − 2 r 1 p 4 k 1 + α 1 h 4 ( 1 + w ∗ + a h 4 ) ( 1 + w ∗ + ( a h 4 + b p 4 ) ) 2 , and a 22 = ( r 2 − γ 2 ) − 2 r 2 h 4 k 2 + α 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.

[ Tr ] | w = w ∗ = 0 ,

[ Det ] | w = w ∗ > 0 ,

∂ ∂w [ Re ( λ 1 , 2 ) ] w = w ∗ ≠ 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 ) − 2 r 1 p 4 k 1 + α 1 h 4 ( 1 + w + a h 4 ) ( 1 + w + ( a h 4 + b p 4 ) ) 2 − 2 r 2 h 4 k 2 + α 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 ) − 2 r 1 p 4 k 1 − 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 + 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, ( 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 .

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 = − ( 2 b 1 b 2 + b 4 ) ± Δ 2 b 1 , w ∗ = v − 1 , b 1 ≠ 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 = ± i 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 = Tr ± i Tr 2 − 4 Det 2 , then ∂ ∂w [ Re ( λ 2 , 3 ) ] w = w ∗ = 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 ≠ 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. 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.

The nullclines of <inline-formula> <mml:math display="inline"> <mml:mi>phn</mml:mi> </mml:math> </inline-formula> system with the parameters listed in <xref ref-type="disp-formula" rid="e95">Eq. (36)</xref>. The global stability of <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>S</mml:mi> <mml:mn>5</mml:mn> </mml:msub> <mml:mo>=</mml:mo> <mml:mrow> <mml:mo stretchy="true">(</mml:mo> <mml:mn>0.21</mml:mn> <mml:mo>,</mml:mo> <mml:mn>0.11</mml:mn> <mml:mo>,</mml:mo> <mml:mn>0.23</mml:mn> <mml:mo stretchy="true">)</mml:mo> </mml:mrow> </mml:math> </inline-formula>.

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 ≤ 0.74 , the solution of the phn system converges to a limit cycle, see Figure 4(b) and Figure 5. It means that for 0.017<w≤0.74, the conditions listed in Theorem 14 are satisfied, and the phn undergoes a Hopf bifurcation in the mentioned region. 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.

The effect of various wind speeds <inline-formula> <mml:math display="inline"> <mml:mi>w</mml:mi> </mml:math> </inline-formula>. Hopf bifurcation at <inline-formula> <mml:math display="inline"> <mml:mi>w</mml:mi> <mml:mo>=</mml:mo> <mml:mn>0.74</mml:mn> </mml:math> </inline-formula>.

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 ≤ 0.03 . This result confirms the occurrence of a transcritical bifurcation at, r 2 = γ 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.

The effect of varying <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>r</mml:mi> <mml:mn>2</mml:mn> </mml:msub> </mml:math> </inline-formula>. Transcritical bifurcation at <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>r</mml:mi> <mml:mn>2</mml:mn> </mml:msub> <mml:mo>=</mml:mo> <mml:mn>0.03</mml:mn> </mml:math> </inline-formula>. The global stability of <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>S</mml:mi> <mml:mn>1</mml:mn> </mml:msub> </mml:math> </inline-formula>.

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 ≤ 0.04 , which means the system faces an occurrence of transcritical bifurcations at r 1 = γ 2 = r 1 TB = 0.04, which confirms the output of Theorem 14. So, phn system loses two of its components for r 1 ≤ 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.

The effect of varying <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>r</mml:mi> <mml:mn>1</mml:mn> </mml:msub> </mml:math> </inline-formula>. The global stability of <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>S</mml:mi> <mml:mn>2</mml:mn> </mml:msub> </mml:math> </inline-formula>.

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.

The effect of varying <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>α</mml:mi> <mml:mn>1</mml:mn> </mml:msub> </mml:math> </inline-formula>. The effect of varying <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>α</mml:mi> <mml:mn>2</mml:mn> </mml:msub> </mml:math> </inline-formula>.

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.

The global stability of <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>S</mml:mi> <mml:mn>3</mml:mn> </mml:msub> </mml:math> </inline-formula> when <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>γ</mml:mi> <mml:mn>1</mml:mn> </mml:msub> <mml:mo>=</mml:mo> <mml:mn>0.57</mml:mn> </mml:math> </inline-formula>.

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.

The effect of varying <inline-formula> <mml:math display="inline"> <mml:mi>β</mml:mi> </mml:math> </inline-formula>.

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.

The sensitivity of the data given in <xref ref-type="disp-formula" rid="e95">Eq. 36</xref> relative to the <inline-formula> <mml:math display="inline"> <mml:msub> <mml:mi>S</mml:mi> <mml:mn>5</mml:mn> </mml:msub> <mml:mo>=</mml:mo> <mml:mrow> <mml:mo stretchy="true">(</mml:mo> <mml:msub> <mml:mi>p</mml:mi> <mml:mn>4</mml:mn> </mml:msub> <mml:mo>,</mml:mo> <mml:msub> <mml:mi>h</mml:mi> <mml:mn>4</mml:mn> </mml:msub> <mml:mo>,</mml:mo> <mml:msub> <mml:mi>n</mml:mi> <mml:mn>4</mml:mn> </mml:msub> <mml:mo stretchy="true">)</mml:mo> </mml:mrow> </mml:math> </inline-formula>. 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. In the context of the existing literature 5, as plant–pollinator mutualism is affected by high levels of global warming, the populations of flowering plants and honey bees are at risk of extinction. In addition, the plant and its pollinator could coexist with limit-cycle behavior if the pollination rate by honey bees is reduced. Additionally, Biswas 6 has shown that plant–pollinator mutualism does not coexist without pollinator behavior modification in the presence of predators, whereas this system can coexist without a predator population. Furthermore, this study emphasizes the impact of high wind speed on the plant–pollinator mutualism model. The study indicates that wind flow could be an important ecological factor, as it has the potential to reduce the frequency of honey bee visits to flowers substantially. Controlling habitat loss is crucial for maintaining the stability and persistence of bee populations, which are directly linked to sustainable honey production. In ecological systems, habitat degradation reduces floral resources and nesting sites, leading to a decline in pollinator abundance. The bifurcation analysis presented in this study shows that appropriate parameter management for blooming plant ( r 1 ), honeybee growth rates ( r 2 ), and wind speed ( w) can mitigate these negative effects by stabilizing population dynamics and preventing critical transitions, such as population collapse or extinction. Appropriate management of the above parameter can maintain the system in a stable regime, ensuring continuous pollination and honey yield. Therefore, strategies aimed at reducing habitat loss not only preserve biodiversity but also contribute to the long-term sustainability of apiculture and agricultural productivity. For future work, we will investigate how pesticide exposure affects the relationship between plants and pollinators.

 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 Dean AM : A simple model of mutualism. Am. Nat. 1983;121(3):409–417. 10.1086/284069 Ollerton J Winfree R Tarrant S : How many flowering plants are pollinated by animals? Oikos. 2011;120(3):321–326. 10.1111/j.1600-0706.2010.18644.x 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. 39122794 10.1038/s41598-024-69537-7 PMC11316071 Fishman MA Hadany L : Plant–pollinator population dynamics. Theor. Popul. Biol. 2010;78(4):270–277. 10.1016/j.tpb.2010.08.002 Jawad S Thirthar AA Nisar KS : The impact of climate change on flowering plants-bees-Vespa orientalis model. Results Control Optim. 2025;20:100583. 10.1016/j.rico.2025.100583 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). 39899584 10.1063/5.0233838 Figueroa LL Grab H Ng WH : Landscape simplification shapes pathogen prevalence in plant‐pollinator networks. Ecol Lett. 2020;23(8):1212–1222. Hakeem E Jawad S Ali AH : How mathematical models might predict desertification from global warming and dust pollutants. MethodsX. 2025;14:103259. 40165851 10.1016/j.mex.2025.103259 PMC11957603 Lever JJ Nes EH van Scheffer M : The sudden collapse of pollinator communities. Ecol. Lett. 2014;17(3):350–359. 24386999 10.1111/ele.12236 Al Nuaimi M Jawad S : Modelling and stability analysis of the competitional ecological model with harvesting. Commun. Math. Biol. Neurosci. 2022;2022:Article-ID. Shalan RN Shireen R Lafta AH : Discrete an SIS model with immigrants and treatment. J. Interdiscip. Math. 2020;24:1201–1206. 10.1080/09720502.2020.1814496 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. 10.24996/ijs.2024.65.5.31 Nezar F Jawad S Winter M : 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. Cresswell JE : A demographic approach to evaluating the impact of stressors on bumble bee colonies. Ecol. Entomol. 2017;42(2):221–229. 10.1111/een.12376 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. 10.1016/j.padiff.2025.101189 Jawad S Roy S Thirthar AA : Deterministic and stochastic risks assessment of excessive CO2 emission on forest biomass under weak allee effect. J. Appl. Math. Comput. 2025;71:9129–9156. 10.1007/s12190-025-02640-8 Thirthar AA Alaoui AL Roy S : Fractional and stochastic dynamics of predator–prey systems: The role of fear and global warming. Eur. Phys. J. B. 2025;98(7):1–21. 10.1140/epjb/s10051-025-00992-5 Ahmed M Jawad S Das D : Impact of dust storms on plant biomass: Model structure and dynamic study. Alex. Eng. J. 2025;126:605–622. 10.1016/j.aej.2025.04.058 Balfour NJ Ratnieks FLW : Wind Alters Plant-Pollinator Community Structure, Bee Foraging Rate & Movements Between Plants. Behav. Ecol. 2025;36:araf067. 40666820 10.1093/beheco/araf067 PMC12260157 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. 19218583 10.1093/aob/mcp035 PMC2701749 Chen F : On a nonlinear nonautonomous predator–prey model with diffusion and distributed delay. J. Comput. Appl. Math. 2005;180(1):33–49. 10.1016/j.cam.2004.10.001 Place CM : Dynamical Systems: Differential Equations, Maps, and Chaotic Behaviour. Routledge;2017. Aamer Z Jawad S Batiha B : Evaluation of the Dynamics of Psychological Panic Factor, Glucose Risk and Estrogen Effects on Breast Cancer Model. Computation. 2024;12(8):160. 10.3390/computation12080160 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. Thirthara AA Jawadb S Shahc K : How does media coverage affect a COVID-19 pandemic model with direct and indirect transmission? 2024. Javaid Y Jawad S Ahmed R : Dynamic complexity of a discretized predator-prey system with Allee effect and herd behaviour. Appl. Math. Sci. Eng. 2024;32(1):2420953. 10.1080/27690911.2024.2420953 Roy S AL-Jaf DS Thirthar AA : 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. 10.1016/j.matcom.2025.11.016 10.5256/f1000research.198543.r480277 Reviewer response for version 2 Seralan Vinoth 1 Referee SRM Institute of Science and Technology, Chennai, Tamil Nadu, India

Competing interests: No competing interests were disclosed.

 12 5 2026 Copyright: © 2026 Seralan V 2026 This is an open access peer review report distributed under the terms of the Creative Commons Attribution Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. recommendation approve

I appreciate the authors for carried out all my comments and suggestions.

Is the work clearly and accurately presented and does it cite the current literature?

Yes

If applicable, is the statistical analysis and its interpretation appropriate?

No

Are all the source data underlying the results available to ensure full reproducibility?

No

Is the study design appropriate and is the work technically sound?

Yes

Are the conclusions drawn adequately supported by the results?

Yes

Are sufficient details of methods and analysis provided to allow replication by others?

Yes

Reviewer Expertise:

Mathematical biology, chaos theory

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

 10.5256/f1000research.189832.r464608 Reviewer response for version 1 Seralan Vinoth 1 Referee SRM Institute of Science and Technology, Chennai, Tamil Nadu, India

Competing interests: No competing interests were disclosed.

 25 3 2026 Copyright: © 2026 Seralan V 2026 This is an open access peer review report distributed under the terms of the Creative Commons Attribution Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. recommendation approve-with-reservations

In this work, the authors studied the dynamics of interactions between blooming plants and honeybee populations, which is crucial for addressing honeybee decline and sustainable ecosystems. The mathematical model on dynamics of a blooming plant, honeybee population, and honey production symbiosis, which is special emphasis on the effect of high-speed wind flow were investigated. Furthermore, the existence of two difference bifurcation phenomena such as the transcritical and Hopf bifurcations have been discussed.   I have few comments as follows:

In abstract, the authors mentioned, Habitat loss, diseases, climate change, pesticides, and predation have all contributed to the decline of pollinator species. Although, it has been discussed in the manuscript. More discussions to be added.

There are some typos, please correct it. For example. blooming plants p, honey honeybees h (honey is repeated). The biological meaning of the parameter w is not mentioned.

The derivation of theoretical bifurcation condition both bifurcations are correct but add a proof write more lines that are numerically satisfied.

In phase portrait Figure 4 (a), it has been shown as chaotic trajectories. But in bifurcation diagram Figure 5, its look like a periodic attractor.

The numerical simulation section provides strong supporting evidence whenever needed, which is a positive aspect. The prescribed-time results need more biologically justified.

The authors suggested to compare their finding by referring some related papers about implementing control of habitat loss in a separate remark. Discuss about sustainable production of honey.

Please summary the conclusions with more biological justification and point out how different the model is from other existing models.

Is the work clearly and accurately presented and does it cite the current literature?

Yes

If applicable, is the statistical analysis and its interpretation appropriate?

No

Are all the source data underlying the results available to ensure full reproducibility?

No

Is the study design appropriate and is the work technically sound?

Yes

Are the conclusions drawn adequately supported by the results?

Yes

Are sufficient details of methods and analysis provided to allow replication by others?

Yes

Reviewer Expertise:

Mathematical biology, chaos theory

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

Jawad Shireen Mathematics, University of Baghdad, Baghdad, Baghdad Governorate, Iraq

Competing interests: No competing interests were disclosed.

 6 4 2026

Thank you very much for giving us the opportunity to revise our manuscript. We really appreciate the helpful comments from the learned reviewers, which helped improve our manuscript. We have revised the manuscript in response to each comment and provided our responses below. We sincerely hope the revised manuscript meets your expectations, and we look forward to hearing from you at your convenience.

In abstract, the authors mentioned, Habitat loss, diseases, climate change, pesticides, and predation have all contributed to the decline of pollinator species. Although, it has been discussed in the manuscript. More discussions to be added.

Response: Thank you for your comment. More discussion of all contributions to the decline of pollinator species has been added. Please see the introduction section. We consider the impact of pesticides in our next study; we mention this at the end of the conclusion section.

There are some typos, please correct it. For example. blooming plants p, honey honeybees h (honey is repeated). The biological meaning of the parameter w is not mentioned.

 Response: Thank you very much for your comment. The repeated word ‘honey’ has been removed. The biological meaning of w is clearly presented in numbers 1-2 of the “Assumptions of the model” section.

The derivation of theoretical bifurcation condition both bifurcations are correct but add a proof write more lines that are numerically satisfied.

 Response: Thank you for your valuable comment. We are glad that the reviewer finds the derivation of the theoretical bifurcation conditions to be correct. In response to your suggestion, we have expanded the numerical simulation section by showing at which specific region the bifurcations could happen.

In phase portrait Figure 4 (a), it has been shown as chaotic trajectories. But in bifurcation diagram Figure 5, its look like a periodic attractor.

 Response: Thank you for your comment. Figure 5 confirms only the accuracy of the Hopf bifurcation for case 0.017

 seralan vinoth

Competing interests: No competing interests were disclosed.

 6 5 2026

I appreciate the authors for carried out all my comments and suggestions.

 10.5256/f1000research.189832.r464611 Reviewer response for version 1 Shabbir Muhammad Sajjad 1 Referee https://orcid.org/0000-0002-5943-6437 University of Engineering and Technology, Lahore, Pakistan

Competing interests: No competing interests were disclosed.

 13 3 2026 Copyright: © 2026 Shabbir MS 2026 This is an open access peer review report distributed under the terms of the Creative Commons Attribution Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. recommendation approve

The article develops a mathematical model describing the interaction between blooming plants, honeybees, and honey production while considering wind speed as an environmental disturbance. The model is formulated using nonlinear ordinary differential equations with a Beddington–DeAngelis functional response. The authors analyze the system’s dynamics through theoretical methods and numerical simulations, examining properties such as stability and persistence. Bifurcation analysis is also performed to study how parameter changes affect system stability. The results indicate that high wind speeds can destabilize the ecosystem and may lead to the extinction of honeybee populations and reduced honey production.

Overall, the manuscript is well organized, mathematically rigorous, and addresses a significant problem in mathematical ecology. I recommend it for indexing provided the following comments are addressed:

Parameter estimation issues: Many of the parameter values used in the simulations appear to be selected for illustrative purposes rather than derived from measured ecological data, which may affect the realism of the simulation results.

Limited environmental factors: The model only considers wind speed as an environmental disturbance, while other important ecological factors such as temperature variation, pesticide exposure, diseases, and habitat loss are not included.

Simplified ecological assumptions: Some biological processes involved in plant–pollinator interactions are simplified in the model, which may overlook complex ecological behaviors occurring in real ecosystems.

High mathematical complexity: Several parts of the analysis are mathematically intensive and may be difficult for readers from ecological or interdisciplinary backgrounds to fully interpret.

Is the work clearly and accurately presented and does it cite the current literature?

Yes

If applicable, is the statistical analysis and its interpretation appropriate?

Not applicable

Are all the source data underlying the results available to ensure full reproducibility?

No source data required

Is the study design appropriate and is the work technically sound?

Yes

Are the conclusions drawn adequately supported by the results?

Yes

Are sufficient details of methods and analysis provided to allow replication by others?

Yes

Reviewer Expertise:

Mathematical biology, Optimal control of epidemiological models, numerical solution of ode

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.