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

Abbreviation

1.   Introduction

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.

2.   Case study in Bangladesh

Bangladesh, with 160 million inhabitants and 147,570 $km^{2}$ 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].

3.   Mathematical model

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:

with initial conditions,

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:

● $r_m$: The intrinsic growth rate of trees.

● $r_n$: The intrinsic reproductive rate of beetles.

● $r$: The critical attack density of beetles.

● $k_m$: The carrying capacity of trees.

● $k_n$: The carrying capacity of beetles.

● $g_k$: The fraction of trees killed by beetles.

● $\alpha$: The defensive rate of trees.

● $\gamma$: 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:

● $Q_m$: 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.

● $P \frac{M}{k_m}$: The probability or frequency of fire, assumed to be related to the density of trees in the forest.

● $P \frac{M}{k_m} Q_m$: The fire strength, which determines additional tree mortality and is influenced by both fire intensity and tree density.

● $Q_n$: 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 ($Q_m = 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 ($Q_n = 0$).

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

In addition, when harvesting and fire are present in the system, they are represented by $E \neq 0$ and $Q_m \neq 0$. 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.

3.1.   Dimensional analysis

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:

Therefore, the dimensionless model is defined as follows:

with initial conditions,

4.   Steady state and stability analysis

4.1.   Presence of both disruptions: Harvesting and fire

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:

with initial conditions,

For equilibrium solutions,

such that

From (4.3) we get,

From (4.4) we have,

If $n = 0$, then from (4.5) we get, $m = \frac{1-q_{s}}{1+q_{m}}.$

Thus, we obtain the boundary equilibrium,

Let the coexistence equilibrium is $E_{*} = (m^{*}, n^{*})$, here $m^{*} > 0$, $n^{*} > 0$, then from (4.5), we get,

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

We can write (4.8) as,

where

Since $\gamma_{1} > 0$, $k_{n} > 0$, $r_{e} > 0$, $q_{m} > 0$, $q_{n} > 0$, $0 < q_{s} < 1$, and $f_{1} < 0$, then $\lim_{n \rightarrow \; +\infty}F(n) = -\infty$. If $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n}\left(1-q_{s}\right)$, then $F(0) > 0$. According to the real continuation method, there is $n^{*} > 0$ such that $F(n^{*}) = 0$. As $k = g_{k}k_{n}$, $0 < g_{k} < 1$, $0 < q_{s} < 1$, then we have $k_{n}-k > 0$, $\left(1-q_{s}\right) < 1$, $r > q_{s}r$, $k_{n} > q_{s}k_{n}$. Consider, $k_{n}-q_{s}k_{n}-k > 0$ and we can write,

Thus, we can conclude that when $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n}\left(1-q_{s}\right)$, the model has at least one coexistence equilibrium $E_{*} = (m^{*}, n^{*})$. We can consider the basic reproduction number,

4.1.1.   Stability analysis

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

Assume that,

The Jacobian matrix is,

which yields

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

The Jacobian matrix at point $E_{0}$ is,

The eigenvalues of the matrix $J_{0}$ can be defined as $\lambda_{0}$, hence, we have $|J_{0}- \lambda_{0}I|$, where, $I$ is the identity matrix. The characteristic equation is,

where

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

Theorem 1. Let $q_{s} > 0$, $k > 0$, ${q_{m} > 0}$, $r_{e} > 0$, $r > 0$, $\alpha_{n} > 0$ such that all of them are positive parameters. Then, we have,

1) If $u_{1}v_{2} > 0$ and $u_{1}+v_{2} < 0$, then the boundary equilibrium $E_{0}$ is locally asymptotically stable.

2) If $u_{1}v_{2} < 0$ or, $u_{1}+v_{2} > 0$, then the boundary equilibrium $E_{0}$ is unstable.

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

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 $E_{0} = \left(\frac{1-q_{s}}{1-q_{m}}, 0\right)$ are:

To ensure stability, the conditions are as follows:

and

Hence, if $u_{1}v_{2} > 0$ and $u_{1}+v_{2} < 0$, the boundary equilibrium $E_{0}$ is locally asymptotically stable.

If the equilibrium $E_{0}$ is unstable, then at least one of the conditions in (4.13) is not satisfied, such that $\text{tr}(J_0) > 0$ or $\text{det}(J_0) < 0$. Now, we have:

Otherwise, we have:

Consequently, if $u_{1}v_{2} < 0$ or $u_{1}+v_{2} > 0$, then the boundary equilibrium $E_{0}$ is unstable.

□

From (4.14) we have,

If

The condition (4.18) is not valid as $q_{s}\ngtr 1$ such that for stable equilibrium, $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right)\ngtr q_{n}\left(1-q_{s}\right)$. Otherwise we can write,

The prerequisite (4.19) is possible for the stability as $q_{s} < 1$ and that is why for the stable equilibrium $E_{0}$, we have $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) < q_{n}\left(1-q_{s}\right)$. So, for stability, the product of $\left(r_{e}-\alpha_{n}\right)$ and $\left(1+q_{m}\right)$ is less than the product of $q_{n}$ and $\left(1-q_{s}\right)$.

Following the relation (4.15), we get,

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

and, we have $q_{s} < 1, \; \; \text{and}\; \; \left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n}\left(1-q_{s}\right).$ The constraint is possible as $q_{s} < 1$ and for the unstable equilibrium $E_{0}$, the required condition is,

The next mathematical inequality prospect, $u_{1} > 0$ and $v_{2} < 0,$ that yields $q_{s}\; > \; 1$, and $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) < q_{n}\left(1-q_{s}\right)$, respectively.

The Eq (4.17) implies

From (4.20), we have for the instability $\frac{\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right)}{\left(1+q_{m}+q_{n}\right)\left(1-q_{s}\right)} > 1$. So, the ratio of $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right)$ and $\left(1+q_{m}+q_{n}\right)\left(1-q_{s}\right)$ is greater than 1.

Theorem 2. Let $q_{s} > 0$, $k > 0$, ${q_{m} > 0}$, $r_{e} > 0$, $r > 0$, $\alpha_{n} > 0$. Then we have,

1) If $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) < q_{n}\left(1-q_{s}\right)$, then the boundary equilibrium $E_{0}$ is locally asymptotically stable.

2) If $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n}\left(1-q_{s}\right)$, then the equilibrium $E_{0}$ is unstable.

Proof. Consider the basic reproduction number,

For stable equilibrium,

For unstable equilibrium,

Thus, if $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) < q_{n}\left(1-q_{s}\right)$, then the boundary equilibrium $E_{0}$ is locally asymptotically stable. If $\left(r_{e}-\alpha_{n}\right)\left(1+q_{m}\right) > q_{n}\left(1-q_{s}\right)$, then the boundary equilibrium $E_{0}$ is unstable. □

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

The Jacobian martix at $E_{*}$ is,

Consider,

Here,

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

while $\text{tr}(J_{*}) = w_{1}+z_{2}, \; \; \text{det}(J_{*}) = w_{1}z_{2}-w_{2}z_{1}.$ For stability analysis, following Theorem 1, we have $\text{tr}(J_{*}) < 0$ and $\text{det}(J_{*}) > 0$.

The equilibrium is stable as long as, $\text{tr}(J_{*}) < 0 \Rightarrow w_{1} < -z_{2}, \; \; \text{and}\; \; \text{det}(J_{*}) > 0 \Rightarrow w_{1} > \frac{w_{2}z_{1}}{z_{2}}.$ It is a remark for the stable situation, $w_{1}$ is less than $-z_{2}$ and $w_{1}$ is greater than the ratio of $w_{2}z_{1}$ and $z_{2}$. To check the unstable equilibrium point $E_{*}$, at least one of the conditions of (4.13) is not satisfied, such that $\text{tr}(J_*) > 0$ or, $\text{det}(J_*) < 0$ and we obtain the following relations.

Hence, we can conclude that for instability, $w_{1}$ is greater than $-z_{2}$ or, $w_{1}$ is less than the ratio of $w_{2}z_{1}$ and $z_{2}$.

Theorem 3. Let any open region $\Omega$ containing equilibrium $E_{*}$ of system (4.1) and $R_{0} > 1$. If $\alpha_{n} < r_e \left(n+\frac{1}{\gamma_{1}}\right)$, then equilibrium $E_{*}$ is globally asympotically stable.

Proof. Assume any open region,

Consider the scalar function, $D = \frac{1}{mn}.$ Next from the model (4.1), we have,

Taking partial derivatives,

Check that,

From (4.27), consider,

Using (4.28), we remark that

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

such that $R_{0} > 1$. Hence the equilibrium $E_{*}$ is globally asymptotically stable. □

4.1.2.   Properties of basic reproduction number $R_{0}$

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

Now,

where, $q_{m} = \frac{Q_{m}P}{r_{m}} \; \; and\; \; q_{s} = \frac{EL}{r_{m}}$.

Consider,

4.1.3.   A numerical example for basic reproduction number $R_{0}$

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:

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

4.2.   Presence of fire disruption

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,

with initial conditions,

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

and, the mathematical equations to establish the next results.

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

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

We can write (4.39) as,

where

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:

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

Consequently, we conclude that if the condition:

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

4.2.1.   Stability analysis

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

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

The Jacobian matrix at point $E_{0}$ is,

such that

The characteristic equation is,

where,

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:

This implies:

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

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

Additionally, the trace condition is given by:

This leads to:

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:

This implies:

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

which confirms instability. Hence, for instability,

Additionally, if

we obtain:

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

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

The Jacobian martix at $E_{*}$ is,

Assume,

where,

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

Here,

For stability,

For instability,

4.3.   Presence of harvesting disruption

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:

with initial conditions,

To determine equilibrium, we set

This leads to system

Solving these equations, we obtain two equilibria:

● The boundary equilibrium,

● The coexistence equilibrium,

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

From the quadratic equation,

where

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^{*})$.

4.3.1.   Stability analysis

Define the Jacobian matrix of the model (4.46),

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

where

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

Here,

For stability, the conditions are:

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

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

If $E_{0}$ is unstable, then:

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

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,

Consider,

Here,

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

Now, for stable equilibrium,

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,

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}}$.

5.   Numerical examples

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).

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$.

5.1.   Computational results in the presence of fire and harvesting disruptions

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,

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.

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

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.

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.

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.

5.2.   Computational results in the presence of fire disruptions

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

which is asymptotically stable under the condition:

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.

The endemic equilibrium of the system is denoted as

We know that if the following condition holds:

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$.

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.

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.

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.

5.3.   Numerical analysis of the model considering harvesting disruptions

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

We establish that if

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$.

The endemic equilibrium of the system is given by

We know that if

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.

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.

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.

6.   How to control beetle attacks in forests

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.

7.   Optimal control strategy

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],

Consider the objective function,

The Hamiltonian function is,

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,

with

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,

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

Hence, we conclude that,

□

For the minimization problem,

The problem is maximization problem as long as

whereas the problem is minimization one, if

7.1.   Mathematical model with pesticide

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

with initial conditions,

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.

By introducing the following non-dimensional variables and parameters,

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

with initial conditions,

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:

The optimal objective function is formulated as:

The optimal control problem is then defined as:

The integrand of the objective function is given by:

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

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

Using Pontryagin's Maximum Principle, we have,

Here,

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

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].$

7.2.   Computational results of pesticide model

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

we determine the value of $\alpha_{n}$. Additionally, using the 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:

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

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.

We assume the conditions:

and assign specific values to achieve accurate results. Here,

Thus, it follows that:

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.

8.   Sensitivity analysis

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

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

with initial condition,

The Jacobian of the system (7.13) is,

For pesticide parameter $w$, we have

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

where

8.1.   Numerical example for sensitivity analysis

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.

9.   Conclusions

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.

Author contributions

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.

Use of Generative-AI tools declaration

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

Acknowledgments

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

Conflict of interest

The authors declare that there is no conflict of interest.