Loading [MathJax]/jax/output/SVG/jax.js
Research article

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

  • Received: 08 November 2023 Revised: 24 December 2023 Accepted: 27 December 2023 Published: 03 January 2024
  • MSC : 34C23, 35K57

  • Three models for the propagation of forest disease are revisited to include the effect of forest fires and disease spread. We study the global stability of the forest-disease model in the absence of forest fires and the spread of disease. When forest fires caused by grass cover are considered, we show that the equilibrium points are locally asymptotically stable. If both forest fires and the spread of disease exist in the second model, then Turing instability can occur. In this case, the system exhibits complex dynamic behavior. To determine the effect of fire on the forest disease model, we obtain the optimal control expression of the key parameter fire factor, and carry out sensitivity analysis. Finally, we use forest biomass data of some provinces in China from 2002 to 2018 for numerical simulation, and the results are in agreement with the theoretical analysis.

    Citation: Xiaoxiao Liu, Chunrui Zhang. Forest model dynamics analysis and optimal control based on disease and fire interactions[J]. AIMS Mathematics, 2024, 9(2): 3174-3194. doi: 10.3934/math.2024154

    Related Papers:

    [1] Sadia Shihab Sinje, Md Kamrujjaman, Rubayyi T. Alqahtani . Forest beetle infestation and its impact on ecosystems: Effects of harvesting practices and fire disruptions. AIMS Mathematics, 2025, 10(4): 9933-9973. doi: 10.3934/math.2025455
    [2] Ahmed Alshehri, Saif Ullah . Optimal control analysis of Monkeypox disease with the impact of environmental transmission. AIMS Mathematics, 2023, 8(7): 16926-16960. doi: 10.3934/math.2023865
    [3] 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
    [4] Ahmed Hamza Osman, Ashraf Osman Ibrahim, Abeer Alsadoon, Ahmad A Alzahrani, Omar Mohammed Barukub, Anas W. Abulfaraj, Nesreen M. Alharbi . Breaking new ground in cardiovascular heart disease Diagnosis K-RFC: An integrated learning approach with K-means clustering and Random Forest classifier. AIMS Mathematics, 2024, 9(4): 8262-8291. doi: 10.3934/math.2024402
    [5] Awatif J. Alqarni . Modeling and numerical simulation of Lumpy skin disease: Optimal control dynamics approach. AIMS Mathematics, 2025, 10(4): 10204-10227. doi: 10.3934/math.2025465
    [6] 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
    [7] Jinji Du, Chuangliang Qin, Yuanxian Hui . Optimal control and analysis of a stochastic SEIR epidemic model with nonlinear incidence and treatment. AIMS Mathematics, 2024, 9(12): 33532-33550. doi: 10.3934/math.20241600
    [8] Ahmed Alshehri, Miled El Hajji . Mathematical study for Zika virus transmission with general incidence rate. AIMS Mathematics, 2022, 7(4): 7117-7142. doi: 10.3934/math.2022397
    [9] Suprawee Lertnaweephorn, Usa Wannasigha Humphries, Amir Khan . Stability analysis and optimal control for leaf brown spot disease of rice. AIMS Mathematics, 2023, 8(4): 9602-9623. doi: 10.3934/math.2023485
    [10] Kai Zhang, Yunpeng Ji, Qiuwei Pan, Yumei Wei, Yong Ye, Hua Liu . Sensitivity analysis and optimal treatment control for a mathematical model of Human Papillomavirus infection. AIMS Mathematics, 2020, 5(3): 2646-2670. doi: 10.3934/math.2020172
  • Three models for the propagation of forest disease are revisited to include the effect of forest fires and disease spread. We study the global stability of the forest-disease model in the absence of forest fires and the spread of disease. When forest fires caused by grass cover are considered, we show that the equilibrium points are locally asymptotically stable. If both forest fires and the spread of disease exist in the second model, then Turing instability can occur. In this case, the system exhibits complex dynamic behavior. To determine the effect of fire on the forest disease model, we obtain the optimal control expression of the key parameter fire factor, and carry out sensitivity analysis. Finally, we use forest biomass data of some provinces in China from 2002 to 2018 for numerical simulation, and the results are in agreement with the theoretical analysis.



    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.

    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, βI+αSI represents nonlinear incidence rate, and μ is the death rate of infected trees (I). We consider the Forest-Disease model

    {dSdt=rS(1SK)βI+αSI,dIdt=βI+αSIμI. (2.1)

    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

    {rS(1SK)βI+αSI=0,βI+αSIμI=0. (2.2)

    Consequently, we can immediately calculate that system (2.1) has only one disease-free equilibrium E0=(S0,0), where S0=K.

    Consider the internal equilibrium point E=(S,I) of (2.1). By using the second equation of system (2.2), when I0, we get S=μ(I+α)β. Substituting into the first equation of system (2.2), define R0=Kβμα, then we obtain

    μrI2+(2μrαKrβ+Kβ2)I+μrα2Krβα=0.

    Because of

    μrα2Kαβr=rα(μαKβ)=rμα2(1Kβμα)=rμα2(1R0)

    we can observe that:

    (ⅰ) If R0<1, then there is no positive root I, and hence, no internal equilibrium point E=(S,I);

    (ⅱ) If R0>1, then there is a internal equilibrium point E=(S,I).

    We present the global stability analysis corresponding to the system (2.1). Consider the disease-free equilibrium point E0=(S0,0):

    Theorem 1. The disease free equilibrium point E0=(S0,0) is global asymptotically stable, if R0<1, where R0=Kβμα.

    Proof. We prove Theorem (1) by constructing a suitable Lyapunov function, where R2+={(S,I)|S0,I0}. Consider V(t):R2+R such that

    V(t)=V1(t)+V2(t),

    where V1(t)=(SS0S0lnSS0),V2(t)=I0. This particular type of Lyapunov function has been considered widely [23]. We obtain V(t)=SS0S0lnSS0+I, and then

    dVdt=SS0S[rS(1SK)βI+αSI]+βI+αSIμI=(SS0)[r(1SS0)βI+αI]+βI+αSIμI=rS0(SS0)2+S0βI+αIμI=rS0(SS0)2+(S0βI+αμ)I<(S0βI+αμ)I.

    Since R0<1, Kβαμ<0. Thus, S0βI+αKβα and we conclude that dVdt<0 along all the trajectories in R2+ and dVdt=0 at E0=(S0,0). Therefore, by LaSalle's theorem, E0 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 R0>1 and βI2α2<min{rK,Sβ2α2}.

    Proof. We consider the following Lyapunov function:

    V=SSSlnSS+IIIlnII.

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

    dVdt=SSS[rS(1SK)βI+αSI]+III(ISβα+IμI)=(SS)[r(1SK)IβI+α]+(II)(βSI+αμ)rK(SS)2+β(SS)(II+αII+α)+β(II)(SI+αSI+α)rK(SS)2+β(II)(SISI)(I+α)(I+α)(SS)2(βI2α2rK)+(II)(βI2α2Sβ2α2).

    We conclude that dVdt<0 along all the trajectories in R2+ and dVdt=0 at E=(S,I). Therefore, by LaSalle's theorem, E is globally asymptotically stable.

    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 1ω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 σ. 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 σ, resulting in a single parameter that incorporates both intensity and frequency into a severity indicator. We consider the Forest-Disease-Fire model as follows:

    {dSdt=rS(1SK)βI+αSIσSωS+1,dIdt=βI+αSIμI. (3.1)

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

    {rS(1SK)βI+αSIσSωS+1=0,βI+αSIμI=0.

    Since I=0, we substitute into the first equation rS(1SK)σSωS+1=0. Then,

    S1=r(Kω1)+r2(1Kω)2+4rKω(rσ)2rω.

    When Kω<1, system (3.1) has only one disease-free equilibrium E1=(S1,0), and thus, the disease-free equilibrium E1=(S1,0) exists.

    We now consider the internal equilibrium point E=(S,I) of (3.1):

    Lemma 1. The internal equilibrium Ef exists if R2=Sfβαμ>1 and S1<μαβ<Sf.

    Proof. The internal equilibrium Ef=(Sf,If) can be determind by

    {If=βμSfα,a0Sf3+a1Sf2+a2Sf+a3=0,

    where a0=rω, a1=rKωr+Kωβ, a2=K(βαμωr+α), a3=Kαμ.

    Define a function f(S)=a0S3+a1S2+a2S+a3. Then, Sf is a solution of f(S)=0. Since a0>0, a3<0, f(S) has at least one positive solution Sf.

    Since If=βμSfα>0, then we need Sf>μαβ. Here, the basic reproduction number R2=Sfβαμ. f(S)=3a0S2+2a1S+a2, and thus, S1 and S2 are solutions of f(S)=0. If Kω>1, then 0<S1<Sf<S2, see Figure 1. According to this, if Kω>1, r<σ, and αμβ>Kω1ω hold, then

    f(αμβ)=αμβ3(Krβ2+rαβμ+Kσβ2Krαβμω+rα2μ2ω)>αμβ3(Krβ2+rαβμ+Kσβ2Krαβμω+rαμ(Kωββ))=αμβ3(Kβ2(σr))>0.
    Figure 1.  Graph of the function f(S).

    Since S1<μαβ<Sf, the internal equilibrium Ef exists.

    Next we will perform a local stability analysis on system (3.1). For E1=(S1,0), the Jacobian matrix is

    J|E1(S1,0)=(rS1K+σωS1(1+S1ω)2S1βα0S1βαμ,)

    and

    =λ2+(rS1KσωS1(1+S1ω)2S1βα+μ)λ+S1(S1βαμ)(rKα+σωα(1+S1ω)2).

    Here, the basic reproduction number R1=S1βαμ<1 is defined.

    Theorem 3. The disease-free equilibrium E1=(S1,0) is locally asymptotically stable if R1<1, Kω<1, and r>σ.

    Proof. Clearly, if Kω<1 and r>σ hold, then

    rK+σω(1+S1ω)2<rK+σω=r+σωKK<r+rωKK=r(ωK1)K<0.

    According to this, then

    λ1+λ2=rS1K+σωS1(1+S1ω)2+S1βαμ<0,

    and

    λ1λ2=S1(S1βαμ)(rKα+σωα(1+S1ω)2)>0.

    We conclude that when the basic reproduction number R1<1, Kω<1, and r>σ, E1 is locally asymptotically stable.

    Theorem 4. If Kω>1, r<σ<Kωr, and αμβ>Kω1ω hold, Ef=(Sf,If) is locally asymptotically stable.

    Proof.

    J|Ef(Sf,If)=(rSfK+Sfσω(1+Sfω)2Sfβα(If+α)2βIfIf+αSfβα(If+α)2μ)

    and

    =λ2+(rSfKSfβα(α+If)2σωSf(1+Sfω)2+μ)λ+μ(αμ(Sfβαμ)Sf2β+(1αμSfβ)Sf(Kσω+r(1+Sfω)2)K(1+Sfω)2).

    Clearly,

    λ1+λ2=rSfK+σωSf(1+Sfω)2+Sfβα(α+If)2μ=rSfKμ+αμ2Sfβ+σωSf(1+Sfω)2<rSfKμ+μ+σSfK2ω=Sf(σKωr)K2ω<0

    and

    λ1λ2=μ(αμ(Sfβαμ)Sf2β+(1αμSfβ)Sf(Kσω+r(1+Sfω)2)K(1+Sfω)2)>μ(αμ(Sfβαμ)Sf2β+(1αμSfβ)SfKω(σ+Kωr)K(1+Sfω)2)>0.

    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,α=37,β=2,σ=0.60137,ω=0.9, and μ=0.6, we obtain positive equilibrium points (Sf=14.5313,If=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 Sf endemic state for fire, and the figure on the right shows the If endemic state for fire.

    Figure 2.  The numerical simulations of system (3.1). (a): component Sf (Locally asymptotically stable); (b): component If (Locally asymptotically stable).

    Remark 1. In Section 2, we analyzed the global stability of the no-fire model. If R0<1, E0=(S0,0) is global asymptotically stable, and if R0>1, the internal equilibrium E=(S,I) exists and is globally asymptotically stable, and we obtain the basic reproduction number R0=Kβμα. In Section 3, we considered the local stability of the model under fire disturbance. If R1=S1βαμ<1, E1 is locally asymptotically stable. If R2=Sfβαμ>1, the internal equilibrium Ef exists and is locally asymptotically stable.

    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.

    Consider the diffusion equations

    {St=DSΔS+rS(1SK)σSωS+1βSII+α,It=DIΔI+βSII+αμI. (4.1)

    In this section, DS and DI represent diffusion coefficients, and Δ=2/S2+2/I2 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+˜S,I+˜I). Then, we can obtain the linearized system as follows:

    {˜St=d1Δ˜S+[r(12SK)σ(1+ωS)2βII+α]˜SSαβ(I+α)2˜I,˜It=d2Δ˜I+βII+α˜S+(Sαβ(I+α)2μ)˜I. (4.2)

    We denote p1=r(12SK)σ(1+ωS)2βII+α,q1=Sαβ(I+α)2,p2=βII+α,q2=Sαβ(I+α)2μ.

    Let d1=εd,d2=d, then we can get ε=d1d2>0. We rewrite Eq (4.2) as

    {˜St=εd˜S+p1˜S+q1˜I,˜It=d˜S+p2˜S+q2˜I. (4.3)

    We can then obtain the Jacobian matrix

    J=(p1εdk2π2q1p2q2k2π2d).

    Then, the characteristic equation of J is

    Dk(λ,ε)=λ2TRkλ+DETk=0,

    with

    TRk=(ε+1)dk2π2+p1+q2,DETk=εd2k4π4dk2π2(p1+q2ε)+p1q2p2q1.

    Now, we assume that

    B0{p1+q2<0,p1q2p2q1>0.

    Lemma 2. Suppose that (B0) 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,

    D0(λ,ε)=λ2TR0λ+DET0=0

    with

    TR0=p1+q2,DET0=p1q2p2q1.

    Both roots of D0(λ,ε)=0 have negative real parts when TR0<0 and DET0>0. Therefore, system (4.1) is asymptotically stable at the positive equilibrium point (S,I). Hence, the lemma is proved.

    Now we consider the existence conditions for Turing instability.

    Assume that

    (B1) 0<ε<ε1, where ε12p2q1+p1q22p22q21p1p2q1q2q22.

    (B2) 0<ε<ε2(d), where ε2(d)p1dπ2q2. ε=ε(d) decreases monotonically in d and intersects with ε=ε1 at the point d=d0. We take εB(d)=min{ε1,ε2(d)}, then

    εB(d)={ε1,0<dd0ε2(d),d>d0.

    Proof. Let x=dk2π2>0, then we rewrite DETk as

    DETk=εx2(p1+q2ε)x+p1q2p2q1.

    Then, we get that the symmetry axis is x=p1+q2ε2ε, and at this time DETk can be taken to a minimum. That is,

    DETkmin=p1q2p2q1(p1+q2ε)24ε.

    Since x>0, we have p1+q2ε>0. If we want DETkmin<0 with the condition p1+q2ε>0, we must satisfy p1q2p2q1(p1+q2ε)24ε<0. We can obtain

    0<ε<2p2q1+p1q22p22q21p1p2q1q2q22,

    when

    p1>0,q1<0,q2<p1,p2>p1q2q1.

    Let ε1=2p2q1+p1q22p22q21p1p2q1q2q22. Then, when 0<ε<ε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 DETk. Then, p1+q2ε2εxp1+q2ε+2ε. k can be taken to the minimum when x takes a value on the symmetry axis, and we get

    k2min=p1+q2ε2ε1dπ2=12dp1+q2εεπ2,

    and

    kmin=12p1+q2εdεπ2.

    By making sure that 12p1+q2εdεπ2>12, we can obtain 0<ε<p1dπ2q2. When p1>0,q1<0,q2<p1,p2>p1q2q1, let ε2(d)=p1dπ2q2, and we can find that ε=ε2(d) decreases monotonically in d and intersects with ε=ε1 at the point d=d0, where

    d0=2p2q1q2+2p1q222q2p1q1(p2q1p1q2)π2(2p2q1+p1q22q2p1q1(p2q1p1q2)).

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

    ε(k,d)=dk2π2p1+p2q1p1q2dk2π2(dk2π2q2),d>dk,

    where dk=p1q2p2q1k2π2p1. Then, DETk=0, when ε=ε(k,d). That is,

    DETk=0ε(k,d)=dk2π2p1+p2q1p1q2dk2π2(dk2π2q2).

    Then

    dεdx=p1x22p1q1x+2p1q2x+p2q1q2p1q22x2(xq2)2,

    and we can get

    d=p2q1+p1q2+p22q21p1p2q1q2p11k2π2.

    Set dM=p2q1+p1q2+p22q21p1p2q1q2p11k2π2, then dεdx=0, when d=dM. At this time, ε(k,d) reaches the maximum value ε1.

    Lemma 3. Suppose that (B0) holds, function ε=ε(k,d) has the following properties: (1) When 0<d<dM(k), ε(k,d) increases monotonically, and when d>dM(k), ε(k,d) decreases monotonically. (2) As for ki<ki+1,kiN,i=1,2,3,, there is only one root dk1,k2(dM(k2),dM(k1)) satisfying ε(k1,d)=ε(k2,d) for d>0. Furthermore,

    ε=ε(d)=ε(k,d),k(dki+1,i+2,ki,i+i),kiN,i=1,2,3,.

    In Figure 3, we characterize a graph of functions ε=ε1,ε=ε2(d), and ε=ε(k,d),d>0,k=1,2,3,. The value of parameters are selected following reference [18], and the value of p1,q1,p2,q2 are calculated from the model (4.2). In Figure 4, we present a graph of the Turing bifurcation line ε=ε(k,d),d>0.

    Figure 3.  The figure of functions ε=ε1,ε=ε2(d), and ε=ε(k,d),d>0,k=1,2,3 in (d,ε) plane. r=0.721,K=90,σ=0.601,α=37.854,β=0.318,ω=0.055,μ=0.02. (S,I)(17.252,236.374), p1=0.0120,q1=0.0028,p2=0.2740,q2=0.0172.
    Figure 4.  The Turing bifurcation line ε=ε(d)(d>0).

    Let DS=0.03,DI=0.41, r=0.721,K=90,σ=0.601,α=37.8537234,β=0.318,ω=0.055, and μ=0.02. We then obtain the positive equilibrium of system (4.2), (S,I)=(17.252,236.374). Through (B1), (B2), we obtain ε1=0.0545,

    εB(d)={0.0545,0<d0.0205,0.0120dπ2+0.0172,d>0.0205,

    and

    ε(k,d)=0.012dk2π20.0005608dk2π2(dk2π2+0.0172),d>dk.

    Then, we can find ε>ε(d), so system (4.2) is asymptotically stable at (S,I), see Figure 5. By setting k=1, we obtain d0,1=+ and d1,2=0.0250246, and select d=d1=0.12(d0,1,d1,2); thus, ε=ε(1,0.12)=0.0011360854. System (4.2) with d=0.12 undergoes a Turing bifurcation near equilibrium (17.252,236.374) at ε=0.0011360854. Figure 6 shows that system (4.2) has a Turing instability near the equilibrium point (17.252,236.374) at this time.

    Figure 5.  The numerical simulations of system (4.2). (a): component S (Locally asymptotically stable); (b): component I (Locally asymptotically stable).
    Figure 6.  The numerical simulations of system (4.2) with ε=0.00111725 and the initial condition at (17.252,236.374). (a): component S (Turing instability); (b): component I (Turing instability).

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

    {dSdt=rS(1SK)βI+αSIσSωS+1,dIdt=βI+αSIμ1IσI, (5.1)

    where μ1=μσ is the natural mortality rate of diseased trees, and σ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.

    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 σ in the system is set to be the time-dependent function σ(t), describing the time-varying fire intensity, and there exists ¯σ such that 0σ(t)¯σ. Within the fixed time [0,T] with T>0, the constraint set reads

    U={σ(t)|0σ(t)¯σ,0tT,σ(t) is Lebesgue measurable}. (5.2)

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

    J(σ)=T0(12I2(t)+12σ2(t))dt,

    with I(0)=I0,σ(0)=σ0. The optimal control problem is rewritten

    J(σ()(t))=minσ(t)UJ(σ(t)), (5.3)

    where the minimization of J(σ(t)) is subject to (5.1) and σU with U from (5.2).

    Using the method in [26], we check the existence of optimal control σ(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 σ>0 such that the integrand L(σ(t)) of the objective functional J is convex and satisfied.

    We have L(σ(t))=12(I2(t)+σ2(t)).

    Theorem 5. For the optimal control problems (5.2) and (5.3), there exists an optimal control σ such that J(σ()(t))=minσ(t)UJ(σ(t)).

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

    H=12(I2+σ2)+λ1[rS(1SK)βI+αSIσSωS+1]+λ2[βI+αSIμ1IσI],

    where λ1,λ2 are adjoint variables satisfying the following adjoint equations:

    {Hλ1|(σ,λ1,λ2,t)=˙λ1,λ1(T)=0,Hλ2|(σ,λ1,λ2,t)=˙λ2,λ2(T)=0.

    The optimality conditions are given by

    Hσ|(σ,λ1,λ2,t)=0.

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

    σ=min{¯σ,max{λ1SωS+1+λ2I,0}}, (5.4)

    where λ1,λ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.

    Hσ|(σ,λ1,λ2,t)=0 gives σλ1SωS+1λ2I=0, where the adjoint variables are satisfied. Furthermore,

    ˙λ1=HS|(σ,λ1,λ2,t)=λ1(r+2rSK+IβI+α+σ(1+Sω)2)λ2IβI+α.˙λ2=HI|(σ,λ1,λ2,t)=I+λ1(SβI+αISβ(I+α)2)+λ2(ISβ(I+α)2SβI+αμσ),

    and the transversality conditions are λ1(T)=0,λ2(T)=0.

    Moreover, since σU, using the lower and upper bounds of ¯σ(t), the optimal σ(t) can be characterized by (5.4). The curves of σ(t), S(t), and I(t) with time for different ω are given in Figure 7.

    Figure 7.  The optimal control σ(t), and the corresponding S(t), I(t). Here, r=0.62117,K=85,α=37,β=2,μ=0.6,¯σ=0.9.

    For the improved model (5.1), in the transient dynamic process, with ω=0, σ=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 (ω=0.054879880). This region has a higher I. When ω=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 (ω=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.

    Table 1.  Annual forest standing tree biomass data for forest fires in each province (unit: m3/ha).
    Year 2002 2003 2008 2013 2018
    Shanghai 10.97706 17.5873 16.9095 27.3642 50.5157
    Inner Mongoria 66.5583 53.7157 49.7467 54.0739 58.3988
    Yunnan 99.7149 89.69645 85.4803 88.4495 93.6613
    Hebei 17.6961 19.7972 20.0179 24.5259 27.3289
    Jilin 111.2575 113.3776 114.6019 120.7763 129.0606
    Jiangsu 18.7235 32.5714 39.9136 39.9136 45.1598

     | Show Table
    DownLoad: CSV

    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.

    Figure 8.  Trends in live forest tree biomass over time at different fire frequencies, with the number of fires per year as units on the x-axis. The grey lines are linear regressions performed within each survey year. Colors represent fire frequency groupings.
    Figure 9.  Longitudinal trends within the plots in the different treatments with respect to fire frequency. Raw data is expressed as points and the solid smoothed line is a generalized additive model fit within the plot. (a) Number of forests in Shanghai. (b) In order, number of forests in Yunnan, Inner Mongoria and Shanxi. (c) In order, number of forests in Guangdong and Hebei. (d) In order, number of forests in Jilin and Henan. (e) Number of forests in Jangsu.

    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.37m3/ha, and through the model, we predicted that Jilin's tree biomass in 2020 would be 132.395m3/ha, which illustrates the reasonableness of the model.

    Figure 10.  Predicted tendency of forest fire in Shanghai, Inner Mongolia, Yunnan, Hebei, Jilin, etc. (a) Number of forests in Shanghai. (b) Number of forests in Inner Mongoria. (c) Number of forests in Yunnan. (d) Number of forests in Hebei. (e) Number of forests in Jilin. (f) Number of forests in Jiangsu. See Table 1 for specific data.

    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 σ 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 ˙y(t)=f(t,y,m), and the absolute sensitivity of variable yi to parameter M:

    Pi(t)=yi(t.M)M,i=1,2.

    We denote the relative sensitivity of a variable to a parameter

    pi(t)=yi(t.M)MMyi,i=1,2

    and the absolute sensitivity equation of parameter Pi

    ˙Pi=fyPi+fm,i=1,2.

    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:

    {˙S=rS(1SK)βI+αSIσSωS+1,˙I=βI+αSIμ1IσI,˙P1=P1(r2rSkβII+ασ(ωS+1)2)P2Sβα(I+α)2SωS+1,˙P2=P1βII+α+P2(Sβα(I+α)2μ1)I. (5.5)

    Assuming the values of time step t=2000 and fire frequency σ=0.6037, the initial values are chosen to be consistent with Section 3, and the sensitivity analysis of σ 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 σ on each component (susceptible and infected trees) in the forest, as shown in Figure 11.

    Figure 11.  Relative sensitivity and sensitivity of σ to S(t) and I(t).

    As shown in Figure 11, the sensitivity P1 and relative sensitivity σP1/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 σ, and when t=0.5, P1=1.606×104<0, σP1S=1.531×107<0. That means that when σ increases by 10%, S will decrease by 1.531×106%, and the sensitivity P2 and relative sensitivity σP2/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 σ, and when t=0.5, P2=3.342<0, σP2S=0.026<0, it means that when σ increases by 10%, I will decrease by 0.260%.

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

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

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

    The authors declare no competing interests.



    [1] N. G. Mcdowell, C. D. Allen, A. T. Kristina, Pervasive shifts in forest dynamics in a changing world, Science, 368 (2020), 964–976. http://dx.doi.org/10.1126/science.aaz9463 doi: 10.1126/science.aaz9463
    [2] R. Seidl, D. Thom, M. Kautz, D. Martin-Benito, Forest disturbances under climate change, Nat. Clim. Chang., 7 (2017), 395–402. http://dx.doi.org/10.1038/NCLIMATE3303 doi: 10.1038/NCLIMATE3303
    [3] J. A. Hicke, C. D. Allen, A. R. Desai, Effects of biotic disturbances on forest carbon cycling in the United States and Canada, Glob. Change Biol., 18 (2012), 7–34. http://dx.doi.org/10.1111/j.1365-2486.2011.02543.x doi: 10.1111/j.1365-2486.2011.02543.x
    [4] I. L. Boyd, P. H. Freersmith, C. A. Gilligan, H. C. J. Godfray, The consequence of tree pests and diseases for ecosystem services, Science, 342 (2013), 823–833. http://dx.doi.org/10.1126/science.1235773 doi: 10.1126/science.1235773
    [5] A. J. Tepley, E. Thomann, T. T. Veblen, Influences of fire-vegetation feedbacks and post-fire recovery rates on forest landscape vulnerability to altered fire regimes, J. Ecol., 106 (2018), 1925–1940. http://dx.doi.org/10.1111/1365-2745.12950 doi: 10.1111/1365-2745.12950
    [6] A. C. Staver, S. Archibald, S. A. Levin, The global extent and determinants of savanna and forest as alternative biome states, Science, 334 (2011), 230–232. http://dx.doi.org/10.1126/science.1210465 doi: 10.1126/science.1210465
    [7] S. I. Higgins, W. J. Bond, W. S. W. Trollope, Fire, resprouting and variability: a recipe for grass-tree coexistence in savanna, J. Ecol., 88 (2000), 213–229. http://dx.doi.org/10.1046/j.1365-2745.2000.00435.x doi: 10.1046/j.1365-2745.2000.00435.x
    [8] W. A. Hoffmann, E. L. Geiger, S. G. Gotsch, D. R. Rossatto, L. C. R. Silva, O. L. Lau, Ecological thresholds at the savanna forest boundary: how plant traits, resources and fire govern the distribution of tropical biomes, Ecol. Lett., 15 (2012), 759–768. http://dx.doi.org/10.1111/j.1461-0248.2012.01789.x doi: 10.1111/j.1461-0248.2012.01789.x
    [9] D. W. Peterson, P. B. Reich, K. J. Wrage, J. Franklin, Plant functional group responses to fire frequency and tree canopy cover gradients in oak savannas and woodlands, J. Veg. Sci., 18 (2007), 3–12. http://dx.doi.org/10.1111/j.1654-1103.2007.tb02510.x doi: 10.1111/j.1654-1103.2007.tb02510.x
    [10] R. K. Meentemeyer, N. J. Cunniffe, A. R. Cook, J. A. Joao, R. D. Hunter, D. M. Rizzo, Epidemiological modeling of invasion in heterogeneous landscapes: spread of sudden oak death in California (1990–2030), Ecosphere, 2 (2011), 1–24. http://dx.doi.org/10.1890/ES10-00192.1 doi: 10.1890/ES10-00192.1
    [11] A. Dobson, M. Crawley, Pathogens and the structure of plant communities, Trends Ecol. Evol., 9 (1994), 393–398. http://dx.doi.org/10.1016/0169-5347(94)90062-0 doi: 10.1016/0169-5347(94)90062-0
    [12] R. M. Anderson, R. May, Infectious Diseases of Humans: Dynamics and Control, Oxford: Oxford University Press, 1991. http://dx.doi.org/10.1002/hep.1840150131
    [13] J. J. Burdon, G. A. Chilvers, Host density as a factor in plant disease ecology, Annu. Rev. Phytopathol, 20 (1982), 143–166. http://dx.doi.org/10.1146/annurev.py.20.090182.001043 doi: 10.1146/annurev.py.20.090182.001043
    [14] L. S. Comita, S. A. Queenborough, S. J. Murphy, Testing predictions of the Janzen-Connell hypothesis: a meta-analysis of experimental evidence for distance-and density-dependent seed and seedling survival, J. Ecol., 102 (2014), 845–856. http://dx.doi.org/10.1111/1365-2745.12232 doi: 10.1111/1365-2745.12232
    [15] J. P. Gerlach, P. B. Reich, K. Puettmann, T. Baker, Species, diversity, and density affect tree seedling mortality from Armillaria root rot, Can. J. For. Res., 27 (1997), 1509–1512. http://dx.doi.org/10.1139/x97-098 doi: 10.1139/x97-098
    [16] D. W. Peterson, P. B. Reich, Prescribed fire in oak savanna: fire frequency effects on stand structure and dynamics, Ecol. Appl., 11 (2001), 914–927. http://dx.doi.org/10.1890/1051-0761(2001)011[0914:PFIOSF]2.0.CO;2 doi: 10.1890/1051-0761(2001)011[0914:PFIOSF]2.0.CO;2
    [17] L. E. Frelich, P. B. Reich, D. W. Peterson, Fire in Upper Midwestern Oak Forest Ecosystems: An Oak Forest Restoration and Management Handbook, Portland: Pacific Northwest Research Station, 2015. http://dx.doi.org/10.2737/PNW-GTR-914
    [18] A. F. A. Pellegrini, A. M. Hein, J. Cavender-Bares, R. A. Montgomery, Disease and fire interact to influence transitions between savanna-forest ecosystems over a multi-decadal experiment, Ecol. Lett., 24 (2021), 1007–1017. http://dx.doi.org/10.1111/ele.13719 doi: 10.1111/ele.13719
    [19] P. Magal, Z. Y. Zhang, A system of state-dependent delay differential equation modelling forest growth Ⅱ: Boundedness of solutions, Nonlinear Anal. Real World Appl., 42 (2018), 334–352. http://dx.doi.org/10.1016/j.nonrwa.2018.01.002 doi: 10.1016/j.nonrwa.2018.01.002
    [20] P. Magal, Z. Y. Zhang, A system of state-dependent delay differential equation modeling forest growth Ⅰ: semiflow properties, J. Evol. Equ., 18 (2018), 1853–1888. http://dx.doi.org/10.1007/s00028-018-0464-0 doi: 10.1007/s00028-018-0464-0
    [21] P. Magal, Z. Y. Zhang, Numerical simulations of a population dynamic model describing parasite destruction in a wild type pine forest, Ecol. Complex., 34 (2018), 147–160. http://dx.doi.org/10.1016/j.ecocom.2017.05.001 doi: 10.1016/j.ecocom.2017.05.001
    [22] P. Magal, Z. Y. Zhang, Competition for light in forest population dynamics: from computer simulator to mathematical model, J. Evol. Equ., 419 (2017), 290–304. http://dx.doi.org/10.1016/j.jtbi.2017.02.025 doi: 10.1016/j.jtbi.2017.02.025
    [23] J. P. Tripathi, S. Tyagi, S. Abbas, Global analysis of a delayed density dependent predator-prey model with Crowley-Martin functional response, Commun. Nonlinear Sci. Numer. Simul., 30 (2015), 45–69. http://dx.doi.org/10.1016/j.cnsns.2015.06.008 doi: 10.1016/j.cnsns.2015.06.008
    [24] H. Chen, C. Zhang, Dynamic analysis of a Leslie Gower type predator prey system with the fear effect and ratio-dependent Holling Ⅲ functional response, Nonlinear Anal. Model Control, 27 (2022), 904–926. http://dx.doi.org/10.15388/namc.2022.27.27932 doi: 10.15388/namc.2022.27.27932
    [25] X. M. Rong, M. Fan, H. P. Zhu, Dynamic modeling and optimal control of cystic echinococcosis, Infect. Dis. Poverty, 10 (2021), 38. http://dx.doi.org/10.1186/s40249-021-00807-6 doi: 10.1186/s40249-021-00807-6
    [26] L. H. Zhou, M. Fan, Q. Hou, Z. Jin, X. D. Sun, Transmission dynamics and optimal control of brucellosis in Inner Mongolia of China, Math. Biosci. Eng., 15 (2018), 543–567. http://dx.doi.org/10.3934/mbe.2018025 doi: 10.3934/mbe.2018025
    [27] R. H. Martin, Logarithmic norms and projections applied to linear differential systems, J. Math. Anal. Appl., 45 (1974), 432–454. http://dx.doi.org/10.1016/0022-247X(74)90084-5 doi: 10.1016/0022-247X(74)90084-5
    [28] X. X. Liu, C. R. Zhang, Stability and optimal control of Tree-Insect Model under forest fire disturbance, Mathematics, 10 (2022), 2563. http://dx.doi.org/10.3390/math10152563 doi: 10.3390/math10152563
  • Reader Comments
  • © 2024 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(1351) PDF downloads(85) Cited by(0)

Figures and Tables

Figures(11)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog