Dynamics of a plant–herbivore model with differential–difference equations

Abstract

This paper studies the behavior of a plant–herbivore model including both differential and difference equations. To analyze global behavior of the model, we consider the solution of the system in a certain subinterval which gives to system of difference equations. The boundedness characters, the periodic nature, both local and global stability conditions of the plant–herbivore system are investigated. Numerical studies indicate that the system exhibits Neimark–Sacker bifurcation for different parameter values in a certain regions.

Keywords:

• plant–herbivore system
• difference equation
• stability
• Neimark–Sacker bifurcation

Mathematics subject classifications:

• 37N25
• 39A28
• 39A30

Public Interest Statement

In this study, we consider a plant–herbivore model which consists of ordinary differential equations. Our aim is to build a better understanding of how both discrete and continuous times affect the dynamic behavior of plant–herbivore interactions. Therefore, we add discrete time to this model and obtain a system of differential equations with piecewise constant arguments which gives system of difference equations. The boundedness characters, the periodic nature, both local and global stability conditions of the system are investigated.

1. Introduction

Classical approaches to modeling plant–herbivore interactions are based on predator–prey system (Caughley & Lawton, 1981; May, 2001). This interaction has been described in much research using discrete and continuous model (Agiza, ELabbasy, EL-Metwally, & Elsadany, 2009; Chattopadhayay, Sarkar, Frıtzsche-Hoballah, Turlıngs, & Bersıer, 2001; Danca, Codreanu, & Bakó, 1997; Das & Sarkar, 2001; Edelstein-Keshet, 1986; Feng, Qiu, Liu, & DeAngelis, 2011; Lebon, Mailleret, Dumont, & Grognard, 2014; Li, 2011; Liu, Feng, Zhu, & DeAngelis, 2008; Mukherjee, Das, & Kesh, 2011; Ortega-Cejas, Fort, & Méndez, 2004; Owen-Smith, 2002; Saha & Bandyopadhyay, 2005; Sui, Fan, Loladze, & Kuang, 2007; Sun, Chakraborty, Liu, Jin, & Anderson, 2014; Zhao, Feng, Zheng, & Cen, 2015). The model of Li (2011) is a system of differential equations with Holling type II functional response where the plant toxin’s influence in herbivores is considered. In study, Mukherjee et al. (2011) have used discrete time model with Holling type II functional response for describing the plant–herbivore interaction.

It is well known that discrete time models governed by difference equations are more appropriate than the continuous time models when the populations have non-overlapping generations. So, a significant number of the study on the mathematical models of plant–herbivore interactions are described by the system of difference equations (Agiza et al., 2009; Danca et al., 1997; Mukherjee et al., 2011; Sui et al., 2007). In addition, working with difference equations instead of differential equations allows us to some advantages. Discrete dynamical models can bring about easier computational methods for the persistence, periodic solutions, boundedness, local and global properties of the dynamical system.

In plant–herbivore interactions, delay differential equations may widely occur due to herbivore damage and deployment of inducible defenses (Das & Sarkar, 2001; Ortega-Cejas et al., 2004; Sun et al., 2014). From this point of view, Sun et al. (2014) and et all have constructed a reaction-diffusion model with delay governed by system of partial differential equations where the effect of time delay on the herbivore cycles is investigated. In addition, the properties of delay differential equations are very close to differential equation with piecewise constant arguments. In Cooke and Györi study (1994), it was pointed out that these equations can be used to get approximate solutions to delay differential equations that include discrete delays. In such biological situations, dynamics of growth and death of populations can be described by differential equations otherwise, difference equations may reflect the interaction of two populations such as competition or predation phenomena (Gurcan, Kartal, Ozturk, & Bozkurt, 2014; Kartal & Gurcan, 2015). In the literature, various types of biological model consisting of differential equations with piecewise constant arguments have been analyzed using the method of reduction to discrete equations (Busenberg & Cooke, 1982; Gopalsamy & Liu, 1998; Gurcan et al., 2014; Kartal & Gurcan, 2015; Liu & Gopalsamy, 1999; Öztürk, Bozkurt, & Gurcan, 2012).

In the present paper, our aim is to build a better understanding of how both discrete and continuous times affect the dynamic behavior of plant–herbivore interactions. So we will reconsider the model (see Chattopadhayay et al., 2001)(1.1) $$\left\{\begin{array}{c}\frac{dx}{dt}=rx(t)\left(1-\frac{x(t)}{K}\right)-\mathit{\alpha }x(t)y(t),\hfill \\ \frac{dy}{dt}=-sy(t)+\mathit{\beta }x(t)y(t),\hfill \end{array}\right.$$(1.1)

as a system of differential equations with piecewise constant arguments such as(1.2) $$\left\{\begin{array}{c}\frac{dx}{dt}=rx(t)\left(1-\frac{x(t)}{K}\right)-\mathit{\alpha }x(t)y(\left[\left[t\right]\right]),\hfill \\ \frac{dy}{dt}=-sy(t)+\mathit{\beta }x\left(\left[\left[t\right]\right]\right)y(t),\hfill \end{array}\right.$$(1.2)

which include both differential and difference equations. In this model, x(t) and y(t) represent the density of plant and herbivore population, respectively, $\left[\left[t\right]\right]$ denotes the integer part of $t\in [0,\infty )$ and all these parameters are positive. The parameter r, K, and $\mathit{\alpha }$ is the intrinsic growth rate, environmental carrying capacity, and specific predation rate of plant species, respectively. s represents the death rate of herbivores and $\mathit{\beta }$ is the conversion factor of herbivores (Chattopadhayay et al., 2001). The logistic term $rx(t)\left(1-\frac{x(t)}{K}\right)$ and the term sy(t) include only a continuous time for the growth of plant and for the death of herbivore, respectively. The predational form $\mathit{\alpha }x(t)y\left(\left[\left[t\right]\right]\right)$ represent the loss of plant population and $\mathit{\beta }x\left(\left[\left[t\right]\right]\right)y(t)$ is conversion factor of herbivores which include both discrete and continuous time for a each populations. So the plant–herbivore interaction is considered in a certain subinterval and is modeled using a system of differential equations with piecewise constant arguments.

2. Local and global stability analysis

System (1.2) can be written an interval $t\in [n,n+1)$ as follows:(2.1) $$\left\{\begin{array}{c}\frac{dx}{dt}-x(t)\left(r-\mathit{\alpha }y(n)\right)=-rk{x}^2(t),\hfill \\ \frac{dy}{y(t)}=(\left(\mathit{\beta }x(n)-s\right))dt,\hfill \end{array}\right.$$(2.1)

where $\frac{1}{K}=k.$

By solving each equations of the system (2.1) and letting $t\to n+1$, we obtain a system of difference equations(2.2) $$\left\{\begin{array}{c}x\left(n+1\right)=\frac{x(n)(r-\mathit{\alpha }y(n))}{(r-\mathit{\alpha }y(n)-rkx(n))e^{-(r-\mathit{\alpha }y(n))}+rkx(n)},\hfill \\ y\left(n+1\right)=y(n)e^{\mathit{\beta }x(n)-s},\hfill \end{array}\right.$$(2.2)

System (2.2) reflects the dynamical behavior of the system of differential equations with piecewise constant arguments. So we will consider the system of difference equation to analyze the global behavior of system (1.2).

The equilibrium points of system (2.2) can be obtained as$$E_0=(0,0),E_1=\left(\frac{1}{k},0\right),E_{\ast }=\left(\frac{s}{\mathit{\beta }},\frac{r}{\mathit{\alpha }}\left(1-\frac{ks}{\mathit{\beta }}\right)\right).$$

We note that the positive equilibrium of the system exists if $\mathit{\beta }>ks$. Now, we will find Jacobian matrix of the system to investigate the dynamic behavior of the model.

Theorem 2.1.

The equilibrium points $E_0$ and $E_1$ are saddle point.

Proof

At the equilibrium point $E_0$, the Jacobian matrix is the form$$J_0=\left(\begin{array}{cc}{e}^r& 0\\ 0& e^{-s}\end{array}\right).$$

The matrix $J_0$ has eigenvalues ${\mathit{\lambda }}_1=e^r,{\mathit{\lambda }}_2=e^{-s}$. Hence ${\mathit{\lambda }}_1>1$ and ${\mathit{\lambda }}_2<1$ and consequently $E_0$ is saddle point. On the other hand, the Jacobian matrix $J_1$ at the point $E_1$ is$$J_1=\left(\begin{array}{cc}{e}^{-r}\hfill & -\frac{\mathit{\alpha }-e^{-r}\mathit{\alpha }}{k}r\hfill \\ 0\hfill & e^{-s+\frac{\mathit{\beta }}{k}}\hfill \end{array}\right)$$

which gives eigenvalues ${\mathit{\lambda }}_1=e^{-r}\text{and}{\mathit{\lambda }}_2=e^{-s+\frac{\mathit{\beta }}{k}}$. Considering the condition $\mathit{\beta }>ks$, we can say that $E_1$ is saddle point.

On the other hand, the Jacobian matrix $J_{\ast }$ at the positive equilibrium point $E_{\ast }$ is$$J_{\ast }=\left(\begin{array}{cc}{e}^{-kr\overline{x}}& \frac{(-1+e^{-kr\overline{x}})\mathit{\alpha }}{kr}\\ \mathit{\beta }\overline{y}& 1\end{array}\right)$$

which yields the following characteristic equation$$p(\mathit{\lambda })={\mathit{\lambda }}^2+\mathit{\lambda }\left(-1-e^{-kr\overline{x}}\right)+e^{-kr\overline{x}}+\frac{(1-e^{-kr\overline{x}})\mathit{\alpha }\mathit{\beta }\overline{y}}{kr}=0.$$

Now, we can apply Schur–Cohn criterion to determine stability conditions of the system with characteristic equation $p(\mathit{\lambda })$.

Theorem 2.2

The positive equilibrium point $E_{\ast }$ of system (2.2) is local asymptotically stable if and only if$$ks<\mathit{\beta }<k+ks.$$

Proof

From the Schur–Cohn criterion, $E_{\ast }$ is local asymptotically stable if and only if

(a)	$$p(1)=1+e^{-kr\overline{x}}+\frac{(1-e^{-kr\overline{x}})\overline{y}}{k}r-1-e^{-kr\overline{x}}>0,$$
(b)	$$p(-1)=1+e^{-kr\overline{x}}+\frac{(1-e^{-kr\overline{x}})\overline{y}}{kr}+1+e^{-kr\overline{x}}>0,$$
(c)	$$D_1^{+}=1+e^{-kr\overline{x}}+\frac{\left(1-e^{-kr\overline{x}}\right)\overline{y}}{kr}>0,$$
(d)	$$D_1^{-}=1-e^{-kr\overline{x}}-\frac{\left(1-e^{-kr\overline{x}}\right)\overline{y}}{kr}>0.$$

The condition (a), (b), and (c) gives the inequalities(2.3) $$p(1)=\frac{\left(1-e^{-kr\overline{x}}\right)\mathit{\alpha }\mathit{\beta }\overline{y}}{kr}>0,$$(2.3) (2.4) $$p(-1)=2+2e^{-kr\overline{x}}+\frac{\left(1-e^{-kr\overline{x}}\right)\mathit{\alpha }\mathit{\beta }\overline{y}}{kr}>0$$(2.4)

and(2.5) $$D_1^{+}=1+e^{-kr\overline{x}}+\frac{\left(1-e^{-kr\overline{x}}\right)\mathit{\alpha }\mathit{\beta }\overline{y}}{kr}>0$$(2.5)

which always hold under the condition $\mathit{\beta }>ks$. From (d), we get$$e^{-kr\overline{x}}+\frac{\left(1-e^{-kr\overline{x}}\right)\mathit{\alpha }\mathit{\beta }\overline{y}}{kr}<1$$

which reveal$$\mathit{\beta }<k+ks.$$

This completes the proof. □

For the parameter values $r=0.2,\mathit{\alpha }=0.6,K=5,\mathit{\beta }=0.01,s=0.02$ and using initial conditions $x(1)=0.2,y(1)=0.22$, it can be seen that the positive equilibrium point $\left(\overline{x},\overline{y}\right)=(2,0.2)$ is local asymptotically stable, where blue and red graphs represent population density of plant and herbivore population, respectively.

Theorem 2.3

Let ${\{x(n),y(n)\}}_{n=-1}^{\infty }$ be a positive solution of system (2.2); then$$x(n)\le \frac{e^r}{k\left(e^r-1\right)}.$$

In addition, if $y(n)<x(n),\text{then}y(n)\le \frac{e^r}{k\left(e^r-1\right)}{e}^{\frac{\mathit{\beta }{e}^r}{k\left(e^r-1\right)}-s}.$

Proof

It can be easily seen that$$x\left(n+1\right)=\frac{x(n)[r-\mathit{\alpha }y(n)]e^{r-\mathit{\alpha }y(n)}}{r-\mathit{\alpha }y(n)+rkx(n)(e^{r-\mathit{\alpha }y(n)}-1)}\le \frac{[r-\mathit{\alpha }y(n)]e^{r-\mathit{\alpha }y(n)}}{\text{rk}(e^{r-\mathit{\alpha }y(n)}-1)}\le \frac{e^r}{k\left(e^r-1\right)}.$$

Also, it can be shown that $y\left(n+1\right)\le \frac{e^r}{k\left(e^r-1\right)}{e}^{\frac{\mathit{\beta }{e}^r}{k\left(e^r-1\right)}-s}$ under the condition $y(n)<x(n)$.

Theorem 2.4

The system has no prime period-two solutions.

Proof

On the contrary, suppose that the system (2.2) has a distinctive prime period-two solutions$$\dots ,\left(w_1,q_1\right),\left(w_2,q_2\right),\left(w_1,q_1\right),\left(w_2,q_2\right)\dots .$$

where $w_1\ne w_2$, and $q_1\ne q_2$, and $w_{\text{i}},q_{\text{i}}$ are positive real numbers for $i\in \{1,2\}$. Then, from system (2.2) one has$$\left\{\begin{array}{c}{w}_1=\frac{w_2[r-\mathit{\alpha }{q}_2]}{[r-\mathit{\alpha }{q}_2-rk{w}_2]e^{-\left[r-\mathit{\alpha }{q}_2\right]}+rk{w}_2},\hfill \\ q_1=q_2{e}^{\mathit{\beta }{w}_2-s}\hfill \\ w_2=\frac{w_1[r-\mathit{\alpha }{q}_1]}{[r-\mathit{\alpha }{q}_1-rk{w}_1]e^{-\left[r-\mathit{\alpha }{q}_1\right]}+rk{w}_1}\hfill \\ q_2=q_1{e}^{\mathit{\beta }{w}_1-s},\hfill \end{array},\right.$$

Since $q_1\ne q_2$, we have $\mathit{\beta }{w}_2-s\ne 0$ and $\mathit{\beta }{w}_1-s\ne 0$. From the second and last equation in the system, we have$$q_1^2=q_2^2{e}^{\mathit{\beta }{w}_2-s-\mathit{\beta }{w}_1+s}.$$

If $q_2$ is written the above equation, we hold$$q_1^2=q_1^2{e}^{\mathit{\beta }{w}_1-s+\mathit{\beta }{w}_2-s}.$$

This equation must satisfy$$\mathit{\beta }{w}_1-s+\mathit{\beta }{w}_2-s=0$$

which is a contradiction $\mathit{\beta }{w}_2-s\ne 0$ and $\mathit{\beta }{w}_1-s\ne 0$.

Theorem 2.5

Let $A_1=r-\mathit{\alpha }y(n)$ and $A_2=\mathit{\beta }x(n)-s$. Suppose that the conditions of Theorem 2.1 hold and

(i)	$$\text{Let}y(n)<\frac{r}{\mathit{\alpha }}\text{and}\overline{x}<\frac{A_1(1+e^{-A_1})}{2{\text{rke}}^{-A_1}}\text{for}x(n)\in \left(0,\frac{2\overline{x}{e}^{-A_1}}{1+e^{-A_1}}\right),$$
(ii)	$$\text{Let}y(n)>\frac{r}{\mathit{\alpha }}\text{and}\overline{x}>\frac{-A_1(1+e^{-A_1})}{2\text{rk}\left(e^{-A_1}-1\right)}\text{for}x(n)\in \left(\frac{2\overline{x}{e}^{-A_1}}{1+e^{-A_1}},2\overline{x}\right),$$
(iii)	$$\text{Let}y(n)>\frac{r}{\mathit{\alpha }}\text{for}x(n)\in \left(2\overline{x},\infty \right),$$
(iv)	$$\text{Let}{A}_2>0\text{fory}(n)\in \left(0,\frac{2\overline{y}}{1+e^{A_2}}\right),$$
(v)	$$\text{Let}{A}_2<0\text{fory}(n)\in \left(\frac{2\overline{y}}{1+e^{A_2}},\infty \right).$$

Then the positive equilibrium point of system (2.2) is global asymptotically stable.

Proof

We define a Lyapunov function as$$V(n)={[q(n)-\overline{q}]}^2,n=0,1,2\dots$$

where $\overline{q}=\left(\overline{x},\overline{y}\right)$ is positive equilibrium point of system (2.2).

The change along the solutions of the system is$$\mathrm{\Delta }V(n)=V\left(n+1\right)-V(n)=\left\{q\left(n+1\right)-q(n)\right\}\left\{q\left(n+1\right)+q(n)-2\overline{q}\right\}.$$

From the first equation in (2.2), we hold;$$\mathrm{\Delta }{V}_1(n)=\left[x\left(n+1\right)-x(n)\right][x\left(n+1\right)+x(n)-2\overline{x}]=x(n)\left[(A_1-rkx(n))(1-e^{-A_1}\right)][A_1\left(x(n)+x(n)e^{-A_1}-2\overline{x}{e}^{-A_1}\right)+rkx(n)(x(n)-2\overline{x})(1-e^{-A_1})].$$

By considering (i), (ii), and (iii), we have $\mathrm{\Delta }{V}_1(n)<0.$ These imply that $\underset{n\to \infty }{lim}x(n)=\overline{x}$. Additionally, we can show that $\mathrm{\Delta }{V}_2(n)<0$ which gives $\underset{n\to \infty }{lim}y(n)=\overline{y}$.

3. Bifurcation analysis

In this section, we investigate existence of stationary bifurcation (fold, transcritical, and pitchfork bifurcation), period doubling bifurcation, and Neimark–Sacker bifurcation for the system (2.2). All of these bifurcations can be analyzed under the set of algebraic condition that is called Schur–Cohn criterion. It is well known that the system may undergo stationary bifurcation if and only if $p(1)=0$, $p(-1)>0$, $D_1^{+}>0$ and $D_1^{-}>0$. On the other hand, inequalities $p(1)>0,$ $p(-1)=0$, $D_1^{+}>0$ and $D_1^{-}>0$ give the conditions of period doubling bifurcation. But considering (2.3) and (2.4), it is easily seen that these conditions do not hold for the system. Therefore, stationary bifurcation and period doubling bifurcation do not exist for the system.

Now, we can investigate the existence of Neimark–Sacker bifurcation for the plant–herbivore model (Hone, Irle, & Thurura, 2010). The algebraic condition of Neimark–Sacker bifurcation can be obtained from the analysis of inequalities $p(1)>0$, $p(-1)>0$, $D_1^{+}>0$ and $D_1^{-}=0$. In local stability analysis, we have already shown that the inequalities $p(1)>0$, $p(-1)>0$, $D_1^{+}>0$ are always exist. Therefore, we will only analyse the equation $D_1^{-}=0$ to determine Neimark–Sacker bifurcation

Figure 1. Stable equilibrium point of the system for $r=0.2,\mathit{\alpha }=0.6,K=5,\mathit{\beta }=0.01,s=0.02,x(1)=0.2\text{and}y(1)=0.22$.

Figure 2. Stable limit cycle for $r=0.2,\mathit{\alpha }=0.6,\mathit{\beta }=0.01,s=0.02$ and $\overline{K}=102$.

condition.

Theorem 3.1

System (2.2) undergoes Neimark–Sacker bifurcation if and only if$$\overline{k}=\frac{\mathit{\beta }}{1+s}.$$

Proof

This result comes from the analysis of $D_1^{-}=0.$

Using the condition of Theorem 3.1 with the parameters given in Figure 1, we have the Neimark–Sacker bifurcation point as $\overline{K}=102$ (Figure 2).

4. Result and discussion

In this paper, dynamics of a discrete-continuous time plant–herbivore model has been investigated. Local and global stability properties of the positive equilibrium point are analyzed. It is interesting to note that when conversion factor of herbivores becomes low then the system converges to a stable situation. On the other hand, we investigate possible bifurcation types for the system and observe that the system exhibits Neimark–Sacker bifurcation. This type of bifurcation has been observed many plant–herbivore models (Liu et al., 2008; Saha & Bandyopadhyay, 2005; Zhao et al., 2015) and shows that periodic or quasi-periodic solutions occurs as a result of limit cycle.

In our manuscript, the parameter K (environmental carrying capacity of plant species) is determined as a bifurcation parameter. When the environmental carrying capacity of plant species reaches to $\overline{K}=102$, the system enters a Neimark–Sacker bifurcation as a result of stable limit cycle (Figure 2). If K exceeds the $\overline{K}$, the system continues oscillatory behavior with growing amplitude (Figure 3). So we can say that the parameter K has a strong effect on the stability of the system so as to control two populations.

Figure 3. Graph of the iteration solution of system (2.2) for $r=0.2,\mathit{\alpha }=0.6,\mathit{\beta }=0.01,s=0.02$ and K = 200.

Additional information

Funding

Funding. The author received no direct funding for this research.

Notes on contributors

S. Kartal

Senol Kartal is an associate professor in Nevsehir Haci Bektas Veli University, Turkey. His PhD thesis is about the population dynamics (Modeling of Tumor Immune System Dynamics Using System of Difference Equations and Its Stability Analysis). His research interests are related to population dynamics, discrete and continuous dynamical system, and bifurcation theory. He has published his research contributions in some internationally renowned journals whose publishers are Elsevier, Taylor & Francis, Wiley, and other journals.