Forest model dynamics analysis and optimal control based on disease and fire interactions

1.   Introduction

Forest dynamics result from the interaction of natural disturbances with the demographic processes of recruitment, growth, and mortality. Forest dynamics are changing due to environmental changes, such as rising temperatures and carbon dioxide, as well as increasing natural disturbances, such as wildfires, pests, droughts, and hurricanes ^[1]. However, disturbance events cannot be considered in isolation because their effects are often compounded. For example, insects are a major cause of tree mortality, and the death of large numbers of trees increases the likelihood of fire ^[2]. Fire and disease are two disturbances that are spatially widespread. They affect ecosystem processes both in isolation and in concert ^[3]. Independently, fire and disease can have a major impact on ecosystems through the death of susceptible trees ^[4], which can lead to a cascade of changes in the ecosystem. When fire and disease occur together, the interactions between them can either dampen the severity of one of the disturbances or amplify it ^[5]. = Here, we examine the impact of two natural disturbances, disease and fire, on forest development using a forest model as a case study. Our objective is to assess how these disturbances can independently or interactively affect the forest ecosystem. Recurrent fire can influence the transition from savanna to forest ^[6], while disease can also have a significant impact on the development of forest ecosystems. The interaction of fire and disease can have unique consequences for forest ecosystems as each disturbance can affect plant growth and mortality through different pathways. Tree survival is much higher in surface fires ^[7], but the growth of fire-sensitive trees can be limited by frequent burning ^[8]. In these ecological zones, grass is the primary fuel for surface fires. Fire intensity decreases as mature trees grow and form a closed canopy to shade out flammable grasses. Tree growth and fire have a negative feedback loop ^[9]. This relationship is further modified by changes in tree species composition ^[10]. Similar to fire, as diseases can be species or strain specific, tree species composition is important in regulating many diseases. Reduced fire frequency may make ecosystems more susceptible to disease if tree species that increase in low fire frequency environments are more susceptible to disease ^[11]. Simply, in the plant domain, the more densely populated a plant is, the more susceptible it is to disease ^[12]. This view has good empirical support in agricultural literature ^[13] as well as in the extensive literature exploring the Jenson Cornell effect in natural ecosystems ^[14]. Many plant pathogens are specific to a genus or species, such as beech bark disease, which affects several conifer species ^[15]. In systems with specific diseases but diverse tree communities, the colonization and growth of resistant species may compensate for the loss of susceptible species over time. For example, in forests dominated by needle oak, other fire-susceptible species (e.g. maple) are rare, which may limit their ability to replenish rapidly when fire is excluded ^[16]. As forest succession continues, coniferous forests are replaced by red oak and maple ^[17].

In ^[18], Adam assumes that fire only affects saplings. However, there are some trees where adults and saplings have the same sensitivity to fire. In response to this question, we improve the model, by considering the relationship between susceptible and infected trees. The number of parameters fitted has been minimized in order to facilitate the exploration of the effects of the parameters on the model. Compared with ^[19,20], we consider both adults and saplings under the synergistic effect of disease and forest fire. Compared with ^[21,22], tree species may diffuse in the direction of regions rich in resources, or may diffuse in areas with few competitors because of the asymmetric competition between species. The effect of spatial diffusion is considered.

The paper is structured as follows: In Section 2, we consider a forest-disease model and investigate the global asymptotic stability of the disease-free equilibrium point and the internal equilibrium point. In Section 3, we build on the Forest-Disease model by considering the disturbance of forest systems by forest fire to develop a forest-disease-fire model. The local stability of the disease-free equilibrium point as well as the internal equilibrium point is investigated. We consider the case where there is seed dispersal in an ecosystem. We add a diffusion term for susceptible and infected trees in Section 4 to explore the Turing instability of the Forest-Disease-Fire model. Numerical simulations of the Turing branch are given. Furthermore, in Section 5, to explore the effect of fire on the model, we improve the model presented in Section 3, consider the optimal control strategy for the improved model, fit the model using data from the State Forestry Administration on forest fires in each province from 2002 to 2018, and perform sensitivity analysis on key parameters in the model.

2.   Global stability the of Forest-Disease model

In this section, based on ^[18], only considering the effect of disease, we develop a Forest-Disease model that assumes the population is divided into two compartments: $S(t)$ represents the susceptible tree, in which all individuals are susceptible to the disease and $I(t)$ represents the infected tree, in which all individuals are infected by the disease and they can transmit the disease to the healthy one. $r$ represents the growth rate of susceptible trees, $K$ is the environmental capacity, $\frac{\beta}{I+\alpha}SI$ represents nonlinear incidence rate, and $\mu$ is the death rate of infected trees $(I)$. We consider the Forest-Disease model

2.1.   The existence of equilibria

In this subsection, we demonstrate the existence of the equilibria. According to the biological meaning, we need all equilibria to be nonnegative. By using system (2.1), the equilibrium satisfies the equations

Consequently, we can immediately calculate that system (2.1) has only one disease-free equilibrium $E^{0} = (S^{0}, 0)$, where $S^{0} = K$.

Consider the internal equilibrium point $E^{*} = (S^{*}, I^{*})$ of (2.1). By using the second equation of system (2.2), when $I^{*}\neq{0}$, we get $S^{*} = \frac{\mu(I^{*}+\alpha)}{\beta}$. Substituting into the first equation of system (2.2), define $R_{0} = \frac{K\beta}{\mu\alpha}$, then we obtain

Because of

we can observe that:

(ⅰ) If $R_{0} < 1$, then there is no positive root $I^{*}$, and hence, no internal equilibrium point $E^{*} = (S^{*}, I^{*})$;

(ⅱ) If $R_{0} > 1$, then there is a internal equilibrium point $E^{*} = (S^{*}, I^{*})$.

2.2.   Global stability of the Forest-Disease model

We present the global stability analysis corresponding to the system (2.1). Consider the disease-free equilibrium point $E^{0} = (S^{0}, 0)$:

Theorem 1. The disease free equilibrium point $E^{0} = (S^{0}, 0)$ is global asymptotically stable, if $R_{0} < 1$, where $R_{0} = \frac{K\beta}{\mu\alpha}$.

Proof. We prove Theorem (1) by constructing a suitable Lyapunov function, where $R^{2}_{+} = \{(S, I)|S\geq{0}, I\geq{0}\}$. Consider $V(t):R^{2}_{+}\longrightarrow{R}$ such that

where $V_{1}(t) = (S-S^{0}-S^{0}\ln{\frac{S}{S^{0}}}), V_{2}(t) = I-0$. This particular type of Lyapunov function has been considered widely ^[23]. We obtain $V(t) = S-S^{0}-S^{0}\ln{\frac{S}{S^{0}}}+I$, and then

Since $R_{0} < 1$, $\frac{K\beta}{\alpha}-\mu < 0$. Thus, $\frac{S^0\beta}{I+\alpha} \leq \frac{K\beta}{\alpha}$ and we conclude that $\frac{dV}{dt} < 0$ along all the trajectories in $R^{2}_{+}$ and $\frac{dV}{dt} = 0$ at $E^{0} = (S^{0}, 0)$. Therefore, by LaSalle's theorem, $E^{0}$ is globally asymptotically stable. □

Next, we present the global stability analysis of the internal equilibrium $E^{*} = (S^{*}, I^{*})$.

Theorem 2. The internal equilibrium $E^{*} = (S^{*}, I^{*})$ exists and is globally asymptotically stable if $R_{0} > 1$ and $\frac{\beta{I^{*}}}{2\alpha^{2}} < \min\{\frac{r}{K}, \frac{S^{*}\beta}{2\alpha^{2}}\}$.

Proof. We consider the following Lyapunov function:

The time derivative of $V$ along the solution of system (2.1) is:

We conclude that $\frac{dV}{dt} < 0$ along all the trajectories in $R^{2}_{+}$ and $\frac{dV}{dt} = 0$ at $E^{*} = (S^{*}, I^{*})$. Therefore, by LaSalle's theorem, $E^{*}$ is globally asymptotically stable. □

3.   Fire based Forest-Disease model

On the basis of system (2.1), we model the fire effect as a factor that kills sensitive trees, with fire-driven mortality as a function of potential intensity, controlled by grass cover, which is implicitly modeled as a decreasing function of the number of sensitive trees $(S)$ in the measurement area $\frac{1}{\omega{S}+1}$. Ecologically, this is because trees shade grassland, but shading is a nonlinear saturation function of tree abundance. Then, both frequency and intensity are related to mortality through the parameter $\sigma$. The functional form of the tree-grass relationship is determined empirically by nonlinear least squares fitting and model selection ^[18]. We model how grass cover changes the fire frequency effect by adding a term $\sigma$, resulting in a single parameter that incorporates both intensity and frequency into a severity indicator. We consider the Forest-Disease-Fire model as follows:

3.1.   The existence of equilibria

By using system (3.1), the equilibrium satisfies the equations

Since $I = 0$, we substitute into the first equation $rS(1-\frac{S}{K})-\sigma\frac{S}{\omega{S}+1} = 0$. Then,

When $K\omega < 1$, system (3.1) has only one disease-free equilibrium $E^{1} = (S^{1}, 0)$, and thus, the disease-free equilibrium $E^{1} = (S^{1}, 0)$ exists.

We now consider the internal equilibrium point $E^{*} = (S^{*}, I^{*})$ of (3.1):

Lemma 1. The internal equilibrium $E^{*}_{f}$ exists if $R_{2} = \frac{{{S^{*}_{f}}\beta }}{\alpha \mu } > 1$ and $S_{1} < \frac{\mu \alpha}{\beta} < S^{*}_{f}$.

Proof. The internal equilibrium $E^{*}_{f} = (S^{*}_{f}, I^{*}_{f})$ can be determind by

where $a_{0} = r\omega$, $a_{1} = {r - K\omega r +K\omega \beta }$, $a_{2} = K\left({\beta - \alpha \mu \omega - r + \alpha } \right)$, $a_{3} = -K\alpha \mu$.

Define a function $f(S) = a_{0}S^{3}+a_{1}S^{2}+a_{2}S+a_{3}$. Then, $S^{*}_{f}$ is a solution of $f(S) = 0$. Since $a_{0} > 0$, $a_{3} < 0$, $f(S)$ has at least one positive solution $S^{*}_{f}$.

Since $I^{*}_{f} = \frac{\beta}{\mu}S^{*}_{f}-\alpha > 0$, then we need $S^{*}_{f} > \frac{\mu \alpha}{\beta}$. Here, the basic reproduction number $R_2 = \frac{{{S^{*}_{f}}\beta }}{\alpha \mu }$. $f'(S) = 3a_{0}S^{2}+2a_{1}S+a_{2}$, and thus, $S_{1}$ and $S_{2}$ are solutions of $f'(S) = 0$. If $K\omega > 1$, then $0 < S_{1} < S^{*}_{f} < S_{2}$, see Figure 1. According to this, if $K\omega > 1$, $r < \sigma$, and $\frac{{\alpha \mu }}{\beta } > \frac{{K\omega - 1}}{\omega }$ hold, then

Since $S_{1} < \frac{\mu \alpha}{\beta} < S^{*}_{f}$, the internal equilibrium $E^{*}_{f}$ exists. □

3.2.   Local stability of the fire based Forest-Disease model

Next we will perform a local stability analysis on system (3.1). For $E^{1} = (S^{1}, 0)$, the Jacobian matrix is

and

Here, the basic reproduction number $R_1 = \frac{{{S^1}\beta }}{\alpha \mu } < 1$ is defined.

Theorem 3. The disease-free equilibrium $E^{1} = (S^{1}, 0)$ is locally asymptotically stable if $R_1 < 1$, $K\omega < 1$, and $r > \sigma$.

Proof. Clearly, if $K\omega < 1$ and $r > \sigma$ hold, then

According to this, then

and

We conclude that when the basic reproduction number $R_1 < 1$, $K\omega < 1$, and $r > \sigma$, $E^{1}$ is locally asymptotically stable. □

Theorem 4. If $K\omega > 1$, $r < \sigma < K\omega{r}$, and $\frac{{\alpha \mu }}{\beta } > \frac{{K\omega - 1}}{\omega }$ hold, $E^{*}_{f} = (S^{*}_{f}, I^{*}_{f})$ is locally asymptotically stable.

Proof.

and

Clearly,

and

□

In order to verify the asymptotic stability of the system, numerical simulation is performed on the model (3.1). Using the parameter values $r = 0.62117, K = 85, \alpha = 37, \beta = 2, \sigma = 0.60137, \omega = 0.9$, and $\mu = 0.6$, we obtain positive equilibrium points $(S^{*}_{f} = 14.5313, I^{*}_{f} = 11.4376)$. These parameters are selected following reference ^[18]. These parameter values satisfy the conditions of the theorem analyzed above. The above analysis shows that the system is asymptotically stable, as shown in Figure 2. The figure on the left shows the $S^{*}_{f}$ endemic state for fire, and the figure on the right shows the $I^{*}_{f}$ endemic state for fire.

Remark 1. In Section 2, we analyzed the global stability of the no-fire model. If $R_{0} < 1$, $E^{0} = (S^{0}, 0)$ is global asymptotically stable, and if $R_{0} > 1$, the internal equilibrium $E^{*} = (S^{*}, I^{*})$ exists and is globally asymptotically stable, and we obtain the basic reproduction number $R_{0} = \frac{K\beta}{\mu\alpha}$. In Section 3, we considered the local stability of the model under fire disturbance. If $R_1 = \frac{{{S^1}\beta }}{\alpha \mu } < 1$, $E^{1}$ is locally asymptotically stable. If $R_{2} = \frac{{{S^{*}_{f}}\beta }}{\alpha \mu } > 1$, the internal equilibrium $E^{*}_{f}$ exists and is locally asymptotically stable.

4.   Fire and diffusion based Forest-Disease model

Species distributions are both static and dynamic, and as forests evolve, communities maintain biodiversity through diffusion. It is therefore of practical importance to study the diffusion of populations of forest tree species. To get closer to the evolutionary process of biodynamic systems, the spatial evolution of the model needs to be further considered. Taking into account the reproduction of tree species, we add diffusion terms for susceptible and infected trees. For model (3.1) under the coupling of the fire, tree species may diffuse in the direction of regions rich in resources, or may diffuse in areas with few competitors because of the asymmetric competition between species, or may diffuse with the influence of external environmental factors such as water flow. Based on the above complex ecological process, it is reasonable to describe the relationship between susceptible and infected trees by using the reaction-diffusion equation.

4.1.   Linear stability analysis

Consider the diffusion equations

In this section, $D_{S}$ and $D_{I}$ represent diffusion coefficients, and $\Delta = {\partial ^2}/\partial {S^2} + {\partial ^2}/\partial {I^2}$ is the Laplace operator in space. We perform a linear stability analysis ^[24] of system (4.1). To simplify the discussion of system (4.1), we transform the positive equilibrium point $(S^{*}, I^{*})$ into $(0, 0)$ by the means of transformation $(S, I) = (S^{*}+\widetilde{S}, I^{*}+\widetilde{I})$. Then, we can obtain the linearized system as follows:

We denote $p_{1} = r(1-\frac{2S^{*}}{K})-\frac{\sigma}{(1+\omega{S^{*}})^{2}}-\frac{\beta{I^{*}}}{I^{*}+\alpha}, q_{1} = -\frac{S^{*}\alpha\beta}{(I^{*}+\alpha)^{2}}, p_{2} = \frac{\beta{I^{*}}}{I^{*}+\alpha}, q_{2} = \frac{S^{*}\alpha\beta}{(I^{*}+\alpha)^{2}}-\mu$.

Let $d_{1} = \varepsilon{d}, d_{2} = d$, then we can get $\varepsilon = \frac{d_{1}}{d_{2}} > 0$. We rewrite Eq (4.2) as

We can then obtain the Jacobian matrix

Then, the characteristic equation of $J$ is

with

Now, we assume that

Lemma 2. Suppose that $(B_{0})$ holds. Then, the ordinary differential equation system corresponding to system (4.1) is locally asymptotically stable at the positive equilibrium point $(S^{*}, I^{*})$.

Proof. We can get the characteristic equation of the ordinary differential equation system when the diffusion term is not added. That is,

with

Both roots of $D_{0}(\lambda, \varepsilon) = 0$ have negative real parts when $TR_{0} < 0$ and $DET_{0} > 0$. Therefore, system (4.1) is asymptotically stable at the positive equilibrium point $(S^{*}, I^{*})$. Hence, the lemma is proved. □

4.2.   Analysis of Turing instability

Now we consider the existence conditions for Turing instability.

Assume that

$(B_{1})$ $0 < \varepsilon < \varepsilon_{1}$, where $\varepsilon_{1}\triangleq\frac{-2p_{2}q_{1}+p_{1}q_{2}-2\sqrt{p^{2}_{2}q^{2}_{1}-p_{1}p_{2}q_{1}q_{2}}}{q^{2}_{2}}$.

$(B_{2})$ $0 < \varepsilon < \varepsilon_{2}(d)$, where $\varepsilon_{2}(d)\triangleq\frac{p_{1}}{d\pi^{2}-q_{2}}$. $\varepsilon = \varepsilon(d)$ decreases monotonically in $d$ and intersects with $\varepsilon = \varepsilon_{1}$ at the point $d = d_{0}$. We take $\varepsilon_{B}(d) = \min\{{\varepsilon_{1}, \varepsilon_{2}(d)}\}$, then

Proof. Let $x = dk^{2}\pi^{2} > 0$, then we rewrite $DET_{k}$ as

Then, we get that the symmetry axis is $x = \frac{p_{1}+q_{2}\varepsilon}{2\varepsilon}$, and at this time $DET_{k}$ can be taken to a minimum. That is,

Since $x > 0$, we have $p_{1}+q_{2}\varepsilon > 0$. If we want $DET_{k_{min}} < 0$ with the condition $p_{1}+q_{2}\varepsilon > 0$, we must satisfy $p_{1}q_{2}-p_{2}q_{1}-\frac{(p_{1}+q_{2}\varepsilon)^{2}}{4\varepsilon} < 0$. We can obtain

when

Let $\varepsilon_{1} = \frac{-2p_{2}q_{1}+p_{1}q_{2}-2\sqrt{p^{2}_{2}q^{2}_{1}-p_{1}p_{2}q_{1}q_{2}}}{q^{2}_{2}}$. Then, when $0 < \varepsilon < \varepsilon_{1}$, it ensures that the symmetry axis is greater than 0, and Turing instability occurs at $(S^{*}, I^{*})$.

Next, we take the value of $x$ between the symmetry axis and the right root of $DET_{k}$. Then, $\frac{p_{1}+q_{2}\varepsilon}{2\varepsilon}\leqslant{x}\leqslant\frac{p_{1}+q_{2}\varepsilon+\sqrt{\triangle}}{2\varepsilon}$. $k$ can be taken to the minimum when $x$ takes a value on the symmetry axis, and we get

and

By making sure that $\sqrt{\frac{1}{2}\cdot\frac{p_{1}+q_{2}\varepsilon}{d\varepsilon\pi^{2}}} > \sqrt{\frac{1}{2}}$, we can obtain $0 < \varepsilon < \frac{p_{1}}{d\pi^{2}-q_{2}}$. When $p_{1} > 0, q_{1} < 0, q_{2} < -p_{1}, p_{2} > \frac{p_{1}q_{2}}{q_{1}}$, let $\varepsilon_{2}(d) = \frac{p_{1}}{d\pi^{2}-q_{2}}$, and we can find that $\varepsilon = \varepsilon_{2}(d)$ decreases monotonically in $d$ and intersects with $\varepsilon = \varepsilon_{1}$ at the point $d = d_{0}$, where

□

The next step is to determine the bounds of the Turing instability. We denote

where $d_{k} = \frac{p_{1}q_{2}-p_{2}q_{1}}{k^{2}\pi^{2}p_{1}}$. Then, $DET_{k} = 0$, when $\varepsilon = \varepsilon_{*}(k, d)$. That is,

Then

and we can get

Set $d_{M} = \frac{-p_{2}q_{1}+p_{1}q_{2}+\sqrt{p^{2}_{2}q^{2}_{1}-p_{1}p_{2}q_{1}q_{2}}}{p_{1}}\cdot\frac{1}{k^{2}\pi^{2}}$, then $\frac{d\varepsilon_{*}}{dx} = 0$, when $d = d_{M}$. At this time, $\varepsilon_{*}(k, d)$ reaches the maximum value $\varepsilon_{1}$.

Lemma 3. Suppose that $(B_{0})$ holds, function $\varepsilon = \varepsilon_{*}(k, d)$ has the following properties: (1) When $0 < d < d_{M}(k)$, $\varepsilon_{*}(k, d)$ increases monotonically, and when $d > d_{M}(k)$, $\varepsilon_{*}(k, d)$ decreases monotonically. (2) As for $k_{i} < k_{i+1}, k_{i}\in{N}, i = 1, 2, 3, \cdots$, there is only one root $d_{k1, k2}\in(d_{M}(k_{2}), d_{M}(k_{1}))$ satisfying $\varepsilon_{*}(k_{1}, d) = \varepsilon_{*}(k_{2}, d)$ for $d > 0$. Furthermore,

In Figure 3, we characterize a graph of functions $\varepsilon = \varepsilon_{1}, \varepsilon = \varepsilon_{2}(d)$, and $\varepsilon = \varepsilon_{*}(k, d), d > 0, k = 1, 2, 3, \cdots$. The value of parameters are selected following reference ^[18], and the value of $p_{1}, q_{1}, p_{2}, q_{2}$ are calculated from the model (4.2). In Figure 4, we present a graph of the Turing bifurcation line $\varepsilon = \varepsilon_{*}(k, d), d > 0$.

4.3.   Numerical simulations of Turing bifurcation

Let $D_{S} = 0.03, D_{I} = 0.41$, $r = 0.721, K = 90, \sigma = 0.601, \alpha = 37.8537234, \beta = 0.318, \omega = 0.055$, and $\mu = 0.02$. We then obtain the positive equilibrium of system (4.2), $(S^{*}, I^{*}) = (17.252,236.374)$. Through ($B_{1}$), ($B_{2}$), we obtain $\varepsilon_{1} = 0.0545$,

and

Then, we can find $\varepsilon > \varepsilon_{*}(d)$, so system (4.2) is asymptotically stable at $(S^{*}, I^{*})$, see Figure 5. By setting $k = 1$, we obtain $d_{0, 1} = +\infty$ and $d_{1, 2} = 0.0250246$, and select $d = d_{1} = 0.12\in(d_{0, 1}, d_{1, 2})$; thus, $\varepsilon_{*} = \varepsilon_{*}(1, 0.12) = 0.0011360854$. System (4.2) with $d = 0.12$ undergoes a Turing bifurcation near equilibrium $(17.252,236.374)$ at $\varepsilon = 0.0011360854$. Figure 6 shows that system (4.2) has a Turing instability near the equilibrium point $(17.252,236.374)$ at this time.

5.   Optimal control strategy and sensitivity analysis

To better describe the effect of fire on model (3.1), we rewrite the model as

where ${\mu _1} = \mu - \sigma$ is the natural mortality rate of diseased trees, and $\sigma{I}$ stands for diseased trees which died due to fire. In forest ecosystems, fire is a potential threat to tree growth, and fire can directly affect forest insect populations and abundance. Therefore, one of the most important tools to reduce fire damage and insect outbreaks is to control fire intensity. Finding the right frequency of fire is an optimization problem such that damage to trees is reduced and beetles can be eliminated. In this section, we analyze the influence of key parameters on the model and use forest area data coupled with the model to construct a more complex forest model in both a fire and insect infested environment in order to gain insight into the effects of the synergistic effects of both fire and pest natural disturbances on the model.

5.1.   Optimal control strategy

Optimization strategies have always played a key role in optimizing ecosystem structure and resource use patterns ^[25] and in identifying appropriate management actions to ensure that the overall benefits of the system are maximized. In the following sections, we will discuss optimal fire frequency management strategies. The $\sigma$ in the system is set to be the time-dependent function $\sigma(t)$, describing the time-varying fire intensity, and there exists $\overline{\sigma}$ such that $0\leq\sigma(t)\leq\overline{\sigma}$. Within the fixed time $[0, T]$ with $T > 0$, the constraint set reads

Since the parameter $\sigma(t)$ is related to mortality, we aim to minimize the frequency and intensity of fire. The quadratic optimal objective function is

with $I(0) = I_{0}, \sigma(0) = \sigma_{0}$. The optimal control problem is rewritten

where the minimization of $J(\sigma(t))$ is subject to (5.1) and $\sigma\in{U}$ with $U$ from (5.2).

Using the method in ^[26], we check the existence of optimal control $\sigma(t)$ by satisfying H5–H8 in the following:

H5: The set of state and control variables are nonempty;

H6: The set $U$ of the control variables is closed and convex;

H7: The right side of each equation in the control problem is continuous with a bounded sum of controls and states above. This can be written as a linear function of $U$ with coefficients that depend on time and state ^[27];

H8: There exists constants $\sigma > 0$ such that the integrand $L(\sigma(t))$ of the objective functional $J$ is convex and satisfied.

We have $L(\sigma(t)) = \frac{1}{2}(I^{2}(t)+\sigma^{2}(t))$.

Theorem 5. For the optimal control problems (5.2) and (5.3), there exists an optimal control $\sigma^{*}$ such that $J^{*}(\sigma^{*}(\cdot)(t)) = \min\limits_{\sigma(t)\in{U}}J(\sigma(t))$.

To obtain the minimum value of (5.3), we give the Hamiltonian function:

where $\lambda_{1}, \lambda_{2}$ are adjoint variables satisfying the following adjoint equations:

The optimality conditions are given by

Theorem 6. The optimal control of (5.2) and (5.3) is given by

where $\lambda_{1}, \lambda_{2}$ are adjoint variables satisfying the second and third Maximum Principle condition.

Proof. According to Pontryagin's Maximum Principle, finding the optimal control of (5.2) and (5.3) is equivalent to minimizing the Hamiltonian function $H$ defined previously.

$\frac{\partial{H}}{\partial{\sigma}}|_{(\sigma^{*}, \lambda_{1}, \lambda_{2}, t)} = 0$ gives $\sigma^{*}-\lambda_{1}\frac{S}{\omega{S}+1}-\lambda_{2}{I} = 0$, where the adjoint variables are satisfied. Furthermore,

and the transversality conditions are $\lambda_{1}(T) = 0, \lambda_{2}(T) = 0$.

Moreover, since $\sigma^{*}\in{U}$, using the lower and upper bounds of $\overline{\sigma}(t)$, the optimal $\sigma^{*}(t)$ can be characterized by (5.4). The curves of $\sigma^{*}(t)$, $S^{*}(t)$, and $I^{*}(t)$ with time for different $\omega$ are given in Figure 7.□

For the improved model (5.1), in the transient dynamic process, with $\omega = 0$, $\sigma = 0$ under optimal control, which in practice can be considered as a forest area without fire, the overall tree cover in the presence of disease is lower than the tree cover in a forest area with low frequency fires ($\omega = 0.054879880$). This region has a higher $I^{*}$. When $\omega = 0.9$, grass cover was higher and fire frequency peaked, limiting the rate of disease spread. Overall tree cover at this time was higher than in diseased woodlands ($\omega = 0$.)

Fire data from some provinces in China are used next. (Shanghai, Jilin, Inner Mongolia, etc.). With these data, parameter values were obtained for several provinces from 2002 to 2018 and are shown in Table 1 (some of the data are given in Table 1), to see more clearly the trends in biomass.

According to the data in Table 1, Figure 8 shows the biomass loss under different fire frequencies in the same year. Figure 9 shows the survival of tree biomass over time under different fire frequencies.

The model is fitted to the actual data, and the results show that the model can predict to some extent the change of forest live wood volume in the next 30 years when the fire frequency is relatively stable, see Figure 10. In 2020, forest tree biomass in Jilin Province reached 138.37${m^3}/ha$, and through the model, we predicted that Jilin's tree biomass in 2020 would be 132.395${m^3}/ha$, which illustrates the reasonableness of the model.

5.2.   Sensitivity analysis

Sensitivity analysis is the most commonly used test in the process of model optimization. Usually, only parameters with large uncertainties can be selected for sensitivity analysis without the need to calculate sensitivity coefficients for each parameter. A study is carried out to analyze the influence of the key parameter $\sigma$ in the sensitivity characterization in order to better understand the effects of fire on forests. The specific steps of the algorithm of the coefficient of sensitivity are follows ^[28]. We consider the function $\dot{y(t)} = f(t, y, m)$, and the absolute sensitivity of variable $y_{i}$ to parameter $M$:

We denote the relative sensitivity of a variable to a parameter

and the absolute sensitivity equation of parameter $P_{i}$

According to the steps of sensitivity analysis, the sensitivity analysis and relative sensitivity analysis are developed for the system, and then the sensitivity equation of the system contains $4$ equations:

Assuming the values of time step $\bigtriangleup{t} = 2000$ and fire frequency $\sigma = 0.6037$, the initial values are chosen to be consistent with Section 3, and the sensitivity analysis of $\sigma$ on susceptible and infected trees was been carried out using the Runge-Kutta integration method of 4 (RK-4) for system (5.1), computed in the MATLAB R2021a environment. The conclusions obtained are able to reflect the degree of influence of $\sigma$ on each component (susceptible and infected trees) in the forest, as shown in Figure 11.

As shown in Figure 11, the sensitivity $P_{1}$ and relative sensitivity $\sigma P_{1}/S(t)$ to susceptible trees $S(t)$ show a trend of less than 0. The results show a decrease in $S(t)$ with the increase of $\sigma$, and when $t = 0.5$, $P_{1} = -1.606\times{10^{-4}} < 0$, $\frac{\sigma{P_{1}}}{S} = -1.531\times{10^{-7}} < 0$. That means that when $\sigma$ increases by $10\%$, $S$ will decrease by $1.531\times{10^{-6}}\%$, and the sensitivity $P_{2}$ and relative sensitivity $\sigma P_{2}/I(t)$ to susceptible trees $I(t)$ show a trend of less than 0. The results show a decrease in $I(t)$ with the increase of $\sigma$, and when $t = 0.5$, $P_{2} = -3.342 < 0$, $\frac{\sigma{P_{2}}}{S} = -0.026 < 0$, it means that when $\sigma$ increases by $10\%$, $I$ will decrease by $0.260\%$.

6.   Discussion

We analyze the process of synergistic interaction between fire and disease in the case of tree species where fire acts on both adults and saplings. A diffusion term is added to this, where tree species may diffuse in the direction of regions rich in resources, or may diffuse in areas with few competitors because of the asymmetric competition between species, or may diffuse with the influence of external environmental factors such as water flow. Based on the above complex ecological process, it is reasonable to describe the relationship between susceptible and infected trees by using the reaction-diffusion equation.

We present mathematical models with generalizations that provide the basis for exploring the effects of fire and insect pests on grasslands and forests. In Section 2, the disease-free equilibrium points and coexistence equilibrium points of the fire-free model are shown. The disease-free equilibrium points and coexistence equilibrium points of the disease model in the fire-free state are globally stable under certain conditions. In Section 3, we factor in fire and grassland cover in the model. Model (3.1) is developed to demonstrate the local stability of the equilibrium point of the model. The dynamical properties of the model in the case of band diffusion are discussed in Section 4. And in Section 5, we discuss the optimal control strategy for fire to find the appropriate frequency of wildfire occurrence. The results show that the biological survival in forest areas without fire but with disease is smaller than that in forest areas where fire occurs and disease is present.

Based on the forest fire data provided by the State Forestry Administration from 2002-2018, the forest development trend after fire was predicted, and the results showed that the fire factor $\sigma$ plays an important role in the evolution of the forest. Further, we use sensitivity analysis to determine the degree of influence of the parameters on the model. Our empirical analysis shows that disease invasion reduces tree biomass, but that fire itself can significantly affect forest biomass in an area, which also takes into account the effect of grassland cover. Our model predictions suggest that the interaction between disease and fire in forested ecoregions jointly affects the distribution of forest ecosystems.

Use of AI tools declaration

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

Acknowledgments

The authors wish to thank the editors and reviewers for their helpful comments.

Conflict of interest

The authors declare no competing interests.