- Research
- Open Access
- Published:
Modelling the population dynamics of brown planthopper, Cyrtorhinus lividipennis and Lycosa pseudoannulata
Advances in Difference Equations volume 2019, Article number: 265 (2019)
Abstract
In this work, a mathematical model is proposed to investigate the population dynamics of brown planthopper, which is a major insect pest of rice, and its natural enemies namely Cyrtorhinus Lividipennis and Lycosa Pseudoannulata. To analyse the model theoretically, the geometric singular perturbation method is employed. Conditions that differentiate dynamic behaviors exhibited by the model are derived. To illustrate our theoretical predictions, computer simulations are also presented showing that the dynamic behavior exhibited by the model correspond to the reported field observations.
Introduction
Rice is considered as the major staple food for over 50% of the population of the world especially for those living in the Asian–Pacific region [1]. It is a rich source of carbohydrates and a range of nutrients. The world’s rice production is mainly produced in the Asian–Pacific region including Thailand [2]. In rice production, the estimates for the averages of the potential losses worldwide is approximately 37, 25 and 13% due to weeds, animal pests and pathogens, respectively [3].
One of the major insect pests of rice is brown planthopper (BPH). It sucks sap from the leaves and lays egg masses in the leaf blade and the leaf sheath to blocks the xylem and phloem [4]. The leaves of the infested rice plants will turn yellow and the rice plants will be dried up and die. The damage is the symptom of hopperburn, it begins in patches and spreads rapidly as BPH moves from dying plants to the others. In addition, virus diseases such as grassy stunt, ragged stunt and wilted stunt can be transmitted from a rice plant to the others by BPH as well [4, 5]. When the outbreak of BPH occurs, the rice yield will be decreased leading to large economic losses.
Although insecticides have been widely used for controlling the pest, BPH has developed their resistance to some major insecticides such as carbamates, organophosphates, neonicotinoids, phenylpyrazoles and pyrethroids [4, 5]. As it is quick, easy to use and cost-effective against insects, chemical control is a popular choice in pest management. However, excessive and irrational use of chemical pesticides could lead to negative effects in the environment such as biodiversity’s reduction and the decrease in population of natural enemies. Alternatively, biological control is a safe and an effective method. Cyrtorhinus Lividipennis and Lycosa Pseudoannulata are two important natural enamies of BPH. Investigating biological control of BPH requires an understanding of the life cycles, fecundity and consumption behavior of BPH, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata.
Cyrtorhinus lividipennis is one of natural enemies of BPH, which mainly preys on eggs and nymphs of BPH [6, 7]. The predatory activity of Cyrtorhinus lividipennis against BPH has been investigated by many researcher and the study indicated that the Cyrtorhinus lividipennis’s preying on BPH’s eggs was an important cause of the decrease in BPH population [6, 8]. A nymph of Cyrtorhinus lividipennis consumes approximately 7.5 BPH’s eggs or 1.4 adult BPH per day for a period of 14 days. An adult Cyrtorhinus lividipennis consume approximately 10.2 BPH’s eggs or 4.7 BPH’s nymphs or 2.4 adult BPH per day for a period of 10 days [8]. However, the population of Cyrtorhinus lividipennis do not reproduce rapidly enough to control an infestation of BPH in a paddy field.
Lycosa pseudoannulata is also a natural enemy of BPH [9, 10]. It plays an important role in controlling BPH populations. It was found that in a 14-day period one Lycosa pseudoannulata could consume approximately 17 BPH’s nymphs or 15–20 adult BPH per day [11]. It inhabits at the lower part of rice plants during the daytime to prey on BPH and moves to the middle and upper sections during the night time to prey on leafhoppers [9, 10, 12].
Therefore, in order to control the population of brown planthoppers which are insect pest of rice efficiently by using two biological control agents, Cyrtorhinus lividipennis and Lycosa pseudoannulata, we should start with studying their population dynamics. In this paper, a mathematical model is proposed to investigate the population dynamics of BPH and its natural enemies, Cyrtorhinus lividipennis and Lycosa pseudoannulata. The analyses of our model will be carried out both theoretically and numerically.
A mathematical model
In what follows, let \(X(t)\) denote the population density of BPH at time t, \(Y(t)\) denote the population density of Cyrtorhinus Lividipennis at time t, and \(Z(t)\) denote the population density of Lycosa Pseudoannulata at time t.
Firstly, the dynamics of the population of BPH is assumed to follow the following equation:
Consider the right-hand side of (1). The first term accounts for the reproduction rate of BPH, while \(a_{1}\) is the intrinsic growth rate of BPH and \(k_{1}\) is the carrying capacity of BPH. The functional responses of Cyrtorhinus Lividipennis and Lycosa Pseudoannulata feeding on BPH were reported to be Holling’s type II [13, 14]. The second and the third terms then account for the death rate of BPH due to the predation of Cyrtorhinus Lividipennis and Lycosa Pseudoannulata, respectively, while α is the attacking rate of Cyrtorhinus Lividipennis to BPH, β is the attacking rate of Lycosa Pseudoannulata to BPH, \(h_{1}\) is the handling time of Cyrtorhinus Lividipennis to BPH, and \(h_{2}\) is the handling time of Lycosa Pseudoannulata to BPH. The last term account for the natural death rate of BPH.
Next, the dynamics of the population of Cyrtorhinus Lividipennis is assumed to follow the following equation:
Consider the right-hand side of (2). The first term on accounts for the reproduction rate of Cyrtorhinus Lividipennis, while \(a_{2}\) and \(k_{2}\) are the intrinsic growth rate and the carrying capacity of Cyrtorhinus Lividipennis, respectively. The functional responses of Cyrtorhinus Lividipennis feeding on BPH were reported to be Holling’s type II [13].The second term then accounts for the reproduction rate of Cyrtorhinus Lividipennis due to the predation on BPH, while \(b_{1}\) is the conversion rate of BPH to Cyrtorhinus Lividipennis. The last term accounts for the natural death rate of Cyrtorhinus Lividipennis.
Finally, the dynamics of the population of Lycosa Pseudoannulata is assumed to follow the following equation:
Consider the right-hand side of (3). The first term accounts for the reproduction rate of Lycosa Pseudoannulata, while \(a_{3}\) and \(k_{3}\) are the intrinsic growth rate and the carrying capacity of Lycosa Pseudoannulata, respectively. The functional responses of Lycosa Pseudoannulata feeding on BPH were reported to be Holling’s type II [14].The second term then accounts for the reproduction rate of Lycosa Pseudoannulata due to the predation on BPH, while \(b_{2}\) and \(h_{2}\) are the conversion rate of BPH to Lycosa Pseudoannulata and the handing time of Lycosa Pseudoannulata to BPH, respectively. The last term accounts for the natural death rate Lycosa Pseudoannulata.
Therefore, our model consists of Eqs. (1)–(3), where all parameters are assumed to be positive.
Model analysis
BPH is a small brownish insect pest. Females start laying eggs which are white and elongated along with the midribs of the leaf sheath and leaf blade. The average number of eggs laid by a female is about 300. Eggs hatch in about 4–8 days into nymph period. The nymph then undergoes five instars for a period of 15–20 days to grow up into adult which can live for 10–20 days. Therefore, the life cycle is completed in about 30–40 days [15].
The green miridbug, Cyrtorhinus lividipennis, is an important natural enemy of BPH and mainly preys on BPH eggs and young nymphs. Adult females lay eggs either singly or in groups within the leaf sheath and the average number of eggs laid by a female is about 93. The incubation period of Cyrtorhinus lividipennis ranges about 6–9 days and then it grow up into nymphal. Nymphal stages undergo four instars for a period of 10–17 days and they develop into adult. The longevity for females range 5–21 days and the males range 7–25 days [16]. Therefore, the total life cycle of green miridbug is about 21–37 days.
The wolf spider, Lycosa pseudoannulata, is one of the predominant spiders in paddy fields and is an important predator of BPH. It has a fork-shaped mark on the back and the abdomen has white markings. This wolf spider can play a major role in keeping down BPH populations. It was found that in a 14-day period each wolf spider killed an average of 17 BPH nymphs/day and it also could kill about 15-20 adult BPH/day [11]. An adult female of Lycosa pseudoannulata lays about 30 eggs. A female wolf spider carries her egg sac under abdomen and after hatching, the spiderlings climb on their mother’s back. Lycosa pseudoannulata inhabits around on water surface and the lower part of rice plants to catch its prey directly without creating web. In 1973, Gavarra et al. [17] reported that Lycosa pseudoannulata egg stage lasted for 59 days and the young spiders took about 170 days to reach maturity. The life cycle of this wolf spider lasted for an average of 116.3 days form egg to egg but the average generation time is about 263.9 days from egg to adult death.
Hence, we assume in what follows that the dynamics of BPH population is fastest while the dynamics of Cyrtorhinus lividipennis population and Lycosa pseudoannulata population are intermediate and slowest, respectively.
To analyse our model of equations (1)–(3) by the geometric singular perturbation method [18,19,20], ε and δ which are the small dimensionless positive parameters will be used to scale the variables and parameters of the system. Letting \(x=X\), \(y=Y\), \(z=Z\), \(c_{1}=a_{1}\), \(c_{2}=\alpha \), \(c_{3} =\beta \), \(c _{4}=\frac{a_{2}}{\varepsilon }\), \(c_{5}=\frac{a_{3}}{\varepsilon \delta }\), \({ \gamma }_{1}=\frac{b_{1}}{\varepsilon }\), \({\gamma }_{2}=\frac{b _{2}}{\varepsilon \delta }\), \(e_{1}=d_{1}\), \(e_{2}=\frac{d_{2}}{\varepsilon }\), \(e_{3}= \frac{d_{3}}{\varepsilon \delta }\), and we obtain
During transitions, when the right-hand sides of Eqs. (4)–(6) are finite and nonzero, \(\vert \dot{y} \vert \) will be of the order ε and \(\vert \dot{z} \vert \) will be of the order εδ represented by the notations \(\dot{x}=O (1 )\), \(\dot{y}=O (\varepsilon )\) and \(\dot{z}=O (\varepsilon \delta )\), respectively, in what follows.
Next, we will show that, for suitable parametric values with the sufficiently small ε and δ, the manifolds \(\{F (x,y,z )=0 \}\), \(\{G (x,y,z )=0 \}\) and \({H (x,y,z )=0}\) are as shown in Figs. 1–3.
Manifold \(\lbrace F=0 \rbrace \)
Since
the manifold \(\{F=0\}\) is composed of the trivial manifold \(x=0\) and the nontrivial manifold
Let us consider the intersection of the nontrivial manifold and the \((x,y)\)-plane. The intersection of the nontrivial manifold and the \((x,y)\)-plane occurs along the curve
and the intersection of the nontrivial manifold and the y-axis occurs at the point for which \(x = 0\) and
Moreover, the intersection of the nontrivial manifold and the x-axis occurs at the point for which \(y = 0\) and
Note that \(x_{1} > 0\) and \(y_{1} > 0\) if and only if
Since
\(T' (x ) = 0\) at the point where
Moreover, \(T'' (x_{2} ) <0\). Therefore, the nontrivial manifold has a relative maximum at the point where \(x = x_{2}\), \(y = T(x_{2}) \equiv y_{2} \) on the \((x,y)\)-plane where we note that \(x_{2} >0\) if and only if
Now, let us consider the intersection of the nontrivial manifold and the \((x,z)\)-plane. The intersection of the nontrivial manifold and the \((x,z)\)-plane occurs along the curve
which intersects the z-axis at the point where \(x = 0\) and
Since
\(U' (x ) = 0\) at the point where
Moreover, \(U'' (x_{3} ) <0\). Therefore, the nontrivial manifold has a relative maximum at the point where \(x = x_{3}\), \(z = U(x_{3}) \equiv z_{3} \) on the \((x,z)\)-plane where we note that \(x_{3} >0\) if and only if
Next, let us consider the intersection of the nontrivial manifold and the \((y,z)\)-plane. The intersection of the nontrivial manifold and the \((y,z)\)-plane occurs along the line
which intersects the y-axis at the point for which \(z=0\) and \(y=y_{1}\). In addition, it intersects the z-axis at the point for which \(y = 0\) and \(z=z_{1}\).
Consider the nontrivial manifold
in the first octant. We observe that \(\frac{\partial R(x,y)}{\partial y}<0\) and
Note that \(\frac{\partial R(x,y)}{\partial x}=0\) along the curve
Manifold \(\lbrace G=0 \rbrace \)
Since
the manifold \(\{G=0\}\) is composed of the trivial manifold \(y = 0\) and the nontrivial manifold
Here, the nontrivial manifold is independent of z and parallel to the z-axis. The intersection of the nontrivial manifold and the y-axis occurs at the point for which \(x=0\) and
Note that \(y_{3}>0\) if and only if
Since \(V^{\prime } (x )=\frac{{\gamma }_{1}c_{2}c_{4}}{{ (e_{2}+ (c_{2}e_{2}h_{1}-{\gamma }_{2}c_{2} )x )}^{2}}>0\), \(V (x )\) is an increasing function of x.
In addition,
Note that \(y_{4}>0\) if
and
The intersection of \(\{F=0 \}\) and \(\{G=0 \}\) is composed of \(\{x=0,y=0\}\), \(\{x=0,y=y_{3}\}\), \(\{y=0,z=U(x)\}\) and \(\{y=V(x),z=R(x,y)\}\). Note that \(\{y=V(x),z=R(x,y)\}\) intersects the \((x,y)\)-plane at the point where \(x=x_{4}\), \(y=T(x_{4})\equiv y_{5}\) and \(z=0\).
Manifold \(\lbrace H=0 \rbrace \)
Since
the manifold \(\{H=0\}\) is composed of the trivial manifold \(z =0\), and the nontrivial manifold
which is independent of y and parallel to the y-axis. The intersection of the nontrivial manifold and the x-axis occurs at the point for which \(z=0\) and
Note that \(x_{5}>0\) if
and
Since
\(W (z )\) is an increasing function of z in the first octant.
In addition,
Note that \({\mathrm{x}}_{\mathrm{6}}>0\) if and only if
Next, the intersection of \(\{F=0 \}\) and \(\{H=0 \}\) is composed of \(\{x=0, z=0\}\), \(\{x=0, z=z_{4} \}\), \(\{y=T(x), z=0\}\) where \(z_{4}=\sqrt{ (\frac{c_{5}}{e_{3}}-k _{3} )}\) and \(\{x=W(z), z=R(x,y)\}\). Moreover, \(\{x=W(z), z=R(x,y) \}\) intersects the \((x,z)\)-plane where \(y=0\) at the point for which \(x=x_{7}\) and \(z=U(x_{7})\equiv {z}_{7}\).
Hence, all six possible equilibrium points of the system of Eqs. (4)–(6) can be listed as follows:
\(S_{0}= (0,0,0 )\),
\(S_{1}= (x_{1},0,0 )= (\frac{ (c_{1}-e_{1} )k _{1}}{c_{1}},0,0 )\),
\(S_{2}= (0,y_{3},0 )= (0,\frac{c_{4}-e_{2}k_{2}}{e _{2}},0 )\),
\(S_{3}= (x_{4},y_{5},0 )\) (the intersection of \(y=V(x)\) and \(z=R(x,y)\) on the \((x,y)\)-plane),
\(S_{4}= (x_{7},0,z_{7} )\) (the intersection of \(x=W(z)\) and \(z=R(x,y)\) on the \((y,z)\)-plane),
and \(S_{5}= (x_{8},y_{8},z_{8} )\) (the intersection of \(x=W(z)\), \(y=V(x)\), and \(z=R(x,y)\)).
The different shapes and positions of the three manifolds \(\{F=0 \}\), \(\{G=0 \}\) and \(\{H=0 \}\) can lead to different dynamic behaviors of the solution of the system of Eqs. (4)–(6). Hence, we identify the possible cases according to the shapes and positions of the three manifolds. Here, only the five possible cases are stated since the other possible cases lead to similar dynamic behaviors occurred in these five cases.
Theorem 1
If the inequalities
hold, the system of Eqs. (4)–(6) will have a periodic solution and the equilibrium points, \(S_{0}\), \(S_{1}\), \(S_{2}\), \(S_{3}\), \(S_{4}\), \(S_{5}\), are all unstable, provided ε and δ are sufficiently small.
Proof
Since the inequalities (33) and (34) hold, the positions and shapes of the manifolds are as shown in Fig. 1. In what follows, we indicate the slow, intermediate and fast transitions by one, two and three arrows, respectively. In this case, there are six equilibrium points, \(S_{0}\), \(S_{1}\), \(S_{2}\), \(S_{3}\), \(S_{4}\) and \(S_{5}\), in the first octant.
Starting from a point \(I=(x_{0},y_{0},z_{0})\), with \(F(x_{0},y_{0},z _{0})\neq 0\) in Fig. 1. Here, the position of the point I is located in front of the nontrivial manifold \(\{F=0\}\) for which \(F < 0\) here. The solution trajectory then moves from the point I to the point J on the nontrivial manifold \(\{F=0\}\) with the fast transition parallel to the x-axis in the direction that x decreases. The solution trajectory then moves along the nontrivial manifold \(\{F=0\}\) from the point J located on the right-hand side of the nontrivial manifold \(\{G=0\}\) where \(G<0\) to the point K on the curve \(\{F=G=0\}\) located in front of the nontrivial manifold \(\{H=0\}\) with intermediate transition parallel to the y-axis in the direction that y decreases. At the point K, \(H>0\) and hence the solution trajectory then moves along the curve \(\{F=G=0\}\) from the point K in the direction that z increases to the point L where the stability of the manifold is lost and the solution trajectory then jumps to the point M on the straight line \(\{x=0, y=y_{3}\}\) on the \((y,z)\)-plane in which \(F=G=0\) with the fast transition parallelled to the x-axis in the direction that x decreases because \(F<0\) here. Here, \(H<0\), the solution trajectory then moves from the point M to the point N along the straight line \(\{x=0, y=y_{3}\}\) on the \((y,z)\)-plane with the slow transition in the direction that z decreases. The point N is located below the nontrivial manifold \(\{F=0\}\), the stability of the manifold is lost again and the solution trajectory then jumps to the point O on the curve \(\{F=G=0\}\) with the fast transition parallel to the x-axis in the direction that x increases because \(F>0\) below nontrivial manifold \(\{F=0\}\). At the point O, \(H>0\) and hence the solution trajectory then moves along the curve \(\{F=G=0\}\) from the point O in the direction that z decreases to the point L where the stability is lost and the solution trajectory then jumps to the point M. The solution trajectory then moves to the point N and O forming a closed cycle \(OLMNO\) and hence a limit cycle occurs.
By starting at other initial point that closes to each equilibrium point, the local stability for each equilibrium point can be determined in a similar manner and the proof is complete. □
Theorem 2
If the inequalities
hold, the solution of the system of Eqs. (4)–(6) tends toward a stable equilibrium point \(S_{3}\) as time passes and the equilibrium points, \(S_{0}\), \(S_{1}\) and \(S_{2}\), are unstable provided ε and δ are sufficiently small.
Proof
Since the inequalities (35) and (36) hold, the positions and shapes of the manifolds are as shown in Fig. 2. In this case, there are four equilibrium points, \(S_{0}\), \(S_{1}\), \(S_{2}\) and \(S_{3}\), in the first octant.
Starting from a generic point \(I=(x_{0},y_{0},z_{0})\) close to \(S_{3}\), with \(F(x_{0},y_{0},z_{0})\neq 0\) in Fig. 2. The position of the point I is located in front of the nontrivial manifold \(\{F=0\}\) for which \(F < 0\) here. The solution trajectory then moves from the point I to the point J on the nontrivial manifold \(\{F=0\}\) with the fast transition parallel to the x-axis in the direction that x decreases. The solution trajectory then moves along the nontrivial manifold \(\{F=0\}\) from the point J located on the right-hand side of the nontrivial manifold \(\{G=0\}\) where \(G<0\) to the point K on the curve \(\{F=G=0\}\) located behind the nontrivial manifold \(\{H=0\}\) with intermediate transition parallel to the y-axis in the direction that y decreases. At the point K, \(H<0\) and hence the solution trajectory then moves along the curve \(\{F=G=0\}\) from the point K in the direction that z decreases to the equilibrium point \(S_{3}\). Therefore, the solution trajectory in this case tends toward the stable equilibrium point \(S_{3}\) as time passes.
By starting at other initial point that closes to each equilibrium point, the local stability for each equilibrium point can be determined in a similar manner and the proof is complete. □
Theorem 3
If the inequalities
hold, the solution of the system of Eqs. (4)–(6) tends toward a stable equilibrium point \(S_{5}\) as time passes and the equilibrium points, \(S_{0}\), \(S_{1}\), \(S_{2}\), \(S_{3}\) and \(S_{4}\), are unstable provided ε and δ are sufficiently small.
Proof
Since the inequalities (37) and (38) hold, the positions and shapes of the manifolds are as shown in Fig. 3. In this case, there are six equilibrium points, \(S_{0}\), \(S_{1}\), \(S_{2}\), \(S_{3}\), \(S_{4}\) and \(S_{5}\), in the first octant.
Starting from a generic point \(I=(x_{0},y_{0},z_{0})\), with \(F(x_{0},y_{0},z_{0})> 0\) in Fig. 3. The solution trajectory then moves from the point I to the point J on the nontrivial manifold \(\{F=0\}\) with the fast transition parallel to the x-axis in the direction that x decreases. The solution trajectory then moves along the nontrivial manifold \(\{F=0\}\) from the point J located on the right-hand side of the nontrivial manifold \(\{G=0\}\) where \(G<0\) to the point K on the curve \(\{F=G=0\}\) located in front of the nontrivial manifold \(\{H=0\}\) with intermediate transition parallel to the y-axis in the direction that y decreases. At the point K, \(H>0\) and hence the solution trajectory then moves along the curve \(\{F=G=0\}\) from the point K in the direction that z increases to the equilibrium point \(S_{5}\) with intermediate transition since \(G>0\) here. Therefore, the solution trajectory in this case tends toward the stable equilibrium point \(S_{5}\) as time passes.
By starting at other initial point that closes to each equilibrium point, the local stability for each equilibrium point can be determined in a similar manner and the proof is complete. □
Theorem 4
If the inequalities
hold, the solution of the system of Eqs. (4)–(6) tends toward a stable equilibrium point \(S_{2}\) as time passes and the equilibrium points, \(S_{0}\) and \(S_{1}\), are unstable provided ε and δ are sufficiently small.
Proof
Since the inequalities (39) and (40) hold, the positions and shapes of the manifolds are as shown in Fig. 4. In this case, there are three equilibrium points, \(S_{0}\), \(S_{1}\) and \(S_{2}\), in the first octant.
Starting from a generic point \(I=(x_{0},y_{0},z_{0})\), with \(F(x_{0},y_{0},z_{0})> 0\) in Fig. 4. The solution trajectory then moves from the point I to the point J on the nontrivial manifold \(\{F=0\}\) with the fast transition parallel to the x-axis in the direction that x increases. The solution trajectory then moves along the nontrivial manifold \(\{F=0\}\) from the point J located on the left-hand side of the nontrivial manifold \(\{G=0\}\) where \(G>0\) to the point K where the stability of the manifold is lost. The solution trajectory then jumps to the point L on the other stable branch \(\{F=0\}\) with the fast transition parallelled to the x-axis in the direction that x decreases because \(F<0\) here. The solution trajectory then moves along the curve \(\{F=H=0\}\) from the point L in the direction that y increases with intermediate transition to the equilibrium point \(S_{2}\) since \(G>0\) here. Therefore, the solution trajectory in this case tends toward the stable equilibrium point \(S_{2}\) as time passes.
By starting at other initial point that closes to each equilibrium point, the local stability for each equilibrium point can be determined in a similar manner and the proof is complete. □
Theorem 5
If the inequalities
hold, the solution of the system of Eqs. (4)–(6) tends toward a stable equilibrium point \(S_{2}\) as time passes and the equilibrium points, \(S_{0}\), \(S_{1}\) and \(S_{4}\), are unstable provided ε and δ are sufficiently small.
Proof
Since the inequalities (41) and (42) hold, the positions and shapes of the manifolds are as shown in Fig. 5. In this case, there are four equilibrium points, \(S_{0}\), \(S_{1}\), \(S_{2}\) and \(S_{4}\), in the first octant.
Starting from a generic point \(I=(x_{0},y_{0},z_{0})\), with \(F(x_{0},y_{0},z_{0})\neq 0\) in Fig. 5. The position of the point I is located in front of the nontrivial manifold \(\{F=0\}\) for which \(F < 0\) here. The solution trajectory then moves from the point I to the point J on the nontrivial manifold \(\{F=0\}\) with the fast transition parallel to the x-axis in the direction that x decreases. The solution trajectory then moves along the nontrivial manifold \(\{F=0\}\) from the point J located on the left-hand side of the nontrivial manifold \(\{G=0\}\) where \(G>0\) to the point K on the curve \(\{F=H=0\}\) with intermediate transition parallel to the y-axis in the direction that y increases. The solution trajectory then moves along the curve \(\{F=H=0\}\) from the point K in the direction that y increases to the point L and then the point M where the stability of the manifold is lost. The solution trajectory then jumps to the point N on the other stable branch \(\{F=H=0\}\) with the fast transition parallelled to the x-axis in the direction that x decreases because \(F<0\) here. The solution trajectory then moves along the curve \(\{F=H=0\}\) from the point N in the direction that y increases with intermediate transition to the equilibrium point \(S_{2}\) since \(G>0\) here. Therefore, the solution trajectory in this case tends toward the stable equilibrium point \(S_{2}\) as time passes.
By starting at other initial point that closes to each equilibrium point, the local stability for each equilibrium point can be determined in a similar manner and the proof is complete. □
Computer simulations
To illustrate our theoretical results, numerical simulations of the system of Eqs. (4)–(6) are carried out and presented in Figs. 6–12. The parametric values of \(c_{2}\), \(C_{3}\), \(h_{1}\) and \(h_{2}\) are obtained from the literature [13, 14, 21] while other parametric values are chosen to satisfy Theorems 1–5.
Figure 6 shows a simulation result of the system of Eqs. (4)–(6) with the parametric values \(c_{1}=0.9~\mathrm{day}^{-1}\), \(c_{2}=0.95~\mathrm{day} ^{-1}\), \(c_{3}=0.239~\mathrm{day}^{-1}\), \(c_{4}=0.91~\mathrm{day}^{-1}\), \(c_{5}=0.49~\mathrm{day}^{-1}\), \(k_{1}=67\), \(k_{2}=0.02\), \(k_{3}=0.9\), \(h_{1}=0.024\), \(h _{2}=0.063\), \(e_{1}=0.00001~\mathrm{day}^{-1}\), \(e_{2}=0.95~\mathrm{day}^{-1}\), \(e_{3}=0.55~\mathrm{day}^{-1}\), \(\gamma _{1}=0.01~\mathrm{day}^{-1}\), \(\gamma _{2}=0.085~\mathrm{day}^{-1}\), \(\epsilon =0.0042\) and \(\delta =0.822\) where \(x(0)=67\), \(y(0)=0.95\) and \(z(0)=0.03\) in which all the conditions in Theorem 1 are satisfied. The solution trajectory projected onto the \((x,y)\)-plane, \((x,z)\)-plane, \((y,z)\)-plane and the time courses of the population densities of BPH, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata populationare are also shown in Fig. 6. We can see that the solution trajectory tends to a limit cycle as predicted in Theorem 1.
Figure 7 shows a simulation result of the system of Eqs. (4)–(6) with the parametric values \(c_{1}=0.9~\mathrm{day}^{-1}\), \(c_{2}=0.95~\mathrm{day} ^{-1}\), \(c_{3}=0.239~\mathrm{day}^{-1}\), \(c_{4}=0.5~\mathrm{day}^{-1}\), \(c_{5}=0.10~\mathrm{day}^{-1}\), \(k_{1}=67.0\), \(k_{2}=0.81\), \(k_{3}=0.9\), \(h_{1}=0.024\), \(h _{2}=0.063\), \(e_{1}=0.0001~\mathrm{day}^{-1}\), \(e_{2}=0.42~\mathrm{day}^{-1}\), \(e_{3}=0.4~\mathrm{day}^{-1}\), \(\gamma _{1}=0.01~\mathrm{day}^{-1}\), \(\gamma _{2}=0.03~\mathrm{day}^{-1}\), \(\epsilon =0.05\) and \(\delta =0.9\) where \(x(0)=67\), \(y(0)=0.95\) and \(z(0)=1\) in which all the conditions in Theorem 2 are satisfied. The solution trajectory projected onto the \((x,y)\)-plane, \((x,z)\)-plane, \((y,z)\)-plane and the time courses of the population densities of BPH, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata populationare are also shown in Fig. 7. We can see that the solution tends to a stable equilibrium as time passes as predicted in Theorem 2.
Figure 8 shows a simulation result of the system of Eqs. (4)–(6) with the parametric values \(c_{1}=0.9~\mathrm{day}^{-1}\), \(c_{2}=0.95~\mathrm{day} ^{-1}\), \(c_{3}=0.239~\mathrm{day}^{-1}\), \(c_{4}=0.8~\mathrm{day}^{-1}\), \(c_{5}=0.445~\mathrm{day}^{-1}\), \(k_{1}=68\), \(k_{2}=0.02\), \(k_{3}=0.54\), \(h_{1}=0.024\), \(h _{2}=0.063\), \(e_{1}=0.0001~\mathrm{day}^{-1}\), \(e_{2}=0.95~\mathrm{day}^{-1}\), \(e_{3}=0.98~\mathrm{day}^{-1}\), \(\gamma _{1}=0.01~\mathrm{day}^{-1}\), \(\gamma _{2}=0.062~\mathrm{day}^{-1}\), \(\epsilon =0.004\) and \(\delta =0.9\) where \(x(0)=67\), \(y(0)=0.95\) and \(z(0)=0.01\), in which all the conditions in Theorem 3 are satisfied. The solution trajectory projected onto the \((x,y)\)-plane, \((x,z)\)-plane, \((y,z)\)-plane and the time courses of the population densities of BPH, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata populationare are also shown in Fig. 8. We can see that the solution tends to a stable equilibrium as time passes as predicted in Theorem 3.
Figure 9 shows a simulation result of the system of Eqs. (4)–(6) with the parametric values \(c_{1}=0.2~\mathrm{day}^{-1}\), \(c_{2}=0.95~\mathrm{day} ^{-1}\), \(c_{3}=0.239~\mathrm{day}^{-1}\), \(c_{4}=0.85~\mathrm{day}^{-1}\), \(c_{5}=0.1~\mathrm{day}^{-1}\), \(k_{1}=67\), \(k_{2}=0.5\), \(k_{3}=0.7\), \(h_{1}=0.024\), \(h _{2}=0.063\), \(e_{1}=0.001~\mathrm{day}^{-1}\), \(e_{2}=0.06~\mathrm{day}^{-1}\), \(e_{3}=0.4~\mathrm{day}^{-1}\), \(\gamma _{1}=0.001~\mathrm{day}^{-1}\), \(\gamma _{2}=0.03~\mathrm{day}^{-1}\), \(\epsilon =0.1\) and \(\delta =0.9\) where \(x(0)=0.1\), \(y(0)=0.1\) and \(z(0)=0.1\) in which all the conditions in Theorem 4 are satisfied. The solution trajectory projected onto the \((x,y)\)-plane, \((x,z)\)-plane, \((y,z)\)-plane and the time courses of the population densities of BPH, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata populationare are also shown in Fig. 9. We can see that the solution tends to a stable equilibrium as time passes as predicted in Theorem 4.
Figure 10 shows a simulation result of the system of Eqs. (4)–(6) with the parametric values \(c_{1}=0.2~\mathrm{day}^{-1}\), \(c_{2}=0.95~\mathrm{day} ^{-1}\), \(c_{3}=0.239~\mathrm{day}^{-1}\), \(c_{4}=0.6~\mathrm{day}^{-1}\), \(c_{5}=0.01~\mathrm{day}^{-1}\), \(k_{1}=67.0\), \(k_{2}=0.51\), \(k_{3}=0.02\), \(h_{1}=0.024\), \(h _{2}=0.063\), \(e_{1}=0.001~\mathrm{day}^{-1}\), \(e_{2}=0.8~\mathrm{day}^{-1}\), \(e_{3}=0.9~\mathrm{day}^{-1}\), \(\gamma _{1}=0.01~\mathrm{day}^{-1}\), \(\gamma _{2}=0.1~\mathrm{day}^{-1}\), \(\epsilon =0.05\) and \(\delta =0.9\) where \(x(0)=67\), \(y(0)=0.1\) and \(z(0)=0.1\) in which all the conditions in Theorem 5 are satisfied. The solution trajectory projected onto the \((x,y)\)-plane, \((x,z)\)-plane, \((y,z)\)-plane and the time courses of the population densities of BPH, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata populationare are also shown in Fig. 10. We can see that the solution tends to a stable equilibrium as time passes as predicted in Theorem 5.
Moreover, we also found that the system of Eqs. (4)–(6) can exhibit a chaotic behavior when the parametric values are \(c_{1}=0.85~\mathrm{day}^{-1}\), \(c_{2}=0.95~\mathrm{day}^{-1}\), \(c_{3}=0.239~\mathrm{day}^{-1}\), \(c_{4}=0.84~\mathrm{day}^{-1}\), \(c_{5}=0.5~\mathrm{day}^{-1}\), \(k_{1}=67.5\), \(k_{2}=0.02\), \(k_{3}=0.9\), \(h_{1}=0.024\), \(h_{2}=0.063\), \(e_{1}=0.0001~\mathrm{day}^{-1}\), \(e _{2}=0.95~\mathrm{day}^{-1}\), \(e_{3}=0.56~\mathrm{day}^{-1}\), \(\gamma _{1}=0.02~\mathrm{day} ^{-1}\), \(\gamma _{2}=0.085~\mathrm{day}^{-1}\), \(\epsilon =0.004\) and \(\delta =0.8\) where \(x(0)=0.52\), \(y(0)=0.8733\) and \(z(0)=0.0987\). The computer simulation is as shown in Fig. 11. The solution trajectory projected onto the \((x,y)\)-plane, \((x,z)\)-plane, \((y,z)\)-plane and the time courses of the population densities of BPH, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata populationare. In addition, Fig. 12 shows that even though the initial values are very slightly different, the solution trajectories will stay close for a period of time, before starting to follow noticeably different paths as time passes.
Conclusion
The developed model for the population dynamics of BPH, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata is analysed theoretically by the geometric singular perturbation technique. The five possible cases leading to different dynamic behaviors are investigated. The conditions for each case to occur are stated in Theorems 1–5. Computer simulations are carried out showing that the results correspond to the theoretical predictions in the five cases.
In addition, we also investigate that our model can demonstrate a chaotic behavior which has been observed in the paddy field when the outbreak of BPH occurs as reported in [22] and hence, our model might be used to investigate the control of BPH by the release of its natural enemies, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata or the use of insecticide and pathogen in the paddy field. Time delays in the development of BPH, Cyrtorhinus Lividipennis and Lycosa Pseudoannulata reported in the literature will be incorporated in our model in our future work so that our model will be more realistic.
References
- 1.
Kole, C.: Cereals and Millets. Genome mapping and molecular breeding in plants, vol. 1. Springer, Berlin (2006)
- 2.
Maclean, J., Hardy, B., Hettel, G.: Chap. 6: rice around the world. In: Hardy, B. (ed.) Rice Almanac, 4th edn., pp. 77–103. International Rice Research Institute, Philippines (2013)
- 3.
Oerke, E.C.: Crop losses to pests. J. Agric. Sci. 144, 13–43 (2006)
- 4.
Zhang, X., Liao, X., Mao, K., Zhang, K., Wan, H., Li, J.: Insecticide resistance monitoring and correlation analysis of insecticides in field populations of the brown planthopper Nilaparvara lugens (Stål) in China 2012–2014. Pestic. Biochem. Physiol. 132, 13–20 (2016)
- 5.
Ramli, N.H., Yusup, S., Kueh, B.W.B., Kamarulzaman, P.S.D., Osman, N., Rahim, M.A., Aziz, R., Mokhtar, S., Ahmad, A.B.: Effectiveness of biopesticides against brown planthopppers (Nilaparvata lugens) in paddy cultivation. Sustain. Chem. Pharm. 8, 16–20 (2018)
- 6.
Katti, G., Pasalu, I.C., Padmakumari, A.P., Padmavathi, C., Jhansilakshmi, V., Krishnaiah, N.V., Bentur, J.S., Prasad, J.S., Rao, Y.K.: Biological control of insect pests of rice. Technical Bulletin, Directorate of Rice Research, Rajendranagar, Hyderabad, AP, India 22, 22 (2007)
- 7.
Sigsgaard, L.: Early season natural control of the brown planthopper nilaparvata lugens: the contribution and interaction of two spider species and a predatory bug. Bull. Entomol. Res. 97, 533–544 (2007)
- 8.
Preetha, G., Stanley, J., Suresh, S., Samiyappan, R.: Risk assessment of insecticides used in rice on miridbug, Cyrtorhinus lividipennis reuter, the important predator of brown planthopper, Nilaparvata lugens (stal.). Chemosphere 80, 498–503 (2010)
- 9.
Khan, A.A., Misra, D.S.: Biology of woft spider Lycosa pseudoannulata boesenberg and strand (Araneae: Lycosidae). Environ. Ecol. 26, 796–799 (2008)
- 10.
Samiayyan, K.: Chap. 15: spider—the generalist super predators in agro-ecosystem. In: Abrol, D.P. (ed.) Integrated Pest Management, 1st edn., pp. 283–310. Academic Press, Coimbatore (2014)
- 11.
Samal, P., Misra, B.C.: Spider: the most effective natural enemies of the brown planthopper in rice. Rice Entomol. Newsl. 3, 31 (1975)
- 12.
Khan, A.A., Misra, D.S.: Study on the life table of wolf spider (Lycosa pseudoannulata boes. et str., Araneae: Lycosidae) in rice ecosystem. Oryza 40, 31–33 (2003)
- 13.
Heong, K., Bleih, S., Lazaro, A.: Predation of Cyrtorhinus lividipennis reuter on eggs of the green leafhopper and brown planthopper in rice. Res. Popul. Ecol. 32, 255–262 (1990)
- 14.
Heong, K., Rubia, E.: Functional response of Lycosa pseudoannulata on brown planthoppers (BPH) and green leafthoppers (GLH). Int. Rice Res. Not. 14(6), 29–30 (1989)
- 15.
Win, S., Muhamad, R., Ahmad, A., Adam, A.: Life table and population parameters of Nilaparvata lugens Stal. (Homoptera: Delphacidae) on rice. Trop. Life Sci. Res. 22(1), 25–35 (2011)
- 16.
Reyes, T., Gabriel, B.: The life history and consumption habits of Cyrtorhinus lividipennis reuter (Hemiptera: Miridae). Philipp. Entomol. 3(2), 79–88 (1975)
- 17.
Gavarra, M., Raros, R.: Studies on the biology of the predatory wolf spider, Lycosa pseudoannulata (Araneae: Lycosidae). Philipp. Entomol. 29(6), 427–444 (1973)
- 18.
Jones, C.: Geometric singular perturbation theory. In: Dynamical Systems, Montecatibi Terme. Lecture Notes in Math. vol. 1609, pp. 44–118 (1994)
- 19.
Kaper, T.J.: An introduction to geometric methods and dynamical systems theory for singular perturbation problems. Analyzing multiscale phenomena using singular perturbation methods. In: Cronin, J., O’Malley, R.E. Jr. (eds.) Proc. Symposia Appl. Math., vol. 56. Am. Math. Soc., Providence (1999)
- 20.
Muratori, S., Rinaldi, S.: Low-and high-frequency oscillations in three-dimensional food chain systems. SIAM J. Appl. Math. 52, 1688–1706 (1992)
- 21.
Jiang, X., Huang, Q., Ling, Y., Chen, Y., Xiao, G., Huang, S., Wu, B., Huang, F., Cai, J., Long, L.: Functional responses of Cyrtorhinus lividipennis to eggs of Nilaparvara lugens are not affected by genetically modified herbicide-tolerant rice. J. Integr. Agric. (2015). https://doi.org/10.1016/S2095-3119(14)60953-9
- 22.
Hu, G., Lu, F., Zhai, B.P., Lu, M.H., Liu, W.C., Zhu, F., Wu, X.W., Chen, G., Zhang, X.X.: Outbreaks of the brown planthopper Nilaparvara lugens (Stål) in the Yangtze River delta: immigration or local reproduction? PLoS ONE 9(2), e88973 (2014). https://doi.org/10.1371/journal.pone.0088973
Funding
We acknowledge the support by the Centre of Excellence in Mathematics, the Commission on Higher Education, Thailand and the Science Achievement Scholarship of Thailand (SAST).
Author information
Affiliations
Contributions
All authors contributed equally to this work. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Manorod, S., Rattanakul, C. Modelling the population dynamics of brown planthopper, Cyrtorhinus lividipennis and Lycosa pseudoannulata. Adv Differ Equ 2019, 265 (2019). https://doi.org/10.1186/s13662-019-2217-y
Received:
Accepted:
Published:
Keywords
- Mathematical model
- Brown planthoppers
- Cyrtorhinus lividipennis
- Lycosa pseudoannulata