Processing math: 17%
Research article Special Issues

Forest beetle infestation and its impact on ecosystems: Effects of harvesting practices and fire disruptions

  • Trees play a vital role in the climate by producing oxygen, supporting ecosystems, and benefiting communities and the environment. They are essential for maintaining ecological balance and overall ecosystem health. However, beetle infestations pose a significant threat to forests, causing severe ecological and economic damage, particularly when aggressive species, such as spruce and mountain pine beetles, heavily attack trees. These infestations can lead to widespread tree mortality and long-term environmental consequences. In this study, we highlight the importance of tree plantations and developed a mathematical model to analyze tree-beetle interactions, incorporating the effects of wildfire and harvesting. We examined equilibrium points, their stability, and the basic reproduction number to understand beetle population dynamics. Additionally, we introduced an optimal control strategy to mitigate beetle attacks, incorporating a pesticide variable into the model to reduce beetle populations. The model was numerically solved to generate visual representations, and optimal control strategies were applied to minimize the impact of beetle infestations using pesticides. Sensitivity analysis was conducted to explore factors influencing beetle reproduction, particularly their preference for trees with larger diameters, thicker bark, and extensive phloem, which enhance brood survival and growth. This study underscores the urgency of implementing effective beetle management strategies to protect and restore forest tree populations, ensuring long-term ecosystem sustainability.

    Citation: Sadia Shihab Sinje, Md Kamrujjaman, Rubayyi T. Alqahtani. Forest beetle infestation and its impact on ecosystems: Effects of harvesting practices and fire disruptions[J]. AIMS Mathematics, 2025, 10(4): 9933-9973. doi: 10.3934/math.2025455

    Related Papers:

    [1] Xiaoxiao Liu, Chunrui Zhang . Forest model dynamics analysis and optimal control based on disease and fire interactions. AIMS Mathematics, 2024, 9(2): 3174-3194. doi: 10.3934/math.2024154
    [2] Rongjie Yu, Hengguo Yu, Min Zhao . Steady states and spatiotemporal dynamics of a diffusive predator-prey system with predator harvesting. AIMS Mathematics, 2024, 9(9): 24058-24088. doi: 10.3934/math.20241170
    [3] Mohammed Alsubhi, Rizwan Ahmed, Ibrahim Alraddadi, Faisal Alsharif, Muhammad Imran . Stability and bifurcation analysis of a discrete-time plant-herbivore model with harvesting effect. AIMS Mathematics, 2024, 9(8): 20014-20042. doi: 10.3934/math.2024976
    [4] Kottakkaran Sooppy Nisar, G Ranjith Kumar, K Ramesh . The study on the complex nature of a predator-prey model with fractional-order derivatives incorporating refuge and nonlinear prey harvesting. AIMS Mathematics, 2024, 9(5): 13492-13507. doi: 10.3934/math.2024657
    [5] Xiaoling Zhou, Chao Yang, Weihua He . The linear $ k $-arboricity of symmetric directed trees. AIMS Mathematics, 2022, 7(2): 1603-1614. doi: 10.3934/math.2022093
    [6] Muhammad Bilal Riaz, Nauman Raza, Jan Martinovic, Abu Bakar, Osman Tunç . Modeling and simulations for the mitigation of atmospheric carbon dioxide through forest management programs. AIMS Mathematics, 2024, 9(8): 22712-22742. doi: 10.3934/math.20241107
    [7] Chengchong Lu, Xinxin Liu, Zhicheng Li . The dynamics and harvesting strategies of a predator-prey system with Allee effect on prey. AIMS Mathematics, 2023, 8(12): 28897-28925. doi: 10.3934/math.20231481
    [8] Devendra Kumar, Jogendra Singh, Dumitru Baleanu . Dynamical and computational analysis of a fractional predator-prey model with an infectious disease and harvesting policy. AIMS Mathematics, 2024, 9(12): 36082-36101. doi: 10.3934/math.20241712
    [9] Ahmed Sedky Eldeeb, Muhammad Ahsan-ul-Haq, Mohamed S. Eliwa . A discrete Ramos-Louzada distribution for asymmetric and over-dispersed data with leptokurtic-shaped: Properties and various estimation techniques with inference. AIMS Mathematics, 2022, 7(2): 1726-1741. doi: 10.3934/math.2022099
    [10] Jie Liu, Qinglong Wang, Xuyang Cao, Ting Yu . Bifurcation and optimal harvesting analysis of a discrete-time predator–prey model with fear and prey refuge effects. AIMS Mathematics, 2024, 9(10): 26283-26306. doi: 10.3934/math.20241281
  • Trees play a vital role in the climate by producing oxygen, supporting ecosystems, and benefiting communities and the environment. They are essential for maintaining ecological balance and overall ecosystem health. However, beetle infestations pose a significant threat to forests, causing severe ecological and economic damage, particularly when aggressive species, such as spruce and mountain pine beetles, heavily attack trees. These infestations can lead to widespread tree mortality and long-term environmental consequences. In this study, we highlight the importance of tree plantations and developed a mathematical model to analyze tree-beetle interactions, incorporating the effects of wildfire and harvesting. We examined equilibrium points, their stability, and the basic reproduction number to understand beetle population dynamics. Additionally, we introduced an optimal control strategy to mitigate beetle attacks, incorporating a pesticide variable into the model to reduce beetle populations. The model was numerically solved to generate visual representations, and optimal control strategies were applied to minimize the impact of beetle infestations using pesticides. Sensitivity analysis was conducted to explore factors influencing beetle reproduction, particularly their preference for trees with larger diameters, thicker bark, and extensive phloem, which enhance brood survival and growth. This study underscores the urgency of implementing effective beetle management strategies to protect and restore forest tree populations, ensuring long-term ecosystem sustainability.



    Terms Description
    IPM Integrated Pest Management.
    GDP Gross Domestic Product.
    CIP Country Investment Plan.

    Beetle infestations are primarily caused by certain species of beetles that burrow into tree bark, creating tunnels and disrupting the tree's vascular system. Different beetle species tend to specialize in particular types of trees. For instance, emerald ash borers predominantly attack ash trees, while pine beetles target pine trees. Bark beetles, a family of insects known for infesting multiple tree species, can cause extensive damage in forests. Their larvae feed on the inner bark (phloem) and outer sapwood, interfering with the tree's ability to transport water and nutrients. Without timely intervention, this damage weakens the tree, making it more vulnerable to infections, environmental stressors, and ultimately, tree death. In cases of severe infestations, beetle outbreaks can lead to widespread tree mortality, particularly in forested areas where beetles rapidly spread from tree to tree. Effective management strategies include tree removal, tree traps, insecticide applications, and cultural practices that promote tree health. Early detection and intervention are crucial in minimizing tree damage and controlling beetle populations. To ensure long-term beetle management and forest health, Integrated Pest Management (IPM) strategies combining multiple control methods may be necessary. As noted by Alexander et al. (1977): "Trees hold profound significance for human beings. The importance of ancient trees is deeply rooted in our collective consciousness; they often symbolize the wholeness of personality in our dreams. The trees we cherish create unique spaces—places of solace, passage, and connection. They have the power to shape social environments in diverse ways".

    Variations in climate affect the population dynamics of bark beetles because temperature affects the number of generations they develop, how long they take to reproduce, how many times they die during the winter, how far they can fly, and how well they can migrate [1]. The European spruce bark beetle primarily targets Norway spruce, potentially affecting over 25% of Europe's current Norway spruce growth population [2]. Drought in recent years has greatly increased the amount of insect-induced tree mortality, which has disastrous consequences for the world's biogeochemical cycles, ecosystem function, atmospheric processes, and the availability of sustainable resources. However, the physiological relationships among insect outbreaks, drought, and tree defenses are poorly understood, which makes it challenging to predict tree mortality in the face of continuous climate change [3]. Over the past 50 years, the European spruce bark beetle, the most commercially significant pest in Europe, has ravaged approximately 150 million square meters of forest. A destructive insect pest common in coniferous forests across Europe [4]. The way that this complex of phloem-feeding insects and their microbial partners partitions the belowground resource depends on the anatomy and physiology of the hosts. The red turpentine beetle, Dendroctonus vallens LeConte, reproduces in both living and dead trees [5]. For a year, a single lodgepole pine is usually used as the brood tree by mountain pine beetles in British Columbia [6]. By increasing the potential for acidification, the harvesting of trees in forests has a severe impact on the richness of animals and birds living there as well as the health of the related aquatic ecosystems [7,8]. A significant loss of habitat for a wide variety of trees and animal species can result from forest fires. If a species cannot adapt or find a suitable home elsewhere, it may experience population decreases or perhaps local extinction. This could throw ecological interactions off balance.

    We can learn about the bifurcation analysis of the spruce budworm model and the population model with harvesting in other studies [9,10]. Decisions about forest management are usually influenced by the mountain pine beetle population's size, power, and range, as shown in [11]. Uncontrolled flames that spread quickly through vegetation are known as wildfires, and they typically occur in grasslands, forests, or other wilderness regions. The impacts of wildfires—both positive and negative—on the environment, economy, and society are explored by Gill et al. [12]. In a heterogeneous environment where functions vary spatially over time, the population dynamics of a single species under the influence of harvesting are analyzed by Kamrujjaman et al. [13]. Harvesting activities can be considered as "experimental manipulations" of entire ecosystems, as they often encompass large areas, including water catchments and wildlife habitats [14]. A simplified Lotka-Volterra model has been developed to study the dynamics of tidal forests and wetlands. The Sundarbans, the largest mangrove forest in Bangladesh, faces significant environmental threats and hazardous events, as discussed in [15,16,17]. Additionally, insights into forest ecosystems governed by an age-structured parabolic-ordinary system can be found in other studies [18].

    In this paper, we examine beetle infestations in Bangladesh, along with pine beetle attacks in forests influenced by two key disruptions: harvesting and fire. We develop a mathematical model to analyze tree-beetle interactions while incorporating the effects of harvesting and fire. Additionally, we introduce a pesticide component into the model as a control measure to mitigate beetle populations. Our objects of the study are:

    ● Conducting a study in Bangladesh for beetle attacks in forests.

    ● Describing the negative effects of beetle attack, harvesting and fire on forest.

    ● Analyzing the steps to control beetles.

    ● The use of pesticide to control beetle outbreak.

    The remainder of this paper is structured as follows: In Section 2, we discuss beetle infestations in forests in Bangladesh. In Section 3, we present the development of a mathematical model for beetle attacks on trees, incorporating the effects of harvesting and fire. In Section 4, we analyze the model by determining equilibrium points, their stability, and the properties of the basic reproduction number. In Section 5, we provide numerical solutions and visual representations for both the beetle-free and endemic stages. In Section 6, we examine strategies for controlling beetle populations. In Section 7, we explore optimal control strategies, including the introduction of a pesticide variable into the model to mitigate beetle outbreaks. In Section 8, we conduct a sensitivity analysis of the model with pesticide intervention. Finally, in Section 9, we summarize our findings and conclusions.

    Bangladesh, with 160 million inhabitants and 147,570 km2 of land, is the eighth most populous country in the world. There is a lot of room for the forestry sector to mitigate the effects of climate change. Gross Domestic Product (GDP) derived from the forest sector, which employed 2% of the labor force. This is the main source of the carbon sink. The way that trees lessen the effects of cyclones and storm surges has already been shown. In terms of providing for their needs, it is also crucial for the locals that reside near woodlands. In addition, the country's Gross Domestic Product (GDP) derived from the forest sector, which employed 2% of the labor force. The substantial amounts of fuel wood, fodder, small timber and poles, thatching grass, medicinal plants, and other forest produce that are extracted both legally and illegally are not included in the GDP list, which gives the impression that the forest sector's contribution to the GDP is low. One of the major barriers to the growth of forest species nurseries and plantations is pest (beetle) and disease. When creating the plan to protect forest from beetle, consideration was given to the national act as well as laws, regulations, policies, and guidelines that dealt with issues related to disease, pest, and chemical control. The plan has been prepared by examining a number of national regulatory mechanisms, such as the Forest Act of 1927, the National Forest Policy of 2016, the Forestry Master Plan of 2017–2036, the Forest Investment Plan of 2017, and the Bangladesh Country Investment Plan (CIP) for forests. The resources listed above comprise the following: Environment, Forestry and Climate Change (2016–2021), Seventh Five Year Plan (2016–2020), National Biodiversity Strategy and Action Plan (2016–2021), National Conservationation Strategy, 2017 (Draft), National Agriculture Policy (2013), National Integrated Pest Management Policy (2002), National Crop and Forest Biotechnology Policy (2012), Pesticide Ordinance (1971), amended Pesticide Rules (1985), Plant Quarantine Act (2011), and the Destructive Organizations Act. The plan included a list of the primary diseases and pests that impact nurseries, plantations, and storage facilities, along with suggestions for integrated pest management (IPM) methods to mitigate their impacts [19].

    A system of ordinary differential equations is employed to illustrate the dynamics of the model, which shares similarities with the Modified Predator-Prey Model with Non-linear Incidence [20,21]. This model, known as the "Beetle-Tree Interaction Model, " considers two major disruptions affecting the forest ecosystem: Harvesting and fire. Harvesting refers to the removal of a certain population of trees at specific intervals, making it crucial to analyze its impact on forest dynamics. On the other hand, forest fires, often triggered by extreme heat, pose a significant threat to trees and can drastically alter the ecosystem.

    The mathematical representation of the model is given by:

    {dMdt=rmM(1MkmgkNr+N)ELMPMkmQmM,t0,dNdt=rnN(1Nkn)α(1bEcQm)knN1+γknNQnPMkmQmN,t0, (3.1)

    with initial conditions,

    M(0)=M0,andN(0)=N0. (3.2)

    The state variables in the model represent:

    M: The density of trees.

    N: The density of mountain pine beetles per tree.

    The parameters governing the dynamics of the system are as follows:

    rm: The intrinsic growth rate of trees.

    rn: The intrinsic reproductive rate of beetles.

    r: The critical attack density of beetles.

    km: The carrying capacity of trees.

    kn: The carrying capacity of beetles.

    gk: The fraction of trees killed by beetles.

    α: The defensive rate of trees.

    γ: The reciprocal of the beetle density scale at which tree defense reaches saturation.

    The harvesting-related parameters include:

    E: The proportionate loss of trees due to harvesting.

    L: The average intensity rate of the tree's response to harvesting.

    b: The impact of harvesting on tree defense against beetles.

    The fire-related parameters are:

    Qm: The proportionate loss of trees due to fire.

    P: A constant representing the average intensity rate of the tree's response to fire, which depends on different thermogenic characteristics.

    PMkm: The probability or frequency of fire, assumed to be related to the density of trees in the forest.

    PMkmQm: The fire strength, which determines additional tree mortality and is influenced by both fire intensity and tree density.

    Qn: The proportionate loss of beetles due to fire, correlated with the beetle density per tree.

    c: The impact of fire on tree defense against beetles.

    This formulation accounts for the interplay between tree population dynamics, beetle infestation, harvesting, and fire, incorporating both ecological and external factors that influence forest health and beetle population stability.

    In the absence of harvesting, i.e., when E=0, there is no proportionate loss of trees due to harvesting. Similarly, if there is no fire, meaning the proportionate loss of trees due to fire is also zero (Qm=0), then the additional mortality of trees caused by fire is eliminated. As a result, when both the proportionate loss of trees due to harvesting and fire is zero, the proportionate loss of beetles due to fire also becomes zero (Qn=0).

    Under these conditions, the original model described by Eq (3.1) simplifies to the following form:

    {dMdt=rmM(1MkmgkNr+N),t0,dNdt=rnN(1Nkn)αknN1+γknN,t0. (3.3)

    In addition, when harvesting and fire are present in the system, they are represented by E0 and Qm0. However, if the effects of harvesting and fire on the tree's defense mechanism are absent, i.e., b=0 and c=0, then these factors do not influence the tree's ability to defend against beetle invasion. Table 1 provides a summary of all variables and parameters of the main model.

    Table 1.  Model parameters and their descriptions.
    Notation Unit Definition
    M tree The density of trees
    N beetles/tree The density of mountain pine beetles per tree
    rm (time)1 The intrinsic growth rate of trees
    rn (time)1 The intrinsic reproductive rate of beetles
    r beetles/tree Critical attack density of beetles
    km tree Carrying capacity of trees
    kn beetles/tree Carrying capacity of beetles
    gk unitless Fraction of trees that are killed by beetles
    α tree.(beetles.time)1 Defensive rate of trees
    γ (beetles/tree)2 Reciprocal of the beetle's density scale at which tree's defense
    reaches saturation
    P (time)1 The average intensity rate of the tree's response to the fire
    Qm Unitless Proportionate loss of trees for fire.
    Qn Unitless Proportionate loss of beetles for fire.
    L (time)1 The average intensity rate of the tree's response to the harvesting.
    E Unitless Proportionate loss of trees for harvesting.

     | Show Table
    DownLoad: CSV

    A dimensionless model is derived by rescaling the variables and parameters of a system to eliminate physical units. This transformation simplifies the mathematical analysis and enhances the model's generality by reducing the number of independent parameters. By minimizing the number of governing parameters, dimensionless models make the system more manageable for analysis and numerical simulations. Since the model is independent of specific units, it facilitates comparisons across-different systems and experimental conditions. Furthermore, dimensionless equations often lead to better-conditioned numerical solutions by avoiding extreme values in computations. For clarity and simplification, we introduce the following dimensionless variables and parameters, resulting in:

    m=Mkm,n=Nkn,re=rnrm,qn=QnQmPrm,qm=QmPrm,qs=ELrm,k=gkkn,γ1=γk2n,τ=trm,αn=α(1bEcQm)knrm,(1+qm)<qn,qn>qm.

    Therefore, the dimensionless model is defined as follows:

    {dmdτ=m(1mknr+knnqsqmm),τ0,dndτ=n[re(1n)αn1+γ1nqnm],τ0, (3.4)

    with initial conditions,

    m(0)=m0,andn(0)=n0,where, 0m1,and0n1. (3.5)

    Beetle attacks occur in trees of any forest. To understand this interaction, we develop a model called the "Beetle-Tree Interaction Model", which describes the relationship between beetles and trees in the forest. The forest ecosystem is affected by two major disturbances: Fire and harvesting. Harvesting refers to the removal of multiple trees at once, while forest fires can be triggered by excessive heat, causing significant harm to trees. When both harvesting and fire impact the forest, the model is modified accordingly, leading to the following mathematical representation:

    {dmdτ=m(1mknr+knnqsqmm),τ0,dndτ=n[re(1n)αn1+γ1nqnm],τ0, (4.1)

    with initial conditions,

    m(0)=m0,andn(0)=n0. (4.2)

    For equilibrium solutions,

    dmdτ=0,anddndτ=0,

    such that

    m(1mknr+knnqsqmm)=0. (4.3)
    n[re(1n)αn1+γ1nqnm]=0. (4.4)

    From (4.3) we get,

    m=0,and(1mknr+knnqsqmm)=0. (4.5)

    From (4.4) we have,

    n=0andre+reγ1nrenren2γ1αnqnmqnγ1mn=0. (4.6)

    If n=0, then from (4.5) we get, m=1qs1+qm.

    Thus, we obtain the boundary equilibrium,

    E0=(1qs1+qm,0).

    Let the coexistence equilibrium is E=(m,n), here m>0, n>0, then from (4.5), we get,

    1mknr+knnqsqmm=0m=rqsr+(knqsknk)n(1+qm)(r+knn). (4.7)

    Putting the value of m from (4.7) in (4.6) we get,

    re+reγ1nrenren2γ1αnqnmqnγ1mn=0re+reγ1nrenren2γ1αnqnrqsr+(knqsknk)n(1+qm)(r+knn)qnγ1(rqsr+(knqsknk)n(1+qm)(r+knn))n=0re(1+qm)(r+knn)+reγ1n(1+qm)(r+knn)ren(1+qm)(r+knn)ren2γ1(1+qm)(r+knn)αn(1+qm)(r+knn)qn(rqsr+(knqsknk)n)qnγ1(rqsr+(knqsknk)n)n=0rer+reknn+reqmr+reqmknn+reγ1nr+reγ1knn2+reγ1nqmr+reγ1qmknn2rernreknn2reqmrnreqmknn2reγ1n2rreγ1knn3reγ1qmn2rreγ1qmknn3αnrαnknnαnqmrαnqmknnqnr+qnqsrqn(knk)n+qnqsknnqnγ1nr+qnqsγ1rnqnγ1(knk)n2+qnqsγ1knn2=0[γ1knre(1+qm)]n3+[(kqnrre(1+qm))γ1+kn((1+qm)(1+γ1)re(1qs)qnγ1)]n2+[(1+qm)[knαn+re(kn+r(1+γ1))]+qn(k(1qs)knrγ1(1qs))]n+r[(reαn)(1+qm)qn(1qs)]=0. (4.8)

    We can write (4.8) as,

    F(n)=f1n3+f2n2+f3n+f4=0, (4.9)

    where

    f1=[γ1knre(1+qm)],f2=[(kqnrre(1+qm))γ1+kn((1+qm)(1+γ1)re(1qs)qnγ1)],f3=[(1+qm)[knαn+re(kn+r(1+γ1))]+qn(k(1qs)knrγ1(1qs))],f4=r[(reαn)(1+qm)qn(1qs)].

    Since γ1>0, kn>0, re>0, qm>0, qn>0, 0<qs<1, and f1<0, then limn+F(n)=. If (reαn)(1+qm)>qn(1qs), then F(0)>0. According to the real continuation method, there is n>0 such that F(n)=0. As k=gkkn, 0<gk<1, 0<qs<1, then we have knk>0, (1qs)<1, r>qsr, kn>qskn. Consider, knqsknk>0 and we can write,

    m=rqsr+(knqsknk)n(1+qm)(r+knn)>0.

    Thus, we can conclude that when (reαn)(1+qm)>qn(1qs), the model has at least one coexistence equilibrium E=(m,n). We can consider the basic reproduction number,

    R0=(reαn)(1+qm)qn(1qs). (4.10)

    First, we would describe local asymptotic and global asymptotic stability. The model is,

    {dmdτ=mm2knmr+knnqsmqmm2,τ0,dndτ=renren2αnn1+γ1nqnmn,τ0. (4.11)

    Assume that,

    g1=mm2knmr+knnqsmqmm2,g2=renren2αnn1+γ1nqnmn.

    The Jacobian matrix is,

    J=(g1mg1ng2mg2n),

    which yields

    J=(12(1+qm)mknr+knnqsknkmn(r+knn)2km(r+knn)qnnre(12n)qnm+αnnγ1(1+γ1n)2αn(1+γ1n)).

    For the equilibrium E0=(1qs1+qm,0):

    The Jacobian matrix at point E0 is,

    J0=(1+qsk(1qs)(1+qm)r0reqn(1qs)(1+qm)αn)=(u1u2v1v2). (4.12)

    The eigenvalues of the matrix J0 can be defined as λ0, hence, we have |J0λ0I|, where, I is the identity matrix. The characteristic equation is,

    λ20tr(J0)λ0+det(J0)=0,

    where

    tr(J0)=u1+v2,det(J0)=u1v2u2v1=u1v2,sincev1=0.

    The following theorem establishes the results of the stability criterion for boundary equilibrium.

    Theorem 1. Let qs>0, k>0, qm>0, re>0, r>0, αn>0 such that all of them are positive parameters. Then, we have,

    1) If u1v2>0 and u1+v2<0, then the boundary equilibrium E0 is locally asymptotically stable.

    2) If u1v2<0 or, u1+v2>0, then the boundary equilibrium E0 is unstable.

    Proof. The eigenvalues, which determine stability, are determined by:

    λ0=tr(J0)±(tr(J0))24det(J0)2.

    If the real components of both eigenvalues are negative, the equilibrium is stable. Therefore, the conditions for the stability of the system (4.1) at the equilibrium E0=(1qs1qm,0) are:

    tr(J0)<0anddet(J0)>0. (4.13)

    To ensure stability, the conditions are as follows:

    det(J0)>0u1v2>0, (4.14)

    and

    tr(J0)<0u1+v2<0. (4.15)

    Hence, if u1v2>0 and u1+v2<0, the boundary equilibrium E0 is locally asymptotically stable.

    If the equilibrium E0 is unstable, then at least one of the conditions in (4.13) is not satisfied, such that tr(J0)>0 or det(J0)<0. Now, we have:

    det(J0)<0u1v2<0. (4.16)

    Otherwise, we have:

    tr(J0)>0u1+v2>0. (4.17)

    Consequently, if u1v2<0 or u1+v2>0, then the boundary equilibrium E0 is unstable.

    From (4.14) we have,

    u1>0andv2>0or, u1<0andv2<0.

    If

    u1>0qs>1andv2>0(reαn)(1+qm)>qn(1qs). (4.18)

    The condition (4.18) is not valid as qs1 such that for stable equilibrium, (reαn)(1+qm)qn(1qs). Otherwise we can write,

    u1<0qs<1andv2<0(reαn)(1+qm)<qn(1qs). (4.19)

    The prerequisite (4.19) is possible for the stability as qs<1 and that is why for the stable equilibrium E0, we have (reαn)(1+qm)<qn(1qs). So, for stability, the product of (reαn) and (1+qm) is less than the product of qn and (1qs).

    Following the relation (4.15), we get,

    u1+v2<01+qs+reqn(1qs)(1+qm)αn<0(reαn)(1+qm)(1+qm+qn)(1qs)<1,

    the inequality ensure the stability. Similarly, Eq (4.16) implies,

    u1<0andv2>0or, u1>0andv2<0,

    and, we have qs<1,and(reαn)(1+qm)>qn(1qs). The constraint is possible as qs<1 and for the unstable equilibrium E0, the required condition is,

    (reαn)(1+qm)>qn(1qs).

    The next mathematical inequality prospect, u1>0 and v2<0, that yields qs>1, and (reαn)(1+qm)<qn(1qs), respectively.

    The Eq (4.17) implies

    u1+v2>0(reαn)(1+qm)(1+qm+qn)(1qs)>1. (4.20)

    From (4.20), we have for the instability (reαn)(1+qm)(1+qm+qn)(1qs)>1. So, the ratio of (reαn)(1+qm) and (1+qm+qn)(1qs) is greater than 1.

    Theorem 2. Let qs>0, k>0, qm>0, re>0, r>0, αn>0. Then we have,

    1) If (reαn)(1+qm)<qn(1qs), then the boundary equilibrium E0 is locally asymptotically stable.

    2) If (reαn)(1+qm)>qn(1qs), then the equilibrium E0 is unstable.

    Proof. Consider the basic reproduction number,

    R0=(reαn)(1+qm)qn(1qs). (4.21)

    For stable equilibrium,

    R0<1(reαn)(1+qm)qn(1qs)<1(reαn)(1+qm)<qn(1qs). (4.22)

    For unstable equilibrium,

    R0>1(reαn)(1+qm)qn(1qs)>1(reαn)(1+qm)>qn(1qs). (4.23)

    Thus, if (reαn)(1+qm)<qn(1qs), then the boundary equilibrium E0 is locally asymptotically stable. If (reαn)(1+qm)>qn(1qs), then the boundary equilibrium E0 is unstable.

    Coexistence equilibrium E=(m,n):

    The Jacobian martix at E is,

    J=(12(1+qm)mknr+knnqsknkmn(r+knn)2km(r+knn)qnnre(12n)qnm+αnnγ1(1+γ1n)2αn(1+γ1n)). (4.24)

    Consider,

    J=(w1w2z1z2). (4.25)

    Here,

    w1=12(1+qm)mknr+knnqs,w2=knkmn(r+knn)2km(r+knn),z1=qnn,z2=re(12n)qnm+αnnγ1(1+γ1n)2αn(1+γ1n).

    The eigenvalues of the matrix J can be defined as λ. The characteristic equation is,

    λ2tr(J)λ+det(J)=0, (4.26)

    while tr(J)=w1+z2,det(J)=w1z2w2z1. For stability analysis, following Theorem 1, we have tr(J)<0 and det(J)>0.

    The equilibrium is stable as long as, tr(J)<0w1<z2,anddet(J)>0w1>w2z1z2. It is a remark for the stable situation, w1 is less than z2 and w1 is greater than the ratio of w2z1 and z2. To check the unstable equilibrium point E, at least one of the conditions of (4.13) is not satisfied, such that tr(J)>0 or, det(J)<0 and we obtain the following relations.

    w1>z2,w1<w2z1z2.

    Hence, we can conclude that for instability, w1 is greater than z2 or, w1 is less than the ratio of w2z1 and z2.

    Theorem 3. Let any open region Ω containing equilibrium E of system (4.1) and R0>1. If αn<re(n+1γ1), then equilibrium E is globally asympotically stable.

    Proof. Assume any open region,

    Ω={(m,n):m>0,n>0}.

    Consider the scalar function, D=1mn. Next from the model (4.1), we have,

    {DY1mn[m(1mknr+knnqsqmm)],DZ1mn[n[re(1n)αn1+γ1nqnm]].

    Taking partial derivatives,

    (DY)m=1n(qm1).(DZ)n=1m(αnγ1(1+γ1n)2re)<1m(αnγ1(1+γ1n)re). (4.27)

    Check that,

    (DY)m<0, where n belongs to Ω.

    From (4.27), consider,

    αnγ1(1+γ1n)<reαn<re(n+1γ1). (4.28)

    Using (4.28), we remark that

    (DZ)n<0, where m belongs to Ω.
    Therefore, (DY)m+(DZ)n<0, where (m,n) belongs to Ω. (4.29)

    According to Dulac Criterion, there is no limit cycle or periodic orbit in the open region Ω. Here,

    (reαn)(1+qm)>qn(1qs),

    such that R0>1. Hence the equilibrium E is globally asymptotically stable.

    Recall the basic reproduction number as defined in (4.21),

    R0=(reαn)(1+qm)qn(1qs). (4.30)

    Now,

    (reαn)(1+qm)=re+reqmαnαnqm=re+reqmα(1bEcQm)knrmα(1bEcQm)knqmrm=re+reqmαknrm+αbEknrm+αcQmknrmαknqmrm+αbEknqmrm+αcQmknqmrm=re+reqmαknrm+αbqsrmknLrm+αcqmrmknPrmαknqmrm+αbqsrmknqmLrm+αcqmrmknqmPrm=(reαknrm+αbqsknL)+(reαknrm+αbqsknL+αcknP)qm+αcknq2mP, (4.31)

    where, qm=QmPrmandqs=ELrm.

    Consider,

    R0=1(reαknrm+αbqsknL)+(reαknrm+αbqsknL+αcknP)qm+αcknq2mPqn(1qs)=1qn(1qs)=(reαknrm+αbqsknL)+(reαknrm+αbqsknL+αcknP)qm+αcknq2mPqn=(reαknrm+αbqsknL)+(reαknrm+αbqsknL+αcknP)qm+αcknq2mP(1qs). (4.32)

    In Figure 1, we can observe the properties of R_{0} . This figure corresponds to Eq (4.32) and represents the relationship between q_{n} and q_{m} . The parameters used in this figure are:

    \begin{align*} r_{e} & = 33.75, & k_{n} & = 1956, & \alpha & = 0.04086, \\ r_{m} & = 0.08, & P & = 0.1, & c & = 0.5. \end{align*}

    The data utilized in this study is obtained from previous studies [20,21,22], which includes Canadian data from the Canadian Forest Service. We consider 0 < q_{s} < 1 and assume certain parameter values such as b = 0.5 , q_{s} = 0.8 , and L = 0.1 to achieve meaningful results.

    Figure 1.  The properties of basic reproduction number R_{0} , where ( r_{e} = 33.75 , k_{n} = 1956 , \alpha = 0.04086 , r_{m} = 0.08 , P = 0.1 , c = 0.5 , b = 0.5 , q_{s} = 0.8 L = 0.1 ).

    Figure 1 illustrates the scenario when R_{0} = 1 . When the basic reproduction number R_{0} is equal to one, each infected individual will, on average, infect exactly one other individual. This typically results in a stable epidemic, where the total number of affected individuals remains constant over time. In epidemiology, R_{0} is a crucial metric for determining a disease's contagiousness and potential for spread. If R_{0} > 1 , it indicates a higher likelihood of exponential disease transmission, potentially leading to an epidemic or outbreak. Conversely, if R_{0} < 1 , the disease is expected to eventually die out since each infected individual infects, on average, fewer than one other person. Furthermore, if R_{0} < 1 , the graph tends to approach the vertical axis, whereas if R_{0} > 1 , it shifts closer to the horizontal axis.

    Every forest has beetle attacks on its trees. To explain the relationship between the forest's beetles and trees, we would create a model. The forest can also be harmed by wildfire. Human activity and natural elements can both start forest fires. In woods, lightning strikes have the potential to start fires, particularly in arid weather. Nearby trees may catch fire due to the heat and ash from volcanic activity. Sometimes, it is necessary to set fire to trees that have been impacted by beetles. When there is fire in forests, the model is,

    \begin{align} \begin{cases} \frac{ dm}{ d\tau} = m\left(1-m-\frac{kn}{r+k_{n}n}-q_{m}m\right), \; \;\; \tau\geq0, \\ \frac{ dn}{ d\tau} = n\left[r_{e}\left(1-n\right)-\frac{\alpha_{n}}{1+\gamma_{1}n}- q_{n}m\right], \; \; \; \tau\geq0, \end{cases} \end{align} (4.33)

    with initial conditions,

    \begin{align} & m(0) = m_{0}, \;\; \text{and}\;\; n(0) = n_0. \end{align} (4.34)

    Following the same procedure as before, we can determine the boundary equilibrium as:

    \begin{equation} E_{0} = \left(\frac{1}{1+q_{m}}, 0 \right), \end{equation} (4.35)

    and, the mathematical equations to establish the next results.

    \begin{align} &r_{e}\left(1-n\right)-\frac{\alpha_{n}}{1+\gamma_{1}n}- q_{n}m = 0 \end{align} (4.36)
    \begin{align} \Rightarrow \;& r_{e}-r_{e}n-\frac{\alpha_{n}}{1+\gamma_{1}n}- q_{n}m = 0 \\ \Rightarrow \;& r_{e}+r_{e}\gamma_{1}n-r_{e}n-r_{e}n^{2}\gamma_{1}-\alpha_{n}-q_{n}m-q_{n}\gamma_{1}mn = 0. \end{align} (4.37)

    Let the coexistence equilibrium is, E_{*} = (m^{*}, n^{*}) , where m^{*} > 0 , n^{*} > 0 . From (4.33), we obtain,

    \begin{align} &1-m-\frac{kn}{r+k_{n}n}- q_{m}m = 0 \\ \Rightarrow \;& \left(r+k_{n}n\right)-kn-m\left(1+q_{m}\right)\left(r+k_{n}n\right) = 0 \\ \therefore\;& m = \frac{r+\left(k_{n}-k\right)n} {\left(1+q_{m}\right)\left(r+k_{n}n\right)}. \end{align} (4.38)

    Putting the value of m from (4.38) in (4.37), we get,

    \begin{align} &r_{e}+r_{e}\gamma_{1}n-r_{e}n-r_{e}n^{2}\gamma_{1}-\alpha_{n}-q_{n}m-q_{n}\gamma_{1}mn = 0 \\ \Rightarrow \;& r_{e}+r_{e}\gamma_{1}n-r_{e}n-r_{e}n^{2}\gamma_{1}-\alpha_{n}-q_{n}\frac{r+\left(k_{n}-k\right)n} {\left(1+q_{m}\right)\left(r+k_{n}n\right)}\\ & -q_{n}\gamma_{1}\left(\frac{r+\left(k_{n}-k\right)n} {\left(1+q_{m}\right)\left(r+k_{n}n\right)}\right)n = 0 \\ \Rightarrow \;& r_{e}\left(1+q_{m}\right)\left(r+k_{n}n\right)+r_{e}\gamma_{1}n\left(1+q_{m}\right)\left(r+k_{n}n\right) \\ &-r_{e}n\left(1+q_{m}\right)\left(r+k_{n}n\right) -r_{e}n^{2}\gamma_{1}\left(1+q_{m}\right)\left(r+k_{n}n\right) -\alpha_{n}\left(1+q_{m}\right)\left(r+k_{n}n\right)\\ &-q_{n}\left(r+\left(k_{n}-k\right)n\right) -q_{n}\gamma_{1}\left(r+\left(k_{n}-k\right)n\right)n = 0. \\ \therefore\;& -\left[\gamma_{1}k_{n}r_{e}\left(1+q_{m}\right)\right]n^{3} \\ &+\left[\left(kq_{n}-rr_{e}\left(1+q_{m}\right)\right)\gamma_{1}+k_{n}\left(\left(1+q_{m}\right)\left(-1+\gamma_{1}\right)r_{e}-q_{n}\gamma_{1}\right)\right]n^{2} \\ &+\left[\left(1+q_{m}\right)\left[-k_{n}\alpha_{n}+r_{e}\left(k_{n}+r\left(-1+\gamma_{1}\right)\right)\right]+q_{n}\left(k-k_{n}-r\gamma_{1}\right)\right]n \\ &+r\left[\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right)-q_{n}\right] = 0. \end{align} (4.39)

    We can write (4.39) as,

    \begin{equation} F(n) = f_{1}n^{3}+f_{2}n^{2}+f_{3}n+f_{4} = 0, \end{equation} (4.40)

    where

    \begin{align*} &f_{1} = -\left[\gamma_{1}k_{n}r_{e}\left(1+q_{m}\right)\right], \\ &f_{2} = \left[\left(kq_{n}-rr_{e}\left(1+q_{m}\right)\right)\gamma_{1}+k_{n}\left(\left(1+q_{m}\right)\left(-1+\gamma_{1}\right)r_{e}-q_{n}\gamma_{1}\right)\right], \\ &f_{3} = \left[\left(1+q_{m}\right)\left[-k_{n}\alpha_{n}+r_{e}\left(k_{n}+r\left(-1+\gamma_{1}\right)\right)\right]+q_{n}\left(k-k_{n}-r\gamma_{1}\right)\right], \\ &f_{4} = r\left[\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right)-q_{n}\right]. \end{align*}

    Here, \gamma_{1} > 0 , k_{n} > 0 , r_{e} > 0 , q_{m} > 0 , q_{n} > 0 , and f_{1} < 0 , then \lim\limits_{n \rightarrow \; +\infty}F(n) = -\infty . If \left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n} , then F(0) > 0 and according to the real continuation method, there is n^{*} > 0, such that F(n^{*}) = 0 . Since k = g_{k} k_{n} and 0 < g_{k} < 1 , it follows that:

    \begin{equation*} k_{n} - k > 0. \end{equation*}

    Thus, we can express m^{*} as:

    \begin{equation} m^{*} = \frac{r + \left(k_{n} - k\right)n^{*}}{\left(1 + q_{m}\right)\left(r + k_{n} n^{*}\right)} > 0. \end{equation} (4.41)

    Consequently, we conclude that if the condition:

    \begin{equation*} \left(r_{e} - \alpha_{n}\right)\left(1 + q_{m}\right) > q_{n}, \end{equation*}

    is satisfied, then the model (4.33) possesses at least one coexistence equilibrium given by:

    \begin{equation*} E_{*} = \left(m^{*}, n^{*}\right). \end{equation*}

    Define the Jocobian matrix for the system (4.33), we obtain,

    \begin{equation} J = \begin{pmatrix} 1-2\left(1+q_{m}\right)m-\frac{kn}{r+k_{n}n}& \frac{k_{n}kmn}{\left(r+k_{n}n\right)^{2}}-\frac{km}{\left(r+k_{n}n\right)} \\ -q_{n}n & r_{e}(1-2n)-q_{n}m+\frac{\alpha_{n}n\gamma_{1}}{\left(1+\gamma_{1}n\right)^{2}}-\frac{\alpha_{n}}{\left(1+\gamma_{1}n\right)}\\ \end{pmatrix}. \end{equation} (4.42)

    For the equilibrium E_{0} = \left(\frac{1}{1+q_{m}}, 0\right) :

    The Jacobian matrix at point E_{0} is,

    \begin{align} & J_{0} = \begin{pmatrix} -1 & -\frac{k}{\left(1+q_{m}\right)r}\\ 0 & r_{e}- \frac{q_{n}}{\left(1+q_{m}\right)}-\alpha_{n} \end{pmatrix} = \begin{pmatrix} e_{1} & e_{2}\\ h_{1} & h_{2}\\ \end{pmatrix}, \end{align} (4.43)

    such that

    \begin{align*} e_{1} = -1, \;\;\; e_{2} = -\frac{k}{\left(1+q_{m}\right)r}, \;\;\; h_{1} = 0, \;\;\; h_{2} = r_{e}- \frac{q_{n}}{\left(1+q_{m}\right)}-\alpha_{n}. \end{align*}

    The characteristic equation is,

    \begin{equation} \lambda_{0}^{2}-\text{tr}(J_{0})\lambda_0+\text{det}(J_{0}) = 0, \nonumber \end{equation}

    where,

    \begin{align} &\text{tr}(J_{0}) = e_{1}+h_{2}, \\ &\text{det}(J_{0}) = e_{1}h_{2}-e_{2}h_{1} = e_{1}h_{2}, \; \; \; \text{since} \;\; h_{1} = 0. \end{align}

    To ensure the stability of the system at equilibrium E_{0} = \left(\frac{1}{1 - q_{m}}, 0\right) , the determinant and trace conditions must be satisfied. The determinant condition is:

    \begin{equation} \text{det}(J_{0}) > 0 \Rightarrow e_{1} h_{2} > 0.\nonumber \end{equation}

    This implies:

    \begin{equation} e_{1} > 0, \quad h_{2} > 0, \quad \text{or} \quad e_{1} < 0, \quad h_{2} < 0. \nonumber \end{equation}

    For e_{1} > 0 and h_{2} > 0 , we obtain:

    \begin{equation} \left(r_{e} - \alpha_{n}\right)\left(1 + q_{m}\right) > q_{n}, \nonumber \end{equation}

    which is not possible since -1 \ngtr 0 . Hence, the stability condition requires:

    \begin{equation} \left(r_{e} - \alpha_{n}\right)\left(1 + q_{m}\right) < q_{n}. \nonumber \end{equation}

    Additionally, the trace condition is given by:

    \begin{equation} \text{tr}(J_{0}) < 0 \Rightarrow e_{1} + h_{2} < 0.\nonumber \end{equation}

    This leads to:

    \begin{equation} \frac{\left(r_{e} - \alpha_{n}\right)\left(1 + q_{m}\right)}{\left(1 + q_{m} + q_{n}\right)} < 1.\nonumber \end{equation}

    Thus, for stability, the ratio of \left(r_{e} - \alpha_{n}\right)\left(1 + q_{m}\right) to \left(1 + q_{m} + q_{n}\right) must be less than 1.

    If equilibrium E_{0} is unstable, then at least one of the conditions in (4.13) is not satisfied, meaning either \text{tr}(J_0) > 0 or \text{det}(J_0) < 0 . The determinant condition for instability is:

    \begin{equation} \text{det}(J_{0}) < 0 \Rightarrow e_{1} h_{2} < 0.\nonumber \end{equation}

    This implies:

    \begin{equation} e_{1} < 0, \quad h_{2} > 0, \quad \text{or} \quad e_{1} > 0, \quad h_{2} < 0. \nonumber \end{equation}

    For e_{1} < 0 and h_{2} > 0 , we obtain

    \begin{equation} \left(r_{e} - \alpha_{n}\right)\left(1 + q_{m}\right) > q_{n}, \nonumber \end{equation}

    which confirms instability. Hence, for instability,

    \begin{equation} \left(r_{e} - \alpha_{n}\right)\left(1 + q_{m}\right) > q_{n}.\nonumber \end{equation}

    Additionally, if

    \begin{equation} \text{tr}(J_{0}) > 0 \Rightarrow e_{1} + h_{2} > 0, \nonumber \end{equation}

    we obtain:

    \begin{equation} \frac{\left(r_{e} - \alpha_{n}\right)\left(1 + q_{m}\right)}{\left(1 + q_{m} + q_{n}\right)} > 1. \nonumber \end{equation}

    Thus, for instability, this ratio must be greater than 1.

    For the equilibrium E_{*} = \left(m^{*}, n^{*}\right) :

    The Jacobian martix at E_{*} is,

    \begin{equation} J_{*} = \begin{pmatrix} 1-2\left(1+q_{m}\right)m^{*}-\frac{kn^{*}}{r+k_{n}n^{*}} & \frac{k_{n}km^{*}n^{*}}{\left(r+k_{n}n^{*}\right)^{2}}-\frac{km^{*}}{\left(r+k_{n}n^{*}\right)} \\ -q_{n}n^{*} & r_{e}(1-2n^{*})-q_{n}m^{*}+\frac{\alpha_{n}n^{*}\gamma_{1}}{\left(1+\gamma_{1}n^{*}\right)^{2}}-\frac{\alpha_{n}}{\left(1+\gamma_{1}n^{*}\right)} \end{pmatrix}. \end{equation} (4.44)

    Assume,

    \begin{equation} J_{*} = \begin{pmatrix} i_{1} & i_{2}\\ j_{1} & j_{2}\\ \end{pmatrix}, \end{equation} (4.45)

    where,

    \begin{align*} &i_{1} = 1-2\left(1+q_{m}\right)m^{*}-\frac{kn^{*}}{r+k_{n}n^{*}}, \;\;\; i_{2} = \frac{k_{n}km^{*}n^{*}}{\left(r+k_{n}n^{*}\right)^{2}}-\frac{km^{*}}{\left(r+k_{n}n^{*}\right)}, \\ &j_{1} = -q_{n}n^{*}, \;\;\; j_{2} = r_{e}(1-2n^{*})-q_{n}m^{*}+\frac{\alpha_{n}n^{*}\gamma_{1}}{\left(1+\gamma_{1}n^{*}\right)^{2}}-\frac{\alpha_{n}}{\left(1+\gamma_{1}n^{*}\right)}. \end{align*}

    The eigenvalues of matrix J_{*} , denoted as \lambda_{*} , satisfy the characteristic equation:

    \begin{equation} \lambda_{*}^{2} - \text{tr}(J_{*}) \lambda_{*} + \text{det}(J_{*}) = 0.\nonumber \end{equation}

    Here,

    \begin{align} &\text{tr}(J_{*}) = i_{1} + j_{2}, \\ &\text{det}(J_{*}) = i_{1}j_{2} - i_{2}j_{1}. \end{align}

    For stability,

    \begin{equation} i_{1} < -j_{2}, \quad i_{1} > \frac{i_{2}j_{1}}{j_{2}}. \nonumber \end{equation}

    For instability,

    \begin{equation} i_{1} > -j_{2}, \quad i_{1} < \frac{i_{2}j_{1}}{j_{2}}.\nonumber \end{equation}

    Trees in any forest are susceptible to beetle attacks. To understand the interaction between beetles and trees, we develop a model that describes their dynamics. Additionally, trees can be affected by harvesting, which may occur for both commercial and non-commercial purposes. Logging or timber extraction is carried out for various reasons, including the removal of beetle-infested trees to control their spread. When tree harvesting is considered, the Beetle-Tree Interaction Model is given by:

    \begin{align} \begin{cases} \frac{ dm}{ d\tau} = m\left(1-m-\frac{kn}{r+k_{n}n}-q_{s}\right), \; \;\; \tau\geq0, \\ \frac{ dn}{ d\tau} = n\left[r_{e}\left(1-n\right)-\frac{\alpha_{n}}{1+\gamma_{1}n}\right], \; \; \; \tau\geq0, \end{cases} \end{align} (4.46)

    with initial conditions,

    \begin{align} & m(0) = m_{0}, \;\; \text{and}\;\; n(0) = n_0. \end{align} (4.47)

    To determine equilibrium, we set

    \frac{dm}{d\tau} = 0, \quad \frac{dn}{d\tau} = 0.

    This leads to system

    \begin{align} &m\left(1 - m - \frac{kn}{r + k_{n} n} - q_{s}\right) = 0, \\ &n\left[r_{e} (1 - n) - \frac{\alpha_{n}}{1 + \gamma_{1} n}\right] = 0. \end{align}

    Solving these equations, we obtain two equilibria:

    ● The boundary equilibrium,

    \begin{equation} E_{0} = \left(1 - q_{s}, 0\right).\nonumber \end{equation}

    ● The coexistence equilibrium,

    \begin{equation} E_{*} = \left(m^{*}, n^{*}\right), \quad m^{*} > 0, \quad n^{*} > 0.\nonumber \end{equation}

    For coexistence equilibrium, solving for m^{*} is

    \begin{equation} m^{*} = \frac{r - q_{s}r + \left(k_{n} - q_{s} k_{n} - k\right) n^{*}}{r + k_{n} n^{*}} > 0.\nonumber \end{equation}

    From the quadratic equation,

    \begin{equation} F(n) = f_{1} n^{2} + f_{2} n + f_{3} = 0, \nonumber \end{equation}

    where

    \begin{align} &f_{1} = -r_{e} \gamma_{1}, \quad f_{2} = r_{e} (\gamma_{1} - 1), \quad f_{3} = r_{e} - \alpha_{n}. \end{align}

    Given r_{e} > \alpha_{n} , there is a positive root n^{*} > 0 ensuring coexistence. Therefore, when r_{e} > \alpha_{n} , the model admits at least one coexistence equilibrium E_{*} = (m^{*}, n^{*}) .

    Define the Jacobian matrix of the model (4.46),

    \begin{equation} J = \begin{pmatrix} 1-2m-\frac{kn}{r+k_{n}n}-q_{s} & \frac{k_{n}kmn}{\left(r+k_{n}n\right)^{2}}-\frac{km}{\left(r+k_{n}n\right)} \\ 0& r_{e}(1-2n)+\frac{\alpha_{n}n\gamma_{1}}{\left(1+\gamma_{1}n\right)^{2}}-\frac{\alpha_{n}}{\left(1+\gamma_{1}n\right)}\\ \end{pmatrix}. \end{equation} (4.48)

    For the equilibrium E_{0} = \left(1-q_{s}, 0\right) : The Jacobian matrix at point E_{0} is,

    \begin{align} & J_{0} = \begin{pmatrix} -1+q_{s} & -k\frac{\left(1-q_{s}\right)}{r}\\ 0 & r_{e}-\alpha_{n} \end{pmatrix} = \begin{pmatrix} a_{1} & a_{2}\\ d_{1} & d_{2}\\ \end{pmatrix}, \end{align} (4.49)

    where

    \begin{align*} &a_{1} = -1+q_{s}, \;\; a_{2} = -k\frac{\left(1-q_{s}\right)}{r}, \;\; d_{1} = 0, \;\; d_{2} = r_{e}- \alpha_{n}. \end{align*}

    The eigenvalues of the matrix J_{0} can be defined as \lambda_{0} . The characteristic equation is,

    \begin{equation} \lambda_{0}^{2}-\text{tr}(J_{0})\lambda_{0}+\text{det}(J_{0}) = 0. \nonumber \end{equation}

    Here,

    \begin{align*} &\text{tr}(J_{0}) = a_{1}+d_{2},\\ &\text{det}(J_{0}) = a_{1}d_{2}-a_{2}d_{1} = a_{1}d_{2}, \; \; \; \text{since} \;\; d_{1} = 0. \end{align*}

    For stability, the conditions are:

    \begin{align} &\det(J_{0}) > 0 \quad \Rightarrow \quad a_{1} d_{2} > 0 \quad \Rightarrow \quad (a_{1} > 0, d_{2} > 0) \text{ or } (a_{1} < 0, d_{2} < 0). \end{align}

    Since q_{s} < 1 , the stable equilibrium condition is r_{e} < \alpha_{n} . Additionally,

    \begin{equation} \frac{(r_{e} - \alpha_{n})}{(1 - q_{s})} < 1, \nonumber \end{equation}

    ensuring stability when the ratio of (r_{e} - \alpha_{n}) to (1 - q_{s}) is less than one.

    If E_{0} is unstable, then:

    \begin{align} &\det(J_{0}) < 0 \quad \Rightarrow \quad a_{1} d_{2} < 0 \quad \Rightarrow \quad (a_{1} < 0, d_{2} > 0) \text{ or } (a_{1} > 0, d_{2} < 0). \end{align}

    For instability, we find that r_{e} > \alpha_{n} is a necessary condition, implying:

    \begin{equation} \frac{(r_{e} - \alpha_{n})}{(1 - q_{s})} > 1. \nonumber \end{equation}

    Thus, instability occurs when the ratio of (r_{e} - \alpha_{n}) to (1 - q_{s}) exceeds one.

    For the equilibrium E_{*} = \left(m^{*}, n^{*}\right) : The Jacobian martix at E_{*} is,

    \begin{equation} J_{*} = \begin{pmatrix} 1-2m^{*}-\frac{kn^{*}}{r+k_{n}n^{*}}-q_{s} & \frac{k_{n}km^{*}n^{*}}{\left(r+k_{n}n^{*}\right)^{2}}-\frac{km^{*}}{\left(r+k_{n}n^{*}\right)} \\ 0 & r_{e}(1-2n^{*})+\frac{\alpha_{n}n^{*}\gamma_{1}}{\left(1+\gamma_{1}n^{*}\right)^{2}}-\frac{\alpha_{n}}{\left(1+\gamma_{1}n^{*}\right)} \end{pmatrix}. \end{equation} (4.50)

    Consider,

    \begin{equation} J_{*} = \begin{pmatrix} x_{1} & x_{2}\\ y_{1} & y_{2}\\ \end{pmatrix}. \end{equation} (4.51)

    Here,

    \begin{align*} &x_{1} = 1-2m^{*}-\frac{kn^{*}}{r+k_{n}n^{*}}-q_{s}, \;\; x_{2} = \frac{k_{n}km^{*}n^{*}}{\left(r+k_{n}n^{*}\right)^{2}}-\frac{km^{*}}{\left(r+k_{n}n^{*}\right)}, \\ &y_{1} = 0, \;\; y_{2} = r_{e}(1-2n^{*})+\frac{\alpha_{n}n^{*}\gamma_{1}}{\left(1+\gamma_{1}n^{*}\right)^{2}}-\frac{\alpha_{n}}{\left(1+\gamma_{1}n^{*}\right)}. \end{align*}

    The eigenvalues of the matrix J_{*} can be defined as \lambda_{*} . The characteristic equation is,

    \begin{equation} \lambda_{*}^{2}-\text{tr}(J_{*})\lambda_{*}+\text{det}(J_{*}) = 0. \nonumber \end{equation}

    Now, for stable equilibrium,

    \begin{align} &\text{tr}(J_{*}) < 0 \quad \text{and} \quad \text{det}(J_{*}) > 0 \\ \Rightarrow \;& x_{1}+y_{2} < 0 \quad \Rightarrow \; x_{1}y_{2}-x_{2}y_{1} > 0 \\ \therefore\;& x_{1} < -y_{2}\quad \therefore\; x_{1} > \frac{x_{2}y_{1}}{y_{2}}. \end{align} (4.52)

    From (4.52), for stability, x_{1} must be less than -y_{2} and greater than the ratio \frac{x_{2}y_{1}}{y_{2}} .

    For unstable equilibrium,

    \begin{align} &\text{tr}(J_{*}) > 0 \quad \text{or} \quad \text{det}(J_{*}) < 0 \\ \Rightarrow \;& x_{1}+y_{2} > 0 \quad \Rightarrow \; x_{1}y_{2}-x_{2}y_{1} < 0 \\ \therefore\;& x_{1} > -y_{2}\quad \therefore\; x_{1} < \frac{x_{2}y_{1}}{y_{2}}. \end{align} (4.53)

    From (4.53), for instability, x_{1} must be greater than -y_{2} or less than the ratio \frac{x_{2}y_{1}}{y_{2}} .

    Our findings indicate that the number of trees generally decreases due to harvesting and fire. Beetles experience two distinct consequences from these disruptions. The disruptions directly reduce the beetle population by causing mortality. However, they weaken the trees' natural defenses, creating a more favorable environment for beetle proliferation. Due to these opposing effects, we aim to numerically analyze two different equilibrium points to understand the overall impact. We consider the following:

    ● Beetle free equilibrium (boundary equilibrium E_{0} ).

    ● Endemic equilibrium (coexistence equilibrium E_{*} ).

    To solve the model, we need different values of the parameters of the model (Table 2).

    Table 2.  Model parameters values.
    Notation Value Source
    \gamma_{1} 156.9398 [20]
    r_{e} 33.75 [20]
    q_{m} [75-100] [20]
    q_{n} 366.75 [20]
    k_{n} 1956 [20]
    \alpha 0.04086 [20]
    r_{m} 0.08 [20]
    P 0.1 [21]
    Q_{m} [0.3-0.6] [21]
    \alpha_{n} [30-849.17295] Estimated [21]
    k 1467 [21]
    r 9.1 [21]
    q_{s} 0.8 Assumed

     | Show Table
    DownLoad: CSV

    Using the relations, \left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) < q_{n}\left(1-q_{s}\right) and \left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n}\left(1-q_{s}\right) , we found the value of \alpha_{n} . We collect data from [20,21,22], which are Canadian data from Canadian Forest Service and consider 0 < q_{s} < 1 .

    In this case, we have solved model (4.1) numerically using different values. The Beetle free equilibrium is E_{0} = \left(\frac{1-q_{s}}{1+q_{m}}, 0\right) . We know that if,

    \begin{align} & \left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) < q_{n}\left(1-q_{s}\right), \end{align} (5.1)

    then the equilibrium E_{0} is locally asymptotically stable. For Eq (5.1), we consider \alpha_{n} = 529.48431 in the absence of a beetle attack in the forest. Thus, the initial condition is n_0 = 0 . The system experiences two disruptions: Fire and harvesting. At the equilibrium point E_{0} , we have n = 0 , where n = \frac{N}{k_{n}} . This implies that the beetle density, N , is zero.

    From Figure 2, it is evident that beetles have no effect on the forest dynamics. Instead, the primary factors influencing the system are harvesting and fire, both of which contribute to a decrease in m . As the parameter q_m increases, the death rate of trees due to fire also increases. Similarly, as q_s increases, the death rate of trees due to harvesting rises. Consequently, an increase in either q_m or q_s leads to a reduction in the value of m . Furthermore, since m is defined as m = \frac{M}{k_{m}} , a decrease in m directly implies a reduction in the tree density M . Thus, higher values of q_m and q_s result in a lower overall tree population.

    Figure 2.  (a) m vs \tau for different values of q_m = 20, 40, 75 , q_s = 0.8 , (b) m vs \tau for different values of q_s = 0.20, 0.40, 0.60, 0.80 , q_m = 75 , where m = \frac{M}{k_{m}} , \tau = tr_{m} , t is in decades, m_{0} = 0.05 , n_{0} = 0 and for the values of \gamma_{1} , r_{e} , r , q_{n} , k , k_{n} the reference table is 2.

    The endemic equilibrium point is E_{*} = \left(m^{*}, n^{*}\right) . We establish that, under the following criteria,

    \begin{align} &\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n}\left(1-q_{s}\right), \end{align} (5.2)

    the equilibrium E_{0} is unstable. When the equilibrium E_{0} is unstable, an endemic equilibrium E_{*} emerges.

    In this scenario, we consider Eq (5.2) with \alpha_{n} = 30 , indicating the presence of a beetle attack in the forest. Additionally, the system is affected by two external disruptions: Fire and harvesting. From Figure 3, we observe that at the endemic equilibrium E_{*} :

    ● In subplot (a), the value of n increases, leading to a rise in the beetle density N . Notably, n converges almost to its maximum value of 1, indicating a significant beetle population.

    ● In subplot (b), due to the combined effects of beetles, harvesting, and fire, the value of m decreases more rapidly compared to the beetle-free scenario in Figure 2. This implies a sharper decline in the density of trees, M , surpassing the reduction observed in Figure 2.

    Thus, the presence of beetles, along with harvesting and fire disturbances, leads to a substantial increase in beetle density and a more pronounced decline in tree density compared to the beetle-free case.

    Figure 3.  (a) n vs \tau , (b) m vs \tau , (c) Trajectories of m and n vs \tau , where n = \frac{N}{k_{n}} , m = \frac{M}{k_{m}} , \tau = tr_{m} , t is in decades, n_{0} = 0.1 , m_{0} = 0.1 , \alpha_{n} = 30 , q_m = 75 , q_s = 0.8 , k = 1467, and the values of \gamma_{1} , r_{e} , r , q_{n} , k_{n} are referred from Table 2.

    From Figure 4, we observe a direct relationship between the parameter k and the number of trees that are destroyed. As the value of k increases, the number of trees that are killed also increases, indicating a stronger impact of the parameter on tree mortality. Furthermore, an increase in k leads to a decrease in m , where m represents the normalized tree density. Since m is given by m = \frac{M}{k_m} , a reduction in m implies a corresponding decline in the total tree density, M . This suggests that as k grows, the overall forest population diminishes significantly. Thus, a higher value of k accelerates tree loss, leading to a more pronounced decline in tree density over time.

    Figure 4.  m vs \tau at m_{0} = 0.01 , n_{0} = 0.1 , k = 200,500, 1200 , where m = \frac{M}{k_{m}} , n = \frac{N}{k_{n}} , \tau = tr_{m} , t is in decades, \alpha_{n} = 30 , q_m = 75 , q_s = 0.8, and the values of \gamma_{1} , r_{e} , r , q_{n} , k_{n} are available in Table 2.

    From Figure 5, we observe the effect of parameter \gamma_1 on the beetle population dynamics. For \gamma_1 = 15 , the curve representing n initially decreases before gradually increasing. When \gamma_1 = 20 , the curve exhibits a slower initial decline compared to \gamma_1 = 15 , but then it increases at a higher rate than for \gamma_1 = 15 . As the value of \gamma_1 continues to rise, the curve of n increases more significantly. This behavior suggests that when \gamma_1 is relatively small, trees exhibit a strong resistance to beetle attacks, leading to an initial suppression of the beetle population. At this stage, the trees' ability to resist beetles is at its peak. However, as \gamma_1 increases further, the trees lose their ability to counteract the beetle infestation effectively. Consequently, the beetle population grows unchecked, leading to a decline in the tree population. Thus, a higher value of \gamma_1 eventually results in a lower tree density and an increased beetle population.

    Figure 5.  n vs \tau at m_{0} = 0.1 , n_{0} = 0.1 , \gamma_1 = 15, 20, 50,150 , k = 1467 , where m = \frac{M}{k_{m}} , n = \frac{N}{k_{n}} , \tau = tr_{m} , t is in decades, \alpha_{n} = 30 , q_m = 75 , q_s = 0.8 and the values of r_{e} , r , q_{n} , k_{n} are available in Table 2.

    In this section, we focus on the computational results of Eq (4.33), specifically considering the impact of fire disruptions on the system. The beetle-free equilibrium is given by

    \begin{align*} E_{0} = \left(\frac{1}{1+q_{m}}, 0\right), \end{align*}

    which is asymptotically stable under the condition:

    \begin{align} \left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) < q_{n}. \end{align} (5.3)

    For the computations in Eq (5.3), we consider the parameter value \alpha_{n} = 849.17295 , representing a scenario where there is no beetle attack in the forest. Consequently, the initial beetle population is zero, i.e., n_0 = 0 . Since at equilibrium point E_0 we have n = 0 , the density of beetles, N , is also zero, as defined by the relation n = \frac{N}{k_{n}} . In this scenario, fire is the only disruption affecting the forest ecosystem. From Figure 6, we observe that in the absence of beetles, fire remains the primary factor influencing tree mortality. As the fire intensity parameter q_m increases, the rate of tree death due to fire also rises. This leads to a decrease in m , which represents the normalized tree density. Since m is defined as m = \frac{M}{k_m} , a reduction in m implies a corresponding decline in the total tree density, M . Thus, our analysis reveals that in a beetle-free environment, fire alone can significantly impact forest density. As the severity of fire disturbances increases, tree density declines at an accelerated rate, contributing to overall forest degradation.

    Figure 6.  m vs \tau at q_m = 20, 40, 75,100 , where m = \frac{M}{k_{m}} , n = \frac{N}{k_{n}} , \tau = tr_{m} , t is in decades, m_{0} = 0.1 , n_{0} = 0 (as there is no beetle attack), and the values of \gamma_{1} , r_{e} , r , q_{n} , k , k_{n} are counted from the Table 2.

    The endemic equilibrium of the system is denoted as

    \begin{align*} E_{*} = \left(m^{*}, n^{*}\right). \end{align*}

    We know that if the following condition holds:

    \begin{align} \left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n}, \end{align} (5.4)

    then the beetle-free equilibrium E_{0} becomes unstable, leading to the existence of an endemic equilibrium E_{*} where both tree and beetle populations coexist.

    For the computational analysis in Eq (5.4), we consider \alpha_{n} = 30 , representing a scenario where beetle infestation is present in the forest, along with fire as a disruptive factor. From Figure 7, we observe that at the endemic equilibrium point E_{*} , the beetle population density, n , increases over time, leading to a rise in the total beetle population, N .

    Figure 7.  (a) n vs \tau (b) m vs \tau , (c) Trajectories of m and n vs \tau , where n = \frac{N}{k_{n}} , m = \frac{M}{k_{m}} , \tau = tr_{m} , t is in decades, n_{0} = 0.1 , m_{0} = 0.1 , \alpha_{n} = 30 , q_m = 100 , k = 1467, and the values of \gamma_{1} , r_{e} , r , q_{n} , k_{n} are used from Table 2.

    The combined effects of beetle infestation and fire disturbances cause a more rapid decline in tree density compared to the beetle-free scenario in Figure 6. Consequently, the normalized tree density, m , decreases at a faster rate than in the absence of beetles. Since tree mortality is primarily driven by beetle attacks and fire in this case, we note a more significant reduction in the total tree population, M , than in Figure 6.

    Unlike previous scenarios, harvesting does not influence tree density in this analysis. As a result, the trees' natural defensive mechanisms against beetles strengthen, leading to a relatively slower increase in the beetle population compared to Figure 3. Similarly, the rate of tree decline is also slower than that observed in Figure 3, where harvesting played a role. Thus, our findings indicate that while beetle infestations and fire contribute significantly to tree loss, the absence of harvesting allows trees to develop a stronger resistance, partially mitigating the severity of beetle population growth and tree mortality.

    From Figure 8, it is evident that as the parameter k increases, the number of trees that are killed also rises. This results in a decrease in the normalized tree density, m . Since M represents the total tree population and is related to m through m = \frac{M}{k_m} , the decline in m implies that the total tree density, M , also decreases correspondingly. Thus, a higher value of k leads to increased tree mortality, reducing the overall forest density.

    Figure 8.  m vs \tau at m_{0} = 0.03 , n_{0} = 0.1 , k = 200,500, 1200 , where m = \frac{M}{k_{m}} , n = \frac{N}{k_{n}} , \tau = tr_{m} , t is in decades, \alpha_{n} = 30 , q_m = 100, and the values of \gamma_{1} , r_{e} , r , q_{n} , k_{n} from the Table 2.

    In Figure 9, we analyze the effect of the resistance parameter \gamma_1 on beetle population dynamics. When \gamma_1 = 10 , the curve representing n (the normalized beetle density) initially declines before gradually increasing. For \gamma_1 = 20 , we observe a similar trend; however, the initial decrease occurs at a slower rate compared to \gamma_1 = 10 , followed by a more pronounced increase. As the value of \gamma_1 continues to rise, the curve of n exhibits a steady upward trend.

    Figure 9.  n vs \tau at m_{0} = 0.1 , n_{0} = 0.1 , \gamma_1 = 10, \; 20, \; 50, \; 150 , k = 1467 , where m = \frac{M}{k_{m}} , n = \frac{N}{k_{n}} , \tau = tr_{m} , t is in decades, \alpha_{n} = 30 , q_m = 100, and the values of r_{e} , r , q_{n} , k_{n} are used from the Table 2.

    This behavior suggests that when \gamma_1 is small, trees have a strong ability to resist beetle infestations, effectively limiting beetle population growth. At this stage, the trees' resistance to beetles is at its peak. However, as \gamma_1 increases further, trees gradually lose their ability to counteract the beetle infestation effectively. As a result, the beetle population expands unchecked, leading to a decline in tree density. Our findings indicate that while lower values of \gamma_1 enhance the trees' resistance against beetles, higher values weaken this defense mechanism, allowing the beetle population to grow while reducing the tree population.

    In this case, we numerically solve model (4.46) using different parameter values to analyze the system's behavior. The beetle-free equilibrium is given by

    E_{0} = \left(1-q_{s}, 0\right).

    We establish that if

    \begin{align} r_{e} < \alpha_{n}, \end{align} (5.5)

    then the equilibrium E_{0} is locally asymptotically stable.

    To explore this scenario, we consider \alpha_{n} = 679.33836 and assume that there is no beetle attack in the forest. This implies an initial beetle population of n_0 = 0 . In this setting, the primary disruption affecting the forest is harvesting. From Figure 10, we observe that harvesting is the sole factor influencing tree density, as beetles have no impact on the system. As the harvesting parameter q_s increases, the tree mortality rate due to harvesting also rises. Consequently, for increasing values of q_s , the normalized tree density m increases at a slow rate, leading to a gradual increase in the total tree density M .

    Figure 10.  m vs \tau at q_s = 0.20, \, 0.40, \; 0.60, \; 0.80 , where m = \frac{M}{k_{m}} , \tau = tr_{m} , t is in decades, m_{0} = 0.1 , n_{0} = 0 (as there is no beetle attack), and the values of \gamma_{1} , r_{e} , r , k , k_{n} are counted from the Table 2.

    The endemic equilibrium of the system is given by

    E_{*} = \left(m^{*}, n^{*}\right).

    We know that if

    \begin{align} r_{e} > \alpha_{n}, \end{align} (5.6)

    then the beetle-free equilibrium E_{0} becomes unstable, and the system admits an endemic equilibrium E_{*} , where both tree and beetle populations coexist.

    For this scenario, we set \alpha_{n} = 30 , enabling beetle infestation to occur in the forest alongside harvesting disruption. From Figure 11, we observe that at the endemic equilibrium point E_{*} , the beetle density, N , increases over time. The beetle population exhibits oscillations, where beetle numbers decrease when trees are harvested but rise again when new trees grow from the harvested stumps. This dynamic interplay leads to fluctuations in beetle density.

    Figure 11.  (a) n vs \tau , (b) n vs \tau in close range, (c) m vs \tau , where n = \frac{N}{k_{n}} , m = \frac{M}{k_{m}} , \tau = tr_{m} , t is in decades, n_{0} = 0.9 , m_{0} = 0.1 , \alpha_{n} = 30 , q_s = 0.8 and the values of \gamma_{1} , r_{e} , r , k , k_{n} are available in Table 2.

    The combined effects of beetle infestation and harvesting lead to a decline in tree density, as observed in Figure 11(b). Since fire is absent in this scenario, its impact on tree mortality is negligible. However, we note that the defensive capacity of trees against beetles is lower in this case compared to Figure 7. As a result, the beetle population grows more rapidly than in Figure 7, while the tree population declines at a faster rate. Moreover, the system takes a longer time \tau to stabilize compared to Figure 3, where both fire and harvesting were present. These results indicate that in the presence of harvesting alone, tree density initially increases at a slow rate but eventually declines as beetles become more prevalent. The competition between beetle growth and tree harvesting leads to oscillatory dynamics in beetle population levels, influencing the overall stability of the forest ecosystem.

    From Figure 12, we observe the following trends:

    ● In subplot (a), as the value of k increases, the number of trees that are killed also increases. Consequently, for higher values of k , the curve representing m exhibits a decreasing trend.

    ● In subplot (b), when the parameter k_n increases, the decline in m occurs at a slower rate, indicating a more gradual decrease in tree density.

    Figure 12.  (a) m vs \tau at m_{0} = 0.01 , n_{0} = 0.1 , k = 450,550, 1250 , k_n = 1956 , (b) m vs \tau at m_{0} = 0.03 , n_{0} = 0.1 , k_n = 2000, 4000, 6000 , k = 1467 , where m = \frac{M}{k_{m}} , n = \frac{N}{k_{n}} , \tau = tr_{m} , t is in decades, \alpha_{n} = 30 , q_s = 0.8 and the values of \gamma_{1} , r_{e} , r , are in Table 2.

    There are many reasons why our trees sustain damage, including an influx of beetles into trees, trees being chopped down for a variety of reasons, and forest fires. As a result, forest's tree population is declining. If a tree is attacked by beetles, the damage gets worse in addition to what is done by fire and harvesting.

    In forests, managing beetles is essential to preserving the ecosystem's health and equilibrium. Significant tree damage by beetles can result in degraded forests and a decline in biodiversity. Here are a few methods that are frequently used in forests to manage beetles:

    Forested regions should be regularly surveyed and monitored in order to identify beetle populations early on. We should identify the beetle species that are present and evaluate their distribution and population dynamics. To precisely map areas of infection, we must use ground surveys. We should use an integrated strategy that incorporates many control techniques, including mechanical, chemical, cultural, and biological control. To reduce our impact on the environment, whenever possible, we have to give preference to non-chemical methods.

    Beetle parasites or native predators can be introduced to assist manage beetle populations. To prevent unforeseen effects on species other than the intended target, this strategy needs to be carefully considered.

    We must use forest management techniques, such as thinning, pruning, and encouraging a diversity of tree species and age classes, to improve tree resilience and decrease insect habitat. In order to reduce the population of adult beetles, we must utilize trapping tools. To stop beetles from infecting healthy trees, diseased trees should be removed. When preventing beetles from accessing susceptible trees, we have to think about erecting physical barriers or applying tree banding methods. When non-chemical approaches fail to manage beetle populations, we can use chemical pesticides as a last resort. It is imperative that we select pesticides that are specifically targeted at that specific species of beetle and have minimal impact on non-target creatures and the environment. All safety precautions and laws pertaining must be observed to the use of pesticides, including the right timing and application techniques. To lessen the chance of beetle spread, we need to make everyone working in the forest aware of how important it is to follow quarantine regulations. We could encourage studies aimed at exploring the biology, behavior, and available management strategies of beetles. It is important to work together with governmental organizations, academic institutions, and other interested parties to exchange information and resources for successful beetle control. Ii is important to inform the public, managers, and landowners in forests about the value of controlling beetles and their part in early detection and prevention. To enable early intervention, we should give instructions on how to recognize beetle symptoms and report possible infestations. In forest ecosystems, we must keep evaluating the feasibility and long-term impacts of beetle management programs. Using these strategies as part of a comprehensive beetle control plan can help forest managers effectively mitigate the consequences of infestations while protecting the resilience and overall health of forest ecosystems.

    In optimal control problems, consider m(t) and n(t) is state variable and w(t) is control variable. The state variables satisfy the ordinary differential equations. The state variable depends on control variable such that [23,24],

    \begin{align} & m'(t) = g\left(t, w(t), m(t), n(t)\right), \end{align} (7.1)
    \begin{align} & n'(t) = s\left(t, w(t), m(t), n(t)\right). \end{align} (7.2)

    Consider the objective function,

    \begin{align} & \max\limits_{w}\int_{t_{0}}^{t_{1}} f\left(t, w(t), m(t), n(t)\right), \;\text{subject to, } \end{align} (7.3)
    \begin{align} & m'(t) = g\left(t, w(t), m(t), n(t)\right) , \end{align} (7.4)
    \begin{align} & n'(t) = s\left(t, w(t), m(t), n(t)\right) , \\ & m(t_{0}) = m_{0} \;\; \text{and}\;\; n(t_{0}) = n_{0}. \end{align} (7.5)

    The Hamiltonian function is,

    \begin{align} & H\left(t, w(t), m(t), n(t), \lambda_{1}(t), \lambda_{2}(t)\right)\\ & = \; f\left(t, w(t), m(t), n(t)\right)+\lambda_{1}(t)g\left(t, w(t), m(t), n(t)\right)+\lambda_{2}(t)s\left(t, w(t), m(t), n(t)\right). \end{align} (7.6)

    We want to maximize the Hamiltonian function H with respect to w at w_{*} , where w_{*} is the optimal value of w . According to Pontryagin's Maximum Principle, the conditions are,

    \begin{align} & \frac{\partial H}{\partial w} = 0\; \text{at}\; w_{*}\; \\ \Rightarrow \; &f_{w}+\lambda_{1}g_{w}+\lambda_{2}s_{w} = 0, \end{align} (7.7)

    with

    \begin{align} &\begin{cases} \lambda_{1}' = -\frac{\partial H}{\partial m} = -\left(f_{m}+\lambda_{1}g_{m}+\lambda_{2}s_{m}\right)\; \text{(adjoint equation)}, \\ \lambda_{2}' = -\frac{\partial H}{\partial n} = -\left(f_{n}+\lambda_{1}g_{n}+\lambda_{2}s_{n}\right)\; \text{(adjoint equation)}, \\ \lambda_{1}(t_{1}) = 0 \; \; \text{and}\;\;\lambda_{2}(t_{1}) = 0 \; \text{(transversality condition)}. \end{cases} \end{align} (7.8)

    Theorem 4. [23] Consider that f\left(t, w(t), m(t), n(t)\right) , g\left(t, w(t), m(t), n(t)\right) and s\left(t, w(t), m(t), n(t)\right) are continuously differentiable functions in four variables and concave in w. Let w_{*} is an optimal control solution for problem, with associated state variables m_{*} and n_{*} , and \lambda_{1} and \lambda_{2} are piecewise differentiable function with \lambda_{1}(t)\geq 0 and \lambda_{2}(t)\geq 0 . Consider t_{0}\leq t \leq t_{1} and for optimal control variable w_{*} , we have,

    \begin{align} &H\left(t, w(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right) \leq H\left(t, w_{*}(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right), \;\mathit{\text{where, }} \end{align} (7.9)
    \begin{align} &H_{w}\left(t, w_{*}(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right) = 0. \end{align} (7.10)

    Proof. For time t_{0}\leq t \leq t_{1} and optimal control variable w_* , we can write,

    \begin{align*} &H\left(t, w_{*}(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right)-H\left(t, w(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right)\\ & = \left[f\left(t, w_{*}(t), m_{*}(t), n_{*}(t)\right)+\lambda_{1}(t)g\left(t, w_{*}(t), m_{*}(t), n_{*}(t)\right)+\lambda_{2}(t)s\left(t, w_{*}(t), m_{*}(t), n_{*}(t)\right)\right] \\ &\;\;- \left[f\left(t, w(t), m_{*}(t), n_{*}(t)\right)+\lambda_{1}(t)g\left(t, w(t), m_{*}(t), n_{*}(t)\right)+\lambda_{2}(t)s\left(t, w(t), m_{*}(t), n_{*}(t)\right)\right]\\ & = \left[f\left(t, w_{*}(t), m_{*}(t), n_{*}(t)\right)-f\left(t, w(t), m_{*}(t), n_{*}(t)\right)\right]\\ &\;\;+\lambda_{1}(t)\left[g\left(t, w_{*}(t), m_{*}(t), n_{*}(t)\right)-g\left(t, w(t), m_{*}(t), n_{*}(t)\right)\right]\\ &\;\;+\lambda_{2}(t)\left[s\left(t, w_{*}(t), m_{*}(t), n_{*}(t)\right)-s\left(t, w(t), m_{*}(t), n_{*}(t)\right)\right]\\ &\geq \left(w_{*}(t)-w(t)\right)f_{w}\left(t, w_{*}(t), m_{*}(t), n_{*}(t)\right) +\lambda_{1}(t)\left(w_{*}(t)-w(t)\right)g_{w}\left(t, w_{*}(t), m_{*}(t), n_{*}(t)\right)\\ &\;\; +\lambda_{2}(t)\left(w_{*}(t)-w(t)\right)s_{w}\left(t, w_{*}(t), m_{*}(t), n_{*}(t)\right)\\ & = \left(w_{*}(t)-w(t)\right)H_{w}\left(t, w_{*}(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right) = 0. \end{align*}

    Hence, we conclude that,

    \begin{align*} H\left(t, w(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right) \leq H\left(t, w_{*}(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right). \end{align*}

    For the minimization problem,

    \begin{align*} H\left(t, w(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right) \geq H\left(t, w_{*}(t), m_{*}(t), n_{*}(t), \lambda_{1}(t), \lambda_{2}(t)\right). \end{align*}

    The problem is maximization problem as long as

    \begin{equation*} \frac{\partial^{2} H}{\partial w^{2}} < 0 \;\; \text{at}\;\; w_{*}, \end{equation*}

    whereas the problem is minimization one, if

    \begin{equation*} \frac{\partial^{2} H}{\partial w^{2}} > 0 \;\; \text{at}\;\; w_{*}. \end{equation*}

    Our objective is to employ pesticides to reduce the number of beetles. The modified model with pesticide is,

    \begin{align} \begin{cases} \frac{ dM}{ dt} = r_{m}M\left(1-\frac{M}{k_{m}}-g_{k}\frac{N}{r+N}\right)-ELM-P\frac{M}{k_{m}}Q_{m}M+U\frac{W}{k_{w}}M\frac{N}{k_n} , \; \;\; t\geq 0, \\ \frac{ dN}{ dt} = r_{n}N\left(1-\frac{N}{k_{n}}\right)-\frac{\alpha\left(1-bE-cQ_{m}\right)k_nN}{1+\gamma k_n N}- Q_{n}P\frac{M}{k_{m}}Q_{m}N-G\frac{W}{k_{w}}N, \; \; \; t\geq 0, \end{cases} \end{align} (7.11)

    with initial conditions,

    \begin{align} M(0) = M_{0}, \;\; \text{and}\;\; N(0) = N_0. \end{align} (7.12)

    The detailed description of the model, along with its variables and parameters, can be found in the main model (3.1). The additional parameters are provided in the following Table 3.

    Table 3.  Model parameters and their descriptions.
    Notation Unit Definition
    M tree The density of trees
    N beetles/tree The density of mountain pine beetles per tree
    r_{m} \text{(time)}^{-1} The intrinsic growth rate of trees
    r_{n} \text{(time)}^{-1} The intrinsic reproductive rate of beetles
    r beetles/tree Critical attack density of beetles
    k_{m} tree Carrying capacity of trees
    k_{n} beetles/tree Carrying capacity of beetles
    g_{k} unitless Fraction of trees that are killed by beetles
    \alpha \text{tree(beetles.time)}^{-1} Defensive rate of trees
    \gamma \text{(beetles/tree)}^{-2} Reciprocal of the beetle's density scale at which tree's defense
    reaches saturation
    P \text{(time)}^{-1} The average intensity rate of tree's response to the fire
    Q_{m} Unitless Proportionate loss of trees for fire.
    Q_{n} Unitless Proportionate loss of beetles for fire.
    L \text{(time)}^{-1} The average intensity rate of tree's response to the harvesting.
    E Unitless Proportionate loss of trees for harvesting.
    U \text{(time)}^{-1} Growth rate of trees due to beetle reduction with the use of pesticide
    W ml Amount of application of pesticide
    k_{w} ml Maximum amount of pesticide can be applied.
    G \text{(time)}^{-1} Death rate of beetles for using pesticide

     | Show Table
    DownLoad: CSV

    By introducing the following non-dimensional variables and parameters,

    \begin{align*} & m = \frac{ M}{ k_{m}}, \; n = \frac{ N}{ k_{n}}, \; r_{e} = \frac{ r_{n}}{ r_{m}}, \; k = g_{k}k_{n}, \; \gamma_{1} = \gamma k_{n}^{2}, \\ & q_{m} = \frac{ Q_{m}P}{ r_{m}}, \; \tau = tr_{m}, \; q_{n} = \frac{ Q_{n}Q_{m}P}{ r_{m}}, \; q_{s} = \frac{ EL}{ r_{m}}, \;\\ & \alpha_{n} = \frac{ \alpha \left(1-bE-cQ_{m}\right)k_{n}}{ r_m}, \; h = \frac{ U}{ r_{m}}, \; e = \frac{ G}{ r_{m}}, \; w = \frac{ W}{ k_{w}}, \; (1+q_{m}) < q_{n}, \; q_{n} > q_{m}. \end{align*}

    We derive the modified dimensionless system of nonlinear differential equations, which is defined as follows:

    \begin{align} \begin{cases} \frac{ dm}{ d\tau} = m\left(1-m-\frac{kn}{r+k_{n}n}-q_{s}-q_{m}m\right)+hwmn , \; \;\; \tau\geq 0, \\ \frac{ dn}{ d\tau} = n\left[r_{e}\left(1-n\right)-\frac{\alpha_{n}}{1+\gamma_{1}n}- q_{n}m\right]-ewn, \; \; \; \tau\geq 0, \end{cases} \end{align} (7.13)

    with initial conditions,

    \begin{align} m(0) = m_{0}, \;\; \text{and}\;\; n(0) = n_0. \end{align} (7.14)

    Our objective is to increase the number of trees by controlling the beetle population through the application of pesticides. Therefore, we aim to maximize the use of pesticides within a given limit.

    Consider that there exists an upper bound \overline{w} such that 0 \leq w(\tau) \leq \overline{w} over a fixed time interval [0, T] , where T > 0 . The control function w is time-dependent, denoted as w = w(\tau) . The constraint on w(\tau) is given by:

    \begin{equation} \varphi = \{w(\tau) \mid 0\leq w(\tau) \leq \overline{w}, \; 0\leq \tau \leq T, \;w(\tau) \text{ is Lebesgue measurable} \}.\nonumber \end{equation}

    The optimal objective function is formulated as:

    \begin{align} A(w(\tau)) = \int_{0}^{T} \left[m(\tau)w(\tau)-m(\tau)^{2}-w(\tau)^{2}\right]d\tau, \quad \text{where} \quad m(0) = m_{0}, \quad n(0) = n_{0}. \end{align} (7.15)

    The optimal control problem is then defined as:

    \begin{equation} A_{*}(w_{*}(\tau)) = \max\limits_{w(\tau) \in \varphi} A(w(\tau)).\nonumber \end{equation}

    The integrand of the objective function is given by:

    \begin{align} V(w(\tau)) = m(\tau)w(\tau) - m(\tau)^{2} - w(\tau)^{2}. \end{align}

    To obtain maximum value of (7.16), the Hamiltonian function is,

    \begin{align} &H = mw-m^{2}-w^{2}, \\ &\;\;+\lambda_{1}\left[m-m^{2}-\frac{knm}{r+k_{n}n}-q_{s}m-q_{m}m^{2}+hwmn\right] +\lambda_{2}\left[r_{e}n-r_{e}n^{2}-\frac{\alpha_{n}n}{1+\gamma_{1}n}- q_{n}mn-ewn\right]. \end{align}

    To find opotimal solutions, the Pontryagin's Maximum Principle is,

    \begin{align} \begin{cases} &\frac{\partial H}{\partial w} = 0, \\ &\lambda'_{1} = -\frac{\partial H}{\partial m}, \;\lambda'_{2} = -\frac{\partial H}{\partial n}, \\ &\lambda_{1}(T) = 0, \; \lambda_{2}(T) = 0. \end{cases} \end{align} (7.16)

    Using Pontryagin's Maximum Principle, we have,

    \begin{align} &\frac{\partial H}{\partial w} = 0 \Rightarrow \; w_{*} = \frac{m+\lambda_{1}hmn-\lambda_{2}en}{2} \;\; \text{at}\;\; w_{*}. \end{align}

    Here,

    \begin{align} \frac{\partial^{2} H}{\partial w^{2}} = -2 < 0. \end{align}

    Thus, it is a maximization problem. We can write,

    \begin{align} &\lambda'_{1} = -\frac{\partial H}{\partial m} = -\left[w-2m+\lambda_{1}\left[1-2m-\frac{kn}{r+k_{n}n}-q_{s}-2q_{m}m+hwn\right]-\lambda_{2}q_{n}n\right], \\ &\lambda'_{2} = -\frac{\partial H}{\partial n} = -\left[\lambda_{1}Y_1+\lambda_{2}Y_2\right], \\ &\lambda_{1}(T) = 0, \; \lambda_{2}(T) = 0, \; \end{align}

    where, Y_1 = \left[\frac{k_{n}kmn}{\left(r+k_{n}n\right)^{2}}-\frac{km}{\left(r+k_{n}n\right)}+hwm\right], \; Y_2 = \left[r_{e}-2r_{e}n+\frac{\alpha_{n}n\gamma_{1}}{\left(1+\gamma_{1}n\right)^{2}}-\frac{\alpha_{n}}{\left(1+\gamma_{1}n\right)}-q_{n}m-ew\right].

    It is now time to demonstrate the computational results of the model (7.13). Using equation

    \begin{equation} \left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n}\left(1-q_{s}\right), \nonumber \end{equation}

    we determine the value of \alpha_{n} . Additionally, using the equation

    \begin{equation} w = \frac{W}{k_{w}}, \nonumber \end{equation}

    we compute the values of w for different values of W . Given W = 3.50, 4.33, 9.40, 13.95, 19.86 , we obtain the corresponding values of w as follows:

    \begin{equation} w = 0.175, \quad 0.2165, \quad 0.47, \quad 0.6975, \quad 0.993.\nonumber \end{equation}

    For plotting the figures, we consider the values of w as:

    \begin{equation} w = 1.75, \quad 2.165, \quad 4.7, \quad 6.975, \quad 9.93.\nonumber \end{equation}

    We collect the values of parameters \gamma_{1} , r_{e} , r , q_{m} , q_{n} , k , and k_{n} from Table 2, and the values of w , h , and e from Table 4. The data were sourced from [20,21,22,25], which are Canadian datasets obtained from the Canadian Forest Service.

    Table 4.  Model parameters values.
    Notation Value Source
    W [3.50-19.86] [25]
    w [1.75-9.93] Estimated [25]
    h 0.4 Assumed
    e 0.6 Assumed
    k_{w} 20 Assumed

     | Show Table
    DownLoad: CSV

    We assume the conditions:

    \begin{equation} 0 < q_{s} < 1, \quad 0 < h < 1, \quad 0 < e < 1, \nonumber \end{equation}

    and assign specific values to achieve accurate results. Here,

    \begin{equation} m = \frac{M}{k_{m}}, \quad n = \frac{N}{k_{n}}, \quad 0 \leq M \leq k_{m}, \quad 0 \leq N \leq k_{n}.\nonumber \end{equation}

    Thus, it follows that:

    \begin{equation} 0 \leq m \leq 1, \quad 0 \leq n \leq 1.\nonumber \end{equation}

    From Figure 13, it can be inferred that using pesticide to kill beetles effectively prevents their population from increasing in the forest. As the value of w gradually increases, the value of n rises slowly, which in turn causes the beetle density N to increase at a slower rate. Therefore, after the pesticide is applied, the beetle population tends to decrease. Consequently, the density of trees also decreases gradually as the beetles' death rate increases.

    Figure 13.  (a) n vs \tau (b) m vs \tau , where n = \frac{N}{k_{n}} , m = \frac{M}{k_{m}} , \tau = tr_{m} , t is in decades, n_{0} = 0.1 , m_{0} = 0.1 , \alpha_{n} = 30 , and Tables 2 and 4 are referred for other parameters values.

    In this section, we compute the sensitivities of the model (7.13). Let us consider the system,

    \begin{align} \begin{cases} & \dot{x_{1}} = f_{1}\left(x_{1}, x_{2}, \cdots, x_{q}\right), \\ & \dot{x_{2}} = f_{2}\left(x_{1}, x_{2}, \cdots, x_{q}\right), \\ &\;\vdots \quad \vdots\\ & \dot{x_{q}} = f_{q}\left(x_{1}, x_{2}, \cdots, x_{q}\right). \end{cases} \end{align} (8.1)

    We use local forward sensitivity analysis for the ordinary differential equations, and the sensitivity equation is,

    \begin{align} & \dot{u_{i}} = \frac{\partial f}{\partial x}u_{i}+\frac{\partial f}{\partial p}, \;\; \text{where}\;\; i = m, n, \; \; \text{and}\;\; p = \text{parameters}, \end{align} (8.2)
    \begin{align} &\dot{u_{i}} = J.u_{i}+F, \;\; \text{where}\;\; u_{i} = u_{i}(\tau), \;\; \dot{u_{i}} = \dot{u_{i}}(\tau) = \frac{ du_{i}}{d\tau}, \;\; J = \frac{\partial f}{\partial x}, \;\; \text{and} \;\; F = \frac{\partial f}{\partial p}, \end{align} (8.3)

    with initial condition,

    \begin{align*} & u_{i}(0) = u_{i0}, \;\; \text{where}\;\; i = m, n. \end{align*}

    The Jacobian of the system (7.13) is,

    \begin{equation} J = \begin{pmatrix} 1-2\left(1+q_{m}\right)m-\frac{kn}{r+k_{n}n}-q_{s}+hwn & \frac{k_{n}kmn}{\left(r+k_{n}n\right)^{2}}-\frac{km}{\left(r+k_{n}n\right)}+hwm\\ -q_{n}n& r_{e}(1-2n)-q_{n}m+\frac{\alpha_{n}n\gamma_{1}}{\left(1+\gamma_{1}n\right)^{2}}-\frac{\alpha_{n}}{\left(1+\gamma_{1}n\right)}-ew\\ \end{pmatrix}. \end{equation} (8.4)

    For pesticide parameter w , we have

    \begin{align} &F = \begin{pmatrix} hmn\\ -en \end{pmatrix}, \;\;\; \;\;\; u_{i} = \begin{pmatrix} u_{m}\\ u_{n} \end{pmatrix}. \end{align} (8.5)

    Substituting the values from (8.4) and (5.1) into (8.3), we get,

    \begin{align*} &\begin{pmatrix} \dot{u_{m}}\\ \dot{u_{n}} \end{pmatrix} = \begin{pmatrix} A_{1} & A_{2} \\ A_{3} & A_{4} \end{pmatrix}\begin{pmatrix} u_{m}\\ u_{n} \end{pmatrix} + \begin{pmatrix} hmn\\ -en \end{pmatrix}, \end{align*}

    where

    \begin{align*} &A_{1} = 1-2\left(1+q_{m}\right)m-\frac{kn}{r+k_{n}n}-q_{s}+hwn, \\ &A_{2} = \frac{k_{n}kmn}{\left(r+k_{n}n\right)^{2}}-\frac{km}{\left(r+k_{n}n\right)}+hwm, \\ &A_{3} = -q_{n}n, \\ &A_{4} = r_{e}(1-2n)-q_{n}m+\frac{\alpha_{n}n\gamma_{1}}{\left(1+\gamma_{1}n\right)^{2}}-\frac{\alpha_{n}}{\left(1+\gamma_{1}n\right)}-ew, \end{align*}
    \begin{align} \therefore\; &\begin{pmatrix} \dot{u_{m}}\\ \dot{u_{n}} \end{pmatrix} = \begin{pmatrix} A_{1}u_{m}+ A_{2}u_{n}+hmn \\ A_{3}u_{m}+ A_{4}u_{n}-en \end{pmatrix}. \end{align}

    It is time to present the results visually. We determine the values of m = 0.0638 and n = 0.646 by solving the model described in equation (7.13). Additionally, we obtain the values of \gamma_1 , r_e , r , q_m , q_n , k , and k_n from Table 2, and the values of w , h , and e from Table 4. Further data is sourced from [20,21,22,25], which provides Canadian data from the Canadian Forest Service. For the purpose of achieving reliable results, we assume certain values as necessary during the analysis.

    From Figure 14, we observe that the amount of u_n decreases more significantly than u_m when pesticides are used.

    Figure 14.  Trajectories of u_{m} and u_{n} vs \tau , where \tau = tr_{m} , t is in decades, u_{m0} = 0.1 , u_{n0} = 0.1 , w = 9.93 , m = 0.0638 , n = 0.646 , \alpha_{n} = 30 , and Tables 2 and 4 are referred for other parameters values.

    We formulated and analyzed a mathematical model to investigate the impact of beetle infestations on forest ecosystems, incorporating key ecological disruptions such as harvesting and wildfires. To understand the system dynamics, we determined the equilibrium and explored their stability conditions. The results were then visually represented through numerical simulations, providing insights into the interactions between tree populations, beetle outbreaks, and environmental stressors. Additionally, we examined the basic reproduction number, a crucial parameter in determining whether beetle populations will persist or decline, and analyzed its ecological significance.

    To mitigate the devastating impact of beetle infestations, we applied an optimal control strategy, aiming to reduce beetle populations while ensuring sustainable forest management. A pesticide variable was incorporated into the mathematical model, offering an additional control mechanism to curb beetle outbreaks. Through sensitivity analysis, we evaluated how different ecological and environmental parameters influence beetle reproduction, infestation severity, and forest resilience.

    Beetle infestations are caused by specific species that bore into tree bark, forming tunnels and damaging the tree's vascular system, which is crucial for water and nutrient transport. Different beetle species tend to specialize in attacking specific tree species, leading to widespread deforestation and ecological imbalance.

    Despite the growing concern over climate-driven beetle outbreaks, our understanding of the physiological interactions between insect infestations, drought stress, and tree defense mechanisms remains incomplete. This knowledge gap hinders our ability to predict long-term forest mortality trends and assess the resilience of tree populations in response to climate change.

    In this study, we examined beetle outbreaks in Bangladesh and pine beetle infestations in forests affected by two major ecological disruptions: Harvesting and fire. We developed a mathematical model that captures the interactions between tree populations, beetle dynamics, and environmental disturbances. The model incorporates:

    ● Harvesting operations, which act as experimental manipulations of forest ecosystems, affecting tree density and beetle habitat availability.

    ● Wildfires, which not only destroy beetle habitats but can also create favorable conditions for infestation by weakening trees.

    ● Pesticide applications, which serve as a control measure to suppress beetle populations and promote tree regeneration.

    By integrating these factors, our model provides a comprehensive framework for understanding beetle infestations and their impact on forest ecology, ecosystem stability, and sustainable resource management. The findings highlight the importance of implementing proactive conservation strategies, including:

    ● Optimal harvesting practices to maintain ecological balance.

    ● Effective pest management interventions to control beetle outbreaks.

    ● Climate adaptation measures to ensure forest sustainability.

    This study emphasizes the need for a holistic approach to beetle infestation control to preserve forest biodiversity, prevent large-scale tree mortality, and support sustainable environmental management.

    Sadia Shihab Sinje: Data curation, formal analysis, methodology, software, writing-original draft; Md Kamrujjaman: Conceptualization, formal analysis, investigation, software, supervision, validation, writing-review & editing; Rubayyi T. Alqahtani: Funding acquisition, resources, validation, writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This work was supported and funded by the Deanship of Scientific Research at Imam Mohammad Ibn Saud Islamic University (IMSIU) (grant number IMSIU-DDRSP2502).

    The authors declare that there is no conflict of interest.



    [1] B. Beckage, L. Gross, W. Platt, Modelling responses of pine savannas to climate change and large-scale disturbance, Appl. Veg. Sci., 9 (2006), 75–82. https://doi.org/10.1111/j.1654-109X.2006.tb00657.x doi: 10.1111/j.1654-109X.2006.tb00657.x
    [2] T. M. Bury, C. T. Bauch, M. Anand, Charting pathways to climate change mitigation in a coupled socio-climate model, PLoS Comput. Biol., 15 (2019), e1007. https://doi.org/10.1371/journal.pcbi.1007000 doi: 10.1371/journal.pcbi.1007000
    [3] J. Huang, M. Kautz, A. M. Trowbridge, A. Hammerbacher, K. F. Raffa, H. D. Adams, et al., Tree defence and bark beetles in a drying world: Carbon partitioning, functioning and modelling, New Phytol., 225 (2020), 26–36. https://doi.org/10.1111/nph.16173 doi: 10.1111/nph.16173
    [4] L. Huo, H. J. Persson, E. Lindberg, Early detection of forest stress from European spruce bark beetle attack, and a new vegetation index: Normalized distance red & SWIR (NDRS), Remote Sens. Environ., 255 (2021), 112240. https://doi.org/10.1016/j.rse.2020.112240 doi: 10.1016/j.rse.2020.112240
    [5] B. H. Aukema, J. Zhu, J. Møller, J. G. Rasmussen, K. F. Raffa, Predisposition to bark beetle attack by root herbivores and associated pathogens: Roles in forest decline, gap formation, and persistence of endemic bark beetle populations, Forest Ecol. Manag., 259 (2010), 374–382. https://doi.org/10.1016/j.foreco.2009.10.032 doi: 10.1016/j.foreco.2009.10.032
    [6] A. C. A. C. McLeod, A review and synthesis of the effects of unsalvaged mountain-pine-beetle-attacked stands on wildlife and implications for forest management, J. Ecosyst. Manage., 7 (2006). https://doi.org/10.22230/jem.2006v7n2a548 doi: 10.22230/jem.2006v7n2a548
    [7] I. Zahan, M. Kamrujjaman, Evolution of dispersal and the analysis of a resource flourished population model with harvesting, Heliyon, 10 (2024), e30737. https://doi.org/10.1016/j.heliyon.2024.e30737 doi: 10.1016/j.heliyon.2024.e30737
    [8] M. M. Adan, M. Kamrujjaman, M. M. Molla, M. Mohebujjaman, C. Buenrostro, Interplay of harvesting and the growth rate for spatially diversified populations and the testing of a decoupled scheme, Math. Biosci. Eng., 20 (2023).
    [9] F. Tasnim, M. Kamrujjaman, Dynamics of Spruce budworms and single species competition models with bifurcation analysis, Biom. Biostat. Int. J., 9 (2020), 217–222. https://doi.org/10.15406/bbij.2020.09.00323 doi: 10.15406/bbij.2020.09.00323
    [10] I. Zahan, M. Kamrujjaman, S. Tanveer, Mathematical study of a resource-based diffusion model with Gilpin–Ayala growth and harvesting, B. Math. Biol., 84 (2022), 120. https://doi.org/10.1007/s11538-022-01074-8 doi: 10.1007/s11538-022-01074-8
    [11] M. A. Wulder, C. C. Dymond, J. C. White, B. Erickson, L. Safranyik, B. Wilson, Detection, mapping, and monitoring of the mountain pine beetle, In: The mountain pine beetle: A synthesis of biology, management, and impacts in lodgepole pine, Natural Resources Canada, 2006, 123154.
    [12] A. M. Gill, S. L. Stephens, G. J. Cary, The worldwide "wildfire" problem, Ecol. Appl., 23 (2013), 438–454. https://doi.org/10.1890/10-2213.1 doi: 10.1890/10-2213.1
    [13] M. Kamrujjaman, K. N. Keya, U. Bulut, M. R. Islam, M. Mohebujjaman, Spatio-temporal solutions of a diffusive directed dynamics model with harvesting, J. Appl. Math. Comput., 69 (2023), 603–630.
    [14] L. T. Bennett, M. A. Adams, Assessment of ecological effects due to forest harvesting: Approaches and statistical issues, J. Appl. Ecol., 41 (2004), 585–598. https://doi.org/10.1111/j.0021-8901.2004.00924.x doi: 10.1111/j.0021-8901.2004.00924.x
    [15] F. Tasnim, M. Kamrujjaman, T. Khan, Structural changes in Mangroves of Sundarban in Bangladesh: Effects of climate change and human disturbances, Model. Earth Syst. Env., 9 (2023), 3553–3566. https://doi.org/10.1007/s40808-023-01699-1 doi: 10.1007/s40808-023-01699-1
    [16] A. Aziz, A. R. Paul, Bangladesh Sundarbans: Present status of the environment and biota, Diversity, 7 (2015), 242–269. https://doi.org/10.3390/d7030242 doi: 10.3390/d7030242
    [17] B. Beckage, L. J. Gross, W. J. Platt, Modelling responses of pine savannas to climate change and large-scale disturbance, Appl. Veg. Sci., 9 (2006), 75–82. https://doi.org/10.1111/j.1654-109X.2006.tb00657.x doi: 10.1111/j.1654-109X.2006.tb00657.x
    [18] F. Tasnim, M. S. Alam, M. Kamrujjaman, Forest dynamics and the analysis of a reaction-diffusion forest model, GANIT J. Bangladesh Math. Soc., 43 (2023), 26–36. https://doi.org/10.3329/ganit.v43i2.70796 doi: 10.3329/ganit.v43i2.70796
    [19] Bangladesh sustainable forest and livelihood (SUFAL), ministry of environment, forest and climate change, Bangladesh, Bangladesh Forest Department, Pest Management Plan, July 2018, Available from: https://bforest.portal.gov.bd/sites/default/files/files/bforest.portal.gov.bd/notices/1d3bd567_72d9_4016_9be8_4de356474014/Pest.
    [20] X. Liu, C. Zhang, Stability and optimal control of tree-insect model under forest fire disturbance, Mathematics, 10 (2022), 2563.
    [21] B. C. Charpentier, M. C. A. Leite, A model for coupling fire and insect outbreak in forests, Ecol. Model., 286 (2014), 26–36. https://doi.org/10.1016/j.ecolmodel.2014.04.008 doi: 10.1016/j.ecolmodel.2014.04.008
    [22] L. Safranyik, B. Wilson, The mountain pine beetle: A synthesis of biology, management and impacts on lodgepole pine, Canadian Forest Service, 2006.
    [23] S. Lenhart, J. T. Workman, Optimal control applied to biological models, Chapman and Hall/CRC, 2007. https://doi.org/10.1201/9781420011418
    [24] L. Zhou, M. Fan, Q. Hou, Z. Jin, X. Sun, Transmission dynamics and optimal control of brucellosis in Inner Mongolia of China, Math. Biosci. Eng., 15 (2017), 543–567. https://doi.org/10.3934/mbe.2018025 doi: 10.3934/mbe.2018025
    [25] E. R. Echegaray, R. A. Cloyd, Effects of reduced-risk pesticides and plant growth regulators on rove beetle (Coleoptera: Staphylinidae) adults, J. Econ. Entomol., 105 (2012), 2097–2106. https://doi.org/10.1603/EC12244 doi: 10.1603/EC12244
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(529) PDF downloads(35) Cited by(0)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog