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

Authors: Leru Zhou; Zhigang Liu; Tiejun Zhou
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

In the paper, we introduce a differential equations model of paddy ecosystems in the fallow season to study the effect of weeds removal from the paddy fields. 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 fields 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.


---

## Body

## 1. Introduction

Paddy field fertility is the premise of high quality, high yield, and low energy consumption of rice. The 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. Therefore, how to use biodiversity in paddy ecosystems to improve the stability and sustainability of agricultural systems has been attracting more and more attention [1–3].There 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–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–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–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–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 saddle-node 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. The 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.The 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 Section3. Then, 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.


## 2. 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. There 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:(1)p˙t=fIHup−hpA−d1p,u˙t=−fIHup+1−rhpA+d1p+d2A,A˙t=bA1−AK+rhpA−d2A,where pt is the weeds biomass per unit area at time t, ut is the inorganic fertilizer content per unit area, and At is the herbivore biomass per unit area. The function fI=fmI/Ik+I is a light effect function. The function Hu=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. The function(2)hp=c1−e−λp,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. The 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 fIHup shows the weeds growth is affected by the light intensity I and the inorganic fertilizer u, d1 is the mortality rate of weeds. There 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. The 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 d2 is the excretion rate of herbivore.To consider the effect of weeds removal from the paddy fields, we introduce the weeds removal termd0p into the first equation of the model (1),(3)p˙t=fIHup−hpA−d0p−d1p,u˙t=−fIHup+1−rhpA+d1p+d2A,A˙t=bA1−AK+rhpA−d2A.The third termd0p on the right of the first equation represents the rate of artificial weeds removal and d0 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¯=fIa, and λ¯=λK. System (3) thereby is in the following form:(4)p˙t=Hup−hpA−d0p−d1p,u˙t=−Hup+1−rhpA+d1p+d2A,A˙t=bA1−A+rhpA−d2A.In system (4), we have removed the bar above some variables or parameters. Note that the parameter a in the function Hu 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, d1>0, and d2>0. The properties of paddy ecosystem without weeds removal have been studied in [29], so we only consider the case d0>0.


## 3. The Existence of Equilibria

Denote the three expressions on the right hand side of system (4) by(5)F1p,u,A=Hup−hpA−d0p−d1p,F2p,u,A=−Hup+1−rhpA+d1p+d2A,F3p,u,A=bA1−A+rhpA−d2A.In order to obtain the equilibrium states, we need to solve the nonlinear equationsFip,u,A=0i=1,2,3. Accordingly, we have F1p,u,A+F2p,u,A+F3p,u,A=0, or(6)d0p=bA1−A.From (6), we know 0≤A≤1.According toF3p,u,A=Ab1−A+rhp−d2=0, we have A=0 or b1−A=d2−rhp.IfA=0, then F1p,u,A=Hu−d0−d1p=0 and F2p,u,A=−Hu+d1p=0; thereby, we obtain that p=0. Hence, there exists one equilibrium of system (4) as follows:(7)p1,u1,A1=0,u1,0,where u1 is an arbitrary nonnegative real constant.If equationb1−A=d2−rhp holds, we have the following:(8)hp=d2−b+bAr.Define(9)dA=1rd2−b+bA,where A∈0,1. Obviously, dA=hp. From hp≥0, we know that dA≥0 when A∈0,1, so we have A≥1−d2/b=A∗. Synthesizing (2) and (8), we obtain the following:(10)e−λp=1−dAc>0.Thus, we havedA<c. It follows that(11)A<1−d2b+rcb.DenoteA∗=min1,1−d2/b+rc/b. Notice that 0≤A≤1, we have the following:(12)max0,A∗≤A<A∗.Therefore, if the variableA locates in the above range, then we get(13)p=−1λln1−dAc.Substitutingp into (6) yields(14)d0ln1−dAc=−λbA1−A.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.Lemma 1.
Ifb>d2, then, equation (14) has at least one positive root A2∈A∗,A∗.Proof.
Fromb>d2 we know A∗>0. Define(15)φA=λbA1−A+d0ln1−dAc.
It is continuous in the intervalsA∗,A∗. Obviously, the value of the function φA at A=A∗ is(16)φA∗=λd21−d2b.
Ifrc>d2, then A∗=1, and we have the following:(17)φA∗=d0ln1−d2rc<0.
Ifrc≤d2, then A∗=1−d2/b+rc/b, and we have the following:(18)limA⟶A∗−0φA=λ1−d2b+rcbd2−rc+d0limA⟶A∗−0ln1−d2−b+bArc=−∞.
Thus, the functionφA must have a zero point A2∈A∗,A∗. Therefore, equation (14) has at least one positive root A2∈A∗,A∗.
We make the following hypothesis on the positive rootA2,(H1).a>d0+d1+d0dA2/b1−A2.
SubstitutingA=A2 into (6) yields(19)p2=bd0A21−A2.
By substitutingA=A2 and p=p2 into F1p,u,A=0, we get(20)Hu=aum+u=d0+d1+A2hp2p2=d0+d1+A2dA2p2.
Then, we obtain(21)u=mb1−A2d0+d1+d0dA2b1−A2a−d0−d1−d0dA2Δ¯¯u2.
Notice thatdA>0 for A∈A∗,A∗, it follows that the numerator of u2 is positive. Denote the denominator of u2 by(22)gA2=b1−A2a−d0−d1−d0dA2.
Obviously, if the assumption (H1) holds, then the denominator(23)gA2=b1−A2a−d0−d1−d0dA2b1−A2>0.
Therefore, we deduce thatu2>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 p1,u1,A1=0,u1,0. If conditions b>d2, and (H1) hold, there also exists at least one positive equilibrium Ep=p2,u2,A2, where u1 is an arbitrary nonegative real constant, A2 is one positive root of equation (14), and p2 and u2 are given in (19) and (21), respectively.The biological significance of conditionb>d2 in Lemma 1 and Theorem 1 is that the intrinsic growth rate b of herbivore should be greater than its excretion rate d2. The use of hypothesis (H1) is to ensure that u2 is positive. Regarding hypothesis (H1), let us make a further explanation. If we set(24)A∗∗=1−d0d2bra−d0−d1+d0,then hypothesis (H1) is equivalent to(25)a>d0+d1,A2<A∗∗.Obvousily,A∗∗>A∗. But the relationship between A∗∗ and A∗ is complex. If d2≤rc, then A∗=1>A∗∗. If d2>rc, then A∗=1−d2/b+rc/b. Furthermore, if d0>a−d1/1+c/d2−rc, then we have A∗>A∗∗. Otherwise, A∗≤A∗∗. Therefore, if conditions(26)d2>rc,d0≤a−d11+c/d2−rc,hold, then we may not need a hypothesis (H1). That is to say, (H1) should be assumed only if d2≤rc or d2>rc and d0>a−d1/1+c/d2−rc.Generally speaking, in order to ensure the normal growth of herbivore, its excretion rated2 should not be greater than its assimilation rate rc of weeds (that is d2≤rc). The case of “d2>rc” may appear in the period when the herbivore is sick. Therefore, under normal circumstances, condition (H1) is necessary to ensure the existence of positive equilibrium.


## 4. The Stability of Equilibria

Now, let us discuss the stabilities of equilibriump∗,u∗,A∗ of system (4). With a coordinate transformation x=p−p∗, y=u−u∗, z=A−A∗, system (4) is converted to(27)x˙t=−d0+d1xt+Hy+u∗x+p∗−Hu∗p∗−hx+p∗z+A∗−hp∗A∗,y˙t=d1xt+d2zt−Hy+u∗x+p∗−Hu∗p∗+1−rhx+p∗z+A∗−hp∗A∗,z˙t=b−2bA∗−d2z−bz2+rhx+p∗z+A∗−hp∗A∗.The linearized system of (27) at the equilibrium 0,0,0 is as follows:(28)x˙t=Hu∗−h′p∗A∗−d0−d1x+p∗H′u∗y−hp∗z,y˙t=d1−Hu∗+1−rh′p∗A∗x−p∗H′u∗y+d2+1−rhp∗z,z˙t=rh′p∗A∗x+b−2bA∗+rhp∗−d2z,where H′u=am/m+u2, and h′p=cλe−λp.The stability of the equilibriump1,u1,A1=0,u1,0 of system (4) is presented in Theorem 2.Theorem 2.
The equilibriump1,u1,A1=0,u1,0 of system (4) is unstable.Proof.
For the equilibriump1,u1,A1=0,u1,0, we have hp1=h0=0 and h′p1=h′0=cλ, so the linearized system (28) becomes as follows:(29)x˙t=Hu∗−d0−d1x,y˙t=d1−Hu∗x+d2z,z˙t=b−d2z.
Its solution under initial conditionx0,y0,z0 is as follows:(30)xt=x0eHu∗−d0−d1t,zt=z0eb−d2t,yt=x0d1−Hu∗Hu∗−d0−d1eHu∗−d0−d1t+z0d2b−d2eb−d2t+y0.
It shows that the equilibrium0,0,0 of system (29) is unstable. Therefore, the equilibrium p1,u1,A1=0,u1,0 of system (4) is also unstable.
Next, we discuss the stability of the positive equilibriumEp. For convenience, we introduce the following symbols to express the coefficients of system (28),(31)a11=Hu2−A2h′p2−d0−d1,a12=p2H′u2,a21=−Hu2+1−rA2h′p2+d1,a23=d2+1−rhp2,a31=rA2h′p2,a33=b−2bA2+rhp2−d2,a13=−hp2,a22=−a12.
Obviously,a21=−a11−a31−d0. At the positive equilibrium Ep, we have b1−A2=d2−rhp2 and hp2=dA2. Thus, we obtain the following:(32)a13=−dA2,a23=dA2+b1−A2,a33=−bA2.
Considering the following facts:(33)H′u2=amm+u22>0,h′p2=cλe−λp2=λc−dA2>0,it is easy to see that a12 and a31 are all positive. As for the positive and negative properties of a11 and a21, we have the following conclusions.Lemma 2.
At the coexistence equilibriumEp, a11>0 and a21<0.Proof.
Let(34)D0A2≜dA2−λp2c−dA2.
From in (13), we have(35)D0A2=dA2+c−dA2ln1−dA2c=c−dA2dA2c−dA2+ln1−dA2c.
Notice thatc>dA2, we get the following result by Lemma 2.1 in [29].(36)−ln1−dA2c<dA2/c1−dA2/c=dA2c−dA2.
It follows thatD0A2>0.
FromF1p2,u2,A2=Hu2p2−A2hp2−d0+d1p2=0, hp2=dA2, we have(37)Hu2−d0−d1=A2dA2p2.
Therefore,(38)a11=A2dA2p2−A2h′p2=A2p2dA2−p2h′p2.
Substituting (33) into the expression of a11 yields(39)a11=A2D0A2p2.
Therefore,a11>0.
From (31) and (33), we obtain that(40)a21=−a11+d0−a31=−a11+d0−rA2h′p2=−a11+d0−rλc−dA2<0.
We introduce the following notations:(41)a1=a12−a11+bA2,a2=bA2a12−a11+a12d0+a31dA2+a12,a3=a12d0bA2+a12a31b2A2−1,and define a quadratic function(42)JA=−2bλA2+λ2rc−2d2+3bA+d0−λrc−d2+b.
It is easy to verify thatlimA⟶±∞JA=−∞.
Ifrc>d2, then(43)JA∗=J1=λrc−d2+d0>0,otherwise(44)JA∗=J1−d2b+rcb=d0>0.
Therefore, the quadratic equationJA=0 has two real roots A−<A∗ and A+>A∗.
According to the properties of polynomialJA, we can determine whether a3 defined in (41) is a positive number.Lemma 3.
IfA2>A−, then a3>0. Otherwise, if A2<A−, then a3<0.Proof.
According to the definition ofa31 in (31), by replacing h′p2 with (33), we get the following:(45)a3=a12bd0A2+rA2h′p22A2−1=a12bA2d0+rλc−dA22A2−1=a12bA2d0+λrc−d2+b−bA22A2−1=a12bA2−2bλA22+λ2rc−2d2+3bA2+d0−λrc−d2+b=a12bA2JA2.
Notice that the following facts hold:JA∗>0; the highest power coefficient of quadratic polynomial JA is a negative number; A− and A+ are two real roots of equation JA=0; and A2<A∗<A+. Therefore, when A∈A−,A+, we get JA>0, and when A<A− or A>A+, we get JA<0.
Thus, considering thatA2<A∗<A+, if A2>A−; then, we have JA2>0, implying a3>0. Otherwise, if A2<A−, then JA2<0, implying a3<0.
Next, we give some conclusions on the stability of the positive equilibriumEp.Theorem 3.
(i)
Suppose thata2>0 and a1a2>a3. If A2>A− holds, then the positive equilibrium p2,u2,A2 is locally asymptotically stable.(ii)
IfA2<A−, then, the positive equilibrium p2,u2,A2 is unstable.Proof.
For the coexistence equilibriumEp=p2,u2,A2, the characteristic equation of the linearized system (28) is as follows:(46)Δλ=λ−a11−a12dA2a11+a31+d0λ+a12−dA2−b1−A2−a310λ+bA2=0.
Expanding the determinant, we get the following:(47)Δλ=λ3+a12−a11+bA2λ2+a31dA2+a12+a12d0+bA2a12−a11λ+a12d0bA2+a12a31b2A2−1=λ3+a1λ2+a2λ+a3=0.(i)
IfA2>A− holds, then we know a3>0 from Lemma 3. Thus we have a1>0 from the assumptions a2>0 and a1a2>a3. Therefore, using the Routh–Hurwitz criterion, the positive equilibrium Ep is locally asymptotically stable.(ii)
IfA2<A−, then a3<0 from Lemma 3. Hence, the characteristic equation (47) has at least one positive real root. Therefore, the positive equilibrium p2,u2,A2 is unstable.
From (41), if a12>a11, then a2>0, and(48)a1a2−a3=bA2a12−a112+b2A22+a31d+a12+a12d0a12−a11+bA2da31+ba12a311−A2>0.
We have the following corollary.Corollary 1.
Suppose thatA2>A− holds. If a12>a11, then, the positive equilibrium Ep is locally asymptotically stable.


## 5. Effects of Weeds Removal and Other Environmental Parameters on Main Factors in Paddy Fields

From Theorem1, the coexistence equilibrium state p2,u2,A2 exists when b>d2 and (H1) hold, where A2 is one of the solutions of equation (14), and p2 and u2 are computed in (19), (21), respectively.From (14), we obtain the following:(49)dd0dA=−λb1−2Aln1−dA/c+b/rcA1−A/1−dA/cln21−dA/c.If1/2≤A≤1, then dd0/dA<0. Otherwise, by Lemma 2.1 in [29], we have the following:(50)dd0dA<−λb2A−1dA/c−dA+b/rA1−A/c−dAln21−dA/c=−λbbA2−2b−d2A+b−d2rc−dAln21−dA/c.Under the conditionb>d2, the discriminant of the quadratic function bA2−2b−d2A+b−d2 is as follows:(51)−2b−d22−4bb−d2=−4d2b−d2<0.Hence,bA2−2b−d2A+b−d2>0. It gives that dd0/dA<0. Thus, we get the following:(52)dAdd0<0.In the coexistence equilibrium state,p2 can still be calculated by (13). It is not difficult to obtain that(53)dp2dd0=bλrc−dAdAdd0<0.Therefore, the biomass of herbivores and weeds decreases monotonously with the increase of the intensity of weeds removal. Obviously, the result that weeds biomassp2 decreases monotonously with the increase of weeds removal intensity d0 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 d1 of weeds, and the half saturation concentration m of inorganic fertilizer.The expression (21) of the inorganic fertilizer content u2 can also be written as follows:(54)u2=md0+d1+d0dA2/b1−A2a−d0+d1+d0dA2/b1−A2.Notice that the herbivores biomassA2 is independent of parameters a, m, and d1, hence u2 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 d1 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 u2 and d0. Because A2 is affected by the parameter d0, we denote d0+d0dA2/b1−A2 by wd0, that is,(55)wd0=d01−1r+d0d2rb1−A2,where d0≥0. Thus, the inorganic fertilizer content u2 can be rewritten as follows:(56)u2=md1+wd0a−d1+wd0.It is easy to obtain that(57)du2dd0=maa−d1−wd02dwdd0,(58)dwdd0=1−1r+d2rb1−A2+d0d2rb1−A22dA2dd0.BecauseA2>A∗=1−d2/b, we have the following:(59)1−1r+d2rb1−A2>0.Therefore,wd0 consists of two parts, one part increases with the increase of d0 and the other part related to A2 decreases with the increase of d0. Hence, there must be a unique d0 such that dw/dd0=0. From this equation, we have the following:(60)dA2dd0=−1−A2d2−1−rb1−A2d0d2.From (14) and (49), we get the following:(61)dd0dA2=−d0λr1−A2c−dA2d0−λr1−2A2c−dA2.Multiplying (60) and (61), we obtain the following:(62)d0=λrc−dA21−2A2+d2λrc−dA2d2−1−rb1−A2.By solving simultaneous equations (14) and (62), we obtain the unique d0∗ that satisfies dw/dd0=0.Notice that(63)d2wdd02=mad2−1−rb1−A22rb1−A22a−d1−wd02>0,so d2w/dd02|d0=d0∗>0. Thus the parameter value d0∗ is a minimum point of wd0. Furthermore, from (57), it is easy to know that the nondifferentiable points of u2 (such that a−d1−wd0=0) are not extreme points. Therefore, the weeds removal intensity d0∗ minimizes the inorganic fertilizer content u2. This shows that with the increase of parameter d0, the content of inorganic fertilizer u2 first has a downward trend, then when d0=d0∗, u2 reaches the minimum, and then u2 will show an upward trend.


## 6. Hopf Bifurcation

In this section, takingm as the Hopf bifurcation parameter, we consider the existence of a Hopf bifurcation of system (4) at the coexistence equilibrium Ep. We always assume that the positive equilibrium Ep exists and A2>A− holds.From equation (14) for determining A2 and expression (19) for calculating p2, we know that A2 and p2 in equilibrium p2,u2,A2 are independent of the half saturation concentration m of inorganic fertilizer, but u2 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 a12, a21, and a22 are also related to m, while a13, a23, a31, and a33 are not related to m. From (39), we can see that the parameter a11 is not related to m either. From (41), the coefficients of the characteristic equation (47) depend on the parameter m. Let us set a1=a1m, a2=a2m and a3=a3m.Lemma3 tells us the constant term a3 of the characteristic polynomial Δλ is positive under the assumption A2>A−. Therefore, a Hopf bifurcation takes place when the real part of a pair of conjugate complex eigenvalues changes sign.We define a function(64)Lm=a1ma2m−a3m,and give the following lemma [33].Lemma 4.
If there exists a critical valuem0 of a parameter m such that a2m0>0, Lm0=0 and L′m0≠0, then system (4) undergoes a Hopf bifurcation at Ep when m passes m0.According to Lemma4, we can know that the Hopf bifurcation critical value m0 of the parameter m is the real root of equation Lm=0 which satisfies the conditions L′m0≠0 and a2m0>0.Becausea3>0 and Lm0=0, the condition a2m0>0 is equivalent to a1m0=a12m0−a11+bA2>0.From (21), we get the following:(65)m+u2=m+mb1−A2d0+d1+d0dA2b1−A2a−d0−d1−d0dA2=mab1−A2gA2.Substituting (65) into (33), by (31) we obtain the following:(66)a12=p2amg2A2mab1−A22=A2g2A2mabd01−A2.It shows thata12 is related to m. From (41), we have the following:(67)Lm=a12−a11+bA2a31dA2+a12+a12d0+bA2a12−a11−a12d0bA2−a12a31b2A2−1=a12−a11+bA2a31+d0+bA2a12+a31dA2−a11bA2−d0bA2+a31b2A2−1a12=a31+d0+bA2a122+b2A22−2a11bA2+a31dA2−a11+b−bA2−a11d0a12+bA2−a11a31dA2−a11bA2.Therefore, the conditionLm0=0 can be rewritten as follows:(68)a31+d0+bA2a122+b2A22−2a11bA2+a31dA2−a11+b−bA2−a11d0−a12+bA2−a11a31dA2−a11bA2=0.Of all the parameters of the equation (68), only a12 is related to the bifurcation parameter m. Therefore, by calculating a12 from equation (68), the bifurcation parameter m can be determined by (66).In addition, from (31), (39), and (66), we obtain that(69)L′m=2a31+d0+bA2a12+b2A22−2a11bA2+a31dA2−a11+b−bA2−a11d0a12′m=−2a31+d0+bA2a12+b2A22−2a11bA2+a31dA2−a11+b−bA2−a11d0×A2g2A2m2abd01−A2.Notice thatgA2>0, the condition L′m0≠0 can be replaced by the following condition:(70)a12≠−b2A22−2a11bA2+a31dA2−a11+b−bA2−a11d02a31+d0+bA2.According to Lemma4, we can get the existence of a Hopf bifurcation of system (4).Theorem 4.
Suppose that conditions (H1),b>d2 and A2>A− hold. If equation (68) has a positive real root a12 satisfying (70) and a12>a11−bA2, then system (4) undergoes a Hopf bifurcation at Ep when m passes m0, where(71)m0=A2g2A2aba12d01−A2.Hence, oncea12 is obtained from equation (68), we can calculate the critical value m0 with (71), which is obtained by (66).According to the conditions of Theorem4, we summarize the process of calculating the critical value m0 of the Hopf bifurcation as follows:(1)
If the conditionb>d2 holds, then A2 is obtained by solving equation (14) and p2 is calculated by (19). If the condition (H1) is validated, then system (4) has a positive equilibrium Ep.(2)
Find the rootA− of the quadratic equation (42) and verify the assumption A2>A−.(3)
Calculatea11 with (39) and a31 with (31).(4)
Solve the quadratic equation (68). If there is a positive real root, then we obtain the parameter a12.(5)
Verify the conditions (70) and a12>a11−bA2. If one of them is not valid, the critical value m0 can not be calculated.(6)
Calculate the critical valuem0 with (71).Remark 1.
If the quadratic equation (68) has two positive real roots for a12, and both of them satisfy all the conditions of Theorem 4, then there are two Hopf bifurcations at the positive equilibrium point Ep.


## 7. An Example

In system (4), take a set of simulation parameters a=40, c=6, r=0.5, λ=6, d0=0.3, d1=0.2, and d2=0.7. We consider two cases of parameter b=1 and b=4.4, which satisfy the condition b>d2.By the process of calculating the critical valuem0 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 condi1 in the table and a=40, it follows that hypothesis (H1) holds, so there is a positive equilibrium point in system (4) by Theorem 1. It is easy to see from the second and seventh rows of Table 1 that A2>A−. From the second line to the fourth line counted backward in the table, we obtain a12>a11−bA2 and a12≠condi2. Therefore, the conditions in Theorem 4 are all satisfied.Table 1
The calculation results of the Hopf bifurcation.
Itemsb=1b=4.4A20.98680.9970p20.04330.0433condi1(1)31.783132.1058A−0.49110.4945a113.89073.9316a3113.696213.8371a12(2)3.61660.0156, 2.7721a11−bA22.9038−0.4553condi2(3)1.40711.3939m00.0202244.318921, 0.024358(1)condi1=d0+d1+d0dA2/b1−A2. (2)a12 is the positive real roots of equation (68). (3)condi2 is the right side of (70).From Table1, we see that the critical value m0≈0.020224 when parameter b=1. Therefore, system (4) undergoes a Hopf bifurcation at Ep when m passes m0≈0.020224. If we let m=0.02<m0≈0.020224, system (4) has a positive equilibrium Ep≈0.0433,0.0774,0.9868. From (41), we get a1≈0.7533, a2≈69.7690>0, and a3≈49.8515. Hence, a1a2−a3≈2.7039>0. According to Theorem 3, the positive equilibrium Ep≈0.0433,0.0774,0.9868 is locally asymptotically stable (see Figure 1). Therefore, if m passes m0≈0.020224 from 0.02, then the positive equilibrium Ep loses its stability, and system (4) has a periodic solution (see Figure 2, where m=0.021).Figure 1
The phase diagram ofu−A of system (4) with a=40, c=6, λ=6, r=0.5, b=1, d0=0.3, d1=0.2, d2=0.7 and m=0.02. It describes the asymptotic stability of the equilibrium 0.0433,0.0774,0.9868.Figure 2
The phase diagram ofu−A of system (4) with a=40, c=6, λ=6, r=0.5, b=1, d0=0.3, d1=0.2, d2=0.7, and m=0.021. It describes the instability of the positive equilibrium and indicates a positive periodic solution exists in the system (4).When the parameterb=4.4, from Table 1, we still know that system (4) has two critical values of the Hopf bifurcation m0≈0.024358 or m0≈4.318921. Therefore, system (4) undergos a Hopf bifurcation at Ep when m passes m0≈0.024358. It also undergoes the other Hopf bifurcation when m passes m0≈4.318921. If we let m=0.02<m0≈0.024358, system (4) has a positive equilibrium Ep≈0.0433,0.0813,0.9970. From (41), we get a1≈3.8315, a2≈64.3037>0, and a3≈208.7775. Thus, we have a1a2−a3≈37.6027>0. Therefore, according to Theorem 3, the positive equilibrium Ep≈0.0433,0.0813,0.9970 is locally asymptotically stable, which means there exists a periodic solution of system (4) when m>m0≈0.024358. If we let m=4.32>m0≈4.318921, system (4) has a positive equilibrium Ep≈0.0433,17.5696,0.9970. From (41), we get a1≈0.4710, a2≈2.0527>0, and a3≈0.9666. Then, we have a1a2−a3≈1.9931×10−4>0. Therefore, according to Theorem 3, the positive equilibrium Ep≈0.0433,17.5696,0.9970 is locally asymptotically stable, implying that there exists a periodic solution of system (4) when m<m0≈4.318921.Combining the above calculation results whenb=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. This is consistent with our calculation in the above.Figure 3
The relationship curve between the maximum real part of the eigenvalue of the characteristic equation (47) and parameter m with a=40, c=6, λ=6, r=0.5,b=4.4,d0=0.3, d1=0.2 and d2=0.7.


## 8. 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, whenb>d2, if the coupling effect coefficient of light and fertilizer a is higher than the sum of weeds removal intensity d0, weeds mortality d1 , and d0dA2/b1−A2 (assumption (H1)), the coexistence equilibrium Ep can be found in the system.Our results are obtained under the condition ofd0>0. When weeds are not removed, i.e., d0=0, the situation is relatively simple. Xiang et al. has considered the case of d0=0 in [29]. They obtained that p2 and u2 as follows:(72)p2=−1λln1−d2rc,u2=mrp2d1+d2rp2a−d1−d2.Obviously, they need conditionsrc>d2 and p2>d2/ra−d1. In Theorem 1, we replaced their condition p2>d2/ra−d1 with (H1) and deleted condition rc>d2. Thus, our results include the case of rc≤d2. Because conditions (H1) and (25) are equivalent, from (25), we know that the intensity of weeds removal d0 should not be too high d0<a−d1.By using the Routh–Hurwitz criterion, we obtain the conditions for local asymptotic stability of the positive equilibriumEp, in which the condition A2>A− is used to guarantee the coefficient a3>0. If d0=0, then we have A2=1 from (6), so the coefficient a3=a12a31b>0 is natural. So there is no requirement for the lower bound of A2 in Theorem 2.5 for the local stability of positive equilibrium in [29]. However, if weeds removal is considered in the system, then whether a3 is positive or negative is affected by the value A2. From (41), we see that A2 should not be too small, otherwise a3 will be negative. This is why we propose the condition A2>A− in Theorem 3.When the coefficientsai of cubic characteristic polynomial (47) are positive, the Hopf bifurcation may be generated when a1a2=a3. According to this principle, we obtain the Hopf bifurcation conditions of system (4) at the positive equilibrium Ep, and give the Hopf bifurcation critical value formula (71) with m being the bifurcation parameter. The half saturation concentration of inorganic fertilizer is a key parameter in Michaelis–Menten uptake kinetics, which is affected by species of weeds. A Hopf bifurcation appears at a positive equilibrium Ep when the parameter m passes through the critical value m0, 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 Section5, 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 d0∗, the soil inorganic fertilizer content is the lowest. The weeds removal intensity is too small to achieve the purpose of weeding. Therefore, the intensity of weeds removal should be greater than the minimum point d0∗ as far as possible. But considering that d0 is restricted by condition (H1), it has an upper bound. For example, we can know from (25) that d0 cannot exceed a−d1 at least. If the minimum point d0∗ is greater than the upper bound, then we should try to use low-intensity weeding. If d0∗ is less than the upper bound, the intensity of weeds removal should be close to the upper bound and not exceed it.



---

*Source: 8716314-2020-11-11.xml*