Catastrophe control of aphid populations model

1.   Introduction

Though the application of catastrophe theory experiences a period of ups and downs ^[1,2,3], catastrophe phenomena such as pest outbreak, draw some ecologists and mathematicians' attention to catastrophe theory ^[4,5,6,7].

In the 50-acre experiment plot near Northwest A and F University, Zhao Huiyan et al. investigated the open system of the wheat aphid populations for 50 days in 1987–1988. The following phenomena are shown ^[8]: $(1)$ The wheat aphid populations shows distinct difference in every growth stages of wheat, for example, the amount starts to increase in wheat joining stage and goes up faster in the filling stage. $(2)$ The amount of the wheat aphid populations shows hysteresis after spraying pesticides; that is, the wheat aphid populations suddenly decrease, and then increase and exceed the amount before. $(3)$ When the density of the wheat aphid populations is low, the individuals don't disturb each other; that is, there is no correlation between the individual growth rate and the corresponding density. However, a relationship appears gradually with the increase of the wheat aphid populations, and there exists a threshold below which the corresponding individual growth rate increases and above which it decreases. So far, Zhao Huiyan et al. have been devoting to the research above. According to data in the investigations, some catastrophe models are established to successfully predict the outbreak of the wheat aphid populations in Guanzhong and Weihe areas of Shanxi Province ^[9]. On the other hand, they employ some catastrophe models to study the ecological mechanism behind them, such as the fold catastrophe model ^[8], the cusp catastrophe model ^[9,10], the swallowtail catastrophe model ^[11,12], the butterfly catastrophe model ^[13] and the elliptic catastrophe model ^[14]. On the basis of the above, an APHIDSim Software is developed to simulate the dynamics of the wheat aphid populations ^[15]. The research above not only realizes the first application of catastrophe theory in ecosystem but also promotes the real combination of ecology and mathematics.

In 1990s, Ma Zhanshan et al. conducted the laboratory experiments of Russian wheat aphid populations in 1994–1995 to investigate the survival, development, reproduction of 800 Russian wheat aphid populations. The following results are got ^{[16,17,18,19]}: $(1)$ Cohort life tables, reproductive heterogeneity tables and reproductive schedule tables are constructed. $(2)$ Some demographic models are proposed to simulate the growth of the populations. $(3)$ The corresponding cusp catastrophe model is obtained by choosing intrinsic growth rate as the state variable of the system, temperature and host plant-growth stage as control variables. The study offers the theoretical support for pest management.

On the basis of Zhao Huiyan and Ma Zhanshan' study, Zhao Lichun et al. apply singular system theory and control theory to wheat field ecosystems, the main results are as followings ^{[20,21,22,23,24]}: $(1)$ For the chemical management of the wheat aphid populations, the cusp catastrophe model with impulsive effects is built and the corresponding controllers are designed not only to minimize the management cost but also to keep the wheat aphid populations at the refuge level. $(2)$ For the swallowtail catastrophe model and the butterfly catastrophe model, the qualitative analysis is studied to explain the outbreak phenomenon of the wheat aphid populations. $(3)$ Based on cusp catastrophe model, the corresponding singular model is established to reveal the catastrophe mechanism of the wheat aphid populations. In a word, the study above not only expands the application field of singular system theory but also offers a new application background for catastrophe control.

It is well known that aphid populations seriously damage wheat by feeding on the different parts of the crop and by transporting some disease causative agents. Hence it is very crucial to manage them in a sustainable way, which belongs to the scope of catastrophe control. In the early 2000s, considering the similarity between the transfer function in automatic control principle and the potential function in catastrophe theory, Sun Rao et al. introduced the catastrophe control into the quantitative research and description of control systems. The catastrophe control is applied in the artificial heart and the integrated navigation system for the first time ^[25]. With development of catastrophe control theory, it has been extensively used in some engineering systems ^{[26,27,28,29,30,31]} and network systems ^[32,33]. But the applications in ecology just starts, many issues need further studying, which are exactly what the paper is going to do.

The organization of the paper is as follows, firstly, a catastrophe model of the wheat aphid populations with delay is proposed and the catastrophe behaviors are analyzed using the qualitative theory of differential equations. Secondly, for biological management such as releasing natural enemy populations, the control model with impulse effects is proposed and the corresponding catastrophe controller is designed to make the wheat field ecosystem in a virtuous cycle. In the mean time, some simulations are carried to prove the results.

2.   Modelling

Considering the following model ^[8]

where $x$ is the density of aphid populations, $a$ is the corresponding intrinsic growth rate, $a, b, c$ are positive constants.

Zhao Huiyan et al. transform the model $(2.1)$ into the following,

where $a_1 = \frac{b}{2c}$, $P = \frac{a}{c}+(\frac{b}{2c})^2$, $a = r-d$ and $r, d$ are the birth rate and the death rate of the wheat aphid populations respectively.

Remark 2.1. In the the model $(2.2)$, $\frac{b}{2c}$ is the threshold below which the individual growth rate of the wheat aphid population increases and above which it decreases.

Remark 2.2. Obviously, the model $(2.2)$ is a fold catastrophe model. The parameter $P$ is not only related to the birth rate and the death rate but also to the threshold. From Figure 1, it is seen that the wheat aphid populations suddenly drops when $P = P_0 = 0$ and breaks out when $P = P_1$, which express the ecological mechanism of the pest sudden drop or outbreak from a view of dynamics. So it is reasonable to take the parameter $P$ as an important index of pest management.

Remark 2.3. Using the bimodality (two stable manifolds $x = 0$ and $x-a = \sqrt{P}$) in the model $(2.2)$, Zhao Huiyan explains the outbreak phenomenon of the wheat aphid populations.

For the predator population, Volterra introduces the delay to describe the effect of feeding predator, the following model is given ^[34]

For the term $\int_{-\infty}^{t}{F(t-\tau)N_1(\tau)d\tau}$ in the model $(2.3)$, N.Macdonald specifies $F(t)$ as the following form

Remark 2.4. For the term $\int_{-\infty}^{t}{F(t-\tau)N_1(\tau)d\tau}$ in the model $(2.3)$, we get the following result

On the other hand, the literature ^[35] shows that the dynamics of the bupleurum aphid populations is closely related to the amount of natural enemies and there exists a delay in the seedling stage.

Considering the above, the delay is introduced into the model $(2.1)$ and the following is obtained

where $\omega > 0$ is a parameter; the meaning of $x, a, b, c$ are the same as the model $(2.1)$.

Remark 2.5. $\omega$ in this model $(2.5)$ is interpreted as the foraging efficiency of the natural enemy populations for the wheat aphid populations For the model $(2.5)$, let $f(t - s) = e^{-a_0(t - s)}$ and $y = a_0\int_{-\infty }^t{e^{-a_0(t-s)}x(s)ds}$, $a_0$ is the saturation coefficient of the natural enemy populations, thus the following is got

Combining the model $(2.5)$ with the model $(2.6)$, the following is obtained

For the model $(2.7)$, let $dt^* = a_0dt, \ \gamma _1 = \frac{1}{{a_0 }} > 0, \ \gamma _2 = \frac{\omega } {{a_0 }} > 0$. For convenience, $dt^*$ is replaced with $dt$, the model is got,

Remark 2.6. For the model $(2.8)$, the coefficient $\gamma _2$ includes the parameter $\omega$ related to the foraging efficiency of natural enemy populations, thus it is seen as the main parameter of the model.

In the following sections, the model is analyzed using the qualitative theory of differential equations.

3.   Model analysis

As the parameter $\gamma_2$ in the model $(2.8)$ is the main parameter related to the natural enemy populations, it should be focused on in the model analysis.

Let

Obviously, $(x_0^*, y_0^*) = (0, 0)$ is an equilibrium of the model $(2.8)$, the others are determined by the following

Let $\Delta = (\gamma_2-b)^2+4ac$, $\Delta > 0$ implies that the model $(2.8)$ has two nonzero equilibriums, that is

Remark 3.1. Let $A_0 = (0, 0), A_1 = (x_1^*, y_1^*), A_2 = (x_2^*, y_2^*)$. Because $x_1^* < 0, y_1^* < 0$, thus equilibrium $A_1$ is meaningless ecologically, the qualitative analysis near the point $A_1$ is omitted.

In the following, the qualitative analysis near the point $A_0$ and $A_2$ are studied, let

where $f(x, y) = \gamma_1x(a+bx-cx^2-\gamma_2y), \ g(x, y) = x-y.$

According to the literature ^[36], the following theorems are obtained.

Theorem 3.1. For the model $(2.8)$,

$(1)$ The equilibrium $A_0$ is a saddle;

$(2)$ If $x_2^* > \frac{b}{{2c}}$ and $2\sqrt {D_2} < 1$, then the equilibrium $A_2$ is a stable node,

where $D_2 = \gamma_1 c(x_2^*)^2+a\gamma_1$ (see Figures 2 and 3.).

Proof. $(1)$ For the equilibrium $A_0$, the corresponding Jacobian determinant is

the characteristic equation is

where $T_0 = a\gamma_1-1, \ D_0 = -a\gamma_1.$

Obviously, $D_0 < 0$ means that the equilibrium $A_0$ is a saddle

$(2)$ Using the same method as the above, the characteristic equation at the equilibrium $A_2$ is got

where $T_2 = -2\gamma_1c(x_2^*)^2+b\gamma_1x_2^*-1, \ D_2 = \gamma_1c(x_2^*)^2+a\gamma_1.$

$D_2 > 0$ means that the equilibrium $A_2$ is a focus or a node, which can be determined by the sign of $Y = T_2^2-4{D_2}$. According to the condition $(2)$, $T_2 < -1$ and $Y > 0$, thus the equilibrium $A_2$ is a node; $T_2 < 0$, thus the equilibrium $A_2$ is stable(see Figure 2). Thus the equilibrium $A_2$ is a stable node (see Figure 3).

Theorem 3.2. For the equilibrium $A_2$ in the model $(2.8)$, the following results are obtained

$(1)$ $x_2^* > \frac{b}{{2c}}$ and $2\sqrt {D_2} > 1$,

$(i)$ if $- 2\sqrt {D_2} < T_2 < -1$, then the equilibrium $A_2$ is a stable focus (see Figures 4(a) and 5(a));

$(ii)$ if $T_2 < -2\sqrt {D_2}$, then the equilibrium $A_2$ is a stable node (see Figures 4(b) and 5(b)).

$(2)$ $x_2^* < \frac{b}{{2c}}$ and $2\sqrt D_2 < 1$,

$(i)$ if $-1 < T_2 < -2\sqrt {D_2}$, then the equilibrium $A_2$ is a stable node (see Figures 6(a) and 7);

$(ii)$ if $-2\sqrt {D_2} < T_2 < 0$, then the equilibrium $A_2$ is a stable focus (see Figures 6(b) and 8);

$(iii)$ if $0 < T_2 < 2\sqrt {D_2}$, then the equilibrium $A_2$ is an unstable focus (see Figures 6(c) and 9);

$(iv)$ if $T_2 > 1$, then the equilibrium $A_2$ is an unstable node (see Figure 6(d) and Figure 10).

$(3)$ $x_2^* < \frac{b}{{2c}}$ and $2\sqrt {D_2} > 1$,

$(i)$ if $-1 < T_2 < 0$, then the equilibrium $A_2$ is a stable focus (see Figures 11(a) and 12(a));

$(ii)$ if $0 < T_2 < 2\sqrt {D_2}$, then the equilibrium $A_2$ is an unstable focus (see Figures 11(b) and 12(b));

$(iii)$ if $2\sqrt {D_2} < T_2$, then the equilibrium $A_2$ is an unstable node (see Figures 11(c) and 12(c)).

Remark 3.2. For the result $(1)$ in Theorem 3.2, the condition $x_2^* > \frac{b}{{2c}}$ shows that the wheat aphid populations corresponding to the equilibrium $A_2$ is greater than its threshold density and stabilizes at the point $A_2$. If the condition $x_2^* > \frac{b}{{2c}}$ is taken as $x_2^* < \frac{b}{{2c}}$, the results are seen in $(3)$ of the Theorem 3.2 (Figure 6).

Remark 3.3. For the result $(2)$ in Theorem 3.2, the followings are found: when $T_2$ changes from $-1$ to $+\infty$, it undergoes two special values: one is $T_2 = 0$, the other is $T_2 = 2\sqrt {D_2}$. For the former, the equilibrium $A_2$ changes from a stable focus to an unstable one and is surrounded by a stable period orbit, thus $T_2 = 0$ is referred to as a catastrophe point at which the wheat aphid populations jump suddenly to the outbreak level (see from Figures 8 and 9). For the latter, the equilibrium $A_2$ changes from the focus to the node, and the stability remains unchanged.

Remark 3.4. For the result $(3)$ in Theorem 3.2, when $T_2$ changes from $-1$ to $+\infty$, there exits two special values: one is $T_2 = 0$, the other is $T_2 = 2\sqrt {D_2}$. For the former, the structure change of the model $(2.8)$ appears; that is, the equilibrium $A_2$ changes from a stable focus to an unstable one and is surrounded by a stable limit cycle (see Figure 10). For the latter, the stability of equilibrium $A_2$ remains unchanged, while the type changes.

For the results $(2)$ and $(3)$ in Theorem 3.2, when $T_2$ goes through $T_2 = 0$, the point $A_2$ of the model $(2.8)$ changes from a stable equilibrium to an unstable one and a stable limit circle appears. The limit circle in the result $(3)$ is located in the domain $\{(x, y)| 0\leq x\leq 4, 0\leq y \leq 1.2 \}$, which indicates that the system is in a virtuous cycle(see Figure 12(b)); while the limit circle in the result $(2)$ undergoes three jumps(see segment $A_0B$, $CD$ and $DA_0$ in Figure 8): The first one means that the amount of the wheat aphid populations changes from $x\approx 0.03$ to $x\approx 9000$ in a short time, which shows that the wheat aphid populations break out when the density of the enemy populations is low. The second one means that the amount of the wheat aphid populations decreases suddenly. The last one means that the migration of the enemy populations results from the lack of food. For the first one, it is very necessary to release additional enemy populations in a short time, the corresponding control is designed in the following section.

4.   Controller design

The aphid populations are typical R-strategy of small body size, short period, rapid multiplication, etc. the breakout phenomena are common. Figure 8(a) shows that $A_0\approx(0.03, 0.02)$ and $B\approx(9000, 0.02)$, which means that the amount of the wheat aphid populations changes dramatically from $x\approx 0.03$ to $x\approx 9000$, that is, the outbreak of the aphid populations. It also indicates that the natural enemy populations in ecological system can't to prevent the outbreak of the aphid populations by self-regulating mechanism, it is necessary to release additional enemy populations in a short time. Impulsive state feedback control may be used to describe the processes, the corresponding model is given

where $\alpha$ is the amount of the natural enemy populations $y(t)$ released, $k$ is the corresponding absorption rate for the wheat aphid populations, $x = h$ is an economic threshold which does damage to wheat.

For the model $(4.1)$, the following remarks are presented:

Remark 4.1. In the model $(4.1)$, the term $\Delta y = \alpha$ means that the additional enemy populations are released; The term $\Delta x = -k\alpha$ represents the amount of the aphid populations decreased by releasing the enemy populations and the control gain $k$ may be regulated to describe the predation effect.

Remark 4.2. The biological interpretation of the model $(12)$ is that the enemy populations should be released to keep the wheat aphid populations at the refuge level as soon as possible when the amount of the wheat aphid populations reaches the level $x = h$.

Remark 4.3. The condition $x_2^* < \frac{b}{{2c}}$ means that the equilibrium $A_2$ is on the left of the line $x = \frac{b}{{2c}}$.

Figure 13 is used to describe the structure of impulse control of the model $(12)$, in which the line $L_1: x = h$ is the impulsive set; the line $L_2: x = h-k\alpha$ is the corresponding image set; the point $A_0$ and the point $A_2$ are the equilibria; the line $\overline{OA}_2$ is Y-isoline; Y-axis and the parabola $\widehat{OA}_2$ are X-isoline; the point $P$ is the intersection of the parabola $\widehat{OA}_2$ and the image set $L_2$.

In the following section, the controller is designed to keep the wheat aphid populations at the refuge level by regulating the amount $\alpha$ of the enemy populations released. And then the qualitative analysis is studied to evaluate the control effect by means of the method in the literature $ ^[39] $, the specific steps are as follows:

Step I  Considering the orbit starting from the point $P$ on the line $L_1$, it jumps onto the point $P^+$ on the line $L_2$, and then moves to the point $P_1$ on the line $L_1$ as time goes. It is worth noting that the orbit becomes vertical as it crosses the parabola $\widehat{OA}_2$ and horizontal as it crosses the line $\overline{OA}_2$, thus the region $M$ encircled by the curve $PP^+P_1P$ is formed (see Figure 14).

The following steps show that all the orbits of the model ${(12)}$ finally go into the region encircled by the curve $PP^+P_1P$ and remain in it forever; that is, the control achieves the purpose for regulating wheat aphid populations.

Step II  If the initial point $Q$ is inside the region encircled by the curve $PP^+P_1P$, the orbit starting from the point hits the point $Q_0$ on the impulsive set $L_1$ and jumps onto the point $Q_0^+$ on the corresponding image set $L_2$, and then it moves on. Such processes may be repeated many times, but the orbit remains in the region forever.

Step III  If the initial point $R$ is outside the region encircled by the curve $PP^+P_1P$, the following cases are considered:

Case 1  The initial point $R$ is located on the left of the region, thus the orbit starting from the point hits the point $R_0$ on the impulsive set $L_1$, jumps onto the point $R_0^+$ on the image set $L_2$, and then goes into the region and remains in it.

Case 2  The initial point $R$ is located on the right of the region, thus the orbit starting from the point goes through the line $L_1$ from the right to the left, and then shows the same dynamics as Case 1.

Step IV  If the initial point $S$ is on the boundary of the region, for example $S = P$, the corresponding orbit shows the same dynamics as Step I.

Remark 4.4. The steps above describes the whole implementing process of controller design, the results show that the orbits starting from the points around the equilibrium $A_2$ eventually enters the region $M$ encircled by the curve $PP^+P_1P$ and remains in it, which indicates that the region is the attractive region of the model $(4.1)$.

Further, the existence of periodical solutions of the model $(4.1)$ may be investigated. In general, the two methods are used to finish the task, one is theory method, for example, the literature ^[37] generalizes some theories of periodical solutions from the differential equations to impulsive differential equations, and these theories are applied in a predator-prey model; The other is geometrical analysis, for example, the literature ^[38] generalizes successor functions and Poincare-Bendixson theorem from the differential equations to impulsive differential equations, and gives the geometrical analysis of a predator-prey model with two state impulses. In the following, the latter is used to study the existence of periodical solutions of the model $(4.1)$.

In Figure 14, Let $P = (x_p, y_p)$, $P_1 = (x_{p_1}, y_{p_1})$, thus the successor function of the point $P$ is

Similarly, let $R = (x_r, y_r)$, the corresponding image point is $R_1 = (x_{r_1}, y_{r_1})$. According to Remark 3.4, we have

Combine $(4.2)$ and $(4.3)$, the existence of periodical solutions of the model $(4.1)$ is proven.

5.   Discussion

In fact, the outbreak of aphid populations is a catastrophe phenomenon, the ecologists usually give qualitative explaination from the ecological views; In the paper, the delay is introduced the catastrophe model to reveal the mechanism behind the outbreak, which provides the mathematical details for the application of catastrophe theory in ecosystems. Due to the rapid multiplication of aphid populations, releasing enemy and spraying become the main measures of aphids management in agricultural production. The paper uses impulsive control to express the process, in which a new method for the theory study of aphid populations is provides. In a sense, the above realizes the combination of ecology, mathematics and catastrophe control theory. On the other hand, the essence of catastrophe control is the criticality control of systems. Literature ^[39] gives the method of identifying the type of critical point and literature ^[40] illustrate the main achievements in the study of critical dynamics in biological systems. Whether they are applied to regulate the catastrophe behavior of the wheat aphid populations is the future work.

Acknowledgments

This work is supported by Natural Sciences Foundation of People's Republic of China (12171004) and Liaoning Provincial Education Department Fund of People's Republic of China (SYSF202010, LJKQ2021179) and Liaoning Provincial Doctoral Research Startup Fund (2021-BS-299).

Conflict of interest

The authors declare no conflicts of interest.