Home Stability and bifurcation analysis of a modified chemostat model
Article Open Access

Stability and bifurcation analysis of a modified chemostat model

  • Mehtap Lafci Büyükkahraman EMAIL logo
Published/Copyright: August 30, 2025
Download Article (PDF)
Become an author with De Gruyter Brill
Author Information Explore this Subject
Demonstratio Mathematica
From the journal Demonstratio Mathematica Volume 58 Issue 1

Abstract

In this work, we derive a discrete-time chemostat model from a continuous-time population model using the forward Euler method. We first determine the model’s fixed points and analyze their local stability. By applying the central manifold theorem and bifurcation theory, we establish the existence of flip and Neimark-Sacker bifurcations, considering the step size as the bifurcation parameter. Theoretical findings are validated through numerical simulations, which also uncover novel dynamic behaviors.

Keywords: chemostat model; stability; flip bifurcation; Neimark-Sacker bifurcation
MSC 2010: 39A05; 39A28; 39A30; 39A33

1 Introduction

The chemostat is one of the most effective artificial models of the natural environment for studying populations due to its resemblance to the open system structure and temporal continuity of the universe, being a dynamic system with ongoing material inputs and outputs. The input and elimination of nutrients are analogs of the continuous cycle of nutrients in nature. The constant occurrence of non-age-specific death, predation, or emigration in nature is akin to the washout of organisms [1]. Additionally, a well-known type of bioreactor for continuous microbial cultivation is the chemostat. Chemostats are widely used in industrial microbiology and biotechnology for continuous culture experiments, allowing precise control over microbial growth conditions [2–5]. Thanks to its popularity, research on several chemostat models is ongoing. In a chemostat model with pulsed input, Zhao et al. [6] and Jiao et al. [7] examined the extinction and persistence of the model. The coexistence of competing species in reaction-diffusion chemostat model has been investigated by Shi et al. [8]. The global dynamics of a cell quota-based model in chemostat have been studied by Alzahrani et al. [9]. With the aid of Fokker-Planck theory, Baratti et al. [10] analytically described the nonlinear stochastic dynamics of a class of two-state bioreactors with isotonic or nonisotonic kinetics. Ruan and Wolkowicz [11] have studied a chemostat model of a single species with a distributed delay. They have demonstrated that a Hopf bifurcation takes place when the average time delay crosses a crucial amount, using this parameter as a bifurcation parameter. Yuan and Han [12] showed that a chemostat-type model with two distributed delays undergoes Hopf bifurcations using the average time delay as a bifurcation parameter. Xu et al. [13] analyzed a food chain chemostat model with a Michaelis-Menten function response and double delays. They have examined the necessary conditions for the Hopf bifurcation of the system at the positive equilibrium using the stability theory of functional differential equations to investigate the conditions for the stability of three equilibria. One can refer previous literature [14–20] and references cited therein for more studies on the chemostat models.

Despite the abundance of studies regarding the topic, there are a few reports about discrete systems in the chemostat up to now, at least to our knowledge. Some of them are as follows. Gage et al. [21] considered a discrete, environmentally coupled, size-specific model of microbial population dynamics in continuous culture. Smith et al. [22,23] proved competitive exclusion for a discrete-time, size-structured, nonlinear matrix model in the chemostat. Arino et al. [24,25] analyzed inhomogeneous, substrate dependent cell division in a time discrete, nonlinear matrix model of size-structured population growth in the chemostat. Zhang et al. [26] demonstrated that inhibitory kinetics can create extraordinarily complex dynamics, including stable equilibria, cycles, and chaos, in contrast to the model with a monotonically increasing uptake function. Jang [27] proposed and analyzed a discrete-time, discrete size-structured model of plasmid-bearing and plasmid-free competition in the chemostat. Amster et al. [28] extended to the nonautonomous framework the results about the dynamics of a discrete and nonlinear matrix model describing the growth of a single microbial population of size structured in an autonomous chemostat, provided a threshold that determines either extinction or persistence of total biomass and establishes a set of sufficient conditions to ensure the existence, uniqueness, and global attractiveness of a ω periodic solution.

In contrast, only a limited number of studies have investigated bifurcation analysis in discrete chemostat models. Al-Basyouni and Khan [29] investigated equilibrium dynamics in a discrete chemostat model using discrete dynamical system theory and bifurcation analysis. Alsulami [30] examined a discrete-time chemostat model using the piecewise constant argument method, identifying fixed points, transcritical and period-doubling bifurcations, and employing feedback and hybrid control techniques to regulate chaos. In addition, numerical simulations support the theoretical findings, demonstrating that the piecewise constant argument method provides better dynamic consistency than the forward Euler method. In [31], we consider the discrete-time simple chemostat model, derived from a continuous-time population model using the forward Euler method, analyzed for its dynamic behavior. Fixed points are identified, and a local stability analysis is performed. Bifurcation theory reveals a flip bifurcation, with the integral step size acting as the bifurcation parameter. A hybrid control approach stabilizes the bifurcation. Numerical simulations confirm the theoretical findings, showcasing novel and interesting dynamics.

The modified chemostat model in [32] can be described as:

(1) d S d t = D ( S F − S ) − ( a + e − b S ) μ max S x S + K S d x d t = μ max S x S + K S − D x ,

where S and x denote the concentrations of substrate and microorganism in a chemostat at time t , respectively. D is the dilution rate, and S F > 0 is the concentration of the feed substrate. The parameter b represents the sensitivity of the cell to the substrate when the cell is in a good environment, while a > 1 and b > 0 are biological constraints. μ max is the maximal growth rate of the microorganism, and K S is the saturation constant.

The interactions between prey, medium predators, and top predators in populations with overlapping generations are usually described using ordinary differential equations (ODEs) since reproduction occurs continuously in these populations [33]. In contrast, discrete-time models represented by difference equations are often more appropriate and realistic when populations have non-overlapping generations [34]. Furthermore, a discretization technique for continuous models is required since, in general, nonlinear continuous systems do not have analytical solutions that can be expressed in terms of a finite representation of elementary functions. Controlling chaos and bifurcations involves designing a controller to modify chaotic behaviors and bifurcation features of a nonlinear system in order to achieve specific desired dynamical behaviors [35].

The motivation for this study arises from the need to analyze chemostat models in a discrete-time framework. Although ODEs are commonly used to describe population dynamics in continuously reproducing populations, discrete-time models based on difference equations are more suitable when generations do not overlap. Additionally, nonlinear continuous systems often lack closed-form analytical solutions, necessitating discretization techniques.

Building on the modified chemostat model from [32], this study employs the forward Euler scheme to derive a discrete-time version of the system. The Euler method approximates derivatives using finite differences, allowing the continuous chemostat model to be transformed into a discrete system. This discrete approach facilitates the analysis of bifurcations and chaotic behavior, which can then be controlled using specialized techniques.

Applying this scheme to the continuous chemostat model (1) yields the following discrete-time system:

(2) S t + 1 = S t + h D ( S F − S t ) − ( a + e − b S t ) μ max S t x t S t + K S x t + 1 = x t + h μ max S t S t + K S − D x t ,

where h > 0 is the step size.

2 Preliminaries

In this section, we present the necessary definitions, lemmas, and theorems required for our study. These concepts are essential for understanding the dynamics of the modeled system and facilitating analytical investigations.

Lemma 2.1

[36] Let F ( λ ) = λ 2 + B λ + C , where B and C are constants. Suppose F ( 1 ) > 0 and λ 1 and λ 2 are two roots of F ( λ ) = 0 . Then,

  1. ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ < 1 if and only if F ( − 1 ) > 0 and C < 1 ,

  2. ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ > 1 if and only if F ( − 1 ) > 0 and C > 1 ,

  3. ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ < 1 if and only if F ( − 1 ) < 0 ,

  4. λ 1 = − 1 and ∣ λ 2 ∣ ≠ 1 if and only if F ( − 1 ) = 0 , B ≠ 0 , 2 ,

  5. λ 1 and λ 2 are a pair of conjugate complex roots and ∣ λ 1 ∣ = ∣ λ 2 ∣ = 1 if and only if B 2 − 4 C < 0 and C = 1 .

Definition 2.2

[37] A fixed point ( x , y ) is called

  1. sink if ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ < 1 , and it is locally asymptotically stable,

  2. source if ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ > 1 , and it is locally unstable,

  3. saddle if ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ < 1 , or ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ > 1 ,

  4. non-hyperbolic if either ∣ λ 1 ∣ = 1 or ∣ λ 2 ∣ = 1 .

Theorem 2.3

[38] (Existence) There exists a C r center manifold for the following system:

(3) x n + 1 = A x n + f ( x n , y n ) y n + 1 = B y n + g ( x n , y n ) , ( x n , y n ) ∈ R c × R s ,

which can be locally represented as a graph as follows:

W c ( 0 ) = { ( x n , y n ) ∈ R c × R s : y n = h ( x n ) , ∣ x n ∣ < δ , h ( 0 ) = 0 , D h ( 0 ) = 0 }

for δ sufficiently small. Moreover, the dynamics of (3)  restricted to the center manifold is, for u n sufficiently small, given by the following c -dimensional map:

(4) u n + 1 = A u n + f ( u n , h ( u n ) ) , u n ∈ R c ,

where f ( 0 , 0 ) = 0 , D f ( 0 , 0 ) = 0 , g ( 0 , 0 ) = 0 , D g ( 0 , 0 ) = 0 , and f : R c × R s → R c , g : R c × R s → R s are C r functions ( r ≥ 2 ) in some neighborhood of the origin, A is a c × c matrix with eigenvalues of modulus one, and B is an s × s matrix with eigenvalues of modulus less than one.

Definition 2.4

[39] The bifurcation associated with the appearance of μ 1 = − 1 is called a flip (or period-doubling) bifurcation, where μ 1 is an eigenvalue of the non-hyperbolic fixed point of the system.

(5) x n + 1 = f ( x n , α ) , x n ∈ R n , α ∈ R ,

where the map f : R n × R → R n is smooth.

Theorem 2.5

[39] Suppose that a one-dimensional system

x n + 1 = f ( x n , α ) , x n ∈ R 1 , α ∈ R 1 ,

with smooth f, has at α = 0 the fixed point x 0 = 0 , and let μ = f x ( 0 , 0 ) = − 1 . Assume that the following nondegeneracy conditions are satisfied:

  1. 1 2 ( f x x ( 0 , 0 ) ) 2 + 1 3 f x x x ( 0 , 0 ) ≠ 0 ,

  2. f x α ( 0 , 0 ) ≠ 0 .

Then, there is a flip bifurcation in the system.

Definition 2.6

[39] The bifurcation corresponding to the presence of μ 1,2 = e ± i θ 0 , 0 < θ 0 < π , is called a Neimark-Sacker bifurcation, where μ 1,2 are the eigenvalues of the non-hyperbolic fixed point of system (5).

Theorem 2.7

[39,40] Suppose a two-dimensional discrete-time system

x n + 1 = f ( x n , α ) , x n ∈ R 2 , α ∈ R 1

has a pair of complex conjugate eigenvalues μ 1,2 = r ( α ) e ± i θ ( α ) at the fixed point ( 0 , 0 ) , where α c is the critical parameter value, r ( α c ) = 1 and θ ( α c ) = θ 0 with smooth f. Let the following conditions be satisfied:

  1. r ′ ( α c ) ≠ 0 ,

  2. e i k θ 0 ≠ 1 for k = 1 , 2, 3, 4 and a ( α c ) = Re ( e i θ 0 c 1 ( α c ) ) ≠ 0 .

Then, there is a neighborhood of origin in which a unique closed curve bifurcates from the origin as α passes through the critical value α c . The sign of a ( α c ) determines the direction of the appearance of the invariant curve in a generic system exhibiting the Neimark-Sacker bifurcation. If a ( α c ) < 0 , then there is a supercritical Neimark-Sacker bifurcation. If a ( α c ) > 0 , then there is a subcritical Neimark-Sacker bifurcation.

3 Existence and analysis of fixed points

This section examines the existence and local stability of all potential fixed points. Finding the fixed points, represented by E i = ( S , x ) , i = 0 , 1 , and their corresponding stability properties, allows for the analysis of the system (2). Setting the right-hand sides of the first equation and second equation in system (2) equal to S and x , respectively, yields the fixed points E 0 = ( S F , 0 ) and

E 1 = D K S μ max − D , S F ( μ max − D ) − D K S ( μ max − D ) a + e − b D K S μ max − D .

Since all biologically meaningful variables are nonnegative, necessary conditions are derived for the existence of each equilibrium point. Given that the input concentration of the nutrient is positive, E 0 is always true. E 1 exists when S F ( μ max − D ) > D K S . For each equilibrium point, similar biological interpretations of the stability conditions can be made. The Jacobian of system (2) at each equilibrium point and its corresponding eigenvalues are used to examine the stability of fixed points [41].

The forward Euler method is conditionally stable, and its stability depends on the size of the time step in relation to the model’s dynamics. A discussion of numerical stability constraints should be included to ensure that the chosen step size does not lead to instability or divergence in the numerical solution. To ensure numerical stability, we analyze the eigenvalues of the Jacobian matrix. For the continuous system, stability is governed by the real parts of the eigenvalues of the continuous system’s Jacobian, J c . However, for the Forward Euler discretization, stability requires the eigenvalues of the discrete system’s Jacobian, J d = I + h J c , where I is the identity matrix, to satisfy ∣ λ d ∣ ≤ 1 . Here, λ d refers to the eigenvalues of the Jacobian matrix of the discrete system. This imposes the constraint h < 2 ∣ λ max ∣ , where λ max is the largest eigenvalue of J c . Choosing a step size larger than this limit leads to numerical instability [44].

We obtain the main results on the dynamic behavior of every fixed point in system (2) as follows.

Theorem 3.1

Assume that D ( S F + K S ) − S F μ max > 0 . E 0 = ( S F , 0 ) is

  1. sink if 0 < h < 2 D ,

  2. source if h > 2 ( S F + K S ) D ( S F + K S ) − S F μ max ,

  3. saddle if 2 D < h < 2 ( S F + K S ) D ( S F + K S ) − S F μ max .

Proof

The Jacobian matrix at fixed point E 0 = ( S F , 0 ) is

(6) J ( E 0 ) = 1 − D h − S F h μ max S F + K S ( a + e − b S F ) 0 1 − D h + S F h μ max S F + K S ,

which has the characteristic polynomial of (6)

F ( λ ) = λ 2 + 2 ( D h − 1 ) ( S F + K S ) − S F h μ max S F + K S λ + ( D h − 1 ) ( ( D h − 1 ) ( S F + K S ) − S F h μ max ) S F + K S .

It is know that D ( S F + K S ) − S F μ max > 0 , so

F ( 1 ) = D h 2 ( D ( S F + K S ) − S F μ max ) S F + K S > 0 .

On the other hand,

F ( − 1 ) = ( D h − 2 ) ( ( D h − 2 ) ( S F + K S ) − S F h μ max ) S F + K S

is expressed as a function of δ .

δ ( h ) = d 2 ( S F + K S ) − d S F μ max S F + K S h 2 + 2 S F μ max − 4 D ( S F + K S ) S F + K S h + 4 ,

and the roots of δ ( h ) = 0 ,

h ¯ 1 = 2 D , h ¯ 2 = 2 ( S F + K S ) D ( S F + K S ) − S F μ max .

It is noted that h ¯ 2 > h ¯ 1 > 0 , because of D ( S F + K S ) − S F μ max > 0 . Hence, it can be found that F ( − 1 ) > 0 if h < h ¯ 1 or h > h ¯ 2 . If h < h ¯ 1 < h ¯ 2 ,

C = ( D h − 1 ) ( ( D h − 1 ) ( S F + K S ) − S F h μ max ) S F + K S < 1 .

If h > h ¯ 2 > h ¯ 1 , C > 1 . On the other hand, if h ¯ 1 < h < h ¯ 2 , it is obtained that F ( − 1 ) < 0 . □

Theorem 3.2

Assume that S F ( μ max − D ) > D K S and

( − D ( S F ( μ max − D ) − D K S ) + S F μ max ( μ max − D ) ) ( μ max − D ) a e b D K S μ max − D + 1 > b D K S μ max ( S F ( μ max − D ) − D K S ) .

The positive equilibrium point

E 1 = D K S μ max − D , S F ( μ max − D ) − D K S ( μ max − D ) a + e − b D K S μ max − D

is

  1. sink fixed point if one of the following conditions holds:

    0 < h < h ¯ 1 when α 2 > 4 α 1 , 0 < h < 1 α 1 when α 2 ≤ 4 α 1 .

  2. source fixed point if one of the following conditions holds:

    h > h ¯ 2 when α 2 > 4 α 1 , h > 1 α 1 when α 2 ≤ 4 α 1 .

  3. saddle fixed point if the following condition holds:

    h ¯ 1 < h < h ¯ 2 when α 2 > 4 α 1 ,

where

(7) h ¯ 1 = α 2 − α 2 2 − 4 α 1 α 2 α 1 α 2 , h ¯ 2 = α 2 + α 2 2 − 4 α 1 α 2 α 1 α 2

and

(8) α 1 = ( D − μ max ) 2 ( D − μ max ) ( D 2 ( S F + K S ) + S F μ max ( μ max − 2 D ) ) D ( D ( S F + K S ) − S F μ max ) − b K S μ max a e b D K S μ max − D + 1 , α 2 = D 2 ( S F + K S ) + S F μ max ( − 2 D + μ max ) K S μ max − b D ( D ( S F + K S ) − S F μ max ) ( D − μ max ) a e b D K S μ max − D + 1 .

Proof

The Jacobian matrix at fixed point

E 1 = D K S μ max − D , S F ( μ max − D ) − D K S ( μ max − D ) a + e − b D K S μ max − D

is

(9) J ( E 1 ) = J 11 J 12 J 21 1 ,

where

(10) J 11 = b D h ( D ( S F + K S ) − S F μ max ) ( D − μ max ) a e b D K S μ max − D + 1 + μ max ( 2 D S F h − S F h μ max + K S ) − D 2 h ( S F + K S ) K S μ max , J 12 = − D h a + e b D K S D − μ max , J 21 = h ( D − μ max ) ( D ( S F + K S ) − S F μ max ) K S μ max a + e b D K S D − μ max .

Hence, the characteristic polynomial of (9) is

(11) F ( λ ) = λ 2 + ( α 2 h − 2 ) λ + α 1 α 2 h 2 − α 2 h + 1 ,

where α 1 and α 2 are in (8). So, we obtain

F ( 1 ) = D h 2 ( D − μ max ) ( D ( S F + K S ) − S F μ max ) K S μ max > 0

because of S F ( μ max − D ) > D K S . On the other hand, it is obtained that

F ( − 1 ) = α 1 α 2 h 2 − 2 α 2 h + 4 .

Denote

(12) γ ( h ) = α 1 α 2 h 2 − 2 α 2 h + 4 .

The function γ ( h ) has two real roots,

h ¯ 1 = α 2 − α 2 2 − 4 α 1 α 2 α 1 α 2 , h ¯ 2 = α 2 + α 2 2 − 4 α 1 α 2 α 1 α 2 .

Three cases can be considered here. Since α 2 > 0 , for the first case, if we assume that α 2 > 4 α 1 , then h ¯ 1 and h ¯ 2 are real and distinct. Hence, we obtain h ¯ 2 > h ¯ 1 > 0 . It is obtained that F ( − 1 ) > 0 when 0 < h < h ¯ 1 or h > h ¯ 2 . If 0 < h < h ¯ 1 < 1 α 1 < h ¯ 2 , it is found that

C = h 2 α 1 α 2 − h α 2 + 1 < 1 .

If h > h ¯ 2 > h ¯ 1 > 1 α 1 , then C > 1 . Finally, if h ¯ 1 < h < h ¯ 2 , it is obtained that F ( − 1 ) < 0 . Second, because of α 2 > 0 , if α 2 = 4 α 1 and third if α 2 < 4 α 1 , then F ( − 1 ) has real roots h ¯ 1 = h ¯ 2 = 1 α 1 and does not have real roots, respectively and for each case, F ( − 1 ) > 0 . Also, If 0 < h < 1 α 1 , then C < 1 or if h > 1 α 1 , then C > 1 . So, from Definition 2.2 and Lemma 2.1, the proof is completed.□

4 Bifurcation analysis

4.1 Flip bifurcation analysis

Lemma 4.1

When α 2 > 4 α 1 , if F ( − 1 ) has the roots h = h ¯ 1 or h = h ¯ 2 and

h ≠ 2 α 2 , 4 α 2 ,

then a flip bifurcation may occur at fixed point E 1 , where h ¯ 1 , h ¯ 2 , and α 2 are in (7) and (8).

Proof

From Lemma 2.1 (4), we know that if h = h ¯ 1 , where h ¯ 1 is the root of F ( − 1 ) = 0 , then the eigenvalues of F ( λ ) in (11) are obtained as

λ 1 = − 1 , λ 2 = 3 − h ¯ 1 α 2 ,

with ∣ λ 2 ∣ ≠ 1 . Therefore, from Definition 2.4, a flip bifurcation may occur when α 2 > 4 α 1 .□

Now, we analyze the existence of flip bifurcation at positive fixed point

E 1 = D K S μ max − D , S F ( μ max − D ) − D K S ( μ max − D ) a + e − b D K S μ max − D ,

utilizing the central manifold and bifurcation theory [38,39]. First, we transform the fixed point E 1 = ( S * , x * ) to the origin and expand the Taylor series to the third order. Thus, system (2) changes into

(13) S t + 1 − S * x t + 1 − x * = J 11 J 12 J 21 1 S t − S * x t − x * + f 1 f 2 ,

where

(14) J 11 = b D h ( D ( S F + K S ) − S F μ max ) ( D − μ max ) a e b D K S μ max − D + 1 + μ max ( 2 D S F h − S F h μ max + K S ) − D 2 h ( S F + K S ) K S μ max , J 12 = − D h a + e b D K S D − μ max , J 21 = h ( D − μ max ) ( D ( S F + K S ) − S F μ max ) K S μ max a + e b D K S D − μ max ,

and

f 1 = b K S μ max ( − D ( b K S + 4 ) μ max + 2 D 2 + 2 μ max 2 ) − 2 ( D − μ max ) 3 a e b D K S μ max − D + 1 2 K S 2 μ max 2 ( D − μ max ) a e b D K S μ max − D + 1 × h ( D ( S F + K S ) − S F μ max ) ( S t − S * ) 2 − a h ( D − μ max ) 2 + h e b D K S D − μ max ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) K S μ max × ( S t − S * ) ( x t − x * ) + h ( D ( S F + K S ) − S F μ max ) 6 K S 3 μ max 3 ( D − μ max ) a e b D K S μ max − D + 1 × b K S μ max ( − 3 D 2 ( b K S + 6 ) μ max + D ( b K S ( b K S + 6 ) + 18 ) μ max 2 − 3 ( b K S + 2 ) μ max 3 + 6 D 3 ) − 6 ( D − μ max ) 4 a e b D K S μ max − D + 1 ( S t − S * ) 3 − 1 2 K S 2 μ max 2 2 a h ( D − μ max ) 3 + h e b D K S D − μ max ( − 2 D 2 ( b K S + 3 ) μ max + D ( b K S ( b K S + 4 ) + 6 ) μ max 2 − 2 ( b K S + 1 ) μ max 3 + 2 D 3 ) ( S t − S * ) 2 ( x t − x * ) + O ( ∣ S t ∣ + ∣ x t ∣ + ∣ h ∣ ) 4 ,

f 2 = h ( D − μ max ) 2 ( D ( S F + K S ) − S F μ max ) K S 2 μ max 2 a + e b D K S D − μ max ( S t − S * ) 2 + h ( D − μ max ) 2 K S μ max ( S t − S * ) ( x t − x * ) + h ( D − μ max ) 3 ( D ( S F + K S ) − S F μ max ) K S 3 μ max 3 a + e b D K S D − μ max ( S t − S * ) 3 + h ( D − μ max ) 3 K S 2 μ max 2 ( S t − S * ) 2 ( x t − x * ) + O ( ∣ S t ∣ + ∣ x t ∣ + ∣ h ∣ ) 4 .

Let h be the bifurcation parameter. Now, we introduce new variables as follows:

(15) u t = S t − S * , v t = x t − x * , h ˜ = h − h ¯ .

Hence,

(16) u t + 1 = S t + 1 − S * , v t + 1 = x t + 1 − x * , h = h ¯ + h ˜ .

If equalities (15) and (16) are substituted in (13), we obtain

(17) u t + 1 v t + 1 = J 11 J 12 J 21 1 u t v t + f 1 ( u t , v t , h ¯ , h ˜ ) f 2 ( u t , v t , h ¯ , h ˜ ) ,

where J 11 , J 12 , J 21 are in (14) for h = h ¯ + h ˜ and

(18) f 1 ( u t , v t , h ¯ , h ˜ ) = b D h ˜ ( D ( S F + K S ) − S F μ max ) ( D − μ max ) a e b D K S μ max − D + 1 + μ max ( 2 D S F h ˜ − S F h ˜ μ max + K S ) − D 2 h ˜ ( S F + K S ) K S μ max u t − D h ˜ a + e b D k D − μ max v t + ( h ¯ + h ˜ ) ( D ( S F + K S ) − S F μ max ) 2 K S 2 μ max 2 ( D − μ max ) a e b D K S μ max − D + 1 × b K S μ max ( − D ( b K S + 4 ) μ max + 2 D 2 + 2 μ max 2 ) − 2 ( D − μ max ) 3 a e b D K S μ max − D + 1 u t 2 − a ( h ¯ + h ˜ ) ( D − μ max ) 2 + ( h ¯ + h ˜ ) e b D K S D − μ max ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) 1 K S μ max u t v t + b K S μ max ( − 3 D 2 ( b K S + 6 ) μ max + D ( b K S ( b K S + 6 ) + 18 ) μ max 2 − 3 ( b K S + 2 ) μ max 3 + 6 D 3 ) − 6 ( D − μ max ) 4 a e b D K S μ max − D + 1 ( h ¯ + h ˜ ) ( D ( S F + K S ) − S F μ max ) 6 K S 3 μ max 3 ( D − μ max ) a e b D K S μ max − D + 1 u t 3 − ( 2 a ( h ¯ + h ˜ ) ( D − μ max ) 3 + ( h ¯ + h ˜ ) e b D K S D − μ max ( − 2 D 2 ( b K S + 3 ) μ max + D ( b K S ( b K S + 4 ) + 6 ) μ max 2 − 2 ( b K S + 1 ) μ max 3 + 2 D 3 ) 1 2 K S 2 μ max 2 u t 2 v t + O ( ∣ u t ∣ + ∣ v t ∣ + ∣ h ¯ ∣ + ∣ h ˜ ∣ ) 4 ,

f 2 ( u t , v t , h ¯ , h ˜ ) = ( b D h ˜ ( D ( S F + K S ) − S F μ max ) ) ( D − μ max ) ( a e b d k μ max − d + 1 ) u t + v t + ( h ¯ + h ˜ ) ( D − μ max ) 2 ( D ( S F + K S ) − S F μ max ) K S 2 μ max 2 a + e b D K S D − μ max u t 2 + ( h ¯ + h ˜ ) ( D − μ max ) 2 K S μ max u t v t + ( h ¯ + h ˜ ) ( D − μ max ) 3 ( D ( S F + K S ) − S F μ max ) K S 3 μ max 3 a + e b D K S D − μ max u t 3 + ( h ¯ + h ˜ ) ( D − μ max ) 3 K S 2 μ max 2 u t 2 v t + O ( ∣ u t ∣ + ∣ v t ∣ + ∣ h ¯ ∣ + ∣ h ˜ ∣ ) 4 .

The equilibrium point of (17) is (0, 0). Therefore, the equilibrium point

E 1 = D K S μ max − D , S F ( μ max − D ) − D K S ( μ max − D ) a + e − b D K S μ max − D

of (2) is transformed to E 1 = ( 0 , 0 ) of (17).

We construct the matrix shown below, where the columns of matrix T consist of the eigenvectors corresponding to the eigenvalues λ 1 = − 1 , λ 2 = 3 − α 2 h ¯ 1 of the Jacobian matrix

J ¯ ( S * , x * ) = J ¯ 11 J ¯ 12 J ¯ 21 1 ,

when h = h ¯ 1 , where

(19) J ¯ 11 = b D h ¯ ( D ( S F + K S ) − S F μ max ) ( D − μ max ) a e b D K S μ max − D + 1 + μ max ( 2 D S F h ¯ − S F h ¯ μ max + K S ) − D 2 h ¯ ( S F + K S ) K S μ max , J ¯ 12 = − D h ¯ a + e b D K S D − μ max , J ¯ 21 = h ¯ ( D − μ max ) ( D ( S F + K S ) − S F μ max ) K S μ max a + e b D K S D − μ max .

We find the T matrix

(20) T = − 2 2 − h ¯ α 2 J ¯ 21 J ¯ 21 .

Next we apply the following transformation:

(21) u t v t = T y t z t .

Hence, (17) can be expressed as

(22) y t + 1 z t + 1 = λ 1 0 0 λ 2 y t z t + g 1 ( y t , z t ) g 2 ( y t , z t ) ,

where

g 1 ( y t , z t ) = ν 1 y t 2 + ν 2 y t z t + ν 3 h ˜ y t + ν 4 h ˜ z t + ν 5 y t 3 + ν 6 h ˜ y t 2 + O ( ∣ y t ∣ + ∣ z t ∣ + ∣ h ˜ ∣ ) 4 , g 2 ( y t , z t ) = ω 1 y t 2 + ω 2 y t z t + ω 3 h ˜ y t + ω 4 h ˜ z t + ω 5 y t 3 + ω 6 h ˜ y t 2 + O ( ∣ y t ∣ + ∣ z t ∣ + ∣ h ˜ ∣ ) 4 ,

and

ν 1 = − 1 K S 2 μ max 2 ( μ max − D ) a + e b D K S D − μ max 2 h ¯ e b D K S D − μ max ( μ max ( D 3 ( 2 J 11 K S ( a J ¯ 21 + 2 b K S + 6 ) + 4 S F J 11 ( b K S + 4 ) − J 12 J ¯ 21 K S ) + μ max ( D 2 ( − J 11 K S ( a J ¯ 21 ( b K S + 6 ) + 2 b K S ( b K S + 4 ) + 12 ) − 2 S F J 11 ( b K S ( b K S + 6 ) + 12 ) + 3 J 12 J ¯ 21 K S ) + D μ max ( J 11 K S ( a J ¯ 21 ( b K S + 6 ) + 4 b K S + 4 ) + 2 S F J 11 ( b K S + 2 ) ( b K S + 4 ) − 3 J 12 J ¯ 21 K S ) + μ max 2 ( J ¯ 21 K S ( J 12 − 2 a J 11 ) − 4 S F ( b J 11 K S + J 11 ) ) ) ) − 4 D 4 J 11 ( S F + K S ) ) + ( a J 11 − J 12 ) ( D − μ max ) 3 ( μ max ( a J ¯ 21 K S + 4 S F ) − 4 D ( S F + K S ) ) + J 11 J ¯ 21 K S μ max ( D − μ max ) e 2 b D K S D − μ max ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) ,

ν 2 = h ¯ 2 K S 3 μ max 3 ( D − μ max ) 2 a + e b D K S D − μ max 2 8 K S μ max ( D − μ max ) ( D ( S F + K S ) − S F μ max ) × a + e b D K S D − μ max ( 2 ( a J 11 − J 12 ) ( D − μ max ) 3 + J 11 e b D K S D − μ max ( − 2 D 2 ( b K S + 3 ) μ max + D ( b K S ( b K S + 4 ) + 6 ) μ max 2 − 2 ( b K S + 1 ) μ max 3 + 2 D 3 ) ) + h ¯ ( a ( D − μ max ) ( D 2 ( S F + K S ) + S F μ max ( μ max − 2 D ) ) + e b D K S D − μ max ( μ max ( S F μ max ( D ( b K S + 3 ) − μ max ) − D 2 ( b S F K S + b K S 2 + 3 S F + K S ) ) + D 3 ( S F + K S ) ) ) × e b D K S D − μ max ( μ max ( D 3 ( 2 J 11 K S ( a J ¯ 21 + 4 b K S + 12 ) + 8 S F J 11 ( b K S + 4 ) − J 12 J ¯ 21 K S ) + μ max ( D 2 ( − J 11 K S ( a J ¯ 21 ( b K S + 6 ) + 4 b K S ( b K S + 4 ) + 24 ) − 4 S F J 11 ( b K S ( b K S + 6 ) + 12 ) + 3 J 12 J ¯ 21 K S ) + D μ max ( J 11 K S ( a J ¯ 21 ( b K S + 6 ) + 8 b K S + 8 ) + 4 S F J 11 ( b K S + 2 ) ( b K S + 4 ) − 3 J 12 J ¯ 21 K S ) + μ max 2 ( J ¯ 21 K S ( J 12 − 2 a J 11 ) − 8 S F ( b J 11 K S + J 11 ) ) ) ) − 8 D 4 J 11 ( S F + K S ) ) + ( a J 11 − J 12 ) ( D − μ max ) 3 ( μ max ( a J ¯ 21 K S + 8 S F ) − 8 D ( S F + K S ) ) + J 11 J ¯ 21 K S μ max ( D − μ max ) e 2 b D K S D − μ max ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) ,

ν 3 = 1 K S μ max ( D − μ max ) a + e b D K S D − μ max 2 e b D K S D − μ max ( − D 2 μ max ( K S ( a J 11 J ¯ 21 + b J 11 K S + b J 12 K S + J 11 ) + b S F K S ( J 11 + J 12 ) + 3 S F J 11 ) + D μ max 2 ( a J 11 J ¯ 21 K S + b S F K S ( J 11 + J 12 ) + 3 S F J 11 ) + D 3 J 11 ( S F + K S ) − S F J 11 μ max 3 ) + a J 11 ( D − μ max ) ( − D μ max ( a J ¯ 21 K S + 4 S F ) + 2 D 2 ( S F + K S ) + 2 S F μ max 2 ) − D J 11 J ¯ 21 K S μ max ( D − μ max ) e 2 b D K S D − μ max ,

ν 4 = 1 K S 2 μ max 2 ( D − μ max ) 2 a + e b D K S D − μ max 2 h ¯ ( a ( D − μ max ) ( D 2 ( S F + K S ) + S F μ max ( μ max − 2 D ) ) + e b D K S D − μ max ( μ max ( S F μ max ( D ( b K S + 3 ) − μ max ) − D 2 ( b S F K S + b K S 2 + 3 S F + K S ) ) + D 3 ( S F + K S ) ) × a J 11 ( D − μ max ) ( D 2 ( S F + K S ) + S F μ max ( μ max − 2 D ) ) + e b D K S D − μ max ( μ max ( S F μ max ( b D K S ( J 11 + J 12 ) + 3 D J 11 − J 11 μ max ) − D 2 ( b S F K S ( J 11 + J 12 ) + K S ( b K S ( J 11 + J 12 ) + J 11 ) + 3 S F J 11 ) ) + D 3 J 11 ( S F + K S ) ) ) − K S μ max ( D − μ max ) a + e b D K S D − μ max 2 e b D K S D − μ max ( − D 2 μ max ( K S ( − a J 11 J ¯ 21 + b J 11 K S + b J 12 K S + J 11 ) + b S F K S ( J 11 + J 12 ) + 3 S F J 11 ) + D μ max 2 ( − a J 11 J ¯ 21 K S + b S F K S ( J 11 + J 12 ) + 3 S F J 11 ) + D 3 J 11 ( S F + K S ) − S F J 11 μ max 3 ) + a J 11 ( D − μ max ) ( μ max ( a D J ¯ 21 K S − 4 D S F + 2 S F μ max ) + 2 D 2 ( S F + K S ) ) + D J 11 J ¯ 21 K S μ max ( D − μ max ) e 2 b D K S D − μ max ,

ν 5 = 1 3 K S 3 μ max 3 ( μ max − D ) a + e b D K S D − μ max 2 h ¯ μ max 2 D 4 3 a 2 J 11 J ¯ 21 K S + e b D K S D − μ max ( 6 J 11 K S ( a J ¯ 21 + b K S + 4 ) + 6 S F J 11 ( b K S + 5 ) − J 12 J ¯ 21 K S ) + 30 a S F J 11 + 24 a J 11 K S − a J 12 J ¯ 21 K S + 3 J 11 J ¯ 21 K S e 2 b D K S D − μ max + 2 S F J 12 K S + 2 J 12 K S 2 + μ max ( μ max ( 2 μ max 2 ( 3 a 2 J 11 J ¯ 21 K S + 6 a S F J 11 − a J 12 J ¯ 21 K S + 2 S F J 12 K S ) + e b D K S D − μ max ( D 2 ( K S ( J 11 ( 3 a J ¯ 21 ( b K S ( b K S + 6 ) + 24 ) + 2 b K S ( b K S ( b K S + 6 ) + 18 ) + 48 ) − 12 J 12 J ¯ 21 ) + 2 S F J 11 ( b K S ( b K S ( b K S + 9 ) + 36 ) + 60 ) ) + μ max ( D ( K S ( 8 J 12 J ¯ 21 − 3 J 11 ( a J ¯ 21 ( b K S ( b K S + 6 ) + 16 ) + 2 b K S ( b K S + 2 ) + 4 ) ) − 2 S F J 11 ( b K S ( b K S ( b K S + 9 ) + 24 ) + 30 ) ) + 2 μ max ( J ¯ 21 K S ( 3 a J 11 ( b K S + 2 ) − J 12 ) + 3 S F J 11 ( b K S ( b K S + 2 ) + 2 ) ) ) ) + 12 D 2 ( K S ( a J 11 ( 3 a J ¯ 21 + 4 ) − a J 12 J ¯ 21 + 2 S F J 12 ) + 10 a S F J 11 + J 12 K S 2 ) − 4 D μ max ( K S ( a J 11 ( 6 a J ¯ 21 + 3 ) − 2 a J 12 J ¯ 21 + 4 S F J 12 ) + 15 a S F J 11 + J 12 K S 2 ) + 3 J 11 J ¯ 21 K S e 2 b D K S D − μ max ( D 2 ( b K S ( b K S + 6 ) + 12 ) − D ( b K S + 2 ) ( b K S + 4 ) μ max + 2 ( b K S + 1 ) μ max 2 ) ) + 2 D 3 ( − e b D K S D − μ max ( K S ( 3 J 11 ( a J ¯ 21 ( b K S + 8 ) + b K S ( b K S + 6 ) + 12 ) − 4 J 12 J ¯ 21 ) + 3 S F J 11 ( b K S ( b K S + 8 ) + 20 ) ) − 4 K S ( a J ¯ 21 ( 3 a J 11 − J 12 ) + 9 a J 11 + 2 S F J 12 ) − 60 a S F J 11 − 3 J 11 J ¯ 21 K S ( b K S + 4 ) e 2 b D K S D − μ max − 6 J 12 K S 2 ) ) − 12 D 5 J 11 ( S F + K S ) a + e b D K S D − μ max ,

(23) ν 6 = e b D K S D − μ max K S 2 μ max 2 ( μ max − D ) a + e b D K S D − μ max ( − D 3 μ max ( 2 J 11 K S ( a J ¯ 21 + 2 b K S + 6 ) + 4 S F J 11 ( b K S + 4 ) − J 12 J ¯ 21 K S ) + D 2 μ max 2 ( K S ( J 11 ( a J ¯ 21 ( b K S + 6 ) + 2 b K S ( b K S + 4 ) + 12 ) − 3 J 12 J ¯ 21 ) + 2 S F J 11 ( b K S ( b K S + 6 ) + 12 ) ) − D μ max 3 ( J 11 K S ( a J ¯ 21 ( b K S + 6 ) + 4 b K S + 4 ) + 2 S F J 11 ( b K S + 2 ) ( b K S + 4 ) − 3 J 12 J ¯ 21 K S ) + μ max 4 ( J ¯ 21 K S ( 2 a J 11 − J 12 ) + 4 S F ( b J 11 K S + J 11 ) ) + 4 D 4 J 11 ( S F + K S ) ) + ( a J 11 − J 12 ) ( D − μ max ) 3 ( 4 D ( S F + K S ) − μ max ( a J ¯ 21 K S + 4 S F ) ) − J 11 J ¯ 21 K S μ max ( D − μ max ) × e 2 b D K S D − μ max ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) ,

ω 1 = − 1 K S 2 μ max 2 ( μ max − D ) a + e b D K S D − μ max h ¯ × ( e b D K S D − μ max ( μ max ( D 3 ( 2 J 21 K S ( a J ¯ 21 + 2 b K S + 6 ) + 4 S F J 21 ( b K S + 4 ) − J ¯ 21 J 22 K S ) + μ max ( D 2 ( − J 21 K S ( a J ¯ 21 ( b K S + 6 ) + 2 b K S ( b K S + 4 ) + 12 ) − 2 S F J 21 ( b K S ( b K S + 6 ) + 12 ) + 3 J ¯ 21 J 22 K S ) + D μ max ( J 21 K S ( a J ¯ 21 ( b K S + 6 ) + 4 b K S + 4 ) + 2 S F J 21 ( b K S + 2 ) ( b K S + 4 ) − 3 J ¯ 21 J 22 K S ) + μ max 2 ( J ¯ 21 K S ( J 22 − 2 a J 21 ) − 4 S F ( b J 21 K S + J 21 ) ) ) ) − 4 D 4 J 21 ( S F + K S ) ) + ( a J 21 − J 22 ) ( D − μ max ) 3 ( μ max ( a J ¯ 21 K S + 4 S F ) − 4 D ( S F + K S ) ) + J ¯ 21 J 21 K S μ max ( D − μ max ) e 2 b D K S D − μ max ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) ) ,

ω 2 = 8 K S μ max ( D − μ max ) h ¯ 2 K S 3 μ max 3 ( D − μ max ) 2 a + e b D K S D − μ max 2 ( D ( S F + K S ) − S F μ max ) × a + e b D K S D − μ max ( 2 ( a J 21 − J 22 ) ( D − μ max ) 3 + J 21 e b D K S D − μ max ( − 2 D 2 ( b K S + 3 ) μ max + D ( b K S ( b K S + 4 ) + 6 ) μ max 2 − 2 ( b K S + 1 ) μ max 3 + 2 D 3 ) ) + h ¯ ( a ( D − μ max ) × ( D 2 ( S F + K S ) + S F μ max ( μ max − 2 D ) ) + e b D K S D − μ max ( μ max ( S F μ max ( D ( b K S + 3 ) − μ max ) − D 2 ( b S F K S + b K S 2 + 3 S F + K S ) ) + D 3 ( S F + K S ) ) )

× ( e b D K S D − μ max ( μ max ( D 3 ( 2 J 21 K S ( a J ¯ 21 + 4 b K S + 12 ) + 8 S F J 21 ( b K S + 4 ) − J ¯ 21 J 22 K S ) + μ max ( D 2 ( − J 21 K S ( a J ¯ 21 ( b K S + 6 ) + 4 b K S ( b K S + 4 ) + 24 ) − 4 S F J 21 ( b K S ( b K S + 6 ) + 12 ) + 3 J ¯ 21 J 22 K S ) + D μ max ( J 21 K S ( a J ¯ 21 ( b K S + 6 ) + 8 b K S + 8 ) + 4 S F J 21 ( b K S + 2 ) ( b K S + 4 ) − 3 J ¯ 21 J 22 K S ) + μ max 2 ( J ¯ 21 K S ( J 22 − 2 a J 21 ) − 8 S F ( b J 21 K S + J 21 ) ) ) ) − 8 D 4 J 21 ( S F + K S ) ) + ( a J 21 − J 22 ) ( D − μ max ) 3 ( μ max ( a J ¯ 21 K S + 8 S F ) − 8 D ( S F + K S ) ) + J ¯ 21 J 21 K S μ max ( D − μ max ) e 2 b D K S D − μ max ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) ) ,

ω 3 = 1 K S μ max ( D − μ max ) a + e b D K S D − μ max 2 e b D K S D − μ max ( − D 2 μ max ( K S ( a J ¯ 21 J 21 + b J 21 K S + b J 22 K S + J 21 ) + b S F K S ( J 21 + J 22 ) + 3 S F J 21 ) + D μ max 2 ( a J ¯ 21 J 21 K S + b S F K S ( J 21 + J 22 ) + 3 S F J 21 ) + D 3 J 21 ( S F + K S ) − S F J 21 μ max 3 ) + a J 21 ( D − μ max ) ( − D μ max ( a J ¯ 21 K S + 4 S F ) + 2 D 2 ( S F + K S ) + 2 S F μ max 2 ) − D J ¯ 21 J 21 K S μ max ( D − μ max ) e 2 b D K S D − μ max ,

ω 4 = 1 K S 2 μ max 2 ( D − μ max ) 2 a + e b D K S D − μ max 2 h ¯ ( a ( D − μ max ) ( D 2 ( S F + K S ) + S F μ max ( μ max − 2 D ) ) + e b D K S D − μ max ( μ max ( S F μ max ( D ( b K S + 3 ) − μ max ) − D 2 ( b S F K S + b K S 2 + 3 S F + K S ) ) + D 3 ( S F + K S ) ) ) ( a J 21 ( D − μ max ) × ( D 2 ( S F + K S ) + S F μ max ( μ max − 2 D ) ) + e b D K S D − μ max ( μ max ( S F μ max ( b D K S ( J 21 + J 22 ) + 3 D J 21 − J 21 μ max ) − D 2 ( b S F K S ( J 21 + J 22 ) + K S ( b K S ( J 21 + J 22 ) + J 21 ) + 3 S F J 21 ) ) + D 3 J 21 ( S F + K S ) ) ) − K S μ max ( D − μ max ) a + e b D K S D − μ max × 2 e b D K S D − μ max ( − D 2 μ max ( K S ( − a J ¯ 21 J 21 + b J 21 K S + b J 22 K S + J 21 ) + b S F K S ( J 21 + J 22 ) + 3 S F J 21 ) + D μ max 2 ( − a J ¯ 21 J 21 K S + b S F K S ( J 21 + J 22 ) + 3 S F J 21 ) + D 3 J 21 ( S F + K S ) − S F J 21 μ max 3 ) + a J 21 ( D − μ max ) ( μ max ( a D J ¯ 21 K S − 4 D S F + 2 S F μ max ) + 2 D 2 ( S F + K S ) ) + D J ¯ 21 J 21 K S μ max ( D − μ max ) e 2 b D K S D − μ max ,

ω 5 = 1 3 K S 3 μ max 3 ( μ max − D ) a + e b D K S D − μ max 2 h ¯ ( μ max ( 2 D 4 ( 3 a 2 J ¯ 21 J 21 K S + e b D K S D − μ max ( 6 J 21 K S ( a J ¯ 21 + b K S + 4 ) + 6 S F J 21 ( b K S + 5 ) − J ¯ 21 J 22 K S ) + 30 a S F J 21 − a J ¯ 21 J 22 K S + 24 a J 21 K S + 3 J ¯ 21 J 21 K S e 2 b D K S D − μ max + 2 S F J 22 K S + 2 J 22 K S 2 ) + μ max ( μ max ( 2 μ max 2 ( 3 a 2 J ¯ 21 J 21 K S + 6 a S F J 21 − a J ¯ 21 J 22 K S + 2 S F J 22 K S ) + e b D K S D − μ max ( D 2 ( K S ( J 21 ( 3 a J ¯ 21 ( b K S ( b K S + 6 ) + 24 ) + 2 b K S ( b K S ( b K S + 6 ) + 18 ) + 48 ) − 12 J ¯ 21 J 22 ) + 2 S F J 21 ( b K S ( b K S ( b K S + 9 ) + 36 ) + 60 ) ) + μ max ( D ( K S ( 8 J ¯ 21 J 22 − 3 J 21 ( a J ¯ 21 ( b K S ( b K S + 6 ) + 16 ) + 2 b K S ( b K S + 2 ) + 4 ) ) − 2 S F J 21 ( b K S ( b K S ( b K S + 9 ) + 24 ) + 30 ) ) + 2 μ max ( J ¯ 21 K S ( 3 a J 21 ( b K S + 2 ) − J 22 ) + 3 S F J 21 ( b K S ( b K S + 2 ) + 2 ) ) ) ) + 12 D 2 ( K S ( a J 21 ( 3 a J ¯ 21 + 4 ) − a J ¯ 21 J 22 + 2 S F J 22 ) + 10 a S F J 21 + J 22 K S 2 ) − 4 D μ max ( K S ( a J 21 ( 6 a J ¯ 21 + 3 ) − 2 a J ¯ 21 J 22 + 4 S F J 22 ) + 15 a S F J 21 + J 22 K S 2 ) + 3 J ¯ 21 J 21 K S e 2 b D K S D − μ max ( D 2 ( b K S ( b K S + 6 ) + 12 ) − D ( b K S + 2 ) ( b K S + 4 ) μ max + 2 ( b K S + 1 ) μ max 2 ) ) + 2 D 3 ( − e b D K S D − μ max ( K S ( 3 J 21 ( a J ¯ 21 ( b K S + 8 ) + b K S ( b K S + 6 ) + 12 ) − 4 J ¯ 21 J 22 ) + 3 S F J 21 ( b K S ( b K S + 8 ) + 20 ) ) − 4 K S ( 3 a J 21 ( a J ¯ 21 + 3 ) − a J ¯ 21 J 22 + 2 S F J 22 ) − 60 a S F J 21 − 3 J ¯ 21 J 21 K S ( b K S + 4 ) e 2 b D K S D − μ max − 6 J 22 K S 2 ) ) ) − 12 D 5 J 21 ( S F + K S ) a + e b D K S D − μ max ) ,

ω 6 = e b D K S D − μ max K S 2 μ max 2 ( μ max − D ) a + e b D K S D − μ max ( − D 3 μ max ( 2 J 21 K S ( a J ¯ 21 + 2 b K S + 6 ) + 4 S F J 21 ( b K S + 4 ) − J ¯ 21 J 22 K S ) + D 2 μ max 2 ( K S ( J 21 ( a J ¯ 21 ( b K S + 6 ) + 2 b K S ( b K S + 4 ) + 12 ) − 3 J ¯ 21 J 22 ) + 2 S F J 21 ( b K S ( b K S + 6 ) + 12 ) ) − D μ max 3 ( J 21 K S ( a J ¯ 21 ( b K S + 6 ) + 4 b K S + 4 ) + 2 S F J 21 ( b K S + 2 ) ( b K S + 4 ) − 3 J ¯ 21 J 22 K S ) + μ max 4 ( J ¯ 21 K S ( 2 a J 21 − J 22 ) + 4 S F ( b J 21 K S + J 21 ) ) + 4 D 4 J 21 ( S F + K S ) ) + ( a J 21 − J 22 ) ( D − μ max ) 3 ( 4 D ( S F + K S ) − μ max ( a J ¯ 21 K S + 4 S F ) ) − J ¯ 21 J 21 K S μ max ( D − μ max ) e 2 b D K S D − μ max × ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) .

Assume that the central manifold W c (0, 0, 0) has the following approximate representation according to the central manifold theorem.

(24) W c ( 0 , 0 , 0 ) = { ( y t , z t , h ˜ ) ∈ R 3 : z t = H ( y t ) = k 1 y t 2 + k 2 h ˜ y t + k 3 h ˜ 2 , k 1 , k 2 , k 3 ∈ R } .

The central manifold given by (24) satisfies the following equation:

N ( H ( y t ) ) = H ( − y t + g 1 ( y t , H ( y t ) ) ) − λ 2 H ( y t ) − g 2 ( y t , H ( y t ) ) = 0 .

Hence,

N ( H ( y t ) ) = γ 1 y t 2 + γ 2 y t 3 + γ 3 h ˜ y t + γ 4 h ˜ y t 2 + γ 5 h ˜ 2 y t + γ 6 h ˜ 2 + γ 7 h ˜ 3 ,

where

γ 1 = k 1 − ω 1 − k 1 λ 2 γ 2 = − k 1 ( 2 ν 1 + ω 2 ) − ω 5 γ 3 = − k 2 − ω 3 − k 2 λ 2 γ 4 = k 2 ( ν 1 − ω 2 ) − k 1 ( 2 ν 3 + ω 4 ) − ω 6 γ 5 = − k 3 ω 2 + k 2 ( ν 3 − ω 4 ) γ 6 = k 3 − k 3 λ 2 γ 7 = − k 3 ω 4 .

By comparing the coefficients in the formula above, we obtain

k 1 = ω 1 1 − λ 2 k 2 = − ω 3 1 + λ 2 k 3 = 0 .

Therefore, system (22) reduces to the following one-dimensional system

(25) y t + 1 = − y t + g 1 ( y t , H ( y t ) ) ,

where

g 1 ( y t , H ( y t ) ) = ν 1 y t 2 + ( k 1 ν 1 + ν 5 ) y t 3 + ( k 2 ν 2 + k 1 ν 4 + ν 6 ) h ˜ y t 2 + k 2 ν 4 h ˜ 2 y t + ν 3 h ˜ y t .

Let us define f ( y , h ˜ ) to represent the right-hand side of system (25) as follows:

f ( y , h ˜ ) = − y t + g 1 ( y t , H ( y t ) ) .

To examine Theorem 2.5 hypotheses, we first compute the partial derivative of f ( y , h ˜ ) with respect to y .

f y ( y , h ˜ ) = − 1 + 2 ν 1 y + 3 ( k 1 ν 1 + ν 5 ) y 2 + 2 ( k 2 ν 2 + k 1 ν 4 + ν 6 ) h ˜ y + k 2 ν 4 h ˜ 2 + ν 3 h ˜ ,

and

f y ( 0 , 0 ) = − 1 , f y y ( 0 , 0 ) = 2 ν 1 , f y y y ( 0 , 0 ) = 6 ( k 1 ν 1 + ν 5 ) .

With these values, the following expression holds:

1 2 ( f y y ( 0 , 0 ) ) 2 + 1 3 ( f y y y ( 0 , 0 ) ) = 2 ( ν 1 2 + k 1 ν 1 + ν 5 ) , f y h ˜ ( 0 , 0 ) = ν 3 .

Thus, we have concluded the following result.

Theorem 4.2

System (2) at the positive fixed point E 1 undergoes flip bifurcation under the following conditions:

  1. α 2 > 4 α 1 ,

  2. ( ν 1 2 + k 1 ν 1 + ν 5 ) ≠ 0 ,

  3. ν 3 ≠ 0 ,

where α 1 , α 2 , and ν 1 , ν 3 , ν 5 are given in (8) and (23).

4.2 Neimark bifurcation analysis

To demonstrate the presence of a Neimark-Sacker bifurcation, we analyze the characteristic polynomial of the Jacobian matrix near the fixed point E 1 as the bifurcation parameter h is varied. The characteristic polynomial (11) yields a pair of complex conjugate eigenvalues of the form

λ 1,2 ( h ) = p ( h ) ± i q ( h ) = r ( h ) e ± i θ ( h ) ,

since B 2 − 4 C < 0 , that is, α 2 < 4 α 1 , where

p ( h ) = 2 − α 2 h 2 , q ( h ) = h α 2 ( 4 α 1 − α 2 ) 2 , r ( h ) = p 2 ( h ) + q 2 ( h ) = α 1 α 2 h 2 − α 2 h + 1 , θ ( h ) = arctan q ( h ) p ( h ) = arctan h α 2 ( 4 α 1 − α 2 ) 2 − α 2 h .

On the other hand, when h = 1 α 1 , that is,

C = h 2 α 1 α 2 − h α 2 + 1 = 1 .

Then, we have

∣ λ 1,2 1 α 1 ∣ = r 1 α 1 = 1 , θ ( h ) = arctan α 2 ( 4 α 1 − α 2 ) 2 α 1 − α 2 = θ 0 .

So, from the fixed point, a Neimark-Sacker bifurcation might develop. We employ the Neimark-Sacker bifurcation theorem [39] to demonstrate the existence of the Neimark-Sacker bifurcation. It is known that critical bifurcation parameter is h c = 1 α 1 and

d ( r ( h ) ) d h ∣ h c = 1 α 1 = α 2 2 > 0

because it is assumed that every parameter is positive. Since ( λ 1,2 ( h c ) ) k = e ± i k θ = 0 ≠ 1 for k = 1 , 2 , 3 , 4 (since α 2 > 0 ), all the conditions of the Neimark-Sacker bifurcation theorem [39] are satisfied. Thus, the Neimark-Sacker bifurcation occurs as h passes through h c = 1 α 1 = h ¯ from left to right. Afterward, we deduce this bifurcation’s normal form. By initially using u t = S t − S * , v t = x t − x * , h ˜ = h − h ¯ , we move E 1 = ( S * , x * ) with S * = D K S μ max − D , x * = S F ( μ max − D ) − D K S ( μ max − D ) a + e − b D K S μ max − D , and the critical bifurcation value h c = 1 α 1 to the origin, respectively. Hence, the system described by (2) transforms as in (17), where J 11 , J 12 , J 21 , f 1 ( u t , v t , h ¯ , h ˜ ) , f 2 ( u t , v t , h ¯ , h ˜ ) are given in (27) and (18). Taking h = 1 α 1 in (9), we obtain

(26) J 1 α 1 = J 11 1 α 1 J 12 1 α 1 J 21 1 α 1 1 = j 11 j 12 j 21 1 ,

where

(27) j 11 = μ max ( D 2 ( 2 S F + k ) − D S F μ max + k ) − D 3 ( S F + K S ) k μ max , j 12 = − D 2 ( D − μ max ) ( D ( S F + K S ) − S F μ max ) a + e b D k D − μ max μ max ( S F ( μ max − 2 D ) − b D K S ( D ( S F + K S ) − S F μ max ) ( D − μ max ) a e b D k μ max − D + 1 ) + D 2 ( S F + K S ) , j 21 = 1 2 D 2 ( S F + K S ) + 2 S F μ max ( μ max − 2 D ) K S μ max − 2 b D ( D ( S F + K S ) − S F μ max ) ( D − μ max ) a e b D k μ max − D + 1 × 2 D ( D − μ max ) ( D ( S F + K S ) − S F μ max ) S F − D k μ max − D D k μ max − D + k 2 a + e − b D k μ max − D .

We now construct a matrix T , whose columns are formed from the real and imaginary components of the eigenvector corresponding to λ 1 = p 1 α 1 − i q 1 α 1 , as follows:

T = − j 12 0 j 11 − p 1 α 1 q 1 α 1 ,

where p 1 α 1 = 1 − α 2 α 1 , q 1 α 1 = α 2 ( 4 α 1 − α 2 ) 2 α 1 . By using the following transform:

u t v t = T y t z t ,

the normal form of system (17) is obtained as follows:

y t + 1 z t + 1 = p 1 α 1 − q 1 α 1 q 1 α 1 p 1 α 1 y t z t + g 1 ( y t , z t ) g 2 ( y t , z t ) ,

where g 1 ( y t , z t ) and g 2 ( y t , z t ) are higher-order terms representing the nonlinearities of the system. These terms are given by

g 1 ( y t , z t ) = ζ 1 y t 2 + ζ 2 y t z t + ζ 3 y t 2 z t + ζ 4 y t 3 , g 2 ( y t , z t ) = ξ 1 y t 2 + ξ 2 y t z t + ξ 3 y t 2 z t + ξ 4 y t 3 ,

where

ζ 1 = 1 4 k 2 μ max 2 ( D − μ max ) h ¯ × 1 α 1 a + e b D K S D − μ max ( μ max ( D 3 ( − 2 e b D K S D − μ max ( 2 α 1 ( a ( j 11 − 1 ) K S + S F j 12 ( b K S + 4 ) + j 12 K S ( b K S + 3 ) ) + a α 2 K S ) − a ( 2 α 1 ( a ( j 11 − 1 ) K S + 8 S F j 12 + 6 j 12 K S ) + a α 2 K S ) + K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) ( − e 2 b D K S D − μ max ) ) + μ max ( e b D K S D − μ max ( D 2 ( 2 α 1 ( K S ( a ( j 11 − 1 ) ( b K S + 6 ) + j 12 ( b K S ( b K S + 4 ) + 6 ) ) + S F j 12 ( b K S ( b K S + 6 ) + 12 ) ) + a α 2 K S ( b K S + 6 ) ) − D μ max ( 2 α 1 ( a ( j 11 − 1 ) K S ( b K S + 6 ) + S F j 12 ( b K S + 2 ) ( b K S + 4 ) + 2 K S ( b j 12 K S + j 12 ) ) + a α 2 K S ( b K S + 6 ) ) + 2 μ max 2 ( 2 α 1 ( a ( j 11 − 1 ) K S + S F ( b j 12 K S + j 12 ) ) + a α 2 K S ) ) + a ( 3 D 2 ( a K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) + 8 α 1 S F j 12 + 4 α 1 j 12 K S ) − D μ max ( 6 a α 1 ( j 11 − 1 ) K S + 3 a α 2 K S + 4 α 1 j 12 ( 4 S F + K S ) ) + μ max 2 ( a K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) + 4 α 1 S F j 12 ) ) + K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) e 2 b D K S D − μ max ( D 2 ( b K S + 3 ) − D ( b K S + 3 ) μ max + μ max 2 ) ) ) ) + 4 D 4 j 12 ( S F + K S ) ,

ζ 2 = − h ¯ α 2 ( 4 α 1 − α 2 ) ( a ( D − μ max ) 2 + e b D K S D − μ max ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) ) 4 α 1 K S μ max ,

ζ 3 = h ¯ 4 α 1 K S 2 μ max 2 ( j 12 α 2 ( 4 α 1 − α 2 ) ( 2 a ( D − μ max ) 3 + e b D K S D − μ max ( − 2 D 2 ( b K S + 3 ) μ max + D ( b K S ( b K S + 4 ) + 6 ) μ max 2 − 2 ( b K S + 1 ) μ max 3 + 2 D 3 ) ) ) , ζ 4 = h ¯ j 12 12 α 1 K S 3 μ max 3 ( D − μ max ) a + e b D K S D − μ max × μ max ( 6 D 4 ( 2 e b D K S D − μ max ( α 1 ( 2 a ( j 11 − 1 ) K S + S F j 12 ( b K S + 5 ) + j 12 K S ( b K S + 4 ) ) + a α 2 K S ) + a ( 2 α 1 ( a ( j 11 − 1 ) K S + 5 S F j 12 + 4 j 12 K S ) + a α 2 K S ) + K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) e 2 b D K S D − μ max ) + μ max ( 6 D 3 ( − e b D K S D − μ max ( α 1 K S ( 2 a ( j 11 − 1 ) ( b K S + 8 ) + j 12 ( b K S ( b K S + 6 ) + 12 ) ) + a α 2 K S ( b K S + 8 ) + α 1 S F j 12 ( b K S ( b K S + 8 ) + 20 ) ) − 4 a ( a K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) + 5 α 1 S F j 12 + 3 α 1 j 12 K S ) + K S ( b K S + 4 ) ( 2 α 1 ( j 11 − 1 ) + α 2 ) ( − e 2 b D K S D − μ max ) ) + μ max ( e b D K S D − μ max ( D 2 ( 2 α 1 ( K S ( 3 a ( j 11 − 1 ) ( b K S ( b K S + 6 ) + 24 ) + j 12 ( b K S ( b K S ( b k + 6 ) + 18 ) + 24 ) ) + S F j 12 ( b K S ( b K S ( b k + 9 ) + 36 ) + 60 ) ) + 3 a α 2 K S ( b K S ( b k + 6 ) + 24 ) ) − D μ max ( 2 α 1 ( 3 K S ( a ( j 11 − 1 ) ( b K S ( b k + 6 ) + 16 ) + j 12 ( b K S ( b K S + 2 ) + 2 ) ) + S F j 12 ( b K S ( b K S ( b K S + 9 ) + 24 ) + 30 ) ) + 3 a α 2 K S ( b K S ( b K S + 6 ) + 16 ) ) + 6 μ max 2 ( 2 a α 1 ( j 11 − 1 ) K S ( b K S + 2 ) + a α 2 K S ( b k + 2 ) + α 1 S F j 12 ( b K S ( b k + 2 ) + 2 ) ) ) + 6 a ( 2 D 2 ( 3 a K S ( 2 α 1 ( j 11 − 1 ) + α 1 ) + 10 α 1 S F j 12 + 4 α 1 j 12 k ) − 2 D μ max ( α 1 K S ( 4 a ( j 11 − 1 ) + j 12 ) + 2 a α 2 k + 5 α 1 S F j 12 ) + μ max 2 ( a K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) + 2 α 1 S F j 12 ) ) + 3 K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) e 2 b D k D − μ max ( D 2 ( b K S ( b k + 6 ) + 12 ) − D ( b K S + 2 ) ( b K S + 4 ) μ max + 2 ( b K S + 1 ) μ max 2 ) ) ) ) − 12 α 1 D 5 j 12 ( S F + K S ) a + e b D k D − μ max ,

ξ 1 = h ¯ 4 α 2 k 2 α 1 ( 4 α 1 − α 2 ) μ max 2 ( D − μ max ) a + e b D k D − μ max × 4 α 1 D 4 j 12 ( S F + k ) ( − a ( 2 α 1 ( j 11 − 1 ) + α 2 ) + ( 2 α 1 ( j 11 − 1 ) + α 2 ) ( − e b D k D − μ max ) + 2 α 1 j 12 ) + D 3 μ max 2 ( 2 α 1 ( j 11 − 1 ) + α 2 ) e b D K S D − μ max ( α 1 ( 2 a ( j 11 − 1 ) K S + 2 S F j 12 ( b K S + 4 ) + j 12 K S ( 2 b K S + 5 ) ) + a α 2 K S ) + ( a ( 2 α 1 ( j 11 − 1 ) + α 2 ) − 2 α 1 j 12 ) ( 2 α 1 ( a ( j 11 − 1 ) K S + 8 S F j 12 + 6 j 12 k ) + a α 2 K S ) + K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) 2 e 2 b D K S D − μ max + D 2 μ max 2 − ( 2 α 1 ( j 11 − 1 ) + α 2 ) e b D k D − μ max ( 2 α 1 ( K S ( a ( j 11 − 1 ) ( b K S + 6 ) + j 12 ( b K S + 1 ) ( b k + 3 ) ) + S F j 12 ( b K S ( b k + 6 ) + 12 ) ) + a α 2 K S ( b K S + 6 ) ) − 3 ( a ( 2 α 1 ( j 11 − 1 ) + α 2 ) − 2 α 1 j 12 ) ( a K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) + 8 α 1 S F j 12 + 4 α 1 j 12 K S ) + K S ( b K S + 3 ) ( 2 α 1 ( j 11 − 1 ) + α 2 ) 2 ( − e 2 b D K S D − μ max ) + D μ max 3 ( 2 α 1 ( j 11 − 1 ) + α 2 ) e b D K S D − μ max ( 2 α 1 ( a ( j 11 − 1 ) K S ( b K S + 6 ) + S F j 12 ( b K S + 2 ) ( b K S + 4 ) + j 12 K S ( 2 b K S − 1 ) ) + a α 2 K S ( b K S + 6 ) ) + ( a ( 2 α 1 ( j 11 − 1 ) + α 2 ) − 2 α 1 j 12 ) ( 6 a α 1 ( j 11 − 1 ) K S + 3 a α 2 K S + 4 α 1 j 12 ( 4 S F + K S ) ) + K S ( b K S + 3 ) ( 2 α 1 ( j 11 − 1 ) + α 2 ) 2 e 2 b D K S D − μ max

+ μ max 4 − 2 ( 2 α 1 ( j 11 − 1 ) + α 2 ) e b D K S D − μ max ( a K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) + α 1 j 12 K S ( 2 b S F − 1 ) + 2 α 1 S F j 12 ) − ( a ( 2 α 1 ( j 11 − 1 ) + α 2 ) − 2 α 1 j 12 ) ( a K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) + 4 α 1 S F j 12 ) + K S ( 2 α 1 ( j 11 − 1 ) + α 2 ) 2 − e 2 b D K S D − μ max ,

ξ 2 = h ¯ 2 K S μ max 1 2 α 1 ( 2 α 1 ( j 11 − 1 ) + α 2 ) ( ( a ( D − μ max ) 2 + e b D K S D − μ max ( − D ( b K S + 2 ) μ max + D 2 + μ max 2 ) ) ) − j 12 ( D − μ max ) 2 ,

ξ 3 = j 12 h ¯ 6 K S 2 μ max 2 3 − α 2 2 α 1 − j 11 + 1 2 a ( D − μ max ) 3 + e b D K S D − μ max ( − 2 D 2 ( b K S + 3 ) μ max + D ( b K S ( b K S + 4 ) + 6 ) μ max 2 − 2 ( b K S + 1 ) μ max 3 + 2 D 3 ) ) + 2 j 12 ( D − μ max ) 3 ) ,

ξ 4 = − j 12 h ¯ 12 α 1 ( 4 α 1 − α 2 ) α 2 a + e b D K S D − μ max K S 3 ( D − μ max ) μ max 3 × − 12 α 1 a + e b D K S D − μ max ( α 1 + 2 α 1 ( j 11 − 1 ) ) j 12 ( S F + K S ) D 5 + 2 ( 3 e 2 b D K S D − μ max K S ( α 2 + 2 α 1 ( j 11 − 1 ) ) 2 + 3 a 2 K S ( α 1 + 2 α 1 ( j 11 − 1 ) ) 2 + 2 a ( 3 e b D K S D − μ max ( α 2 + 2 α 1 ( j 11 − 1 ) ) K S + α 1 j 12 ( 15 S F + 11 K S ) ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) + 2 α 1 e b D K S D − μ max j 12 ( 3 S F ( b K S + 5 ) + K S ( 3 b K S + 11 ) ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) + 4 α 1 2 j 12 2 K S ( S F + K S ) ) μ max D 4 + 2 ( − 12 a 2 K S ( α 1 + 2 α 1 ( j 11 − 1 ) ) 2 − 3 e 2 b D K S D − μ max K S ( b K S + 4 ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) 2 − 4 a α 1 j 12 ( 15 S F + 7 K S ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) − e b D K S D − μ max ( 3 a α 2 K S ( b K S + 8 ) + 3 α 1 S F j 12 ( b K S ( b K S + 8 ) + 20 ) + α 1 K S ( 6 a ( j 11 − 1 ) ( b K S + 8 ) + j 12 ( 3 b K S ( b K S + 6 ) + 28 ) ) ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) − 4 α 1 2 j 12 2 K S ( 4 S F + 3 K S ) ) μ max 2 D 3 + ( 3 e 2 b D K S D − μ max K S ( b K S ( b K S + 6 ) + 12 ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) 2 + e b D K S D − μ max ( 3 a α 2 K S ( b K S ( b K S + 6 ) + 24 ) + 2 α 1 ( S F j 12 ( b K S ( b K S ( b K S + 9 ) + 36 ) + 60 ) + K S ( 3 a ( j 11 − 1 ) ( b K S ( b K S + 6 ) + 24 ) + j 12 ( b K S ( b K S ( b K S + 6 ) + 18 ) + 12 ) ) ) ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) + 12 ( 3 a 2 K S ( α 2 + 2 α 1 ( j 11 − 1 ) ) 2 + 2 a α 1 j 12 ( 5 S F + K S ) ( α 1 + 2 α 1 ( j 11 − 1 ) ) + 2 α 1 2 j 12 2 K S ( 2 S F + K S ) ) ) μ max 3 D 2

− ( 24 a 2 K S ( α 2 + 2 α 1 ( j 11 − 1 ) ) 2 + 3 e 2 b D K S D − μ max K S ( b K S + 2 ) ( b K S + 4 ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) 2 + 4 a α 1 j 12 ( 15 S F − K S ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) + e b D K S D − μ max ( 3 a α 2 K S ( b K S ( b K S + 6 ) + 16 ) + 2 α 1 ( K S ( j 12 ( 3 b K S ( b K S + 2 ) − 2 ) + 3 a ( j 11 − 1 ) ( b K S ( b K S + 6 ) + 16 ) ) + S F j 12 ( b K S ( b K S ( b K S + 9 ) + 24 ) + 30 ) ) ) × ( α 2 + 2 α 1 ( j 11 − 1 ) ) + 8 α 1 2 j 12 2 K S ( 4 S F + K S ) ) μ max 4 D + 2 ( 3 a 2 K S ( α 2 + 2 α 1 ( j 11 − 1 ) ) 2 + 3 e 2 b D K S D − μ max K S ( b K S + 1 ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) 2 + a ( 2 α 1 j 12 ( 3 S F − K S ) + 3 e b D K S D − μ max × ( α 2 + 2 α 1 ( j 11 − 1 ) ) K S ( b K S + 2 ) ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) + α 1 e b D K S D − μ max j 12 ( 3 S F ( b K S ( b K S + 2 ) + 2 ) − 2 K S ) ( α 2 + 2 α 1 ( j 11 − 1 ) ) + 4 α 1 2 S F j 12 2 K S ) μ max 5 ,

where α 1 , α 2 , and j 11 , j 12 are as in (8) and (27), respectively.

To determine the direction of the Neimark-Sacker bifurcation, we use the formula [42] given by:

(28) a 1 α 1 = − R e ( 1 − 2 λ ) λ ¯ 2 1 − λ η 20 η 11 − 1 2 ∣ η 11 ∣ 2 − ∣ η 02 ∣ 2 + R e ( λ ¯ η 21 ) ,

where λ = λ 1 , λ ¯ = λ 2 , and

η 20 = 1 8 [ ( g y y 1 − g z z 1 + 2 g y z 2 ) + i ( g y y 2 − g z z 2 − 2 g y z 1 ) ] , η 11 = 1 4 [ ( g y y 1 + g z z 1 ) + i ( g y y 2 + g z z 2 ) ] , η 02 = 1 8 [ ( g y y 1 − g z z 1 − 2 g y z 2 ) + i ( g y y 2 − g z z 2 + 2 g y z 1 ) ] , η 21 = 1 16 [ ( g y y y 1 + g y z z 1 + g y y z 2 + g z z z 2 ) + i ( g y y y 2 + g y z z 2 − g y y z 1 − g z z z 1 ) ] ,

which are calculated at the point ( 0 , 0 , 1 α 1 ) .

As a result, we can give the following theorem.

Theorem 4.3

System (2) undergoes a Neimark-Sacker bifurcation at the fixed point E 1 = ( S * , x * ) if α 2 < 4 α 1 and a 1 α 1 ≠ 0 , this bifurcation is supercritical (subcritical) if a 1 α 1 < 0 ( a 1 α 1 > 0 ).

5 Numerical simulations

The goal is to demonstrate the validity of the key findings in this section. We perform numerical calculations, including phase diagrams and bifurcation diagrams, for system (2). The values of the parameters used are chosen from [43], and a , b are estimated.

Example 5.1

By using the values, the parameters ( D , S F , μ max , K S , a , b ) = ( 6 × 1 0 − 2 , 1 0 − 4 , 0.81 , 3 × 1 0 − 4 , 1.5 , 0.5 ) and h = 1 , the dynamical behavior of system (2) can be shown to change. Because h < 11.3685 , that is, h = 1 , the positive fixed point E 1 = ( 0.000024 , 0.0000304001 ) is a sink, and if h is chosen in the interval 11.3685 < h < 33.3331 , it is a saddle, otherwise it is source for t ∈ [ 0 , 1000 ] in Figure 1(a). Furthermore, setting the same parameter values with the same initial condition ( S 0 , x 0 ) = ( 1 0 − 5 , 1 0 − 5 ) , the phase portrait diagram for t ∈ [ 0 , 1000 ] is obtained as shown in Figure 1(b). The blue curve goes to the equilibrium point E 1 , which is shown with the pink star by starting from the initial condition with the black star since E 1 is stable. The above results support Theorem 3.2. Also, it is known that flip bifurcation occurs at E 1 = ( 0.000024 , 0.0000304001 ) for h = 11.3685 obtained in Theorem 4.2. For t ∈ [ 0 , 1000 ] and h ∈ [ 0 , 17.5 ] , the image is shown for S t and x t , respectively, in Figure 2 (a) and (b).

These results underscore the importance of selecting an appropriate time step h to ensure stable dynamics. For the forward Euler method, stability requires h < 2 ∣ λ max ∣ . Here λ max = 0.20083 , so h < 2 0.20083 ≈ 9.9588 . For stability, the time step h must satisfy h < 9.9588 . Choosing h larger than this value may lead to instability, causing the numerical solution to diverge from the expected behavior. On the other hand, selecting a smaller value for h ensures greater numerical accuracy, but this comes at the expense of increased computational cost. Thus, the stability of the system is intricately tied to the step size: choosing a step size that is too large compromises stability, while a smaller step size ensures more precise numerical results at the expense of higher computational demand. Balancing stability and accuracy is key when implementing numerical simulations of such systems.

Figure 1 (a) Stability analysis of the equilibrium point E 1 {E}_{1} = (0.000024, 0.0000304001) for system (2), showing convergence behavior around E 1 {E}_{1} . (b) Phase portrait in the ( S t , x t ) \left({S}_{t},{x}_{t}) -plane illustrating the trajectory of the system under parameters D = 6 × 1 0 − 2 D=6\times 1{0}^{-2} , S F = 1 0 − 4 {S}_{F}=1{0}^{-4} , μ max = 0.81 {\mu }_{\max }=0.81 , K S = 3 × 1 0 − 4 {K}_{S}=3\times 1{0}^{-4} , a = 1.5 a=1.5 , b = 0.5 b=0.5 , h = 1 h=1 , with initial condition ( S 0 , x 0 ) = ( 1 0 − 5 , 1 0 − 5 ) \left({S}_{0},{x}_{0})=\left(1{0}^{-5},1{0}^{-5}) .
Figure 1

(a) Stability analysis of the equilibrium point E 1 = (0.000024, 0.0000304001) for system (2), showing convergence behavior around E 1 . (b) Phase portrait in the ( S t , x t ) -plane illustrating the trajectory of the system under parameters D = 6 × 1 0 − 2 , S F = 1 0 − 4 , μ max = 0.81 , K S = 3 × 1 0 − 4 , a = 1.5 , b = 0.5 , h = 1 , with initial condition ( S 0 , x 0 ) = ( 1 0 − 5 , 1 0 − 5 ) .

Figure 2 Bifurcation diagram of system (2) with respect to the parameter h ∈ [ 0 , 17.5 ] h\in \left[0,17.5] , illustrating qualitative changes in system dynamics. The parameters are fixed as D = 6 × 1 0 − 2 D=6\times 1{0}^{-2} , S F = 1 0 − 4 {S}_{F}=1{0}^{-4} , μ max = 0.81 {\mu }_{\max }=0.81 , K S = 3 × 1 0 − 4 {K}_{S}=3\times 1{0}^{-4} , a = 1.5 a=1.5 , b = 0.5 b=0.5 , with the initial condition ( S 0 , x 0 ) = ( 1 0 − 5 , 1 0 − 5 ) \left({S}_{0},{x}_{0})=\left(1{0}^{-5},1{0}^{-5}) .
Figure 2

Bifurcation diagram of system (2) with respect to the parameter h ∈ [ 0 , 17.5 ] , illustrating qualitative changes in system dynamics. The parameters are fixed as D = 6 × 1 0 − 2 , S F = 1 0 − 4 , μ max = 0.81 , K S = 3 × 1 0 − 4 , a = 1.5 , b = 0.5 , with the initial condition ( S 0 , x 0 ) = ( 1 0 − 5 , 1 0 − 5 ) .

Example 5.2

Let us consider system (2) with the parameters ( D , S F , μ max , K S , a , b ) = ( 6 × 1 0 − 2 , 1 0 − 4 , 0.68 , 3.1 × 1 0 − 4 , 1.5 , 15 × 1 0 3 ) and the initial conditions ( S ( 0 ) , x ( 0 ) ) = ( 1 0 − 5 , 1 0 − 5 ) . When h is treated as a bifurcation parameter, we observe that α 2 − 4 α 1 = − 0.0125747 < 0 , indicating that the Jacobian matrix has a pair of complex conjugate eigenvalues. The critical bifurcation value is h c = 1 α 1 = 22.0471 . If we choose h = 21.5 < h c = 22.0471 , then the unique positive fixed point E 1 = ( 0.00003 , 0.0000327466 ) is sink as shown in Figure 3 (a) and phase portrait as shown in Figure 3(b). Otherwise, E 1 becomes the source and system (2) undergoes Neimark-Sacker bifurcation as in Figure 4(a) and (b) for S t and x t , respectively. On the other hand, using equation (28), we can obtain a 1 α 1 = 4.13005 × 1 0 10 > 0 , which shows that the Neimark-Sacker bifurcation is subcritical.

From the stability perspective, for the forward Euler method, the step size h must satisfy the stability condition h < 2 ∣ λ max ∣ . Given that λ max = 0.15167 , we obtain h < 2 0.15167 ≈ 13.1864 . Therefore, to ensure numerical stability, the time step must be less than 13.1864. If the time step is chosen larger than this threshold, instability may arise, and the numerical results may become inaccurate. Conversely, selecting a smaller step size improves the accuracy of the results, though it increases the computational cost. Thus, the choice of h is crucial for balancing both stability and accuracy in numerical simulations.

Figure 3 (a) Stability analysis of the equilibrium point E 1 = ( 0.00003 , 0.0000327466 ) {E}_{1}=\left(0.00003,0.0000327466) of system (2), showing local convergence behavior near E 1 {E}_{1} . (b) Phase portrait in the ( S t , x t ) \left({S}_{t},{x}_{t}) -plane for t ∈ [ 0 , 1 , 000 ] t\in \left[0,1,000] , illustrating the system’s trajectory under the parameter set: D = 6 × 1 0 − 2 D=6\times 1{0}^{-2} , S F = 1 0 − 4 {S}_{F}=1{0}^{-4} , μ max = 0.68 {\mu }_{\max }=0.68 , K S = 3.1 × 1 0 − 4 {K}_{S}=3.1\times 1{0}^{-4} , a = 1.5 a=1.5 , b = 1.5 × 1 0 4 b=1.5\times 1{0}^{4} , h = 21.5 h=21.5 , with initial condition ( S 0 , x 0 ) = ( 1 0 − 5 , 1 0 − 5 ) \left({S}_{0},{x}_{0})=\left(1{0}^{-5},1{0}^{-5}) .
Figure 3

(a) Stability analysis of the equilibrium point E 1 = ( 0.00003 , 0.0000327466 ) of system (2), showing local convergence behavior near E 1 . (b) Phase portrait in the ( S t , x t ) -plane for t ∈ [ 0 , 1 , 000 ] , illustrating the system’s trajectory under the parameter set: D = 6 × 1 0 − 2 , S F = 1 0 − 4 , μ max = 0.68 , K S = 3.1 × 1 0 − 4 , a = 1.5 , b = 1.5 × 1 0 4 , h = 21.5 , with initial condition ( S 0 , x 0 ) = ( 1 0 − 5 , 1 0 − 5 ) .

Figure 4 Bifurcation diagram of system (2) with respect to the parameter h ∈ [ 21.8 , 22.3 ] h\in \left[21.8,22.3] , illustrating dynamic transitions within this range. The parameters are set as D = 6 × 1 0 − 2 D=6\times 1{0}^{-2} , S F = 1 0 − 4 {S}_{F}=1{0}^{-4} , μ max = 0.68 {\mu }_{\max }=0.68 , K S = 3.1 × 1 0 − 4 {K}_{S}=3.1\times 1{0}^{-4} , a = 1.5 a=1.5 , b = 1.5 × 1 0 4 b=1.5\times 1{0}^{4} , with the initial condition ( S 0 , x 0 ) = ( 1 0 − 5 , 1 0 − 5 ) \left({S}_{0},{x}_{0})=\left(1{0}^{-5},1{0}^{-5}) .
Figure 4

Bifurcation diagram of system (2) with respect to the parameter h ∈ [ 21.8 , 22.3 ] , illustrating dynamic transitions within this range. The parameters are set as D = 6 × 1 0 − 2 , S F = 1 0 − 4 , μ max = 0.68 , K S = 3.1 × 1 0 − 4 , a = 1.5 , b = 1.5 × 1 0 4 , with the initial condition ( S 0 , x 0 ) = ( 1 0 − 5 , 1 0 − 5 ) .

6 Conclusion

In this study, we successfully constructed a discrete-time chemostat model by discretizing a continuous-time population model using the forward Euler method. Our analysis of the model’s fixed points and their stability provided insight into the dynamic behavior of the system. By applying the central manifold theorem and bifurcation theory, we demonstrated that the system exhibits flip and Neimark-Sacker bifurcations as the step size is varied, indicating the presence of complex dynamics and potential for periodic behavior. The numerical simulations validated our theoretical results and highlighted novel dynamic phenomena that could be of interest in understanding ecological and biological systems modeled by chemostats. The findings emphasize the importance of step size in influencing the stability and dynamics of discrete-time models and suggest that careful consideration of discretization methods is essential in studying population dynamics and other systems subject to similar modeling approaches. Future work could explore the impact of other discretization methods or further investigate the effects of varying parameters in more complex models, enhancing the understanding of bifurcations and stability in discrete-time systems.

Acknowledgments

The author sincerely thanks anonymous referees for their careful review of the manuscript and their valuable comments.

  1. Funding information: The author states no funding information.

  2. Author contribution: The author confirms the sole responsibility for the conception of the study, presented results and manuscript preparation.

  3. Conflict of interest: The author declares that there is no conflict of interest.

  4. Ethical approval: This study does not involve human participants or animals; therefore, ethical approval is not required.

  5. Data availability statement: Data sharing is not applicable to this article as no data sets were generated or analyzed during the current study.

References

[1] H. L. Smith and P. Waltman, The Theory of the Chemostat: Dynamics of Microbial Competition, Cambridge University Press, Cambridge, 1995. 10.1017/CBO9780511530043Search in Google Scholar

[2] A. Novick and L. Szilard, Description of the chemostat, Science 112 (1950), no. 2920, 715–716, DOI: https://doi.org/10.1126/science.112.2920.715. 10.1126/science.112.2920.715Search in Google Scholar PubMed

[3] D. S. Kompala, D. Ramkrishna, N. B. Jansen, and G. T. Tsao, Investigation of bacterial growth on mixed substrates: experimental evaluation of cybernetic models, Biotechnol. Bioeng. 28 (1986), no. 7, 1044–1055. 10.1002/bit.260280715Search in Google Scholar PubMed

[4] P. M. Doran, Bioprocess Engineering Principles, Academic Press, San Diego, 1995. 10.1016/B978-012220855-3/50013-4Search in Google Scholar

[5] M. Chi and W. Zhao, Dynamical analysis of two-microorganism and single nutrient stochastic chemostat model with Mnod-Haldane response function, Complexity 2019 (2019), no. 1, 8719067, DOI: https://doi.org/10.1155/2019/8719067. 10.1155/2019/8719067Search in Google Scholar

[6] Z. Zhao, L. Chen, and X. Song, Extinction and permanence of chemostat model with pulsed input in a polluted environment, Commun. Nonlinear Sci. Numer. Simul. 14 (2009), no. 4, 1737–1745, DOI: https://doi.org/10.1016/j.cnsns.2008.01.009. 10.1016/j.cnsns.2008.01.009Search in Google Scholar

[7] J. Jiao, X. Yang, L. Chen, and S. Cai, Effect of delayed response in growth on the dynamics of a chemostat model with impulsive input, Chaos Solitons Fractals 42 (2009), no. 4, 2280–2287, DOI: https://doi.org/10.1016/j.chaos.2009.03.132. 10.1016/j.chaos.2009.03.132Search in Google Scholar

[8] J. Shi, Y. Wu, and X. Zou, Coexistence of competing species for intermediate dispersal rates in a reaction-diffusion chemostat model, J. Dynam. Differential Equations 32 (2020), 1085–1112, DOI: https://doi.org/10.1007/s10884-019-09763-0. 10.1007/s10884-019-09763-0Search in Google Scholar

[9] E. O. Alzahrani, M. M. El-Dessoky, and P. Dogra, Global dynamics of a cell quota-based model of light-dependent algae growth in a chemostat, Commun. Nonlinear Sci. Numer. Simul. 90 (2020), 105295, DOI: https://doi.org/10.1016/j.cnsns.2020.105295. 10.1016/j.cnsns.2020.105295Search in Google Scholar

[10] R. Baratti, J. Alvarez, S. Tronci, M. Grosso, and A. Schaum, Characterization with Fokker-Planck theory of the nonlinear stochastic dynamics of a class of two-state continuous bioreactors, J. Process Control 102 (2021), 66–84, DOI: https://doi.org/10.1016/j.jprocont.2021.04.004. 10.1016/j.jprocont.2021.04.004Search in Google Scholar

[11] S. Ruan and G. S. K. Wolkowicz, Bifurcation analysis of a chemostat model with a distributed delay, J. Math. Anal. Appl. 204 (1996), no. 3, 786–812, DOI: https://doi.org/10.1006/jmaa.1996.0468. 10.1006/jmaa.1996.0468Search in Google Scholar

[12] S. Yuan and M. Han, Bifurcation analysis of a chemostat model with two distributed delays, Chaos Solitons Fractals 20 (2004), no. 5, 995–1004, DOI: https://doi.org/10.1016/j.chaos.2003.09.048. 10.1016/j.chaos.2003.09.048Search in Google Scholar

[13] X. Xu, Y. Qiu, X. Chen, H. Zhang, Z. Liang, and B. Tian, Bifurcation analysis of a food chain chemostat model with Michaelis-Menten functional response and double delays, AIMS Math. 7 (2022), no. 7, 12154–12176, DOI: https://doi.org/10.3934/math.2022676. 10.3934/math.2022676Search in Google Scholar

[14] T. Bayen, J. Harmand, and M. Sebbah, Time-optimal control of concentration changes in the chemostat with one single species, Appl. Math. Model. 50 (2017), 257–278, DOI: https://doi.org/10.1016/j.apm.2017.05.037. 10.1016/j.apm.2017.05.037Search in Google Scholar

[15] P. De Leenheer, B. Li, and H. L. Smith, Competition in the chemostat: some remarks, Can. Appl. Math. Q. 11 (2003), no. 3, 229–248. Search in Google Scholar

[16] S. B. Hsu, T. K. Luo, and P. Waltman, Competition between plasmid-bearing and plasmid-free organisms in a chemostat with an inhibitor, J. Math. Biol. 34 (1995), no. 2, 225–238, DOI: https://doi.org/10.1007/bf00178774. 10.1007/BF00178774Search in Google Scholar

[17] W. Li, J. Ji, L. Huang, and Y. Zhang, Complex dynamics and impulsive control of a chemostat model under the ratio threshold policy, Chaos Solitons Fractals 167 (2023), 113077, DOI: https://doi.org/10.1016/j.chaos.2022.113077. 10.1016/j.chaos.2022.113077Search in Google Scholar

[18] A. Schaum, S. Tronci, R. Baratti, and J. Alvarez, On the dynamics and robustness of the chemostat with multiplicative noise, IFAC-Pap. 54 (2021), no. 3, 342–347, DOI: https://doi.org/10.1016/j.ifacol.2021.08.265. 10.1016/j.ifacol.2021.08.265Search in Google Scholar

[19] G. S. K. Wolkowicz, H. Xia, and S. Ruan, Competition in the chemostat: a distributed delay model and its global asymptotic behavior, SIAM J. Appl. Math. 57 (1997), no. 5, 1281–1310, DOI: https://doi.org/10.1137/s0036139995289842. 10.1137/S0036139995289842Search in Google Scholar

[20] S. Yuan and T. Zhang, Dynamics of a plasmid chemostat model with periodic nutrient input and delayed nutrient recycling, Nonlinear Anal. Real World Appl. 13 (2012), no. 5, 2104–2119, DOI: https://doi.org/10.1016/j.nonrwa.2012.01.006. 10.1016/j.nonrwa.2012.01.006Search in Google Scholar

[21] T. B. Gage, F. M. Williams, and J. B. Horton Division synchrony and the dynamics of microbial populations: A size-specific model, Theor. Popul. Biol. 26 (1984), no. 3, 296–314, DOI: https://doi.org/10.1016/0040-5809(84)90035-2. 10.1016/0040-5809(84)90035-2Search in Google Scholar PubMed

[22] H. L. Smith, A discrete, size-structured model of microbial growth and competition in the chemostat, J. Math. Biol. 34 (1996), no. 7, 734–754, DOI: https://doi.org/10.1007/bf00161517. 10.1007/BF00161517Search in Google Scholar

[23] H. L. Smith and X. Q. Zhao, Competitive exclusion in a discrete-time, size-structured chemostat model, Discret. Contin. Dyn. Syst. - B 1 (2001), no. 2, 183–192, DOI: https://doi.org/10.3934/dcdsb.2001.1.183. 10.3934/dcdsb.2001.1.183Search in Google Scholar

[24] J. Arino, J. L. Gouzé, and A. Sciandra, A discrete, size-structured model of phytoplankton growth in the chemostat, PhD Thesis, INRIA 2000. Search in Google Scholar

[25] J. Arino, J. L. Gouzé, and A. Sciandra, A discrete, size-structured model of phytoplankton growth in the chemostat: Introduction of inhomogeneous cell division size, J. Math. Biol. 45 (2002), no. 4, 313–336, DOI: https://doi.org/10.1007/s002850200160. 10.1007/s002850200160Search in Google Scholar PubMed

[26] D. Zhang, X. Cai, and L. Wang, Complex dynamics in a discrete-time size-structured chemostat model with inhibitory kinetics, Discrete Contin. Dyn. Syst. Ser. B 24 (2019), no. 7, 3439–3451, DOI: https://doi.org/10.3934/dcdsb.2018327. 10.3934/dcdsb.2018327Search in Google Scholar

[27] S. R. J. Jang, A discrete, size-structured chemostat model of plasmid-bearing and plasmid-free competition, J. Difference Equ. Appl. 11 (2005), no. 7, 619–633, DOI: https://doi.org/10.1080/10236190412331334509. 10.1080/10236190412331334509Search in Google Scholar

[28] P. Amster, G. Robledo, and D. Sepulveda, Dynamics of a discrete size-structured chemostat with variable nutrient supply, Commun. Nonlinear Sci. Numer. Simul. 132 (2024), 107904, DOI: https://doi.org/10.2139/ssrn.4569230. 10.1016/j.cnsns.2024.107904Search in Google Scholar

[29] K. S. N. Al-Basyouni and A. Q. Khan, Bifurcation analysis of a discrete-time chemostat model, Math. Probl. Eng. 2023 (2023), 7518261, DOI: https://doi.org/10.1155/2023/7518261. 10.1155/2023/7518261Search in Google Scholar

[30] I. M. Alsulami, On the stability, chaos and bifurcation analysis of a discrete-time chemostat model using the piecewise constant argument method, AIMS Math. 9 (2024), no. 12, 33861–33878, DOI: https://doi.org/10.3934/math.20241615. 10.3934/math.20241615Search in Google Scholar

[31] M. L. Büyükkahraman, Bifurcation analysis and chaos control of simple chemostat model with discrete time, Ann. Math. Sci. Appl. 10 (2025), no. 1, 39–59, DOI: https://doi.org/10.4310/amsa.250305060812. 10.4310/AMSA.250305060812Search in Google Scholar

[32] J. Yang, Y. Tan, and R. A. Cheke, Complex dynamics of an impulsive chemostat model, Int. J. Bifurcation Chaos 29 (2019), no. 8, 1950101, DOI: https://doi.org/10.1142/S0218127419501013. 10.1142/S0218127419501013Search in Google Scholar

[33] R. Ma, Y. Bai, and F. Wang, Dynamical behavior analysis of a two-dimensional discrete predator-prey model with prey refuge and fear factor, J. Appl. Anal. Comput. 10 (2020), no. 4, 1683–1697, DOI: https://doi.org/10.11948/20190426. 10.11948/20190426Search in Google Scholar

[34] Q. Zhou, F. Chen, and S. Lin, Complex dynamics analysis of a discrete amensalism system with a cover for the first species, Axioms 11 (2022), no. 8, 365, DOI: https://doi.org/10.3390/axioms11080365. 10.3390/axioms11080365Search in Google Scholar

[35] Q. Din, Complexity and chaos control in a discrete-time prey-predator model, Commun. Nonlinear Sci. Numer. Simul. 49 (2017), 113–134, DOI: https://doi.org/10.1016/j.cnsns.2017.01.025. 10.1016/j.cnsns.2017.01.025Search in Google Scholar

[36] Z. Hu, Z. Teng, and L. Zhang, Stability and bifurcation analysis of a discrete predator-prey model with nonmonotonic functional response, Nonlinear Anal. Real World Appl. 12 (2011), no. 4, 2356–2377, DOI: https://doi.org/10.1016/j.nonrwa.2011.02.009. 10.1016/j.nonrwa.2011.02.009Search in Google Scholar

[37] X. Liu and D. Xiao, Complex dynamic behaviors of a discrete-time predator-prey system, Chaos Solitons Fractals 32 (2007), no. 1, 80–94, DOI: https://doi.org/10.1016/j.chaos.2005.10.081. 10.1016/j.chaos.2005.10.081Search in Google Scholar

[38] S. Wiggins, S. Wiggins, and M. Golubitsky, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer, New York, 2003. Search in Google Scholar

[39] Y. A. Kuznetsov, I. A. Kuznetsov, and Y. Kuznetsov, Elements of Applied Bifurcation Theory, Springer, New York, 1998. Search in Google Scholar

[40] P. Baydemir, H. Merdan, E. Karaoglu, and G. Sucu, Complex dynamics of a discrete-time prey-predator system with Leslie type: stability, bifurcation analyses and chaos, Int. J. Bifurcation Chaos 30 (2020), no. 10, 2050149, DOI: https://doi.org/10.1142/s0218127420501497. 10.1142/S0218127420501497Search in Google Scholar

[41] S. P. Otto and T. Day, A Biologist’s Guide to Mathematical Modeling in Ecology and Evolution, Princeton University Press, New Jersey, 2011. 10.2307/j.ctvcm4hndSearch in Google Scholar

[42] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, vol. 42, Springer Science & Business Media, New York, 2013. Search in Google Scholar

[43] S. R. Hansen and S. P. Hubbell, Single-nutrient microbial competition: qualitative agreement between experimental and theoretically forecast outcomes, Science 207 (1980), no. 4438, 1491–1493, DOI: https://doi.org/10.1126/science.6767274. 10.1126/science.6767274Search in Google Scholar PubMed

[44] R. L. Burden and J. D. Faires, Numerical Analysis, 9th International Edition, Brooks/Cole, Cencag Learning, 2011. Search in Google Scholar
Received: 2024-10-09
Revised: 2025-06-08
Accepted: 2025-06-25
Published Online: 2025-08-30

© 2025 the author(s), published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

  1. A new perspective on generalized Laguerre polynomials
  2. Research Articles
  3. On approximation by Stancu variant of Bernstein-Durrmeyer-type operators in movable compact disks
  4. Circular n,m-rung orthopair fuzzy sets and their applications in multicriteria decision-making
  5. Grand Triebel-Lizorkin-Morrey spaces
  6. Coefficient estimates and Fekete-Szegö problem for some classes of univalent functions generalized to a complex order
  7. Proofs of two conjectures involving sums of normalized Narayana numbers
  8. On the Laguerre polynomial approximation errors and lower type of entire functions of irregular growth
  9. New convolutions for the Hartley integral transform imbedded in the Banach algebras and convolution-type integral equations
  10. Some inequalities for rational function with prescribed poles and restricted zeros
  11. Lucas difference sequence spaces defined by Orlicz function in 2-normed spaces
  12. Evaluating the efficacy of fuzzy Bayesian networks for financial risk assessment
  13. Fixed point results for contractions of polynomial type
  14. Estimation for spatial semi-functional partial linear regression model with missing response at random
  15. Investigating the controllability of differential systems with nonlinear fractional delays, characterized by the order 0 < η ≤ 1 < ζ ≤ 2
  16. New forms of bilateral inequalities for K-g-frames
  17. Rate of pole detection using Padé approximants to polynomial expansions
  18. Existence results for nonhomogeneous Choquard equation involving p-biharmonic operator and critical growth
  19. Note on the shape-preservation of a new class of Kantorovich-type operators via divided differences
  20. Geršhgorin-type theorems for Z1-eigenvalues of tensors with applications
  21. New topologies derived from the old one via operators
  22. Blow up solutions for two-dimensional semilinear elliptic problem of Liouville type with nonlinear gradient terms
  23. Infinitely many normalized solutions for Schrödinger equations with local sublinear nonlinearity
  24. Nonparametric expectile shortfall regression for functional data
  25. Advancing analytical solutions: Novel wave insights and methodologies for beta fractional Kuralay-II equations
  26. A generalized p-Laplacian problem with parameters
  27. A study of solutions for several classes of systems of complex nonlinear partial differential difference equations in ℂ2
  28. Towards finding equalities involving mixed products of the Moore-Penrose and group inverses by matrix rank methodology
  29. ω ℒ -biprojective and ω ¯ -contractible Banach algebras
  30. Coefficient functionals for Sakaguchi-type-Starlike functions subordinated to the three-leaf function
  31. Solutions of several general quadratic partial differential-difference equations in ℂ2
  32. Inequalities for the generalized trigonometric functions with respect to weighted power mean
  33. Optimization of Lagrange problem with higher-order differential inclusion and special boundary-value conditions
  34. Hankel determinants for q-starlike functions connected with q-sine function
  35. System of partial differential hemivariational inequalities involving nonlocal boundary conditions
  36. A new family of multivalent functions defined by certain forms of the quantum integral operator
  37. A matrix approach to compare BLUEs under a linear regression model and its two competing restricted models with applications
  38. Weighted composition operators on bicomplex Lorentz spaces with their characterization and properties
  39. Behavior of spatial curves under different transformations in Euclidean 4-space
  40. Commutators for the maximal and sharp functions with weighted Lipschitz functions on weighted Morrey spaces
  41. A new kind of Durrmeyer-Stancu-type operators
  42. A study of generalized Mittag-Leffler-type function of arbitrary order
  43. On the approximation of Kantorovich-type Szàsz-Charlier operators
  44. Split quaternion Fourier transforms for two-dimensional real invariant field
  45. Review Article
  46. Characterization generalized derivations of tensor products of nonassociative algebras
  47. Special Issue on Differential Equations and Numerical Analysis - Part II
  48. Existence and optimal control of Hilfer fractional evolution equations
  49. Persistence of a unique periodic wave train in convecting shallow water fluid
  50. Existence results for critical growth Kohn-Laplace equations with jumping nonlinearities
  51. Monotonicity and oscillation for fractional differential equations with Riemann-Liouville derivatives
  52. Nontrivial solutions for a generalized poly-Laplacian system on finite graphs
  53. Stability and bifurcation analysis of a modified chemostat model
  54. Special Issue on Nonlinear Evolution Equations and Their Applications - Part II
  55. Analytic solutions of a generalized complex multi-dimensional system with fractional order
  56. Extraction of soliton solutions and Painlevé test for fractional Peyrard-Bishop DNA model
  57. Special Issue on Recent Developments in Fixed-Point Theory and Applications - Part II
  58. Some fixed point results with the vector degree of nondensifiability in generalized Banach spaces and application on coupled Caputo fractional delay differential equations
  59. On the sum form functional equation related to diversity index
  60. Special Issue on International E-Conference on Mathematical and Statistical Sciences - Part II
  61. Simpson, midpoint, and trapezoid-type inequalities for multiplicatively s-convex functions
  62. Converses of nabla Pachpatte-type dynamic inequalities on arbitrary time scales
  63. Special Issue on Blow-up Phenomena in Nonlinear Equations of Mathematical Physics - Part II
  64. Energy decay of a coupled system involving a biharmonic Schrödinger equation with an internal fractional damping
  65. Special Issue on Some Integral Inequalities, Integral Equations, and Applications - Part II
  66. Nonlinear heat equation with viscoelastic term: Global existence and blowup in finite time
  67. New Jensen's bounds for HA-convex mappings with applications to Shannon entropy
  68. Special Issue on Approximation Theory and Special Functions 2024 conference
  69. Ulam-type stability for Caputo-type fractional delay differential equations
  70. Faster approximation to multivariate functions by combined Bernstein-Taylor operators
  71. (λ, ψ)-Bernstein-Kantorovich operators
  72. Some special functions and cylindrical diffusion equation on α-time scale
  73. (q, p)-Mixing Bloch maps
  74. Orthogonalizing q-Bernoulli polynomials
  75. On better approximation order for the max-product Meyer-König and Zeller operator
  76. Moment-based approximation for a renewal reward process with generalized gamma-distributed interference of chance
  77. A note on linear compositions of the Mellin convolution operators in the weighted Mellin-Lebesgue spaces
  78. Special Issue on Variational Methods and Nonlinear PDEs
  79. A note on mean field type equations
  80. Ground states for fractional Kirchhoff double-phase problem with logarithmic nonlinearity
  81. Solution of nonlinear Langevin equations involving Hilfer-Hadamard fractional order derivatives and variable coefficients
https://doi.org/10.1515/dema-2025-0160
Keywords for this article
chemostat model; stability; flip bifurcation; Neimark-Sacker bifurcation
Creative Commons
BY 4.0

Articles in the same Issue

  1. A new perspective on generalized Laguerre polynomials
  2. Research Articles
  3. On approximation by Stancu variant of Bernstein-Durrmeyer-type operators in movable compact disks
  4. Circular n,m-rung orthopair fuzzy sets and their applications in multicriteria decision-making
  5. Grand Triebel-Lizorkin-Morrey spaces
  6. Coefficient estimates and Fekete-Szegö problem for some classes of univalent functions generalized to a complex order
  7. Proofs of two conjectures involving sums of normalized Narayana numbers
  8. On the Laguerre polynomial approximation errors and lower type of entire functions of irregular growth
  9. New convolutions for the Hartley integral transform imbedded in the Banach algebras and convolution-type integral equations
  10. Some inequalities for rational function with prescribed poles and restricted zeros
  11. Lucas difference sequence spaces defined by Orlicz function in 2-normed spaces
  12. Evaluating the efficacy of fuzzy Bayesian networks for financial risk assessment
  13. Fixed point results for contractions of polynomial type
  14. Estimation for spatial semi-functional partial linear regression model with missing response at random
  15. Investigating the controllability of differential systems with nonlinear fractional delays, characterized by the order 0 < η ≤ 1 < ζ ≤ 2
  16. New forms of bilateral inequalities for K-g-frames
  17. Rate of pole detection using Padé approximants to polynomial expansions
  18. Existence results for nonhomogeneous Choquard equation involving p-biharmonic operator and critical growth
  19. Note on the shape-preservation of a new class of Kantorovich-type operators via divided differences
  20. Geršhgorin-type theorems for Z1-eigenvalues of tensors with applications
  21. New topologies derived from the old one via operators
  22. Blow up solutions for two-dimensional semilinear elliptic problem of Liouville type with nonlinear gradient terms
  23. Infinitely many normalized solutions for Schrödinger equations with local sublinear nonlinearity
  24. Nonparametric expectile shortfall regression for functional data
  25. Advancing analytical solutions: Novel wave insights and methodologies for beta fractional Kuralay-II equations
  26. A generalized p-Laplacian problem with parameters
  27. A study of solutions for several classes of systems of complex nonlinear partial differential difference equations in ℂ2
  28. Towards finding equalities involving mixed products of the Moore-Penrose and group inverses by matrix rank methodology
  29. ω ℒ -biprojective and ω ¯ -contractible Banach algebras
  30. Coefficient functionals for Sakaguchi-type-Starlike functions subordinated to the three-leaf function
  31. Solutions of several general quadratic partial differential-difference equations in ℂ2
  32. Inequalities for the generalized trigonometric functions with respect to weighted power mean
  33. Optimization of Lagrange problem with higher-order differential inclusion and special boundary-value conditions
  34. Hankel determinants for q-starlike functions connected with q-sine function
  35. System of partial differential hemivariational inequalities involving nonlocal boundary conditions
  36. A new family of multivalent functions defined by certain forms of the quantum integral operator
  37. A matrix approach to compare BLUEs under a linear regression model and its two competing restricted models with applications
  38. Weighted composition operators on bicomplex Lorentz spaces with their characterization and properties
  39. Behavior of spatial curves under different transformations in Euclidean 4-space
  40. Commutators for the maximal and sharp functions with weighted Lipschitz functions on weighted Morrey spaces
  41. A new kind of Durrmeyer-Stancu-type operators
  42. A study of generalized Mittag-Leffler-type function of arbitrary order
  43. On the approximation of Kantorovich-type Szàsz-Charlier operators
  44. Split quaternion Fourier transforms for two-dimensional real invariant field
  45. Review Article
  46. Characterization generalized derivations of tensor products of nonassociative algebras
  47. Special Issue on Differential Equations and Numerical Analysis - Part II
  48. Existence and optimal control of Hilfer fractional evolution equations
  49. Persistence of a unique periodic wave train in convecting shallow water fluid
  50. Existence results for critical growth Kohn-Laplace equations with jumping nonlinearities
  51. Monotonicity and oscillation for fractional differential equations with Riemann-Liouville derivatives
  52. Nontrivial solutions for a generalized poly-Laplacian system on finite graphs
  53. Stability and bifurcation analysis of a modified chemostat model
  54. Special Issue on Nonlinear Evolution Equations and Their Applications - Part II
  55. Analytic solutions of a generalized complex multi-dimensional system with fractional order
  56. Extraction of soliton solutions and Painlevé test for fractional Peyrard-Bishop DNA model
  57. Special Issue on Recent Developments in Fixed-Point Theory and Applications - Part II
  58. Some fixed point results with the vector degree of nondensifiability in generalized Banach spaces and application on coupled Caputo fractional delay differential equations
  59. On the sum form functional equation related to diversity index
  60. Special Issue on International E-Conference on Mathematical and Statistical Sciences - Part II
  61. Simpson, midpoint, and trapezoid-type inequalities for multiplicatively s-convex functions
  62. Converses of nabla Pachpatte-type dynamic inequalities on arbitrary time scales
  63. Special Issue on Blow-up Phenomena in Nonlinear Equations of Mathematical Physics - Part II
  64. Energy decay of a coupled system involving a biharmonic Schrödinger equation with an internal fractional damping
  65. Special Issue on Some Integral Inequalities, Integral Equations, and Applications - Part II
  66. Nonlinear heat equation with viscoelastic term: Global existence and blowup in finite time
  67. New Jensen's bounds for HA-convex mappings with applications to Shannon entropy
  68. Special Issue on Approximation Theory and Special Functions 2024 conference
  69. Ulam-type stability for Caputo-type fractional delay differential equations
  70. Faster approximation to multivariate functions by combined Bernstein-Taylor operators
  71. (λ, ψ)-Bernstein-Kantorovich operators
  72. Some special functions and cylindrical diffusion equation on α-time scale
  73. (q, p)-Mixing Bloch maps
  74. Orthogonalizing q-Bernoulli polynomials
  75. On better approximation order for the max-product Meyer-König and Zeller operator
  76. Moment-based approximation for a renewal reward process with generalized gamma-distributed interference of chance
  77. A note on linear compositions of the Mellin convolution operators in the weighted Mellin-Lebesgue spaces
  78. Special Issue on Variational Methods and Nonlinear PDEs
  79. A note on mean field type equations
  80. Ground states for fractional Kirchhoff double-phase problem with logarithmic nonlinearity
  81. Solution of nonlinear Langevin equations involving Hilfer-Hadamard fractional order derivatives and variable coefficients
Sign up now to receive a 20% welcome discount
Institutional Access
How does access work?
Have an idea on how to improve our website?
Please write us.
© 2025 De Gruyter Brill
Downloaded on 20.9.2025 from https://www.degruyterbrill.com/document/doi/10.1515/dema-2025-0160/html
Scroll to top button