Diffusion Dynamics and Impact of Noise on a Discrete-Time Ratio-Dependent Model: An Analytical and Numerical Approach

Abstract

The paper deals with the dynamical behavior of a discrete-time ratio-dependent predator–prey system. The predator dependence is one of the main concerns of the system. The stability analysis of this 2-dimensional map was carried out analytically. Numerical simulation was carried out to verify the analytical results. We analyzed some specific features that could arise in discrete system. Basin of attraction was found for the endemic equilibrium state. We extended the numerical simulation for the maximal Lyapunov exponent. The presence of positive Lyapunov exponent indicated chaotic behavior of the map. The sensitive dependence on initial condition is one of the criteria for a discrete system. We showed that the system is sensitive on the initial conditions. We also carried out the analysis of diffusion and impact of noise.

1. Introduction

The pioneering work of [1,2] on population biology is the threshold created in today’s explosion of the field of population biology. The area has become vast, enriching the field of epidemiology and medicine. The framework of the Lotka–Volterra model is the two differential equations with a simple proportionality between the prey consumption and predation production. The proportional elements have been many functions and the modeling has become more realistic. The proportional functions are the functional response or trophic interaction functions. The Holling type I II III [3] are old to represent a predator–prey system. The prey dependence [4] trophic interactions in the population model have reigned for long time, even today [5]. Such a category of population models cannot produce the situations of biological control. If the trophic function depends on the single variable $N/P$ (N, the prey, P, the predator), then the essential properties of predator dependence are rendered, called ratio-dependence [4]. The work on this ratio-dependence model has mainly focused on continuous models. Such a model only gives some bifurcation and limit cycle as the population evolution. But the discrete time predator–prey models have richer dynamics than the continuous counterpart. The discrete model can explain the basin of attraction, fractal dimension and the chaotic behavior of an attractor and the sensitive dependence of initial conditions. Danca demonstrated that the discrete-time predator–prey model with Holling type I functional response exhibits a chaotic behavior [6]. More realistic results are obtained in the works of [7,8]. The discrete time ratio-dependent predator–prey explains the chaotic dynamics of the model system. The study sheds more light on the complex dynamical behavior. In this chapter, we investigated the stability behavior and different features of a discrete type ratio-dependent predator–prey system.

In this paper, we investigated the stability behavior and different features of a discrete type ratio-dependent predator–prey system. In Section 3, we studied the local stability analysis by using Jury’s conditions. In Section 4, we performed and presented diffusion dynamics. Section 5 deals with noise impact on the stochastic system. In Section 6, we presented computer simulations with notable significances as maximal Lyapunov exponents responsible for chaos of the system and existence of orbit due to sensitivity. In Section 7, the entire compilation of analysis with results is presented as concluding remarks.

2. The Discrete-Time Model Equations

We begin with the continuous-time predator–prey system. In most predator–prey systems, the predator is to search for food and also compete for food, and the more realistic functional response is a function which is dependent on both prey and predator densities. The Michaelis-Menten type functional response introduced by [4] is given by

$$g(N/P)=\frac{aN}{P+ahN},$$

(1)

where N, P are, respectively, the prey and predator densities, a is the maximum prey consumption rate, and h is the predator handling time. This ratio-dependent functional response is strongly supported by numerous fields, laboratory experiments, and ecological literatures [9,10,11,12]. With this functional response, the general type ratio-dependent predator–prey system due to [9] is

$$\begin{array}{c}\hfill \dot{N}=rNf\left(N\right)-g(N/P)P,\\ \hfill \dot{P}=-dP+cg(N/P)P.\end{array}$$

(2)

Here, $rNf\left(N\right)$ is the prey growth rate in the absence of predator, d is the death rate of the predator in absence of prey. Here, we consider the logistic growth rate of the prey given by

$$rNf\left(N\right)=rN\left(1-\frac{N}{K}\right).$$

(3)

So, the general ratio-dependent predator–prey model becomes

$$\begin{array}{c}\hfill \dot{N}=rN\left(1-\frac{N}{K}\right)-\frac{aN}{P+ahN}P,\\ \hfill \dot{P}=-dP+c\frac{aN}{P+ahN}P,\end{array}$$

(4)

where r, K are, respectively, the intrinsic growth rate and environmental carrying capacity of the prey in the absence of predator.

The number of independent parameters can be reduced by using the dimensionless quantities for prey N, predator P and time t as follows: $N=Kx$, $P=Kahy$, and $t=t^{{}^{\prime }}/r$. Then, the system becomes

$$\begin{array}{c}\hfill \frac{dx}{dt}=x(1-x)-\alpha \frac{xy}{x+y},\\ \hfill \frac{dy}{dt}=-\beta y+\gamma \frac{xy}{x+y},\end{array}$$

(5)

where $\alpha =a/r$, $\beta =d/r$, and $\gamma =c/rh$. It is available in any standard literature. Thus, $\alpha$ is the consumption ability, $\beta$ is the death rate of the predator, and $\gamma$ is the predator growing ability.

The model system (5) is the standard ratio-dependent predator–prey model studied by many authors [5]. In this work, we consider the discrete-time predator–prey model analogous to its continuous counterpart (5). In the continuous time model (5), we replace the derivatives with the divided differences

$$\frac{dx}{dt}=\frac{x(t+\Delta t)-x\left(t\right)}{\Delta t},\frac{dy}{dt}=\frac{y(t+\Delta t)-y\left(t\right)}{\Delta t},$$

(6)

and then, taking $\Delta t=1$, the above system of equations reduces to

$$\begin{array}{c}\hfill x_{n+1}=x_n+x_n(1-x_n)-\alpha \frac{x_n{y}_n}{x_n+y_n},\\ \hfill y_{n+1}=y_n-\beta y_n+\gamma \frac{x_n{y}_n}{x_n+y_n},\end{array}$$

(7)

where $x_n=x\left(t\right)$, $y_n=y\left(t\right)$. In this study, the system of Equation (7) represents the main equations of interest. Thus, the continuous system which represents a flow of two predator and prey is now a 2-dimensional map of two species from the discrete time n to $n+1$. If this map is denoted by T, then $T(x_n,y_n)=(x_{n+1},y_{n+1})$ governed by the set of Equation (7). This map will determine the population evolution from the initial state $t=0$ to a state $t=n$, so that $(x_n,y_n)=T^n(x_0,y_0)$. This evolution will move through different complex dynamical phenomena.

3. Steady State Analysis

We studied the stability behavior of the two-dimensional map governed by the set in Equation (7). The fixed points of the map are given by $T(x_n,y_n)=(x_n,y_n)$. So, we have the following three fixed points: (i) Extinction of both predator and prey population: $E_0=(0,0)$, (ii) extinction of predator population: $E_A=(1,0)$, and (iii) the coexistence of both the population predator and prey: $E_I(x_I,y_I)$, where $x_I=1-\frac{\alpha }{\gamma }(\gamma -\beta )$ and $y_I=\frac{\gamma -\beta }{\beta }{x}_I$.

The first two fixed points, $E_0$ and $E_A$, are always biologically feasible. The feasibility of the fixed point $E_I$ demands $\alpha <\frac{\gamma }{\gamma -\beta }$ with $\beta <\gamma$.

To study the stability behavior of the system near the fixed points, we linearize the model equations about the fixed points. If $(X_n,Y_n)$ be the small perturbations corresponding to the fixed point $(x_{*},y_{*})$, then $x_n=x_{*}+X_n$ and $y_n=y_{*}+Y_n$. The linearized system of the above system of equations is

$$\left(\begin{array}{c}{X}_{n+1}\\ Y_{n+1}\end{array}\right)=\mathbf{J}\left(E(x_{*},y_{*})\right)\left(\begin{array}{c}{X}_n\\ Y_n\end{array}\right),$$

(8)

where $\mathbf{J}\left(E(x_{*},y_{*})\right)=\left(\begin{array}{cc}2(1-x_{*})-\frac{\alpha y_{*}^2}{{(x_{*}+y_{*})}^2}& -\frac{\alpha x_{*}^2}{{(x_{*}+y_{*})}^2}\\ \frac{\gamma y_{*}^2}{{(x_{*}+y_{*})}^2}& 1-\beta +\frac{\gamma x_{*}^2}{{(x_{*}+y_{*})}^2}\end{array}\right)$.

We seek a solution of the form

$$\left(\begin{array}{c}{X}_n\\ Y_n\end{array}\right)=\left(\begin{array}{c}A\\ B\end{array}\right){\lambda }^n$$

(9)

to Equation (8), where A and B are arbitrary constants. So, the characteristic equation will determine the value of $\lambda$ which, in turn, will help in giving the dynamical behavior of the fixed point. The characteristic equation of the system of Equation (8) is

$$f\left(\lambda \right)\equiv {\lambda }^2-P\lambda +Q=0,$$

(10)

where $P=Trace\mathbf{J}\left(E(x_{*},y_{*})\right)$, and $Q=\left|\mathbf{J}\right(E(x_{*},y_{*})\left)\right|$. Let ${\lambda }_1$ and ${\lambda }_2$ be characteristic roots of Equation (10) and also suppose that $f(-1)>0$. Then, the equilibrium point $E(x_{*},y_{*})$ is

(i)

stable if $|{\lambda }_1|<1$, $|{\lambda }_2|<1$, that is, if $f\left(1\right)>0$ and $Q<1$.

(ii)

unstable if $|{\lambda }_1|>1$, $|{\lambda }_2|>1$, that is, if $f\left(1\right)>0$ and $Q>1$.

(iii)

saddle if $|{\lambda }_1|<1$, $|{\lambda }_2|>1$ (or $|{\lambda }_1|>1$, $|{\lambda }_2|<1$), that is, if $f\left(1\right)<0$. We now analyze the stability behavior at several fixed points.

3.1. Dynamical Behavior of $E_0$

As the functional response is not defined at $(0,0)$, the fixed point $E_0(0,0)$ is most complicated for the ratio-dependent model system. We choose a neighboring point $(ϵ,ϵ)$ in the first quadrant corresponding to the fixed point $E_0(0,0)$. Then, the Jacobian matrix $J\left(E_0\right)$ is given by

$$\mathbf{J}\left(E_0\right)=li{m}_{ϵ\to 0+}\left(\begin{array}{cc}2(1-ϵ)-\frac{\alpha {ϵ}^2}{{(ϵ+ϵ)}^2}& -\frac{\alpha {ϵ}^2}{{(ϵ+ϵ)}^2}\\ \frac{\gamma {ϵ}^2}{{(ϵ+ϵ)}^2}& 1-\beta +\frac{\gamma {ϵ}^2}{{(ϵ+ϵ)}^2}\end{array}\right)=\left(\begin{array}{cc}2-\frac{\alpha }{4}& -\frac{\alpha }{4}\\ \frac{\gamma }{4}& 1-\beta +\frac{\gamma }{4}\end{array}\right).$$

Then, the linearized system is

$$\begin{array}{c}\hfill x_{n+1}=\left(2-\frac{\alpha }{4}\right)x_n-\frac{\alpha }{4}{y}_n\\ \hfill y_{n+1}=\frac{\gamma }{4}{x}_n+\left(1-\beta +\frac{\gamma }{4}\right)y_n.\end{array}$$

(11)

The characteristic equation corresponding to the fixed point $E_0(0,0)$ is

$${\lambda }^2-\left(3-\beta -\frac{\alpha -\gamma }{4}\right)\lambda +\left(2-\frac{\alpha }{4}\right)\left(1-\beta +\frac{\gamma }{4}\right)+\frac{\alpha \gamma }{16}=0.$$

(12)

Using the above criteria for dynamical behavior we list the conditions as follows: the fixed point $E_0(0,0)$ is

(i)

stable if (a) $\beta <1$ and $\alpha >max\left\{4-\frac{\gamma }{\beta },12-\frac{3\gamma +8}{\beta },\frac{4(1-2\beta )+\gamma }{1-\beta }\right\}$, or, (b) $\beta >1$ and $max\left\{4-\frac{\gamma }{\beta },12-\frac{3\gamma +8}{\beta }\right\}<\alpha <\left\{\frac{4(2\beta -1)-\gamma }{\beta -1}\right\}$,

(ii)

unstable if (a) $\beta <1$ and $max\left\{4-\frac{\gamma }{\beta },12-\frac{3\gamma +8}{\beta }\right\}<\alpha <\left\{\frac{4(1-2\beta )+\gamma }{1-\beta }\right\}$, or, (b) $\beta >1$ and $\alpha >max\left\{4-\frac{\gamma }{\beta },12-\frac{3\gamma +8}{\beta },\frac{4(2\beta -1)-\gamma }{\beta -1}\right\}$,

(iii)

saddle point if $12-\frac{3\gamma +8}{\beta }<\alpha <4-\frac{\gamma }{\beta }$.

3.2. Dynamical Behavior of $E_A$

The Jacobian matrix of system (8) at $E_A(1,0)$ is

$$\mathbf{J}\left(E_A\right)=\left(\begin{array}{cc}0& -\alpha \\ 0& 1-\beta +\gamma \end{array}\right).$$

The corresponding linearized system is given by

$$\begin{array}{c}\hfill x_{n+1}=-\alpha y_n,y_{n+1}=\left(1-\beta +\gamma \right)y_n.\end{array}$$

(13)

So, $E_A$ is stable if $\beta \in (\gamma ,2+\gamma )$ or a saddle point if $\beta$ does not belong to $(\gamma ,2+\gamma )$.

3.3. Dynamical Behavior of $E_I$

The Jacobian matrix of System (8) corresponding to $E_I$ is

$$\mathbf{J}\left(E_I\right)=\left(\begin{array}{cc}\alpha \left(\frac{\gamma -\beta }{\gamma }\right)\left(2-\frac{\gamma -\beta }{\gamma }\right)& -\frac{\alpha {\beta }^2}{{\gamma }^2}\\ \frac{{(\gamma -\beta )}^2}{\gamma }& 1-\beta +\frac{{\beta }^2}{\gamma }\end{array}\right).$$

The linearized system corresponding to the interior equilibrium point is

$$\begin{array}{c}\hfill x_{n+1}=\alpha \left(\frac{\gamma -\beta }{\gamma }\right)\left(2-\frac{\gamma -\beta }{\gamma }\right)x_n-\frac{\alpha {\beta }^2}{{\gamma }^2}{y}_n\\ \hfill y_{n+1}=\frac{{(\gamma -\beta )}^2}{\gamma }{x}_n+\left(1-\beta +\frac{{\beta }^2}{\gamma }\right)y_n.\end{array}$$

(14)

The characteristic equation of the above system is

$$\begin{array}{c}{\lambda }^2-\left[\frac{\alpha (\gamma -\beta )}{\gamma }\left(2-\frac{\gamma -\beta }{\gamma }\right)+1-\beta +\frac{{\beta }^2}{\gamma }\right]\lambda \\ +\frac{\alpha (\gamma -\beta )}{\gamma }\left[\left(2-\frac{\gamma -\beta }{\gamma }\right)\left(1-\beta +\frac{{\beta }^2}{\gamma }\right)+\frac{\beta (\gamma -\beta )}{{\gamma }^2}\right]=0.\end{array}$$

(15)

From the characteristic Equation (15), we can conclude that the fixed point $E_I$ is

(i)

stable if $\alpha$, $\beta$, and $\gamma$ lie in the parametric space

$$V^S=\left\{(\alpha ,\beta ,\gamma ):\left(1+\frac{\alpha (\gamma -\beta )}{\gamma }(2-\frac{\gamma -\beta }{\gamma })\right)\left(2-\beta +\frac{{\beta }^2}{\gamma }\right)+\frac{\alpha {(\gamma -\beta )}^2{\beta }^2}{{\gamma }^3}>0,\alpha <\frac{{\gamma }^3}{D}\right\}.$$

If mathematical observations drawn here interpret the biological conclusion as if the consumption ability is less than the predator growing ability, then the model system attains its stability, and its dynamics explored more clearly to pick some interesting remarks:

(ii)

unstable if $\alpha$, $\beta$ and $\gamma$ lie in the parametric space

$$V^U=\left\{(\alpha ,\beta ,\gamma ):\left(1+\frac{\alpha (\gamma -\beta )}{\gamma }(2-\frac{\gamma -\beta }{\gamma })\right)\left(2-\beta +\frac{{\beta }^2}{\gamma }\right)+\frac{\alpha {(\gamma -\beta )}^2{\beta }^2}{{\gamma }^3}>0,\alpha >\frac{{\gamma }^3}{D}\right\},$$

where $D=(\gamma -\beta )[(\beta +\gamma )(1-\beta +\frac{{\beta }^2}{\gamma })+\beta (\gamma -\beta )]$.

If mathematical observations drawn here interpret the biological conclusion as if the consumption ability is greater than the predator growing ability, then the model system reaches unstable and, further, there is a chance to attain its stability by additional driven forces or parameters with an appropriate controller/adapter to exhibit its dynamics as effectively and interestingly, which is an advanced wide study.

(iii)

saddle if $\alpha >\frac{\gamma }{\gamma -\beta }$ which is a contradiction for the existence of $E_I$.

If mathematical observations drawn here interpret the biological conclusion as if the difference between the predator growing ability and its death rate is greater than the ratio of predator growing ability to consumption ability, then the model system reaches saddle point, which is an uncertainty where it is neither stable nor unstable and no further dynamics are exhibited by the model system.

These are the conditions that the system exhibits different dynamical behavior.

4. Diffusive Structure and Its Dynamic Forces

In this segment, we deliberated the exceptional influences of transmission of the ideal structure of prey predator systems in environmental science, modeled by diffusion equations. Although the dispersal system is a relatively simple model for the raid of prey species by predators in a spatial domain, the solutions exhibit an extensive spectrum of ecologically pertinent behavior. Spatiotemporal dynamics includes chaos and target patterns [13,14,15]. The study of such spatiotemporal dynamics is an intensive area of research, and there are still many unanswered questions concerning these solution types [15,16,17]. Constructing a simple structure, which consists of prey-predator system with constant harvesting rates and temporal affects, in this segment, we next examine the steadiness of the proposed model under the external dynamical forces as our aim in the next segment. The actual dynamics of the species spread is a result of the interplay between diffusion and deterministic factors. To study the effect of diffusion of ecological population on the model system, let us consider the diffusive equation system as

$$\frac{dx}{dt}=x(1-x)-\alpha \frac{xy}{x+y}+D_1{x}_{uu},$$

(16)

$$\frac{dy}{dt}=-\beta y+\gamma \frac{xy}{x+y}+D_2{y}_{uu},$$

(17)

where $x=x\left(u,t\right),y=y\left(u,t\right)$, u is a space variable, and $x(u,0)>0$; $y(u,0)>0;$ for $u\in \left[0,\mathrm{R}\right].$ The trivial fluctuation edge conditions are specified by ${\left[x_u\right]}_{u=0,R}=0;{\left[y_u\right]}_{u=0,R}=0.$

Now, let us consider the ideal (16) and (17) underneath trivial fluctuations edge ailments. To analyze the role of transmission on this ideal, we deliberate the linear ideal of the structure (16) and (17) about the interior steady state $E_3(x^{*},y^{*})$ as given by

$$\frac{dX}{dt}=-x^{*}X+D_1(-p^{2}X),$$

(18)

$$\frac{dY}{dt}=D_2(-p^{2}Y),$$

(19)

by putting $x=x^{*}+X$; $y=y^{*}+Y$ and assuming the solutions of Equations (16) and (17) in the form

$$X={\alpha }_1{e}^{\lambda t}cospu,Y={\alpha }_2{e}^{\lambda t}cospu,$$

(20)

where p is the wave numeral of perturbation, $\lambda$ is the frequency numeral, and ${\alpha }_i,i=1,2,3$ are the amplitudes. The characteristic equation of (16) and (17) using (20) is

$${\mu }^2+A_1\mu +B_1=0,$$

(21)

where $A_1=x^{*}+p^2(D_1+D_2)$; $B_1=p^4{D}_1{D}_2+D_2{x}^{*}{p}^2$.

Now, our main aim is to find the ailments for diffusive unsteadiness of model system (16) and (17), for this, rewrite $B_1$ as $G\left(p^2\right)$, where $G\left(p^2\right)=D_1{D}_2{\left(p^2\right)}^2+D_2{x}^{*}{p}^2$.

The system (16) and (17) is unstable if one of the above roots of the Equation (21) is optimistic. An essential ailment for a solution to be optimistic is that $D_1{D}_2{\left(p^2\right)}^2+D_2{x}^{*}{p}^2>0$, which implies that

$$p^2>-(x^{*}/D_1).$$

(22)

Since the wave number p is real number, then the above statement is achievable, if $x^{*}>0$, which is always achieved, because by nature $x^{*}$ is real number. The sufficient condition for optimitive of one of the solutions of the Equation (22) is $G\left(p^2\right)<0$. Since $G\left(p^2\right)$ is an expression in $p^2$ where p is the wave number, non-zero positive quantity, the sufficient condition reduces to

$$\Delta <(-p^2),$$

(23)

where $\Delta =D_1/D_2$. Thus, the diffusion of the prey-predator populations drive the ecological system into unstable oscillation when (22) and (23) are satisfied. According to Routh–Hurwitz principle, the essential and adequate ailments for local steadiness of $E(x^{*},y^{*})$ are $A_1>0$; $B_1>0.$

Theorem 1.

The point $E(x^{*},y^{*})$ is locally asymptotically stable in the attendance of transmission if $x^{*}+p^2(D_1+D_2)>0$ and $p^4{D}_1{D}_2+D_2{x}^{*}{p}^2>0$.

The above statement is followed immediately by R–H criteria.

Theorem 2.

(i) The system in the absence of spatiotemporal attributes at the inner steady state $E(x^{*},y^{*})$ attains steadiness, then the corresponding uniform steady state of the model system (16) and (17) in the presence of spatiotemporal attributes also attains steadiness.

(ii) If the inner steady state $E_3(x^{*},y^{*})$ of the non-spatial heterogeneity system is unstable, then the respective steady state of the spatiotemporal model system (16) and (17) under initial and boundary settings and attain steadiness by increasing or decreasing the spatiotemporal attributes suitably.

Proof. 

Let us define the function $V_l\left(t\right)={\int }_0^{R}V\left(x,y\right)du$, where $V\left(x,y\right)$ is defined in Stability analysis section. Differentiating $V_l$ w.r.t. t along the solutions of the diffusive model (16) and (17), we get,

$$V_l^{\prime }\left(t\right)={\int }_0^R\left(V_x{x}_t+V_y{y}_t\right)du=I_R+I_D,$$

(24)

where

$$I_R={\int }_0^R{V}^{\prime }\left(t\right)dx,I_D={\int }_0^R\left(D_1{V}_x{x}_{uu}+D_2{V}_y{y}_{uu}\right)du.$$

Using the analysis in [14], we get

$I_D=-D_1{\int }_0^R{V}_{xx}{\left(x_u\right)}^{2}du-D_2{\int }_0^R{V}_{yy}{\left(y_u\right)}^{2}du$

$$=-D_1{\int }_0^R(x^{*}/x^2){\left(x_u\right)}^{2}du-D_2{\int }_0^R(y^{*}/y^2){\left(y_u\right)}^{2}du.$$

(25)

From (24) and (25), it is observed that if $I_R<0$, then $V_l^{{}^{\prime }}\left(t\right)$ is negative. If $I_R>0$, then it is clearly showing if there is an increase in the spatiotemporal attributes $D_1$ and $D_2$ adequately huge numeral, $V_l^{{}^{\prime }}\left(t\right)$ as $-ve$. Henceforth are the succeeding portion of the theorem grasps. □

5. Impact of Environmental Noise on the Model System

Ecological systems are characteristically forced by a number of drivers such as climate and natural disturbances that are not constant in time but fluctuate. With the exception of processes dominated by deterministic oscillations, a significant part of environmental variability is random because of the uncertainty intrinsic in weather patterns, climate fluctuations, and episodic disturbances like earth quakes, landslides, fires, insect outbreaks, etc. The recurrence of random drivers in bio-geophysical processes motivates the study of how a stochastic environment may affect and determine the dynamics of natural systems. We now begin introducing noise on the proposed model (5) to analyze the role of random environmental fluctuations on stability. The random fluctuations make the parameters of the model to oscillate about their average values. We consider such randomness to the model (5) by incorporating additive white noises. The white noise perturbation included will change any parameter $\nu$ of the model as $\nu +{\alpha }_{ii}{\psi }_i\left(t\right)$, where ${\alpha }_{ii}$ is the amplitude of the noise and ${\psi }_i\left(t\right)$ is a Gaussian white noise process at time t, but the deterministic and stochastic models have same equilibria which will also now fluctuate about their mean states. By considering the randomly fluctuating driving forces in the form of additive noise to the model (5), we get the following stochastic model

$$\frac{dx}{dt}=x(1-x)-\alpha \frac{xy}{x+y}+{\alpha }_{11}{\psi }_1\left(t\right),$$

(26)

$$\frac{dy}{dt}=-\beta y+\gamma \frac{xy}{x+y}+{\alpha }_{22}{\psi }_2\left(t\right),$$

(27)

where ${\alpha }_{11},{\alpha }_{22}$ are real constants, and $\psi \left(t\right)=\left[{\psi }_1\left(t\right),{\psi }_2\left(t\right)\right]$ is a two dimensional Gaussian white noise process agreeable $E\left[{\psi }_i\left(t\right)\right]=0;i=1,2$; $E\left[{\psi }_i\left(t\right){\psi }_j\left(t^{\prime }\right)\right]={\delta }_{ij}\delta \left(t-t^{\prime }\right);i,j=1,2$, where ${\delta }_{ij}$ is the Kronecker symbol; $\delta$ is the delta-Dirac function.

In this analysis, we focus on the dynamics of the model (26) and (27) at the interior equilibrium point $E\left(x^{*},y^{*}\right)$ only according to the method introduced by Nisbet and Gurney [18] and Carletti [19]. Let $x\left(t\right)=u_1\left(t\right)+S^{*};$$y\left(t\right)=u_2\left(t\right)+P^{*}$ and by considering only the consequence of linear stochastic perturbations.

Hence, the model (26) and (27) reduces to the following linear system

$$u_1^{{}^{\prime }}\left(t\right)=-u_1\left(t\right)S^{*}+{\alpha }_{11}{\psi }_1\left(t\right),$$

(28)

$$u_2^{{}^{\prime }}\left(t\right)={\alpha }_{22}{\psi }_2\left(t\right).$$

(29)

Taking the Fourier transform of (28) and (29), we get,

$${\alpha }_{11}{\tilde{\psi }}_1\left(\omega \right)=\left(i\omega \right){\tilde{u}}_1\left(\omega \right)+S^{*}{\tilde{u}}_1\left(\omega \right),$$

(30)

$${\alpha }_{22}{\tilde{\psi }}_2\left(\omega \right)=\left(i\omega \right){\tilde{u}}_2\left(\omega \right).$$

(31)

The matrix form of (30) and (31) is

$$M\left(\omega \right)\tilde{u}\left(\omega \right)=\tilde{\psi }\left(\omega \right),$$

(32)

where $M\left(\omega \right)=\left(\begin{array}{cc}A\left(\omega \right)& B\left(\omega \right)\\ C\left(\omega \right)& D\left(\omega \right)\end{array}\right)$; $\tilde{u}\left(\omega \right)=\left[\begin{array}{c}{\tilde{u}}_1\left(\omega \right)\\ {\tilde{u}}_2\left(\omega \right)\end{array}\right]$; $\tilde{\psi }\left(\omega \right)=\left[\begin{array}{c}{\alpha }_1{\tilde{\psi }}_1\left(\omega \right)\hfill \\ {\alpha }_2{\tilde{\psi }}_2\left(\omega \right)\hfill \end{array}\right]$;

$$A\left(\omega \right)=i\omega +S^{*};B\left(\omega \right)=0;C\left(\omega \right)=0;D\left(\omega \right)=i\omega .$$

(33)

Hence, the solution of (32) is given by $\tilde{u}\left(\omega \right)=K\left(\omega \right)\tilde{\psi }\left(\omega \right)$, where

$$K\left(\omega \right)={\left[M\left(\omega \right)\right]}^{-1}.$$

(34)

The solution components of (34) are given by

$${\tilde{u}}_i\left(\omega \right)=\sum _{j=1}^2{K}_{ij}\left(\omega \right){\tilde{\psi }}_j\left(\omega \right);i=1,2.$$

(35)

The spectrum of $u_i,i=1,2$ are given by $S_{u_i}\left(\omega \right)={\sum }_{j=1}^2{\alpha }_j{\left|K_{ij}\left(\omega \right)\right|}^2;i=1,2$.

Hence, the intensities of fluctuations in the variable $u_i,i=1,2$ are given by

$${\sigma }_{u_i}^2=\frac{1}{2\pi }\sum _{j=1}^2{\int }_{-\infty }^{\infty }{\alpha }_j{\left|K_{ij}\left(\omega \right)\right|}^{2}d\omega ;i=1,2.$$

From (35), we obtain

${\sigma }_{u_1}^2=\frac{1}{2\pi }\left\{{\int }_{-\infty }^{\infty }{\alpha }_{11}{\left|\frac{D\left(\omega \right)}{\left|M\left(\omega \right)\right|}\right|}^{2}d\omega +{\int }_{-\infty }^{\infty }{\alpha }_{22}{\left|\frac{B\left(\omega \right)}{\left|M\left(\omega \right)\right|}\right|}^{2}d\omega \right\}$

$${\sigma }_{u_2}^2=\frac{1}{2\pi }\left\{{\int }_{-\infty }^{\infty }{\alpha }_1{\left|\frac{A\left(\omega \right)}{\left|M\left(\omega \right)\right|}\right|}^{2}d\omega +{\int }_{-\infty }^{\infty }{\alpha }_2{\left|\frac{C\left(\omega \right)}{\left|M\left(\omega \right)\right|}\right|}^{2}d\omega \right\},$$

where $\left|M\left(\omega \right)\right|=R\left(\omega \right)+iI\left(\omega \right)$; $R\left(\omega \right)=-{\omega }^2$; $I\left(\omega \right)=\omega S^{*}$.

If we consider the noise effect on any one of the species, which is with either ${\alpha }_{11}=0$ or ${\alpha }_{22}=0$ then we have ${\sigma }_{u_1}^2=0$; ${\sigma }_{u_2}^2=0$ if ${\alpha }_{11}=0$; ${\sigma }_{u_1}^2=\frac{{\alpha }_{11}}{2\pi }{\int }_{-\infty }^{\infty }\frac{{\left(i\omega \right)}^2}{R^2\left(\omega \right)+I^2\left(\omega \right)}d\omega$; ${\sigma }_{u_2}^2=\frac{{\alpha }_{11}}{2\pi }{\int }_{-\infty }^{\infty }\frac{{(i\omega +S^{*})}^2}{R^2\left(\omega \right)+I^2\left(\omega \right)}d\omega$ if ${\alpha }_{22}=0$.

The population variances point out the stability of population for smaller values of mean square fluctuations, while the larger values of population variances indicate the instability of the populations.

6. Computer Simulations

We want to verify the analytical results obtained in the previous section for the stability of System (7). We numerically calculate the maximal Lyapunov exponent and show the chaotic behavior of the system for the presence of the positive Lyapunov exponent of the system. And also we want to analyze the sensitive dependence of the system trajectories on initial conditions.

At first, we verify that the system moves to a stable state for a particular values of the parameters. We take the parameters value $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$ satisfying the stability criteria listed in the previous section. Then, the 2-dimensional map (7) exhibits a stable movement towards the endemic equilibrium point $E_I(x_I,y_I)=(0.0776,0.3290)$ (see Figure 1a).

The eigenvalues of the system are complex numbers for the same parameter value, and its absolute value is $\left|\lambda \right|=0.9704<1$. This fact also suggests that the system is stable about the endemic equilibrium point $E_I$. At this stage, we want to state that, for the adjoining Figure 1a, the phase-portraits are obtained from several initial states. All the initial populations converge to a single state (point) $E_I$, the attractor. The basin of attractor is the set of all initial points which converge to $E_I$. Here, the basin of attraction is $B=\left\{\right(0.08,0.30),(0.07,0.34),(0.07,0.28),(0.075,0.36),(0.09,0.34),(0.065,0.28),(0.085,0.38\left)\right\}$, and the corresponding attractor is $E_I(0.0776,0.3290)$.

If we increase the consumption ability $\alpha$ of the predator, the system remains stable asymptotically to $E_I$. There is a value of $\alpha =\widehat{\alpha }$ at which the system exhibits a stable periodic oscillation. We take $\alpha =1.2105(>\widehat{\alpha })$ when the system exhibits a periodic oscillation as illustrated in Figure 1b. We have taken here several initial states from the basin of attraction B, and we have obtained the phase-portraits in Figure 1b. Then the absolute value of the eigenvalue is $\left|\lambda \right|=1.0000$. So, at $\alpha =\widehat{\alpha }$ the system has a Hopf type bifurcation where the stable solution becomes periodic. We now tune up the consumption ability of the predator. At $\alpha =1.2223$ the periodic oscillation losses its stability, and the phase-portrait reduces to an invariant closed curve as illustrated in Figure 2 around $E_I(x_I,y_I)=(0.0776,0.3290)$. Then, the absolute value of the eigenvalue is $\left|\lambda \right|=1.0048>1$.

6.1. Order of Chaos by Lyapunov Exponent

We now calculate the maximal Lyapunov exponent of the system, which can give information on chaotic behavior of the system of the 2-dimensional map. A positive Lyapunov exponent indicates that the system exhibits a chaos.

We first give here the mathematical definition of the Lyapunov exponent. We consider the 2-dimensional mapping (7). A point ${\mathbf{X}}_0$ is called an n-cycle point if

$$T^n\left({\mathbf{X}}_0\right)\equiv ToToTo\dots oT\left({\mathbf{X}}_0\right)={\mathbf{X}}_0.$$

The periodic point is stable if the whole orbit $\{{\mathbf{X}}_0,{\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_{n-1}\}$ is stable. For an n-cycle,

$$\mid \frac{d{T}^n}{d\mathbf{X}}\left({\mathbf{X}}_0\right)\mid =\prod _{j=0}^{n-1}\mid \frac{dT}{d\mathbf{X}}\left({\mathbf{X}}_0\right)\mid .$$

(36)

The nth root of this quantity,

$$\Lambda ={\left(\prod _{j=0}^{n-1}\mid \frac{dT}{d\mathbf{X}}\left(X_0\right)\mid \right)}^{1/n},$$

(37)

is the geometric mean rate of stretching along the entire orbit, and this is defined as characteristic multiplier. If for large n the logarithm of the above exists, then it is called the Lyapunov exponent based at ${\mathbf{X}}_0$. We have calculated the maximal Lyapunov exponent, in which its positive value is the indicator of chaos.

We calculate the maximal Lyapunov exponents for several values of $\alpha$. We have observed that the Lyapunov exponent is negative for some values of $\alpha$, and it becomes positive for some other values of $\alpha$. So, chaos occurs in the system. We have shown the different values of the maximal Lyapunov exponent for different values of the consumption ability $\alpha$ of the predator in Figure 3.

6.2. Existence of Orbit Due to Sensitivity

We analyzed the sensitive dependence of the 2-dimensional map (7). We considered two neighboring initial conditions $(x_0,y_0)$ and $(x_0+\delta ,y_0)$ and also the initial conditions $(x_0,y_0)$ and $(x_0,y_0+\delta )$, where $\delta$ is very small quantity. This consideration gives us two neighboring orbits in which behavior after some iterations are observed.

First, we take two neighboring trajectories with initial state $(0.02,0.085)$ and $(0.021,0.085)$. We considered a small change in the x-coordinate of quantity $0.001$ only. The time series evolution of the two initial trajectories is shown in Figure 4a for prey and Figure 4b for predator. Initially, the two trajectories overlap, but after some iterations, they show a clear distinction. Secondly, we took two neighboring trajectories with initial state $(0.02,0.085)$ and $(0.02,0.086)$. We considered a small change in the y-coordinate of quantity $0.001$ only. The time series evolution of the two initial trajectories is shown in Figure 5a for prey and Figure 5b for predator. Initially, the two trajectories overlap, but after some iterations, they show a clear distinction.

Thus, we can find a sensitive dependence on initial conditions in the 2-dimensional map. This implies that the long term prediction about the system is not possible. Even in short term prediction of the populations, the small error in the initial conditions can magnify and prediction becomes worthless.

Furthermore, the analytical findings observed in diffusion dynamics and noise process through numerical simulations as follows:

Example 1.

For the parameters $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$ and diffusion coefficients are $D_1=0.00001$, $D_2=0.00002$. See Figure 6.

Example 2.

For the parameters $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$ and diffusion coefficients are $D_1=0.001$, $D_2=0.002$. See Figure 7.

Example 3.

For the parameters $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$ and diffusion coefficients are $D_1=1$, $D_2=2$. See Figure 8.

Example 4.

For the parameters $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$. See Figure 9.

Example 5.

For the parameters $\alpha =1.2105$, $\beta =1.21$, $\gamma =1.20$. See Figure 10.

Example 6.

For the parameters $\alpha =2.2223$, $\beta =2.21$, $\gamma =2.20$. See Figure 11.

Example 7.

For the parameters $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$. See Figure 12.

Example 8.

For the parameters $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$, ${\alpha }_{11}=1$, ${\alpha }_{22}=1$. See Figure 13.

Example 9.

For the parameters $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$, ${\alpha }_{11}=5$, ${\alpha }_{22}=30$. See Figure 14.

Example 10.

For the parameters $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$, ${\alpha }_{11}=5$, ${\alpha }_{22}=60$. See Figure 15.

7. Concluding Remarks

To study the discrete type ratio-dependent model, we obtained several results. We considered the Michaelis–Menten type functional response. In finding the fixed point, we saw that the endemic fixed point exists when the consumption ability $\alpha$ of the predator precedes some critical value ${\alpha }_1=\gamma /(\gamma -\beta )$ with $\beta <\gamma$.

We saw that the origin $E_0(0,0)$ is stable if the consumption ability $\alpha$ of the predator exceeds some critical value ${\alpha }_2=max\left\{4-\frac{\gamma }{\beta },12-\frac{3\gamma +8}{\beta },\frac{4(1-2\beta )+\gamma }{1-\beta }\right\}$ when the death rate $\beta$ of the predator precedes 1. When the death rate $\beta$ of the predator exceeds 1, the stability criteria demands the consumption ability $\alpha$ to belong in $\left(max\left\{4-\frac{\gamma }{\beta },12-\frac{3\gamma +8}{\beta }\right\},\frac{4(2\beta -1)-\gamma }{\beta -1}\right)$. According to the mathematical analysis, the axial fixed point $E_A$ is stable if the death rate $\beta$ of the predator lies between $\gamma$ and $2+\gamma$. Otherwise, it is a saddle point. The endemic equilibrium state $E_I$ is stable if $(\alpha ,\beta ,\gamma )\in V^S$ and unstable if $(\alpha ,\beta ,\gamma )\in V^U$. Here, $V^S$ and $V^U$ are the regions in the positive octant of $\alpha \beta \gamma$-space.

The 2-dimensional discrete map exhibits a chaotic behavior as the maximal Lyapunov exponents are positive for some values of $\alpha$. In the continuous time model, such a irregular behavior is not possible.

The dependence of the time series of two species on the initial conditions is evident from the numerical analysis of the model. This type of dependence is not possible in continuous-time model, and the population have no possibility of evolving two different final asymptotic states. For a chaotic trajectory, the nearby trajectories begin to diverge.

The discrete-time model has totally different behavior than the continuous-time model which presents in Figure 9, Figure 10, Figure 11 and Figure 12. For a low or reasonable death rate in predator system attains stable with attractive oscillations, which is observed in Figure 12 (with set of parameters $\alpha =1.14$, $\beta =0.21$, $\gamma =1.10$). For increased death rates in predator, system dynamics and its stability is shown in Figure 10 and Figure 11 at various carrying capability and death rates in predator.

Analytical findings and results in diffusion process are quite interesting in graphical view. Diffusion dynamics are explored highly in simulation. Figure 6a,b are the graphical outcomes in the process of Diffusion with the coefficients $D_1=0.00001$, $D_2=0.00002$ at the interior equilibrium point. Figure 7a,b are the graphical outcomes in the process of Diffusion with the coefficients $D_1=0.001$, $D_2=0.002$ at the interior equilibrium point. Figure 8a,b are the graphical outcomes in the process of Diffusion with the coefficients $D_1=1$, $D_2=2$ at the interior equilibrium point. In all these figures, prey exhibits rich dynamics. Figure 6a,b show clearly, at very low value of diffusion coefficient, the system attains stable rapidly. Figure 7a,b show clearly, at low value of diffusion coefficient, the system attains stable moderately. Figure 8a,b show clearly, at high value of diffusion coefficient, the system attains stable little lately. Diffusion process and graphical views are playing a different and attractive role on the proposed model, as well as on any other real world complex scenario.

The additional or dynamical forces, like white noise, exhibit interesting results in the process of stochasticity. Stochastic study on any model enhances its stability performance. The perturbation techniques used in the stochastic process help us to strengthen the steadiness of the proposed model. Numerical simulation helps with the rich dynamics of the stochastic study. The analytical findings in the stochastic process are explored attractively through computer simulations in Figure 13, Figure 14 and Figure 15. It is very clear, in graphs at high amplitudes of prey and predator, the system is highly oscillatory, and at low amplitudes of prey and predator (with stochastic parameters ${\alpha }_{11}=1$, ${\alpha }_{22}=1$), in Figure 13, the system attains stablility with low oscillations. Figure 14 (with stochastic parameters ${\alpha }_{11}=5$, ${\alpha }_{22}=30$) and Figure 15 (with stochastic parameters ${\alpha }_{11}=5$, ${\alpha }_{22}=60$) show that the system attains stablility with high oscillations. Thus, all the analytical findings in different segments are observed, supported, and highly elevated in numerical simulations.