Analysis of dynamic properties on forest restoration-population pressure model

1.   Introduction

Forest is the main part of the whole terrestrial ecosystem. The generation and distribution of interspecies models play an important role in ecological protection and biochemical response. One of the typical models is forest restoration model ^{[1,2,3,4,5,6]}. Many mathematicians and biologists have been paying close attention to these kinds of models since the 1960s.

Okubo ^[7] defined the process of diffusion as a set of particles, which may be a molecule or a living, through the individual's random movement in space and time. The diffusion models had been widely used in the fields of chemistry, physics and ecology ^[8]. Diffusion models had also been successfully applied to describe land cover change. Jesse ^[9] and Svirezhev ^[10] indicated that when combining reaction terms, these diffusion models became particularly suitable for simulating land use and land cover change because they took into account the dynamics of spatial structure and land cover class. Acevedo et al. ^[11] described a quantitative model of forest restoration based on diffusion logistic growth (DLG), which provided a new method for study of forest restoration. While considering spatial dependence, the relationship between continuous space and time of forest restoration is studied. Holmes et al. ^[12] pointed out that application of DLG model in forest restoration assumed land-cover change could be described as traveling waves diffusing from source to outside at constant rate. It is assumed that tree colonization can be regarded as a small scale forest restoration process and spread irregularly in space.

Recently, land-use and land-cover change (LUCC) had a wide range of impacts on biodiversity, climate and ecosystem services. In Vitousek, Nunes et al. and Houet et al. ^[13,14,15], the authors believed that LUCC was the most vital human change on earth. Pereira et al. ^[16] said it would straightly affect biodiversity. In ^[17,18], the authors thought the climate was also affected by LUCC. Lambin et al. ^[19] noted that although LUCC described usual land cover change, most studies payed a lot of attentions to large-scale deforestation processes used by humans, such as agriculture. Forest restoration is becoming more and more common in agricultural intensification processes, especially in the tropics ^[20,21,22].

Human activity is one of the main factors leading to change of natural environment. The contradiction between population and forest resources is becoming more and more prominent, so it is important to establish forest restoration-population pressure model to consider the balance between population growth and forest resources.

The objective of this paper is to analyse dynamic properties of forest restoration-population pressure model. In the second section, we establish the population pressure reaction diffusion model under Dirichlet boundary conditions according to the idea ^[23,24]. And we study the local stability of positive equilibrium of the model in the third section. In the fourth section, the existence conditions of Turing instability and Hopf bifurcation are given in Theorem 4.1 and 4.2, respectively. We obtain two kinds of steady-state bifurcation types, that is, simple bifurcation and double bifurcation. At the same time, we select some values for the parameters, and get numerical simulations to support the results in section $4$.

2.   Construction of forest restoration-population pressure model

In the classical applications of partial differential equation models to population ecology, organisms are assumed to have Brownian random dispersion, the rate of which is invariant in time and space

where $u\left(x, y, t\right)$ is the density of organisms at spatial coordinates $x$, $y$, and time $t$, and $D$ is the diffusion coefficient that measures dispersal rate.

Holmes et al. ^[12] pointed out that the diffusion logistic growth (DLG) model is an extension of two-dimensional Fisher Eq ^[25]. DLG consists of two parts: Brownian random dispersion and logistic population growth. This model can be represented by a partial differential equation:

where $u$ stands for forest radio, $t$ for time, ${D_u}$ for diffusion rate, $x$ and $y$ for spatial position, ${r_u}$ for natural growth rate, ${K_u}$ for bearing capacity.

In order to complete the model, we must specify the initial values and boundary conditions. The forest ratio in region $R = \left\{ {(x, y)\left| {0 \le x \le {L_x}, 0 \le y \le {L_y}} \right.} \right\}$ is defined as initial value at time $t = 0$, where ${L_x}$ and ${L_y}$ represent landscape range of the $x$ and $y$ axes. At the same time, we apply Dirichlet boundary $u = 0$ on the boundary of $R$, $\partial R$:

The current application of the model includes three parameters: diffusion rate (${D_u}$), endogenous growth rate (${r_u}$) and carrying capacity (${K_u}$). The diffusion rate describes rate of change of land cover radio. In the case of land cover change, intrinsic growth rate describes the rate of change of land cover category. Proportion of this land cover category can be increased to the maximum, that is, carrying capacity (${K_u} = 1$).

With regard to influence of population on forest resources, we use predator-prey systems for reference and describe the relationship between population and forest as Holling Ⅱ functional responses ^[26]. According to ^[27] and based on Eq (2.1), forest restoration model is considered under population pressure, and the following model is established:

where $v$ is the population density, $a$ is the half-saturation constant, $\alpha = c \beta$, $c$ is the conversion coefficiency of the population. The term ${{\beta u}}/\left({{1 + au}}\right)$ denotes the functional response of the population, all parameters are assumed to be positive.

We specify the value of solution of the differential equation at the boundary, then for all $t \ge 0$, the system complies with homogeneous Dirichlet boundary conditions

where $u$ and $v$ are state variables, $\left({{u_0}, {v_0}} \right)$ is a uniform steady-state solution which is independent of variables, $t$, $x$, $y$ and satisfies

For convenience of calculation, we have recorded ${D_u}$, ${D_N}$, ${r_u}$, ${r_v}$ and $1/{K_u}$, $1/{K_v}$ as ${d_1}$, ${d_2}$, ${r_1}$, ${r_2}$ and $P$, $Q$, respectively. Therefore, we can obtain

3.   Linear stability analysis

Consider

from the first equation of Eq (3.1), we can obtain

Substituting $v$ in the second equation of Eq (3.1) with Eq (3.2), we obtain

where

In the following discussion, we suppose that Eq (2.3) has a positive equilibrium $\left({{u_*}, {v_*}} \right)$.

In order to simplify discussions, we transform homogeneous state solution $\left({{u_*}, {v_*}} \right)$ into $\left({0, 0} \right)$ by the means of transformation $\left({u, v} \right) = \left({{u_*} + \tilde u, {v_*} + \tilde v} \right)$. We obtain

here and below, we call

Rewrite Eq (3.4) as follows

Moreover, we explicitly incorporate length $l$ into equations by transforming $x = l\tilde x$ and $y = l\tilde y$, which changes the domain $\Omega = \left[{0, l} \right] \times \left[{0, l} \right]$ to the unit square $\tilde \Omega = \left[{0, 1} \right] \times \left[{0, 1} \right]$, then Eq (2.3) become

with

For the sake of simplify and calculation, we express $\tilde \Omega$, $\tilde x$, $\tilde y$, and $\tilde u$, $\tilde v$ as $\Omega$, $x$, $y$, and $u$, $v$.

Let

and

we adapt Eq (3.5) as an operator equation

here $U: = \left({u, v} \right)$, and map $\Phi :X \to Y$ is defined as

By differentiating $U$ on homogeneous equilibrium solution ${U_*} \equiv \left({{u_*}, {v_*}} \right) \equiv \left({0, 0} \right)$, we have the linearization of mapping $\Phi$ into $L$

The stability of $U$ can be analyzed by the solution of variational problem

we decompose X in direct sum

and $L$ maps ${X_{m, n}}$ to itself. Therefore,

The limitation of $L$ in subspace ${X_{m, n}}$ is a matrix of $2 \times 2$

Therefore,

4.   Analysis of Turing instability and Hopf bifurcation

In this section, we mainly study the existence of Turing instability and Hopf bifurcation. We define a wave number $k$ with ${k^2} = {{m^2} + {n^2}} \in {\mathbb {N}}$, and denote $\theta = {d_2}/{d_1}$, ${d_1} = d$, then the characteristic equation of ${A_{m, n}}$ is as follows:

with

We assume that

$\left({{{\rm{A}}_0}} \right) - {{b_1}{b_2}}/{{a_2}} < {a_1} < - {a_2}.$

$\left({{{\rm{A}}_1}} \right) 0 < \theta < {\theta _1}, {\theta _1} \buildrel \Delta \over = {\frac{a_1 a_2+2 b_1 b_2}{a_1^2}-2 \sqrt{\frac{a_1 a_2 b_1 b_2+b_1^2 b_2^2}{a_1^4}}}$.

Under the assumption $\left({{{\rm{A}}_0}} \right)$, all eigenvalues of ${\Lambda_0}\left({\lambda, \theta} \right)$ have negative real parts, and ${\rm{Trace}}\left(k \right) < 0$, for $k \in {\mathbb {N}}$.

Lemma 4.1. Suppose that $\left({{{\rm{A}}_1}} \right)$ holds, then $\mathop {min}\limits_{k \in {\mathbb {R}_ + }} {\rm{Det}}\left(k \right) < 0$.

Proof. Let $x = {{d{k^2}{\pi ^2}}}/{{{l^2}}} > 0$, then we rewrite ${\rm{Det}}\left(k \right)$ as

From function (4.2), when $x = \left({{{a_2} + {a_1}\theta}}\right)/{{2\theta }}$, ${\rm{Det}}\left(k \right)$ can be taken to a minimum, that is, ${\rm{Det}}{\left(k \right)_{min}} = {a_1}{a_2} + {b_1}{b_2} - {{{{\left({{a_2} + {a_1}\theta } \right)}^2}}}/{{4{\theta }}}$. Since $x > 0$, we have ${a_2} + {a_1}\theta > 0$. If we want ${\rm{Det}}{\left(k \right)_{min}} < 0$ with the condition ${a_2} + {a_1}\theta > 0$, $\theta$ must satisfy $0 < \theta < \left({a_1 a_2+2 b_1 b_2}\right)/{a_1^2}-2 \sqrt{\left({a_1 a_2 b_1 b_2+b_1^2 b_2^2}\right)/{a_1^4}}$. Hence, the lemma is proved.

Let $k_{min}^2$ be the minimal point of function ${\rm{Det}}\left(k \right)$ on ${k^2} \in {\mathbb {R}_ + }$, then

We make the following assumption to ensure $k_{min} > 0$

$\left({{{\rm{A}}_2}} \right) 0 < \theta < {\theta _2}\left(d \right), {\theta _2}\left(d \right) \buildrel \Delta \over = {\frac{a_2 l^2}{\pi ^2 d - a_1 l^2}}$.

$\theta = {\theta _2}\left(d \right)$ decreases monotonically in $d$ and intersects with $\theta = {\theta _1}$ at the point $d = {d_0}$. Let ${\theta _{\rm{A}}}\left(d \right) = \mathop {min}\limits_{d > 0} \left\{ {{\theta _1}, {\theta _2}\left(d \right)} \right\}$, then

Hence, we have the following lemma.

Lemma 4.2. Suppose that $\left({{{\rm{A}}_0}} \right)$ holds, then assumptions $\left({{{\rm{A}}_1}} \right)$ and $\left({{{\rm{A}}_2}} \right)$ hold if and only if $0 < \theta < {\theta _{\rm{A}}}\left(d \right)$, $d > 0$.

4.1.   Existence of Turing instability

4.1.1.   Stability analysis of Turing instability

Denote

where ${d_k} = {\frac{{\left({a_1}{a_2}+ {b_1}{b_2} \right){l^2}}}{{{a_2} {k^2}{\pi ^2}}}}$, then ${\rm{Det}}\left(k \right) = 0$ when $\theta = {\theta _{\rm{T}}}\left({k, d} \right)$.

Let ${d_M}\left(k \right)$ be the point at which monotonicity of function changes, that is, function $\theta = {\theta _{\rm{T}}}\left({k, d} \right)$ increases monotonically if ${d_k} < d < {d_M}\left(k \right)$, and $\theta = {\theta _{\rm{T}}}\left({k, d} \right)$ decreases monotonically if ${d_M}\left(k \right) < d < + \infty$. Hence, $\theta = {\theta _{\rm{T}}}\left({k, d} \right)$ can take the maximum value $\theta_1$ at ${d_M}\left(k \right)$.

Lemma 4.3. Suppose that $\left({{{\rm{A}}_0}} \right)$ holds, function $\theta = {\theta _{\rm{T}}}\left({k, d} \right)$ has the following properties.

(i) As for $k_i < k_{i+1}$, $k_i \in \mathbb {N}, i = 1, 2, 3\cdots$, there is only one root ${d_{k_1, k_2}} \in \left({{d_M}\left({k_2} \right), {d_M}\left(k_1 \right)} \right)$ satisfies ${\theta _{\rm{T}}}\left({k_1, d} \right) = {\theta _{\rm{T}}}\left({k_2, d} \right)$ for $d > 0$. Furthermore,

(ii) Define ${d_{0, k_1}} = + \infty$, and

Then

Moreover, ${\theta _{\rm{T}}}\left(d \right) = {\theta _{\rm{A}}}\left(d \right)$ if and only if $d = {d_M}\left(k \right)$, $k \in \mathbb {N}$.

We display the relation between $\theta = {\theta _1}$, $\theta = {\theta _2}\left(d \right)$ and $\theta = {\theta _{\rm{T}}}\left({k, d} \right)$ in Figure 1, where $d > 0$, $k = 5, 10, 13 \cdots$.

Theorem 4.1. Suppose that $\left({{{\rm{A}}_0}} \right)$ holds.

(1) For any given ${k_1} \in \mathbb {N}$, when $\theta = {\theta _{\rm{T}}}\left({{k_1}, d} \right)$, system (3.5) occurs Turing bifurcation at $\left({{u_*}, {v_*}} \right)$.

(2) $\theta = {\theta _{\rm{T}}}\left(d \right)$, $d > 0$ is the critical curve of Turing instability.

(i) If $\theta > {\theta _{\rm{T}}}\left(d \right)$, $d > 0$, system (3.5) is asymptotically stable at $\left({{u_*}, {v_*}} \right)$.

(ii) If $0 < \theta < {\theta _{\rm{T}}}\left(d \right)$, $d > 0$, Turing instability occurs in system (3.5) at $\left({{u_*}, {v_*}} \right)$.

Proof. When $\theta = {\theta _{\rm{T}}}\left({{k_1}, d} \right)$, we have ${\rm{Det}}\left({{k_1}} \right) = 0$. Then Eq (4.1) becomes

The Eq (4.5) has a zero root, and the other root of ${\Lambda_{{k_1}}}\left({\lambda, \theta } \right)$ has negative real part. That is, Eq (3.5) occurs Turing bifurcation at $\left({{u_*}, {v_*}} \right)$ when $\theta = {\theta _{\rm{T}}}\left({{k_1}, d} \right)$.

When $\theta > {\theta _{\rm{T}}}\left(d \right)$, $d > 0$, ${\rm{Det}}{\left({{k_1}} \right)} > 0$. $\left({{{\rm{A}}_0}} \right)$ ensures that ${\rm{Trace}}\left({{k_1}} \right) < 0$, then all roots of ${\Lambda_{{k_1}}}\left({\lambda, \theta } \right)$ have negative real parts. On the other hand, when $0 < \theta < {\theta _{\rm{T}}}\left(d \right)$, $d > 0$, ${\rm{Det}}{\left({{k_1}} \right)} < 0$, so system (3.5) is Turing instability.

Remark 4.1. When $\theta > {\theta _{\rm{T}}}\left(d \right)$, $d > 0$, system (3.5) occurs steady-state bifurcation at positive equilibrium $\left({{u_*}, {v_*}} \right)$.

(1) If $m = n$, the bifurcation type is a simple steady-state bifurcation.

(2) If $m \ne n$, the bifurcation type is a double steady-state bifurcation.

4.1.2.   Numerical simulations of Turing bifurcation

Let ${r_1} = 1$, ${r_2} = 0.01$, $P = 10$, $Q = 20$, $a = 0.01$, $\alpha = 8$, $\beta = 11$, $l = 3$, we obtain positive equilibrium of system (2.3) is

and ${a_1} = -0.0996$, ${a_2} = 0.0573$, ${b_1} = 0.11$, ${b_2} = 0.6544$. Through $\left({{{\rm{A}}_1}} \right)$, $\left({{{\rm{A}}_2}} \right)$, Eqs (4.3, 4.4), we obtain ${\theta_1} = 0.0119$,

and

By setting ${k_1} = 5$, ${m_1} = 3$, ${n_1} = 4$, we obtain ${d_{0, 5}} = +\infty$ and ${d_{5, 10}} = 0.0535$. Select $d = {d_1} = 0.06 \in \left[{{d_{5, 10}}, {d_{0, 5}}} \right)$, thus ${\theta _{\rm{T}}} = {\theta _{\rm{T}}}\left({5, 0.06} \right) = 0.0097$. Equation (3.5) with $d = 0.06$ undergoes Turing bifurcation near equilibrium $\left({0.01, 0.0818} \right)$ at $\theta = 0.0097$. Since ${m_1} \ne {n_1}$, that is, the bifurcation type of Eq (3.5) is double steady-state bifurcation, and we derive the following mixed patterns.

4.2.   Existence of Hopf bifurcation

4.2.1.   Stability analysis of Hopf bifurcation

We assume that

${\left({{{\rm{A}}_0}} \right)}^{'} {a_1} > max\{ -{a_2}, - {{b_1}{b_2}}/{{a_2}}\}.$

Denote

where ${d_k} = {\frac{\left({a_1} + {a_2}\right){l^2}}{{k^2}{\pi ^2}}}$, then ${\rm{Trace}}\left(k \right) = 0$ when $\theta = {\theta _{\rm{H}}}\left({k, d} \right)$.

$\theta = {\theta _{\rm{H}}}\left({k, d} \right)$ decreases monotonically in $d$ and intersects with $\theta = {\theta _1}$ at the point $d = {d_{\rm{H}}}$. On the basis of Lemma 4.1, it is known that $\mathop {min}\limits_{k \in {\mathbb {R}_ + }} {\rm{Det}}\left(k \right) < 0$ if $0 < \theta < {\theta _1}$, then we obtain the following theorem.

Theorem 4.2. Suppose that ${\left({{{\rm{A}}_0}} \right)}^{'}$ holds. For any given ${k_1} \in \mathbb {N}$, when $\theta = {\theta _{\rm{H}}}\left({{k_1}, d} \right) > {\theta_1}$, $0 < d < {d_{\rm{H}}}$, system (3.5) occurs Hopf bifurcation at positive equilibrium $\left({{u_*}, {v_*}} \right)$.

Proof. $\theta = {\theta _{\rm{H}}}\left({{k_1}, d} \right)$, we have ${\rm{Trace}}\left({{k_1}} \right) = 0$, then Eq (4.1) becomes

Since $\theta > {\theta_1}$, through assumption ${\left({{{\rm{A}}_0}} \right)}^{'}$, ${\rm{Det}}\left({{k_1}} \right) > 0$, we know that Eq (4.7) has a pair of pure imaginary roots. The above satisfies conditions for the generation of Hopf bifurcation.

4.2.2.   Numerical simulations of Hopf bifurcation

Let ${r_1} = 1$, ${r_2} = 0.01$, $P = 1$, $Q = 10$, $a = 0.01$, $\alpha = 8$, $\beta = 8$, $l = 3$, we can solve positive equilibrium of Eq (2.3) to obtain

We then obtain ${a_1} = -0.0059$, ${a_2} = 0.009$, ${b_1} = 0.024$, ${b_2} = 1$. Through assumption $\left({{{\rm{A}}_1}} \right)$ and (4.6), we have ${\theta_1} = 0.0008$, and

By setting ${m_3} = 5$, ${n_3} = 12$, we obtain ${d_{\rm{H}}} = 1.6714\times{10^{-5}}$. We choose $d = 1.5206\times{10^{-5}}$ which is satisfied $d < {d_{\rm{H}}}$, and obtain ${\theta_{\rm{H}}} = 0.1$. The relationship between $\theta = {\theta_1}$ and $\theta = {\theta_{\rm{H}}}\left({k, d} \right)$ can be reflected by Figure 5.

Hopf bifurcation occurs at $\theta = 0.1$, we can obtain hyperhexagonal patterns in Figures 6 and 7.

Moreover, from Figures 8 and 9, we observe that $u\left({x, t} \right)$ and $v\left({x, t} \right)$ are periodic in relation to $t$, that is, there is a stable bifurcation periodic solution near the positive equilibrium of Eq (2.3).

5.   Conclusions

DLG system to model forest recovery as the spread of land-cover classes in continuous space and time. However, there is only one variable $u$ in DLG system. In this paper, we establish a reaction diffusion model based on DLG system and Holling Ⅱ functional responses, we can use this model to reflect the variable relationship between forest restoration and population pressure. When the population pressure increases, the density of forest restoration will decrease. On the contrary, when the population pressure decreases, the density of forest restoration will increase and the situation of forest restoration will be improved. We also consider the effect of population pressure on forest restoration, through Eq (2.3).

Under Dirichlet boundary conditions, we obtain a positive equilibrium. According to linear stability analysis of positive equilibrium solution of Eq (3.5) and discussion of eigenvalue attributes of the matrix, we obtain the conditions for Turing instability and Hopf bifurcation. Furthermore, the steady-state bifurcation is divided into simple steady-state bifurcation and double steady-state bifurcation. In addition, we verify the existence of Turing instability and Hopf bifurcation by numerical simulation, and obtain different kinds of patterns.

Our work is to further study the diffusion logistic growth (DLG) model, which helps ease the relationship between forest restoration and population pressure.

Acknowledgments

This research is supported by the Fundamental Research Funds for the Central Universities (2572017BB03) and Science and Technology Innovation Talent Research Fund of Harbin (2017RAQXJ108). The authors wish to express their gratitude to the editors and the reviewers for the helpful comments.

Conflict of Interests

The authors declare there is no conflict of interest