Qualitative Analysis of the Effect of Weeds Removal in Paddy Ecosystems in Fallow Season

In the paper, we introduce a diﬀerential equations model of paddy ecosystems in the fallow season to study the eﬀect of weeds removal from the paddy ﬁelds. We found that there is an unstable equilibrium of the extinction of weeds and herbivores in the system. When the intensity of weeds removal meets certain conditions and the intrinsic growth rate of herbivores is higher than their excretion rate, there is a coexistence equilibrium state in the system. By linearizing the system and using the Routh–Hurwitz criterion, we obtained the local asymptotically stable conditions of the coexistence equilibrium state. The critical value formula of the Hopf bifurcation is presented too. The model demonstrates that weeds removal from paddy ﬁelds could largely reduce the weeds biomass in the equilibrium state, but it also decreases the herbivore biomass, which probably reduces the content of inorganic fertilizer in the soil. We found a particular intensity of weeds removal that could result in the minimum content of inorganic fertilizer, suggesting weeds removal should be kept away from this intensity.


Introduction
Paddy field fertility is the premise of high quality, high yield, and low energy consumption of rice. e method of increasing rice yield by applying a large amount of chemical fertilizers and pesticides has threatened the sustainable development of grain and the ecological security. erefore, how to use biodiversity in paddy ecosystems to improve the stability and sustainability of agricultural systems has been attracting more and more attention [1][2][3].
ere are many factors affecting the fertility of a paddy ecosystem, such as soil nutrient availability, light, moisture, weeds, insects, microorganisms, and so on. Not only the number of these factors is big but also the relations among them are complicated. In addition, there are many human disturbances in a paddy ecosystem. Because of these problems, the current research focuses on the experiments and data analyses of paddy ecosystems but rarely on the dynamics of the systems by establishing mathematical models. For example, through experiments, it is found that putting ducks in rice fields could improve soil fertility, control weeds, and reduce rice diseases [4][5][6][7]. Our interest is in using models of differential equations to study the complex nonlinear relationships in paddy ecosystems. For the forest and aquatic ecosystems, some matured mathematical models have been established [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. For example, Hofmann and Ambler developed a system consisting of ten coupled ordinary differential equations to describe the interactions among nitrate, ammonium, two phytoplankton size components, five copepods, and a debris pool on the continental shelf outside the southeastern United States [13]. For an aquatic ecosystem, it is very important to explore the reproductive mechanism of phytoplankton. Because periodic solutions in mathematics can represent reproduction, many scholars have studied the existences of periodic solutions, almost periodic solutions, and Hopf bifurcations of some aquatic ecosystems [8,10,12,14,[17][18][19]27], and some scholars have studied the effect of time delay on nutrient metabolism kinetics [17].
In recent years, we have established some mathematical models of paddy ecosystems from different aspects [28][29][30][31][32]. We built a differential equations model of a paddy ecosystem in fallow season in [28], and studied the interaction between weeds and soil nutrient availability with the discovery of the existence of the stable node, unstable saddle point, or saddlenode in the system. Xiang et al. found that adding the factor of domesticated animals in a paddy ecosystem in the fallow season could improve the level of nutrient availability in the soil [29], and some new features such as a Hopf bifurcation appear in the system.
In addition to farming animals in paddy fields, weeds removal is also a common activity of farmland management. Farmers gather weeds from paddy fields to reduce the consumption of soil nutrients and to feed animals too. Removing weeds from the paddy fields has some unavoidable impacts on the paddy ecosystem. We are interested in how the factors in a paddy ecosystem are affected by this activity and what principles should be followed to measure it. e main purpose of this paper is to analyze the effect of weeds removal in a paddy ecosystem in the fallow season by establishing a differential equations model and find a strategy of weeds removal. e rest of the paper is organized as follows. In the next section, the mathematical model of the paddy ecosystem with weeds removal is introduced. We present the conditions of existence of equilibria in Section 3. en, we consider the stability of equilibria in Section 4. In Section 5, we analyze the effect of weeds removal on the main factors in the paddy ecosystem. We study the Hopf bifurcation from the coexistence equilibrium in Section 6. Some numerical simulations are carried out in Section 7. Finally, we give our conclusion.

Modeling of the Paddy Ecosystem in the Fallow Season with Weeds Removal
We only consider the effect of weeds removal on three main components of a paddy ecosystem in fallow season: weeds, inorganic fertilizer, and herbivores. ere is an obvious nutrient-dependent consumption ring among the three: weeds are used as food for farmed animals in paddy fields, weeds are nourished by the effective nutrients in the soil, and the effective nutrients can be transformed from animal manure. In [29], Xiang et al. proposed the following differential equation model of a paddy ecosystem in the fallow season to analyze the interactions among the three main components: where p(t) is the weeds biomass per unit area at time t, u(t) is the inorganic fertilizer content per unit area, and A(t) is the herbivore biomass per unit area. e function f(I) � (f m I/(I k + I)) is a light effect function. e function H(u) � (au/(m + u)) is used as Michaelis-Menten uptake kinetics of the general form, where m is the half saturation concentration of inorganic fertilizer and a is the maximum uptake rate of inorganic fertilizer u. e function is the rate of grazing weeds per herbivore, where c is the maximum rate of grazing weeds by a herbivore and λ is a constant affecting the grazing rate. e three terms on the right of the first equation in system (1), respectively, represent the growth of weeds, the herbivore grazing, and the natural mortality of weeds, which the first term f(I)H(u)p shows the weeds growth is affected by the light intensity I and the inorganic fertilizer u, d 1 is the mortality rate of weeds. ere are four terms on the right hand side of the second equation in system (1), the first term is the consumption of inorganic fertilizer by the growth of weeds, the second term represents the inorganic fertilizer converted from some weeds that have been wasted by the herbivore, the third term is the inorganic fertilizers transformed from dead weeds, and the last term is herbivore excrements. e three terms on the right hand side of the last equation in system (1) represent in turn the growth of herbivore, the herbivore assimilating of weeds and herbivore excretion, where b is the intrinsic growth rate of herbivore and K is the largest environmental capacity of herbivore, r is the assimilating ratio of herbivore, and d 2 is the excretion rate of herbivore.
To consider the effect of weeds removal from the paddy fields, we introduce the weeds removal term d 0 p into the first equation of the model (1), e third term d 0 p on the right of the first equation represents the rate of artificial weeds removal and d 0 is called the weeds removal intensity.
To reduce the system parameters, we make the following transformations, p � (p/K), u � (u/K), A � (A/K), m � (m/K), a � f(I)a, and λ � λK. System (3) thereby is in the following form: In system (4), we have removed the bar above some variables or parameters. Note that the parameter a in the function H(u) represents both the light effect and the uptake rate of inorganic fertilizer. Here, we call the parameter a the coupling effect coefficient of light and fertilizer.
In the sense of ecology, the parameters in system (4) are nonnegative, and 0 < r < 1, d 1 > 0, and d 2 > 0. e properties of paddy ecosystem without weeds removal have been studied in [29], so we only consider the case d 0 > 0.

The Existence of Equilibria
Denote the three expressions on the right hand side of system (4) by In order to obtain the equilibrium states, we need to solve the nonlinear equations From (6), we know 0 ≤ A ≤ 1.
Hence, there exists one equilibrium of system (4) as follows: where u 1 is an arbitrary nonnegative real constant.
holds, we have the following: (2) and (8), we obtain the following: us, we have d(A) < c. It follows that Denote A * � min 1, 1 − (d 2 /b + (rc/b)}. Notice that 0 ≤ A ≤ 1, we have the following: erefore, if the variable A locates in the above range, then we get Substituting p into (6) yields Formula (14) is regarded as an equation of an unknown quantity A. For the existence of the roots of equation (14), we give the following Lemma 1. (14) has at least one positive root A 2 ∈ (A * , A * ).
If rc > d 2 , then A * � 1, and we have the following: , and we have the following: us, the function φ(A) must have a zero point erefore, equation (14) has at least one positive root A 2 ∈ (A * , A * ).
We make the following hypothesis on the positive root en, we obtain Discrete Dynamics in Nature and Society , it follows that the numerator of u 2 is positive. Denote the denominator of u 2 by Obviously, if the assumption (H1) holds, then the denominator erefore, we deduce that u 2 > 0. Based on the above analysis, we obtain the following conclusions about the existence of equilibria in system (4). □ Theorem 1. System (4) has weeds and herbivore extinction equilibrium and (H1) hold, there also exists at least one positive equilibrium E p � (p 2 , u 2 , A 2 ), where u 1 is an arbitrary nonegative real constant, A 2 is one positive root of equation (14), and p 2 and u 2 are given in (19) and (21), respectively. e biological significance of condition b > d 2 in Lemma 1 and eorem 1 is that the intrinsic growth rate b of herbivore should be greater than its excretion rate d 2 . e use of hypothesis (H1) is to ensure that u 2 is positive. Regarding hypothesis (H1), let us make a further explanation. If we set then hypothesis (H1) is equivalent to Obvousily, A * * > A * . But the relationship between A * * and A * is complex.
hold, then we may not need a hypothesis (H1). at is to say, (H1) should be assumed only if Generally speaking, in order to ensure the normal growth of herbivore, its excretion rate d 2 should not be greater than its assimilation rate rc of weeds (that is d 2 ≤ rc). e case of "d 2 > rc" may appear in the period when the herbivore is sick. erefore, under normal circumstances, condition (H1) is necessary to ensure the existence of positive equilibrium.
Next, we discuss the stability of the positive equilibrium E p . For convenience, we introduce the following symbols to express the coefficients of system (28), . us, we obtain the following: Considering the following facts: it is easy to see that a 12 and a 31 are all positive. As for the positive and negative properties of a 11 and a 21 , we have the following conclusions.

Proof. Let
From in (13), we have Notice that c > d(A 2 ), we get the following result by Lemma 2.1 in [29]. erefore, Substituting (33) into the expression of a 11 yields erefore, a 11 > 0. From (31) and (33), we obtain that We introduce the following notations: and define a quadratic function It is easy to verify that lim A⟶±∞ J(A) � − ∞.
If rc > d 2 , then otherwise Proof. According to the definition of a 31 in (31), by replacing h ′ (p 2 ) with (33), we get the following: Notice that the following facts hold: J(A * ) > 0; the highest power coefficient of quadratic polynomial J(A) is a negative number; A − and A + are two real roots of equation Next, we give some conclusions on the stability of the positive equilibrium E p . is unstable.
Proof. For the coexistence equilibrium E p � (p 2 , u 2 , A 2 ), the characteristic equation of the linearized system (28) is as follows: Expanding the determinant, we get the following: (i) If A 2 > A − holds, then we know a 3 > 0 from Lemma 3. us we have a 1 > 0 from the assumptions a 2 > 0 and a 1 a 2 > a 3 . erefore, using the Routh-Hurwitz criterion, the positive equilibrium E p is locally asymptotically stable. (ii) If A 2 < A − , then a 3 < 0 from Lemma 3. Hence, the characteristic equation (47) has at least one positive real root. erefore, the positive equilibrium (p 2 , u 2 , A 2 ) is unstable.

(48)
We have the following corollary. □ Corollary 1. Suppose that A 2 > A − holds. If a 12 > a 11 , then, the positive equilibrium E p is locally asymptotically stable.

Effects of Weeds Removal and Other Environmental Parameters on Main Factors in Paddy Fields
From eorem 1, the coexistence equilibrium state (p 2 , u 2 , A 2 ) exists when b > d 2 and (H1) hold, where A 2 is one of the solutions of equation (14), and p 2 and u 2 are computed in (19), (21), respectively.
From (14), we obtain the following: If (1/2) ≤ A ≤ 1, then (dd 0 /dA) < 0. Otherwise, by Lemma 2.1 in [29], we have the following: Under the condition b > d 2 , the discriminant of the quadratic function bA 2 It gives that (dd 0 /dA) < 0. us, we get the following: In the coexistence equilibrium state, p 2 can still be calculated by (13). It is not difficult to obtain that dp 2 erefore, the biomass of herbivores and weeds decreases monotonously with the increase of the intensity of 6 Discrete Dynamics in Nature and Society weeds removal. Obviously, the result that weeds biomass p 2 decreases monotonously with the increase of weeds removal intensity d 0 is consistent with the conventional understanding, but interestingly, it has nothing to do with the coupling effect coefficient a of light and fertilizer, the mortality rate d 1 of weeds, and the half saturation concentration m of inorganic fertilizer. e expression (21) of the inorganic fertilizer content u 2 can also be written as follows: Notice that the herbivores biomass A 2 is independent of parameters a, m, and d 1 , hence u 2 is proportional to the half saturation concentration m of inorganic fertilizer. If the efficiency of photosynthesis of weeds or absorbing inorganic fertilizers increases (increases a), the content of inorganic fertilizer at a steady state will decrease. If the mortality rate d 1 of weeds is increasing, the content of inorganic fertilizer could be increasing. However, the relationships between inorganic fertilizer biomass and other parameters are not so simple. We only consider the relationship between u 2 and d 0 . Because A 2 is affected by the parameter d 0 , we denote where d 0 ≥ 0. us, the inorganic fertilizer content u 2 can be rewritten as follows: It is easy to obtain that we have the following: erefore, w(d 0 ) consists of two parts, one part increases with the increase of d 0 and the other part related to A 2 decreases with the increase of d 0 . Hence, there must be a unique d 0 such that (dw/dd 0 ) � 0. From this equation, we have the following: From (14) and (49), we get the following: Multiplying (60) and (61), we obtain the following: (62) By solving simultaneous equations (14) and (62), we obtain the unique d * 0 that satisfies (dw/dd 0 ) � 0. Notice that us the parameter value d * 0 is a minimum point of w(d 0 ). Furthermore, from (57), it is easy to know that the nondifferentiable points of u 2 (such that a − d 1 − w(d 0 ) � 0) are not extreme points. erefore, the weeds removal intensity d * 0 minimizes the inorganic fertilizer content u 2 . is shows that with the increase of parameter d 0 , the content of inorganic fertilizer u 2 first has a downward trend, then when d 0 � d * 0 , u 2 reaches the minimum, and then u 2 will show an upward trend.

Hopf Bifurcation
In this section, taking m as the Hopf bifurcation parameter, we consider the existence of a Hopf bifurcation of system (4) at the coexistence equilibrium E p . We always assume that the positive equilibrium E p exists and A 2 > A − holds.
From equation (14) for determining A 2 and expression (19) for calculating p 2 , we know that A 2 and p 2 in equilibrium (p 2 , u 2 , A 2 ) are independent of the half saturation concentration m of inorganic fertilizer, but u 2 is related to m. We still find that the root A − of equation (42) is also independent of m. According to (31), we can see that the parameters a 12 , a 21 , and a 22 are also related to m, while a 13 , a 23 , a 31 , and a 33 are not related to m. From (39), we can see that the parameter a 11 is not related to m either. From (41), the coefficients of the characteristic equation (47) depend on the parameter m. Let us set a 1 � a 1 (m), a 2 � a 2 (m) and a 3 � a 3 (m).
Lemma 3 tells us the constant term a 3 of the characteristic polynomial Δ(λ) is positive under the assumption A 2 > A − . erefore, a Hopf bifurcation takes place when the real part of a pair of conjugate complex eigenvalues changes sign.
From (21), we get the following: (65) Substituting (65) into (33), by (31) we obtain the following: It shows that a 12 is related to m. From (41), we have the following: (67) erefore, the condition L(m 0 ) � 0 can be rewritten as follows: Of all the parameters of the equation (68), only a 12 is related to the bifurcation parameter m. erefore, by calculating a 12 from equation (68), the bifurcation parameter m can be determined by (66).
In addition, from (31), (39), and (66), we obtain that Notice that g(A 2 ) > 0, the condition L ′ (m 0 ) ≠ 0 can be replaced by the following condition: According to Lemma 4, we can get the existence of a Hopf bifurcation of system (4).
Hence, once a 12 is obtained from equation (68), we can calculate the critical value m 0 with (71), which is obtained by (66).
According to the conditions of eorem 4, we summarize the process of calculating the critical value m 0 of the Hopf bifurcation as follows: (1) If the condition b > d 2 holds, then A 2 is obtained by solving equation (14) and p 2 is calculated by (19). If the condition (H1) is validated, then system (4) has a positive equilibrium E p .
By the process of calculating the critical value m 0 of a Hopf bifurcation in Section 6, we obtain the calculation results of the Hopf bifurcation, which are listed in Table 1. From the line of condi 1 in the table and a � 40, it follows that hypothesis (H1) holds, so there is a positive equilibrium point in system (4) by eorem 1. It is easy to see from the second and seventh rows of Table 1 that A 2 > A − . From the second line to the fourth line counted backward in the table, we obtain a 12 > a 11 − bA 2 and a 12 ≠ cond i 2 . erefore, the conditions in eorem 4 are all satisfied.
When the parameter b � 4.4, from   (1) cond   Combining the above calculation results when b � 4.4, we see that when the parameter m passes through the critical values 0.024358 and 4.318921 from left to right, the positive equilibrium of system (4) undergoes a process of change from stability to instability and then to stability. Figure 3 depicts the relationship between the maximum real part of the eigenvalue of the characteristic equation (47) with b � 4.4 and the parameter m. It can be seen from Figure 3 that the maximum real part of the eigenvalue is positive in the interval [0.024358, 4.318921], indicating that system (4) is unstable in this range, while the maximum real part of the system outside the interval is negative, which indicates that system (4) is stable in this range. is is consistent with our calculation in the above.

Conclusion
By establishing a differential equation model of paddy ecosystem in the fallow season with weeds removal and analyzing its stability and Hopf bifurcation, we found the interaction among weeds, inorganic fertilizer, and herbivores in this system.
In a paddy ecosystem in the fallow season, because of the human management activities of weeds removal and animal farming, the extinction of weeds and herbivores may occur. However, the extinction state of weeds and herbivores is unstable. From (30), we know that as long as the growth rate of herbivores is higher than their excretion rate, the extinction of herbivores in paddy fields can be avoided.
Furthermore, when b > d 2 , if the coupling effect coefficient of light and fertilizer a is higher than the sum of weeds removal intensity d 0 , weeds mortality d 1 , and (d 0 d(A 2 )/(b(1 − A 2 ))) (assumption (H1)), the coexistence equilibrium E p can be found in the system.
Our results are obtained under the condition of d 0 > 0. When weeds are not removed, i.e., d 0 � 0, the situation is relatively simple. Xiang et al. has considered the case of d 0 � 0 in [29]. ey obtained that p 2 and u 2 as follows: Obviously, they need conditions rc > d 2 and p 2 > d 2 /(r(a − d 1 )). In eorem 1, we replaced their condition p 2 > d 2 /(r(a − d 1 )) with (H1) and deleted condition rc > d 2 . us, our results include the case of rc ≤ d 2 . Because conditions (H1) and (25) are equivalent, from (25), we know that the intensity of weeds removal d 0 should not be too high By using the Routh-Hurwitz criterion, we obtain the conditions for local asymptotic stability of the positive equilibrium E p , in which the condition A 2 > A − is used to guarantee the coefficient a 3 > 0. If d 0 � 0, then we have A 2 � 1 from (6), so the coefficient a 3 � a 12 a 31 b > 0 is natural. So there is no requirement for the lower bound of A 2 in eorem 2.5 for the local stability of positive equilibrium in [29]. However, if weeds removal is considered in the system, then whether a 3 is positive or negative is affected by the value A 2 . From (41), we see that A 2 should not be too small, otherwise a 3 will be negative. is is why we propose the condition A 2 > A − in eorem 3.
When the coefficients a i of cubic characteristic polynomial (47) are positive, the Hopf bifurcation may be generated when a 1 a 2 � a 3 . According to this principle, we obtain the Hopf bifurcation conditions of system (4) at the positive equilibrium E p , and give the Hopf bifurcation critical value formula (71) with m being the bifurcation parameter. e half saturation concentration of inorganic fertilizer is a key parameter in Michaelis-Menten uptake 10 Discrete Dynamics in Nature and Society kinetics, which is affected by species of weeds. A Hopf bifurcation appears at a positive equilibrium E p when the parameter m passes through the critical value m 0 , which indicates that different kinds of weeds in paddy field may lead to different dynamic properties of the paddy ecosystem. It is necessary to know the half saturation concentration of inorganic fertilizer for each weeds, so as to grasp the stability of the paddy ecosystem in advance. From our numerical simulation results, the paddy ecosystem with high half saturation concentration of weeds is more prone to Hopf bifurcation. From the analysis in Section 5, we also see that weeds removal can reduce the biomass of weeds in an equilibrium state. However, this activity also reduces the biomass of herbivores in paddy fields. If the intensity of weeds removal is d * 0 , the soil inorganic fertilizer content is the lowest. e weeds removal intensity is too small to achieve the purpose of weeding. erefore, the intensity of weeds removal should be greater than the minimum point d * 0 as far as possible. But considering that d 0 is restricted by condition (H1), it has an upper bound. For example, we can know from (25) that d 0 cannot exceed a − d 1 at least. If the minimum point d * 0 is greater than the upper bound, then we should try to use lowintensity weeding. If d * 0 is less than the upper bound, the intensity of weeds removal should be close to the upper bound and not exceed it.

Data Availability
No data were used to support this study. In our research, we use numerical simulation to verify our main results. e system parameters used in the simulation have been included within this paper.

Conflicts of Interest
All authors declare no conflicts of interest in this paper.