Stability Switching in Lotka-Volterra and Ricker-Type Predator-Prey Systems with Arbitrary Step Size

Abstract

Dynamical properties of numerically approximated discrete systems may become inconsistent with those of the corresponding continuous-time system. We present a qualitative analysis of the dynamical properties of two-species Lotka-Volterra and Ricker-type predator-prey systems under discrete and continuous settings. By creating an arbitrary time discretisation, we obtain stability conditions that preserve the characteristics of continuous-time models and their numerically approximated systems. Here, we show that even small changes to some of the model parameters may alter the system dynamics unless an appropriate time discretisation is chosen to return similar dynamical behaviour to what is observed in the corresponding continuous-time system. We also found similar dynamical properties of the Ricker-type predator-prey systems under certain conditions. Our results demonstrate the need for preliminary analysis to identify which dynamical properties of approximated discretised systems agree or disagree with the corresponding continuous-time systems.

1. Introduction

Population density variations of interacting species in ecosystems are often modelled as discrete-time systems using difference equations (such as the Ricker model [1]). Our study is motivated by the dynamical properties of discrete systems with an arbitrary step size. We consider a simple Ricker-type predator-prey system in discrete form, with a unit time step, given by

$$\begin{array}{c}\begin{array}{cc}\hfill N(t+1)& =N\left(t\right)e^{r(1-\frac{N\left(t\right)}{K})-\alpha P\left(t\right)}\hfill \\ \hfill P(t+1)& =P\left(t\right)e^{\alpha \gamma N\left(t\right)-c},\hfill \end{array}\end{array}$$

(1)

where $N\left(t\right)$ and $P\left(t\right)$ denote the prey and predator population at time t, respectively; r is the rate of prey population increase; K is the prey carrying capacity, $\alpha$ is the predator attack rate; $\gamma$ is the conversion rate of eaten prey to sustenance for the predators; and c is the predator starvation rate in the absence of prey [2,3]. The parameters $r,K,\alpha ,\gamma$ and c are real and positive constants. The system (1) is used to calculate the annual population densities after a time step of one year [2]. We extend the Ricker-type model with arbitrary constant step size h as

$$\begin{array}{c}\begin{array}{cc}\hfill N(t+h)& \hfill =N\left(t\right)\left(1+h\left(e^{X\left(t\right)}-1\right)\right)\hfill \\ \hfill P(t+h)& \hfill =P\left(t\right)\left(1+h\left(e^{Y\left(t\right)}-1\right)\right),\hfill \end{array}\end{array}$$

(2)

where

$$\begin{array}{c}\begin{array}{cc}\hfill X\left(t\right)& =r\left(1-\frac{N\left(t\right)}{K}\right)-\alpha P\left(t\right)\hfill \\ \hfill Y\left(t\right)& =\alpha \gamma N\left(t\right)-c.\hfill \end{array}\end{array}$$

(3)

System (2) is the generalised version of (1) that considers the unit increments as a parameter such that system (2) is the same as system (1) if $h=1$.

In this paper, we also consider models such as Lotka-Volterra

$$\begin{array}{c}\begin{array}{cc}\hfill N(t+h)& \hfill =N\left(t\right)\left(1+hX\left(t\right)\right)\hfill \\ \hfill P(t+h)& \hfill =P\left(t\right)\left(1+hY\left(t\right)\right),\hfill \end{array}\end{array}$$

(4)

where h is the discrete step size. The step size h plays a critical role by permitting the users to choose an appropriate time discretisation for each model. We observe that these generalized discrete systems incorporating an arbitrary fixed step size are forward Euler’s approximations of the respective continuous-time systems.

Beyond the step-size selections, the robustness of the model parameters is critical; however, it is often a source of uncertainty in models based on real data. A slight variation of the model parameters may change the equilibria and directly affect the system stability and robustness of solutions of the system [4,5,6]. If the parameters are estimated from data, lack of information and the inability to collect sufficient real-world data in ecological systems can lead to an imprecise set of model parameters [7]. Therefore, investigating a suitable set of parameters that agrees with the selection of stable or unstable dynamics is essential when constructing population models—especially in approximating continuous-time systems [8] (see Acknowledgements for further clarification).

Following the idea of parametrising the step size, ref. [9] derived the stability properties for a continuous-time Lotka-Volterra-type predator-prey system with scaled model parameters and showed unstable and stable population dynamics for derived conditions in step-size selections (see [10,11] for similar studies for Lotka-Volterra-type predator-prey models). To the best of our knowledge, no study has investigated the dynamic inconsistency under arbitrary step size for the Ricker-type ordinary differential equation (ODE) predator-prey model [12], as given by

$$\begin{array}{c}\begin{array}{cc}\hfill N^{\prime }\left(t\right)& \hfill =N\left(t\right)\left(e^{X\left(t\right)}-1\right)\hfill \\ \hfill P^{\prime }\left(t\right)& \hfill =P\left(t\right)\left(e^{Y\left(t\right)}-1\right).\hfill \end{array}\end{array}$$

(5)

Hence, we study the required conditions for stable and unstable population dynamics of ODE system (5) and their discretised system (2) with generalised step size. Consequently, the results identify similar or different dynamical properties of approximated discrete systems compared to the corresponding continuous-time model.

We perform a comparable study on qualitative analysis of the stability properties of a commonly used continuous-time population model, a logistic growth Lotka-Volterra-type predator-prey system [13]

$$\begin{array}{c}\begin{array}{cc}\hfill N^{\prime }\left(t\right)& \hfill =N\left(t\right)X\left(t\right)\hfill \\ \hfill P^{\prime }\left(t\right)& \hfill =P\left(t\right)Y\left(t\right).\hfill \end{array}\end{array}$$

(6)

Here, the prey population is influenced by the prey natural growth, prey restricted growth in terms of prey carrying capacity and prey death caused by predator attacks. The predator abundance is governed by population growth due to predation and natural death. We indicate that the generalised Lotka-Volterra model (4) is the approximated discrete system of the continuous-time ODE model (6).

The dynamical properties of predator-prey systems of continuous-time models have been studied extensively in the literature (see the examples at [14,15,16]), but the factors that can perturb the dynamical properties (e.g., stabilising or destabilising the system) have not been fully analysed for discrete approximations of continuous-time systems.

We fill the gap of identifying the constraints that deliver similar or different dynamical properties of two continuous-time models and their approximated discrete system in terms of the arbitrary step size under parameter space. More precisely, we derive the stability properties of discrete systems associated with corresponding continuous-time models such that the interconnection of these two systems is purely observable.

The comparison of both approximated (discrete) and actual (continuous) systems reveals the conditions that must be satisfied by the time discretisation. Our study also highlights the importance of the changes in (some) model parameters that direct the system to have stable or unstable population dynamics. Therefore, our work contributes to choosing the best-suited step size for a particular discrete-time system and helps to understand the stabilising and destabilising factors of the continuous system and approximated discretised system under parameter space.

This paper is organised as follows. In Section 2, a qualitative study to construct stability constraints is performed through a Jacobian analysis of the discrete systems followed by that of the continuous systems. In Section 3, we derive the dynamical properties of discrete systems connected to the continuous-time systems and investigate the conditions that lead the system solutions to become stable in terms of selecting suitable sets of model parameters and determining the effects of the step size in time discretisation. We then demonstrate population dynamics to justify our theoretical findings through numerical simulations in Section 4.

2. Equilibrium Stability and Dynamics

To understand the interconnection of dynamical properties in generic discretised systems with continuous-time systems under an arbitrary step size, we investigate the equilibrium points and the possible stable states of all Ricker-type and Lotka-Volterra systems mentioned in this paper. We first examine the system dynamics at the fixed points of the generalised discrete systems (2) and (4).

In order to develop stability conditions for discrete systems, we first consider the fixed point iteration formulation of a nonlinear map:

$$\chi (t+h)=F\left(\chi \right(t\left)\right),$$

(7)

where $\chi \left(t\right)\in {\mathbb{R}}^m$, and F has a Lipschitz condition. A point $\eta$ is said to be a fixed point satisfying $\eta =F\left(\eta \right)$, where F is a map, such that $F:I\to I$, and I is a region in ${\mathbb{R}}^m$. It is proven that the fixed point $\eta$ is asymptotically stable if there exists a norm such that

$$\left|\right|J_F\left|\right|<1,$$

(8)

where $J_F$ is the Jacobian matrix of the discrete system evaluated at $\eta$ [17]. Since we consider two dimensional systems, we have the following characterisations at $\eta$ in terms of the two eigenvalues ${\lambda }_1$ and ${\lambda }_2$ of the Jacobian matrix of the iterated map:

(i)

$\left|{\lambda }_i\right|<1,i=1,2$; a sink, locally asymptotically stable.

(ii)

$\left|{\lambda }_i\right|>1,i=1,2$; a source.

(iii)

One of $\left|{\lambda }_i\right|>1$ and other $\left|{\lambda }_i\right|<1,i=1,2$; a saddle.

(iv)

One of $\left|{\lambda }_i\right|=1$ and other $\left|{\lambda }_i\right|\ne 1,i=1,2$; non-hyperbolic.

Note that our investigation is based on this eigenvalue classification. Furthermore, we can find flip, saddle and Hopf bifurcations when we traverse the boundaries of different stability domains.

The equilibrium points for the Ricker-type and Lotka-Volterra discrete models can be obtained by solving

$$\begin{array}{c}\begin{array}{cc}\hfill N(t+h)& \hfill =N\left(t\right)\hfill \\ \hfill P(t+h)& \hfill =P\left(t\right)\hfill \end{array}\end{array}$$

(9)

for N and P. For the generalised Ricker-type discrete model, (9) implies $X\left(t\right)=Y\left(t\right)=0$, and this also holds for the discrete Lotka-Volterra model, (4). The fixed points of both the discrete formulations of the generalised Ricker-type discrete model and the discrete Lotka-Volterra model ((2) and (4), respectively) are the same, namely $E_1\equiv (0,0)$, $E_2\equiv (K,0)$ and $E_3\equiv \left(\frac{c}{\alpha \gamma },\frac{r}{\alpha }\left(1-\frac{c}{K\alpha \gamma }\right)\right)$. Further, these are the same equilibrium points obtained for the continuous time models since the same condition $X\left(t\right)=Y\left(t\right)=0$ was satisfied, and we use the same notation. We denote the Jacobians of discrete maps in (2) and (4) by ${\widehat{J}}_R$ and ${\widehat{J}}_{LV}$, respectively. Then, some analysis gives, after omitting the dependence on t,

$${\widehat{J}}_R=\left[\begin{array}{cc}1+h\left(e^X-1-\frac{r}{K}N{e}^X\right)& -\alpha Nh{e}^X\\ \alpha \gamma Ph{e}^Y& 1+h\left(e^Y-1\right)\end{array}\right]$$

and

$${\widehat{J}}_{LV}=\left[\begin{array}{cc}1+h\left(X-\frac{r}{K}N\right)& -\alpha hN\\ \alpha \gamma hP& 1+hY\end{array}\right].$$

By using Jacobians, we derive the connection of the stability properties of the discrete systems to their respective ODE systems in the next section.

We then classify the equilibria of the continuous-time models (5) (Ricker-type) and (6) (Lotka-Volterra) based on the stability since unstable and stable equilibria behave differently in population dynamics. Both systems (5) and (6) have fixed points, namely $E_1\equiv (0,0)$, $E_2\equiv (K,0)$ and $E_3\equiv \left(\frac{c}{\alpha \gamma },\frac{r}{\alpha }\left(1-\frac{c}{K\alpha \gamma }\right)\right)$. Then, the Jacobian matrix of (5) is

$$J_R=\left[\begin{array}{cc}{e}^X-1-\frac{r}{K}N{e}^X& -\alpha N{e}^X\\ \alpha \gamma P{e}^Y& e^Y-1\end{array}\right]$$

and the Jacobian matrix for (6) is

$$J_{LV}=\left[\begin{array}{cc}X-\frac{r}{K}N& -\alpha N\\ \alpha \gamma P& Y\end{array}\right].$$

We observe similar stability properties for both ODE models even though the Jacobian matrices are different. The Jacobian matrices evaluated at $E_1=(0,0)$ are

$$J_R{\mid }_{E_1}=\left[\begin{array}{cc}{e}^r-1& 0\\ 0& e^{-c}-1\end{array}\right]\mathrm{and}{J}_{LV}{\mid }_{E_1}=\left[\begin{array}{cc}r& 0\\ 0& -c\end{array}\right].$$

Thus, the eigenvalues at $E_1$ are

$$\begin{array}{c}\begin{array}{cc}\hfill {\lambda }_R{\mid }_{E_1}& =\left\{e^r-1,e^{-c}-1\right\}\hfill \\ \hfill {\lambda }_{LV}{\mid }_{E_1}& =\left\{r,-c\right\}\hfill \end{array}\end{array}$$

(10)

where ${\lambda }_R$ and ${\lambda }_{LV}$ are the eigenvalues for system (5) and (6), respectively. $E_1$ is unstable for both models regardless of any choices of parameter values since it is a saddle point such that one eigenvalue is positive and one is negative. This means that the population never returns to $E_1$ after a small deviation of the population variation.

The Jacobian matrices evaluated at $E_2=(K,0)$ are

$$J_R{\mid }_{E_2}=\left[\begin{array}{cc}-r& -K\alpha \\ 0& e^{\alpha \gamma K-c}-1\end{array}\right]\mathrm{and}{J}_{LV}{\mid }_{E_2}=\left[\begin{array}{cc}-r& -K\alpha \\ 0& \alpha \gamma K-c\end{array}\right]$$

and so

$$\begin{array}{c}\begin{array}{cc}\hfill {\lambda }_R{\mid }_{E_2}& =\left\{-r,e^{\theta }-1\right\}\hfill \\ \hfill {\lambda }_{LV}{\mid }_{E_2}& =\left\{-r,\theta \right\}\hfill \end{array}\end{array}$$

(11)

where $\theta =\alpha \gamma K-c$. The $\theta$ value indicates a condition to determine the properties of the eigenvalues. Therefore, investigations for stability properties are presented in terms of $\theta$ where appropriate. Then, $E_2$ is asymptotically stable for both models if the parameters satisfy $\theta <0$ (i.e., $\alpha \gamma K<c$) since both eigenvalues are then real and negative. This means that, if there is a small deviation of the population densities away from $E_2$, the prey and predators can return to prey-carrying capacity ($N=K$) and no predators, respectively.

This implies that the prey population can be sustained by reaching its maximum carrying capacity without the presence of predators even if there are a small number of predators present in the ecosystem or if the prey population experiences high mortality. Moreover, this case shows the predator-prey existence under low predator populations where their ability to survive through food availability is determined by $K<\frac{c}{\alpha \gamma }$. On the other hand, for both models, $E_2$ is an unstable saddle point if $\theta >0$ and if $\theta =0$ $E_2$ is a non-hyperbolic point where the system dynamics depend on the nonlinear terms of the model equations since it cannot be predicted from the eigenvalue analysis of the Jacobian matrix.

The Jacobian matrices evaluated at $E_3=\left(\frac{c}{\alpha \gamma },\frac{r}{\alpha }\left(1\frac{c}{K\alpha \gamma }\right)\right)$ are

$$J_R{\mid }_{E_3}=J_{LV}{\mid }_{E_3}=\left[\begin{array}{cc}-\frac{rc}{\alpha \gamma K}& -\frac{c}{\gamma }\\ r\gamma \left(1-\frac{c}{\alpha \gamma K}\right)& 0\end{array}\right],$$

and since these Jacobian matrices are the same for both ODE models, we obtain identical results for stability analysis. The eigenvalues calculated for $E_3$ satisfy the characteristic polynomial

$${\lambda }^2+T\lambda +D=0$$

(12)

where $T=\frac{rc}{\alpha \gamma K}$ and $D=c\left(r-T\right)=\theta T$. Thus,

$${\lambda }_R{\mid }_{E_3}={\lambda }_{LV}{\mid }_{E_3}=\left\{{\lambda }_1,{\lambda }_2\right\}$$

(13)

where ${\lambda }_j=\frac{-T\pm \sqrt{T^2-4D}}{2}=\frac{-T\pm \sqrt{T\left(T-4\theta \right)}}{2}$, $j=\{1,2\}$. Note that, if $T\in (0,4\theta )$, the imaginary component for the continuous-time models is $\frac{\sqrt{T(T-4\theta )}}{2}i$, and the larger the imaginary component is, the more oscillatory are the dynamics (note that the maximum imaginary component occurs when $T=2\theta$). The system stability status at a particular equilibrium point can be observed by looking at the sign of the eigenvalues and whether they are real or complex. We can conclude that $E_3$ is asymptotically stable if

$$\theta >0,$$

(14)

with oscillatory dynamics if $T\in (0,4\theta )$, and $E_3$ is an unstable saddle point if $\theta <0$. $E_3$ is the only non-trivial equilibrium point that has positive populations for both prey and predators. Note that, if $\theta =0$, we can have a non-hyperbolic property.

3. Deriving Connections of Dynamical Properties in Discrete and Continuous Systems

The connections between discrete systems and ODE systems are derived under arbitrary step size. Stability analysis for the discrete-time systems is simplified using derivations from the respective ODE systems where necessary. These stability constraints are observed through a Jacobian analysis, and the factors that could affect the system dynamics are analysed through variations of model parameters and by defining a particular range of step size.

We observe that

$${\widehat{J}}_R=I+h{J}_R$$

and

$${\widehat{J}}_{LV}=I+h{J}_{LV}.$$

Hence, $\lambda ({\widehat{J}}_R)=1+h\lambda \left(J_R\right)$ and $\lambda ({\widehat{J}}_{LV})=1+h\lambda \left(J_{LV}\right)$, where the eigenvalues of $J_R$ and $J_{LV}$ are given in (10), (11) and (13). Thus, from the eigenvalue analysis, we have locally asymptotic stability for discrete mappings if $|1+h\lambda (J_R\left)\right|<1$ and similarly for $\lambda \left(J_{LV}\right)$. In the case that $\lambda \left(J_R\right),\lambda \left(J_{LV}\right)<0$, this leads to $-h\lambda \left(J_R\right)<2$ and $-h\lambda \left(J_{LV}\right)<2$. Therefore, depending on the step size and the eigenvalues of the continuous-time system, new conditions exist in discretised systems for stable population dynamics.

From the eigenvalue classification, $E_1$ is not asymptotically stable, but there is a saddle point if $h(1-e^{-c})<2$ or $0<hc<2$ for Ricker-type and Lotka-Volterra discrete maps, respectively. For $\theta <0$, $E_2$ is asymptotically stable for the Ricker-type discrete model, if $h<\left\{\frac{2}{r},\frac{2}{1-e^{\theta }}\right\}$, and for the discrete Lotka-Volterra model, if $h<\left\{\frac{2}{r},\frac{-2}{\theta }\right\}$. Furthermore, if $\theta =0$ and $h\ne \frac{2}{r}$, $E_2$ is non-hyperbolic for both models. For $\theta <0$, $E_2$ is a saddle point for the Ricker-type discrete model if one of the following conditions holds

(i)

$\theta <-r\mathrm{and}\frac{2}{r}<h<\frac{2}{1-e^{\theta }},$

(ii)

$\theta >-r\mathrm{and}\frac{2}{r}>h>\frac{2}{1-e^{\theta }}$.

For $\theta <0$, $E_2$ is a saddle point for the Lotka-Volterra discrete model if one of the following conditions holds

(i)

$\theta <ln(1-r)\mathrm{and}\frac{2}{r}<h<\frac{-2}{\theta },$

(ii)

$\theta >ln(1-r)\mathrm{and}\frac{2}{r}>h>\frac{-2}{\theta }.$

In the case of the fixed point $E_3$, we can classify the stability associated with the two cases if the eigenvalues of $J_R$ and $J_{LV}$ are given as ${\lambda }_1$ and ${\lambda }_2$. From (12), the eigenvalues satisfy ${\lambda }^2+T\lambda +D=0$, where $T=\frac{rc}{\alpha \gamma K}$ and $D=c(r-T)=\theta T$. For $\theta >0$,

(i)

If ${\lambda }_1$ and ${\lambda }_2$ are real and negative where ${\lambda }_1={\lambda }_2$, then $E_3$ is asymptotically stable if

$$h<\frac{T}{4}.$$

This happens only if $T=4\theta$.

(ii)

If both ${\lambda }_1$ and ${\lambda }_2$ are real and negative where ${\lambda }_1\ne {\lambda }_2$, then $E_3$ is asymptotically stable if the step size satisfies

$$h<\left\{-\frac{2}{{\lambda }_1},-\frac{2}{{\lambda }_2}\right\}.$$

This happens if $T-4\theta >0$.

(iii)

If both ${\lambda }_1$ and ${\lambda }_2$ are complex conjugate eigenvalues (say $a\pm ib$), then, from the above, $a^2+b^2=D=\theta T$ and $T=-2a$. This occurs when $T-4\theta <0$. With these complex eigenvalues, the population dynamics lead to oscillations with time. Then, the bound for the step size is

$$h<\left\{\frac{-2a}{a^2+b^2}\right\}=\frac{T}{D}=\frac{1}{\theta }.$$

(15)

This can only happen if $0<1+hT(h\theta -1)$. Note that, if $h=\frac{1}{\theta }$, then both eigenvalues have magnitude one.

For $\theta <0$, say ${\lambda }_1>0$ and ${\lambda }_2<0$, then $E_3$ is non-hyperbolic if $h=-\frac{2}{{\lambda }_2}$. For $\theta =0$, $E_3$ is non-hyperbolic if $h\ne \frac{2}{T}$.

Overall, the stability constraints are different in continuous-time models and their corresponding discrete systems. A summary of the stability analysis of all eigenvalues for the Ricker-type and Lotka-Volterra discrete and continuous-time models are given in

Table 1. Stability status for discrete and continuous Ricker-type (RK) and Lotka-Volterra (LV) models, where $E_1\equiv (0,0)$, $E_2\equiv (K,0),E_3\equiv \left(\frac{c}{\alpha \gamma },\frac{r}{\alpha }\left(1-\frac{c}{K\alpha \gamma }\right)\right)$, and ${\lambda }_1,{\lambda }_2$ are eigenvalues of $E_3$ calculated from (13).

Stability	Model	${\mathit{E}}_{\mathbf{1}}$		${\mathit{E}}_{\mathbf{2}}$
		Discrete	Continuous	Discrete	Continuous
Asym. stable	RK	-	-	if $\theta <0$ and $h<\{\frac{2}{r},\frac{2}{1-e^{\theta }}\}$	if $\theta <0$
	LV	-	-	if $\theta <0$ and $h<\{\frac{2}{r},-\frac{2}{\theta }\}$	as above
Non- hyperbolic	RK	if $h=\frac{2}{1-e^{-c}}$	-	if $\theta =0,h\ne \frac{2}{r}$ if $\theta <0,h=\frac{2}{1-e^{\theta }},h\ne \frac{2}{r}$ if $h=\frac{2}{r},\theta \ne 0,\theta \ne ln(1-r)$	if $\theta =0$
	LV	if $h=\frac{2}{c}$	-	if $\theta =0,h\ne \frac{2}{r}$ if $\theta <0,h=-\frac{2}{\theta },\theta \ne 0,\theta \ne -r$ if $h=\frac{2}{r},\theta \ne 0,\theta \ne -r$	as above
Saddle	RK	if $h<\frac{2}{1-e^{-c}}$	always a saddle point	if $\theta <0,\theta <-r,\frac{2}{r}<h<\frac{2}{1-e^{\theta }}$ if $\theta <0,\theta >-r,\frac{2}{r}>h>\frac{2}{1-e^{\theta }}$ if $\theta >0,h<\{\frac{2}{r},\frac{2}{1-e^{\theta }}\}$	if $\theta >0$
	LV	if $h<\frac{2}{c}$	always a saddle point	if $\theta <0,\theta <ln(1-r),\frac{2}{r}<h<-\frac{2}{\theta }$ if $\theta <0,\theta >ln(1-r),\frac{2}{r}>h>-\frac{2}{\theta }$	as above
Stability	Model	${\mathit{E}}_{\mathbf{3}}$
		Discrete			Continuous
Asym. stable	RK	if $\theta >0$, $T=4\theta ,0<h<\frac{T}{4}$ if $\theta >0,T-4\theta >0,h<\{-\frac{2}{{\lambda }_1},-\frac{2}{{\lambda }_2}\}$ if $\theta >0,T-4\theta <0,h<\frac{1}{\theta },0<1+hT(h\theta -1)$			if $\theta >0$
	LV	as above			as above
Non- hyperbolic	RK	if $\theta >0,T-4\theta >0,h=-\frac{2}{{\lambda }_i},h\ne -\frac{2}{{\lambda }_j},i\ne j,i,j=\{1,2\}$ if $\theta <0,h=-\frac{2}{{\lambda }_2},{\lambda }_2<0,{\lambda }_1>0$ if $\theta =0,h\ne \frac{2}{T}$			if $\theta =0$
	LV	as above			as above
Saddle	RK	if $\theta >0,T-4\theta >0,-\frac{2}{{\lambda }_i}<h<-\frac{2}{{\lambda }_j},i\ne j,i,j=\{1,2\}$			if $\theta <0$
	LV	as above			as above

. Stability criteria evaluated at equilibrium point $E_3$ are similar for both Lotka-Volterra and Ricker-type models. At equilibrium point $E_2$, the stability conditions are similar for both models under continuous-time setting only and different otherwise. We consider the case (iii) when $\theta >0$ for further simulations since it is a stable spiral where the population returns to a steady state, $E_3$.

Thus, additional constraints are required for stable population dynamics in approximated discrete systems compared with those in the continuous-time system (which is $\theta >0$). Therefore, the model dynamics of approximated discrete systems depend on the selected step size and (some) model parameters. We use (15) to observe the stability with different step sizes and ranges of model parameters in the next section.

4. Numerical Results

We present a numerical simulation study for the Lotka-Volterra and Ricker-type discrete systems to illustrate the theoretical findings discussed in Section 3. The stability condition (15), where the discrete systems become stable at $E_3$, are demonstrated for some parameter ranges and step sizes. The impacts of model parameter variations on the system that change the model dynamics are then investigated through a few examples. This numerical study clearly demonstrates the changes to the population over time in the long-term scale.

We observed that some parameters impact the system dynamics and can stabilise or destabilise the populations if the time discretisation is fixed. In our discrete-time models, model stability is governed by the parameters $\alpha ,\gamma ,K,c$ and r according to the derivation of (15) given that step size h is fixed. We first observed the stability condition (15) for discrete systems at unit step size $h=1$ by assigning the parameter values described in [2], where $\alpha =0.05,\gamma =0.01,K=2500,c=0.2$ and $r=0.5$. Then, $\theta =\alpha \gamma K-c=\frac{21}{20}$ and $T=\frac{rc}{\alpha \gamma K}=\frac{2}{25}$ (note that we use these parameters for the numerical simulation study of this paper; some model parameters are subject to change depending on the analysis, and changes to the model parameters are given in each figure caption). Then, (15) specifies stability for step size

$$h<\frac{20}{21},$$

where $\theta >0,T-4\theta <0$ and $1+hT(h\theta -1)>0.$ In this case, since the eigenvalues are complex, the size of the imaginary component of the eigenvalue evaluated at equilibrium point $E_3$ is $\frac{\sqrt{206}}{50}$. For this choice of parameters with step size $h=1$, the system is not asymptotically stable at $E_3$ since annual discretisation does not satisfy $h<\frac{20}{21}$. Thus, the system dynamics diverge while oscillating around the equilibrium point $E_3$ (blue curves in Figure 1). If $\alpha =0.048$ with the other parameters the same and a one-year step size $h=1$, then the system may not become asymptotically stable since $\theta =1,T-4\theta <0$ and $1+hT(h\theta -1)>0$.

In this case, the condition for stability in Equation (15) is violated: $h=\frac{1}{\theta }=1$. The populations seem to be converging to $E_3$ but oscillates around $E_3$ (red curves in Figure 1). If $\alpha =0.04$ and $h=1$ with the other parameters the same, then $\theta =\alpha \gamma K-c=\frac{1}{20}>0,T-4\theta <0$ and $1+hT(h\theta -1)>0$; thus, $h<\frac{1}{\theta }$, and the condition (15) is true for this case. Therefore, the populations converge to $E_3=(500,10)$; see the black curves in Figure 1. Thus, there is considerable sensitivity to the choice of parameters and, hence, to the step size, which controls the dynamics of the model.

Moreover, if we choose $c=0.25+ϵ$, rather than $c=0.2$ as previously, where $ϵ$ is a small positive value and the other parameters are as in [2], then

$$\theta =\alpha \gamma K-c=\frac{25}{20}-\frac{5}{20}-ϵ=1-ϵ,$$

and so (15) gives

$$h<\frac{1}{1-ϵ}.$$

Then, for the step size $h=1$, the condition (15) is true such that the system generates populations that converge to the fixed points (except for $ϵ=0$). Therefore, if the step size is fixed, a small deviation of a model parameter or specific parameters that define model stability can make dramatic changes in the population dynamics.

Returning to the case of $c=0.2$, we varied c by $c+\zeta$ for the step size $h=1$ where $\zeta$ is a small positive constant value. The numerical results of population densities over time in Figure 2 confirm the stability conditions developed for discrete models as in Equation (15). For small values of $\zeta$ near zero, the populations are oscillatory diverging; see the blue curves in Figure 2. Special dynamics are observed for $\zeta =0.05$ (red curves in Figure 2), in which case, $c=0.25$ and corresponds to the case of $\zeta =0$ studied previously.

Fixed point convergence is observed for $\zeta =\{0.1,0.2\}$, which satisfies the stability condition in (15); see the black curves in Figure 2. The population dynamics according to the small changes of the parameter c shows how critical the sensitivity of the parameters is. This is valid for other parameters that affect the stability of the discrete model—in our case, $\alpha ,\gamma ,K,c$ and r. Moreover, the numerical simulations of the discrete models support the developed analytical results in discrete systems.

Figure 3 shows the fixed-point convergence region (yellow) for $\beta \in [0,8\times {10}^{-4}]$, where $\beta =\alpha \gamma$. Therefore, parameter c has different boundaries for converging populations and preserving positive predator populations as demonstrated in the red and black curves of Figure 3, respectively.

Special dynamics of system (1) are observed when $\alpha =0.048$, which is on the upper bound of the fixed-point converging region (the red curves in Figure 1). The black curves in Figure 1 represent the converging predator and prey densities to their fixed points as the $\alpha$ is chosen from the yellow region of Figure 3. If $\alpha$ is chosen from the oscillatory divergence region, the predator-prey densities tend to oscillate continuously and increase in amplitude; see the blue curves in Figure 1.

As a guide to selecting a suitable step size h with the required stability property, the impacts of variable step sizes are investigated (see Figure 4). The fixed-point convergence region becomes larger for small step sizes. Figure 4 indicates that larger step sizes are more likely to show oscillatory divergence.

Even though the model structures of the discrete Ricker-type and Lotka-Volterra models are different, the derived stability conditions for these models are similar in certain cases. We examine the frequency plots of the population dynamics for both models with different step sizes. The red and blue curves in Figure 5 provide evidence of the special behaviour when $\theta =\alpha \gamma K-c=21/20$ (with $h=\frac{1}{\theta }$). These curves have a similar pattern that converges to a fixed point at a later time. The continuously converging pattern of these curves shows similar characteristics of the fixed-point convergence at the beginning, and, from the long-time observations, the populations never approach fixed points. Therefore, this is an exceptional case that behaves as an upper bound for the fixed-point convergence region.

5. Discussion and Conclusions

We investigated the stability of the population dynamics of Ricker-type and Lotka-Volterra discrete and continuous-time population models. Based on the Jacobian analysis, important constraints were generated to identify the stability regions of the discrete systems. The generalised discrete systems can be viewed as the Euler numerical scheme of the approximate solution of the continuous-time models.

The discretised solution may or may not have the same properties as observed in the original continuous-time model. Therefore, a qualitative analysis of the model dynamics in both the original continuous-time model and the relevant discretised system is essential when modelling to inform ecological decisions. We showed that the two models have similar conditions to stabilise the system depending on the constraints as shown in Table 1. This novel work increases the understanding of the similar behaviours of two structurally different predator-prey models, in a nonlinear setting.

The derived stability properties for continuous-time and approximated discrete-time solutions are different at each equilibrium point. Therefore, choosing a suitable time discretisation depends on the stability criteria that are determined by the dynamical properties. Our results are outlined in Table 1, which includes the results for asymptotic stability, non-hyperbolic stability and saddle points. We assume the populations in real ecosystems coexist at, or move to, a stable state, which is $E_3$ with asymptotic stability. Therefore, selecting a suitable time discretisation is essential to approximate the respective continuous-time system if the model parameters are given.

We highlighted the significance of understanding the choices of model parameters that impact the stability of the system. Some parameters need to be verified carefully due to the strong effect on system dynamics, such as the parameters determined by $\theta$. Small changes in these parameters lead to large deviations in the population count when the step size of the time discretisation is fixed (e.g., one year). A prior analysis on the impact of selecting a suitable time discretisation that is relevant to selected parameter sets is essential for better performance of the models. We also numerically showed that the small changes in parameter values and step sizes stabilise or destabilise the system and form different dynamics in population densities.

Our results on the Ricker-type discrete model are consistent with the stabilisation that was found in [3]. However, the discrete model stability analysis that shows the significance of dynamical properties over the parameter space goes beyond previous studies on investigating the stabilising and destabilising factors defined in [18,19]. Populations become more stable and converge to the equilibrium point if the step size is small. This generates a time discretisation (into small time intervals) where the discrete systems behave more closely to the continuous systems.

Certainly, to preserve the characteristics in numerical simulations, this idea supports the small step size recommendation in Euler’s scheme, since it is a first-order method [20]. These results lead to selecting suitable values for the step size in terms of preserving the required stability states when approximating discrete systems to respective continuous systems. On the other hand, if the step size is too small, then, on larger systems, the simulations may take a very long time to run, especially over large time intervals. We note that there are non-standard discretisation methods of nonlinear ODEs to preserve the dynamic consistency regardless of the selection of the step size [21,22].

This is not a complete study of the dynamics of discrete mapping. For example, from (7), we note $\chi (t+2h)=F\left(\chi \right(t+h\left)\right)=F\left(F\right(\chi \left(t\right)\left)\right)$. Hence, we can study the fixed points of $F\circ F$, and the ensuing dynamics would lead to period-two dynamics depending on the nature of the map F. The two models will have different dynamics in this regard. The theories on limit cycles and bifurcation analysis [23,24,25] were not considered; however, this can be seen through our numerical simulations. Of course, this can be generalised to period dynamics of any integer order and potentially lead to chaotic dynamics (of iterations of the logistic map). Moreover, parameter-identification methods can be applied to estimate parameters within a chaotic system [26].

The demonstrated stability analysis can be applied to other forms of two-species continuous-time and discrete-time population models in returning more complex dynamical systems, such as controlling species, functional responses, time delay [27,28] and the Allee effect. Our work can be extended to study the dynamics of three or more species systems [29] and to understand the stabilising and destabilising factors before obtaining the model outcomes.

Finally, this work has implications in uncertainty quantification. In this setting, populations of models (with the same structure but different parameter sets) are constructed based on these models satisfying a set of common outputs. If there is sensitivity of the dynamics to the parameters, as is the case here, then this can make the process of uncertainty quantification also sensitive.

In conclusion, the stability criteria for the continuous-time models depend only on the model parameters defined by $\theta$; however, when the continuous-time models are discretised through numerical approximations, the stability criteria depend on the step size along with (some) model parameters. Therefore, the dynamical properties of the original ODE systems are different from the respective discretised systems with the same parameter values unless a suitable step size is defined.

To obtain a better understanding of the model dynamics, the discretised systems should be thoroughly investigated under an eigenvalue analysis to identify suitable step sizes that agree with the model parameters. In real predator-prey systems, prior knowledge of the system behaviour enhances the understanding of the factors that stabilise or destabilise the populations over time. If the parameters are estimated from a data set, this type of theoretical study provides a detailed analysis for selecting a suitable numerical simulation method that discretises the time component with step size.