Behavior of solutions of a discrete population model with mutualistic interaction

Sibi C. Babu; D. S. Dilip; Smitha Mary Mathew

doi:10.1515/cmb-2023-0121

Article Open Access

Behavior of solutions of a discrete population model with mutualistic interaction

Sibi C. Babu , D. S. Dilip and Smitha Mary Mathew

Published/Copyright: February 29, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Computational and Mathematical Biophysics Volume 12 Issue 1

Abstract

We focus on the stability analysis of two types of discrete dynamic models: a discrete dynamic equation and a discrete dynamics system consisting of two equations with mutualistic interaction given by

x n + 1 = a + b x n λ − ( x n − 1 + x n − k ) c + x n − 1 + x n − k

and

x n + 1 = a 1 + b 1 y n λ − ( y n − 1 + x n − k ) c 1 + y n − 1 + x n − k , y n + 1 = a 2 + b 2 x n λ − ( x n − 1 + y n − k ) c 2 + x n − 1 + y n − k ,

respectively, where k ∈ { 2 , 3 , … } , the constants a , a 1 , ≥ 0 , b , b 1 > 0 , and c , c 1 ≥ 0 be the initial densities, finite rate of increase and the limiting constant associated with the density of the species, respectively, and a 2 ≥ 0 , b 2 > 0 , and c 2 ≥ 0 be the initial densities, finite rate of increase, and the limiting constant associated with the density of the mutually interacting species, respectively. k , a positive integer represents a time delay in the system and λ ≥ 1 shows a decay factor based on the sum of two past time step population densities. Our main objective is to understand the impact of mutualistic interactions on the stability of discrete dynamic systems. To illustrate the boundedness and stability of these models, we also provide animated plots and bifurcation diagrams.

Keywords: boundedness; difference equations; persistence; asymptotic behavior

MSC 2010: 39A22

1 Introduction

Mathematical models of real-life problems often take the form of differential equations, difference equations, functional equations, and integro-difference equations (see [1–3,13,16–18,31,32,34,35]). Models using difference equations are widely used in various fields of science, particularly in the study of population dynamics. These models are typically represented as equations or systems and often include exponential terms. Examples of models with exponential terms can be found in previous studies [5,7–9,14,16,19,23–30,33].

In [12,21,22,28], authors worked on boundedness, persistence, and global stability properties of some difference equations which contain exponential terms. The periodic and chaotic nature of difference equations are discussed in [11,15,20].

The Ricker model given by N t + 1 = N t e r 1 − N t k has been widely utilized to study the dynamics of discrete populations, but its applicability is often constrained by its simplicity. To address this limitation, researchers have proposed modifications and extensions to the traditional Ricker equation to better capture the nuances of real-world ecological systems. For example, in [4,6,10], authors modified the equations defined in [12,21,22] by generalizing the exponent base, replacing e with an arbitrary constant, and studied the boundedness, persistence, and stability properties. In this article, we analyze the proposed generalized model, exploring its behavior under different parameter regimes. Through this approach, we aim to contribute to the ongoing discourse surrounding population dynamics, offering a flexible model that can be adapted to a wide range of ecological scenarios. We focus on the stability analysis of two types of discrete dynamic models: a discrete dynamic equation and a discrete dynamics system consisting of two equations with mutualistic interaction.

(1) x n + 1 = a + b x n λ − ( x n − 1 + x n − k ) c + x n − 1 + x n − k

and

(2) x n + 1 = a 1 + b 1 y n λ − ( y n − 1 + x n − k ) c 1 + y n − 1 + x n − k , y n + 1 = a 2 + b 2 x n λ − ( x n − 1 + y n − k ) c 2 + x n − 1 + y n − k ,

In the first model, the density of the species depends on the previous, second previous, and ( k + 1 ) th previous generations’ densities. In the second model, the density of the species depends on its ( k + 1 ) th previous generation’s density and is influenced by the previous and second previous generations’ densities of another species. This mutual influence is assumed to work in both directions, with the second species also depending on its ( k + 1 ) th previous generation’s density and the previous and second previous generations’ densities of the first species. Mutualistic interactions, where two species benefit from their interaction, is an important aspect of many dynamic systems, particularly in biology. Understanding the stability of such systems is crucial for predicting their behavior and controlling their outcomes.

Equation (1) represents the density of a species at the ( n + 1 ) th stage, denoted as x n + 1 , which depends on its previous, second previous, and ( k + 1 ) th previous generations’ densities. Similarly, (2) represents the density of a species at the ( n + 1 ) th stage, denoted as x n + 1 , which depends on its ( k + 1 ) th previous generation’s density and is influenced by previous and second previous generations’ densities of another species’, say y n + 1 . Also, that y n + 1 mutually depends on its ( k + 1 ) th previous generation’s density, and the previous, second previous generations of x n + 1 .

We observe the conditions under which the density of species is bounded, persistent, convergent, and asymptotic. Moreover, using python, we animate both plots and bifurcation diagrams of (1) and (2) to illustrate the boundedness and stability results.

2 Boundedness and persistence

Here, we observe how the upper bound varies as we study the boundedness, persistence, and invariants of (1) and (2).

Theorem 2.1

Suppose

(3) p = b λ − 2 a c + 2 a < 1 .

Then the positive solution { x n } of (1)is bounded and persists.
I = a , a 1 − p is an invariant set of (1).
Let ε > 0 be an arbitrary number and { x n } be an arbitrary solution of (1). Consider the interval J = a , a + ε 1 − p . Then ∃ n 0 ∈ N such that x n ∈ J for all n ≥ n 0 .

Proof

From (1), we have
(4) x n ≥ a for n = 1 , 2 , … .
Furthermore,
(5) x n + 1 ≤ a + b x n λ − 2 a c + 2 a = a + p x n for n > k .
Consider the non-homogeneous equation
(6) u n + 1 = a + p u n , ∀ n > k .
An arbitrary solution { u n } of (6) for n = k + 2 , k + 3 , … is given by
(7) u n = u k + 1 p n + a 1 − p .
Thus, the relations (3), (6), and (7) give that { u n } is bounded. Let the solution { u n } be such that u 1 = x 1 , u 2 = x 2 , … , u k + 1 = x k + 1 . Thus, from (5) and (6) we obtain
(8) x n ≤ u n , n = 1 , 2 , … .
Then from (4) and (8) we find that solution of (1) is bounded.
Suppose x n > 0 be a solution of (1). Let x − k , x − k + 1 , … , x 0 ∈ a , a 1 − p .

Then from (1) we obtain
a ≤ x 1 = a + b x 0 λ − ( x − 1 + x − k ) c + ( x − 1 + x − k ) ≤ a + b x 0 λ − 2 a c + 2 a = a + p x 0 ≤ a + a p 1 − p = a 1 − p .
Therefore, x 1 ∈ a , a 1 − p .
Working inductively we have x n ∈ a , a 1 − p for n = 1 , 2 , 3 , … .
Let L = limsup x n and l = liminf x n . Then,
L ≤ a + b L λ − 2 l c + 2 l . = a + p L
Thus, a ≤ L ≤ a 1 − p .
Then for all ε > 0 , we can pick out an n 0 ∈ N such that
x n < L + ε = a 1 − p + ε < a + ε 1 − p .
Therefore, x n ∈ a , a + ε 1 − p ∀ n ≥ n 0 .□

For λ = 3 , a = 2 , b = 400 , and for the initial values from I , Animation 1 shows the plot frames of (1) when k = 2 and k = 5 , and for various values of c . We see that the solution is bounded as stated in Theorem 2.1.

In the next theorem, we examine how the upper bound varies due to the influence of { y n } .

Theorem 2.2

Consider the difference system (2) such that

p 1 = b 1 λ − ( a 1 + a 2 ) c 1 + a 1 + a 2 and p 2 = b 2 λ − ( a 1 + a 2 ) c 2 + a 1 + a 2 .

Let

(9) p 1 p 2 < 1 .

Then the positive solution { ( x n , y n ) } of (2)is bounded and persists.
a 1 , P 1 − p 1 p 2 × a 2 , Q 1 − p 1 p 2 is an invariant set for (2), where P = a 1 + p 1 a 2 and Q = a 2 + p 2 a 1 .
Let ε > 0 be arbitrary and consider the intervals
(10) J 1 = a 1 , P + ε 1 − p 1 p 2 and J 2 = a 2 , Q + ε 1 − p 1 p 2 .
Then ∃ an n 1 ∈ N such that x n and y n will be in J 1 and J 2 , respectively, for all n ≥ n 1 .

Proof

From (2), we see that, for n = 1 , 2 , … ,
(11) x n ≥ a 1 , y n ≥ a 2 .
So every solution of (2) persists.

Also (2) implies
(12) x n + 1 ≤ a 1 + p 1 y n ,

(13) y n + 1 ≤ a 2 + p 2 x n for n > k .
Using (12) and (13) and for n > k ,
(14) x n + 1 ≤ P + p 1 p 2 x n − 1 , y n + 1 ≤ Q + p 1 p 2 y n − 1 ,
where P = a 1 + p 1 a 2 and Q = a 2 + p 2 a 1 .

Consider the system
(15) u n + 1 = P + p 1 p 2 u n − 1 , v n + 1 = Q + p 1 p 2 v n − 1 , for n = k + 1 , k + 2 , …

Let { ( u n , v n ) } be a solution of (15) such that
(16) u 1 = x 1 , u 2 = x 2 , … , u k + 1 = x k + 1 and v 1 = y 1 , v 2 = y 2 , … , v k + 1 = y k + 1 .
Using (15) and (16), we obtain u k + 2 > 0 , v k + 2 > 0 .

Working inductively it follows that u n > 0 , v n > 0 for n = k + 2 , k + 3 , … .

Moreover, from (15) for n = k + 2 , k + 3 , …
(17) u n = r 1 ( p 1 p 2 ) n ∕ 2 + r 2 ( − 1 ) n ( p 1 p 2 ) n ∕ 2 + P 1 − p 1 p 2 , v n = s 1 ( p 1 p 2 ) n ∕ 2 + s 2 ( − 1 ) n ( p 1 p 2 ) n ∕ 2 + Q 1 − p 1 p 2 ,
where r 1 , r 2 are constants defined by x k , x k + 1 and s 1 , s 2 are constants defined by y k , y k + 1 .

Using (14) and (17) we can prove by induction that
(18) x n ≤ u n , y n ≤ v n for n ≥ 1 .
Then from (11), (17), and (18), we see that the solution of (2) is bounded.
Let
(19) I 1 = a 1 , P 1 − p 1 p 2 and I 2 = a 2 , Q 1 − p 1 p 2 .
Suppose ( x n , y n ) > 0 be a solution of (2) and
(20) x − k , x − k + 1 , … , x 0 ∈ I 1 , y − k , y − k + 1 , … , y 0 ∈ I 2 .
Then from (2) and (20) we have
a 1 ≤ x 1 ≤ a 1 + b 1 λ − ( y − 1 + x − k ) c 1 + y − 1 + x − k y 0 ≤ a 1 + b 1 λ − ( a 1 + a 2 ) c 1 + a 1 + a 2 Q 1 − p 1 p 2 = a 1 + p 1 Q 1 − p 1 p 2 = a 1 + p 1 ( a 2 + p 2 a 1 ) 1 − p 1 p 2 = P 1 − p 1 p 2 .
Thus, x 1 ∈ I 1 . Similarly, we can show that y 1 ∈ I 2 . Working inductively, we obtain, x n ∈ I 1 , y n ∈ I 2 for n = 1 , 2 , 3 , …
From Theorem 2.2(a), we obtain
(21) 0 < l 1 = liminf n → ∞ x n , 0 < l 2 = liminf n → ∞ y n , L 1 = limsup n → ∞ x n < ∞ , L 2 = limsup n → ∞ y n < ∞ .
It follows from (2) and (21)
L 1 ≤ a 1 + b 1 L 2 λ − ( l 1 + l 2 ) c 1 + l 1 + l 2 , l 1 ≥ a 1 + b 1 l 2 λ − ( L 1 + L 2 ) c 1 + L 1 + L 2 , L 2 ≤ a 2 + b 2 L 1 λ − ( l 1 + l 2 ) c 2 + l 1 + l 2 , l 2 ≥ a 2 + b 2 l 1 λ − ( L 1 + L 2 ) c 2 + L 1 + L 2 .
Thus,
L 1 ≤ a 1 + p 1 a 2 + p 1 p 2 L 1 .
This implies that
a 1 ≤ L 1 ≤ P 1 − p 1 p 2 .
Similarly, we have
a 2 ≤ L 2 ≤ Q 1 − p 1 p 2 .
Thus, ∃ an n 1 ∈ N such that x n ∈ J 1 and y n ∈ J 2 for all n ≥ n 1 .□

For λ = 3 , a 1 = 3.2 , a 2 = 4.4 , b 1 = 20000.4, b 2 = 40000.6 , c 2 = 10.1 , we see from Animation 2 that the plot frames of (2) when k = 2 and k = 5 , and for various values of c 1 , the solution is bounded as stated in Theorem 2.2.

3 Global behavior

We provide conditions for the solutions of (1) and (2) to be locally and globally asymptotically stable. The equilibrium points of (1) are the roots of

(22) x ¯ = a + b x ¯ λ − 2 x ¯ c + 2 x ¯ .

Set g ( x ) = a + b x λ − 2 x c + 2 x − x .

Then, we obtain

(23) g ′ ( x ) = b λ 2 x ( c + 2 x ) − 2 b x ln λ λ 2 x ( c + 2 x ) − 2 b x λ 2 x ( c + 2 x ) 2 − 1 ,

Also,

(24) lim x → ∞ g ( x ) = − ∞ , g ( 0 ) = a > 0 .

Assume z be a root of g ( x ) = 0 . Then

g ( z ) = 0 ⇒ a + b z λ − 2 z c + 2 z − z = 0 ⇒ z ≥ a and z − a = b z λ − 2 z c + 2 z .

Then (23) becomes

g ′ ( z ) = − a z − 2 ( z − a ) ln λ − 2 ( z − a ) c + 2 z < 0 .

Therefore, ∃ an ε > 0 such that g ′ ( x ) < 0 , ∀ x ∈ ( z − ε , z + ε ) . Hence, g is decreasing here. Suppose g has solutions greater than the z and let z 1 be the smallest solution of g such that z 1 > z . There exists an ε 1 such that g is decreasing in ( z 1 − ε 1 , z 1 + ε 1 ) . Since g ( z + ε ) < 0 , g ( z 1 − ε 1 ) > 0 , and g is continuous, there exists a root in ( z + ε , z 1 − ε 1 ) , which contradicts our assumption. Similarly, g has no solution in ( a , z ) . Therefore, the solution of g ( x ) = 0 is unique.

The following lemma describes the convergence of the positive solutions of (1).

Lemma 3.1

Suppose p < 1 and

(25) b < ( 2 a + c ) ( 1 − p ) ( 1 − p + 2 a ln λ ) .

Then, as n → ∞ , the solution { x n } of (1) approaches the unique equilibrium.

Proof

Let f : I 3 → I be a continuous function defined by

f ( x , y , z ) = a + b x λ − ( y + z ) c + y + z ,

where I = a , a 1 − p .

Let m , M ∈ I and

(26) m = a + b m λ − 2 M c + 2 M ,

(27) M = a + b M λ − 2 m c + 2 m .

From (26) and (27), we have

(28) ( 2 a + c ) ( M − m ) = b M λ − 2 m − b m λ − 2 M = b λ − 2 ( M + m ) [ M λ 2 M − m λ 2 m ] .

Also, there exists a ξ , m ≤ ξ ≤ M such that

(29) M λ 2 M − m λ 2 m = λ 2 ξ [ 1 + 2 ξ ln λ ] [ M − m ] .

Hence, (28) implies

(30) ( 2 a + c ) ∣ M − m ∣ = b λ − 2 ( M + m − ξ ) [ 1 + 2 ξ ln λ ] ∣ M − m ∣ ≤ b 1 + 2 a ln λ 1 − p ∣ M − m ∣ ∣ M − m ∣ ≤ b ( 2 a + c ) 1 + 2 a ln λ 1 − p ∣ M − m ∣ .

Therefore from (25) and (30), we obtain M = m .

From Theorem 1.15 of [15], we obtain { x n } that converges to x ¯ .□

Lemma 3.3 shows that, for system (2), a unique equilibrium solution exists and it is a global attractor. The following theorem is a modification of Theorem 1.16 of [15].

Theorem 3.2

Let f : [ a , b ] × [ c , d ] × [ c , d ] → [ a , b ] and g : [ a , b ] × [ a , b ] × [ c , d ] → [ c , d ] be continuous real-valued functions. Consider the difference equations

(31) x n + 1 = f ( x n − k , y n , y n − 1 ) , y n + 1 = g ( x n , x n − 1 , y n − k ) , n = 0 , 1 , …

with initial condition ( x − k , y − k ) , ( x − k + 1 , y − k + 1 ) , … , ( x 0 , y 0 ) ∈ [ a , b ] × [ c , d ] . Suppose the following statements are true:

f ( u , v , w ) is non-increasing in u , w , and is non-decreasing in v .
g ( u , v , w ) is non-decreasing in u , and is non-increasing in v and w .
If ( m , M , r , R ) ∈ [ a , b ] 2 × [ c , d ] 2 is a solution of the system of equations
m = f ( M , r , R ) , r = g ( m , M , R ) , M = f ( m , R , r ) , R = g ( M , m , r ) ,
then m = M and r = R .
Then, there is only one equilibrium point ( x ¯ , y ¯ ) for system (31), and every solution of system (31) converges to ( x ¯ , y ¯ ) as n approaches ∞ .

Proof

The proof is similar to Theorem 1.16 of [15].□

Lemma 3.3

Consider system (2) such that relation (9) holds true. Suppose that

(32) [ ( a 1 + c 1 ) ( 1 − p 1 p 2 ) + Q ] [ ( a 2 + c 2 ) ( 1 − p 1 p 2 ) + P ] b 1 b 2 ( 1 − p 1 p 2 ) 2 ( a 1 ln λ + 1 ) ( a 2 ln λ + 1 ) λ − 2 ( P + Q 1 − p 1 p 2 ) < 1 .

Then the positive solution of system (2) converges to the equilibrium point ( x ¯ , y ¯ ) .

Proof

Let f : a 1 , P 1 − p 1 p 2 × a 2 , Q 1 − p 1 p 2 × a 2 , Q 1 − p 1 p 2 → a 1 , P 1 − p 1 p 2 and g : a 1 , P 1 − p 1 p 2 × a 1 , P 1 − p 1 p 2 × a 2 , Q 1 − p 1 p 2 → a 2 , Q 1 − p 1 p 2 be two continuous functions such that

f ( x , y , z ) = a 1 + b 1 y λ − ( x + z ) c 1 + x + z and g ( x , y , z ) = a 2 + b 2 x λ − ( y + z ) c 2 + y + z .

Let m , M ∈ a 1 , P 1 − p 1 p 2 and r , R ∈ a 2 , Q 1 − p 1 p 2 such that

(33) M = a 1 + b 1 R λ − ( m + r ) c 1 + m + r , m = a 1 + b 1 r λ − ( M + R ) c 1 + M + R ,

(34) R = a 2 + b 2 M λ − ( m + r ) c 2 + m + r , r = a 2 + b 2 m λ − ( M + R ) c 2 + M + R ,

(35) m ≤ M and r ≤ R .

From (33) we obtain

(36) M c 1 + M m + M r − a 1 c 1 − a 1 m − a 1 r = b 1 R λ m + r .

(37) m c 1 + M m + m R − a 1 c 1 − a 1 M − a 1 R = b 1 r λ M + R .

Subtracting (37) from (36) and using (35) we have

(38) ( a 1 + c 1 + r ) [ M − m ] − ( m − a 1 ) [ R − r ] ≥ b 1 λ − m λ R + r [ R λ R − r λ r ] .

Furthermore, for some θ 1 , r < θ 1 < R we obtain

(39) R λ R − r λ r = λ θ 1 ( θ 1 ln λ + 1 ) [ R − r ] .

Substituting (39) in (38) and since m ≥ a 1 we have

( a 1 + c 1 + r ) [ M − m ] ≥ b 1 λ − m λ R + r λ θ 1 ( θ 1 ln λ + 1 ) [ R − r ] . ≥ b 1 λ − ( m + R ) ( θ 1 ln λ + 1 ) [ R − r ] .

This implies

a 1 + c 1 + Q 1 − p 1 p 2 [ M − m ] ≥ b 1 λ − ( m + R ) ( a 2 ln λ + 1 ) [ R − r ] .

Thus,

(40) [ R − r ] ≤ [ ( a 1 + c 1 ) ( 1 − p 1 p 2 ) + Q ] b 1 ( 1 − p 1 p 2 ) ( a 2 ln λ + 1 ) λ − ( P + Q 1 − p 1 p 2 ) ≤ [ M − m ] .

From (34) we obtain

(41) R c 2 + R m + R r − a 2 c 2 − a 2 m − a 2 r = b 2 M λ m + r .

(42) r c 2 + r M + r R − a 2 c 2 − a 2 M − a 2 R = b 2 m λ M + R .

Subtracting (42) from (41) we obtain

( a 2 + c 2 ) [ R − r ] + [ R m − r M ] + a 2 [ M − m ] = b 2 M λ − r λ m − b 2 m λ − R λ M .

Using (35) we have

(43) ( a 2 + c 2 + m ) [ R − r ] − ( r − a 2 ) [ M − m ] ≥ b 2 λ − r λ M + m [ M λ M − m λ m ] .

For some θ 2 , m < θ 2 < M , we have

(44) M λ M − m λ m = λ θ 2 ( θ 2 ln λ + 1 ) [ M − m ] .

Substituting (44) in (43) and since r ≥ a 2 we have

( a 2 + c 2 + m ) [ R − r ] ≥ b 2 λ − r λ M + m λ θ 2 ( θ 2 ln λ + 1 ) [ M − m ] ≥ b 2 λ − ( M + r ) ( θ 2 ln λ + 1 ) [ M − m ] .

This implies

a 2 + c 2 + P 1 − p 1 p 2 [ R − r ] ≥ b 2 λ − ( M + r ) ( a 1 ln λ + 1 ) [ M − m ] ,

Thus,

(45) [ M − m ] ≤ [ ( a 2 + c 2 ) ( 1 − p 1 p 2 ) + P ] b 2 ( 1 − p 1 p 2 ) ( a 1 ln λ + 1 ) λ − ( P + Q 1 − p 1 p 2 ) ≤ [ R − r ] .

Then (40) and (45) gives

(46) [ M − m ] ≤ [ ( a 1 + c 1 ) ( 1 − p 1 p 2 ) + Q ] [ ( a 2 + c 2 ) ( 1 − p 1 p 2 ) + P ] b 1 b 2 ( 1 − p 1 p 2 ) 2 ( a 1 ln λ + 1 ) ( a 2 ln λ + 1 ) λ − 2 ( P + Q 1 − p 1 p 2 ) [ M − m ]

and

(47) [ R − r ] ≤ [ ( a 1 + c 1 ) ( 1 − p 1 p 2 ) + Q ] [ ( a 2 + c 2 ) ( 1 − p 1 p 2 ) + P ] b 1 b 2 ( 1 − p 1 p 2 ) 2 ( a 1 ln λ + 1 ) ( a 2 ln λ + 1 ) λ − 2 P + Q 1 − p 1 p 2 [ R − r ] .

Therefore, from (32), (46), and (47), we have M = m and R = r . Consequently, every solution of (2) tends to a unique positive equilibrium ( x ¯ , y ¯ ) .□

In the following theorems, we derive the globally asymptotically stability conditions for the equilibrium solutions of (1) and (2).

Theorem 3.4

Let (25) hold. Suppose that

(48) b ≤ 1 2 ( c + 2 a ) λ 2 a ,

and

(49) b < ( c + 2 a ) λ 2 a 1 − 2 a ln λ + 1 c + 2 a .

Then x ¯ of (1) is globally asymptotically stable.

Proof

(50) μ k + 1 − α μ k − β μ k − 1 − β = 0

is the characteristic equation corresponding to (1), where α = b λ − 2 x ¯ c + 2 x ¯ and β = − b x ¯ λ − 2 x ¯ [ 1 + ( c + 2 x ¯ ) ln λ ] ( c + 2 x ¯ ) 2 .

Since x ¯ satisfies (22) we have

(51) a < x ¯ ≤ a 1 − p .

Then by using (48), (49), and (51) we obtain

(52) ∣ α ∣ + 2 ∣ β ∣ = b λ − 2 x ¯ c + 2 x ¯ + 2 b x ¯ λ − 2 x ¯ [ 1 + ( c + 2 x ¯ ) ln λ ] ( c + 2 x ¯ ) 2 ≤ b λ − 2 a c + 2 a + 2 a b λ − 2 a ( c + 2 a ) ( 1 − p ) ln λ + 1 c + 2 a ≤ p + 2 a ln λ + 1 c + 2 a . Since, 0 ≤ p ≤ 0.5 , we have p 1 − p ≤ 1 < 1 .

Therefore, from (52) and from Theorem 1.6 of [15], the equilibrium point x ¯ is locally asymptotically stable. By Lemma 3.1, we conclude that x ¯ is asymptotically stable globally.□

To illustrate the global asymptotic stability stated in Theorem 3.1, we take λ = 3 , a = 0.03 and we animate bifurcation diagrams of (1) for various values of c between 1 and 16 with bifurcation parameter b . The black and red vertical lines in the animation correspond to the points b = 1 2 ( c + 2 a ) λ 2 a and b = ( c + 2 a ) λ 2 a 1 − 2 a ln λ + 1 c + 2 a , respectively. We can see in Animation 3, that for values of b between 0 and the vertical line (black or red) near 0, (1) is globally asymptotically stable.

Theorem 3.5

Suppose conditions (9) and (32) hold. Also suppose that

(53) b 1 ≤ ( c 1 + a 1 + a 2 ) ( c 2 + a 1 + a 2 ) λ 2 ( a 1 + a 2 ) 2 b 2

and

(54) b 1 < ( c 1 + a 1 + a 2 ) λ ( a 1 + a 2 ) p 2 + a 2 ln λ + 1 c 2 + 4 a 2 2 ln λ + 1 c 1 ln λ + 1 c 2 + 2 a 2 ln λ + 1 c 1 − 1 × 1 − ( a 1 + a 2 + a 1 p 2 ) ln λ + 1 c 1 − ( a 1 + a 2 + 2 a 1 p 2 ) ln λ + 1 c 2 − ( 6 a 1 a 2 + 4 p 2 a 1 2 ) ln λ + 1 c 1 ln λ + 1 c 2 ,

then ( x ¯ , y ¯ ) of (2) is globally asymptotically stable.

Proof

The linearized system of (2) is

(55) x n + 1 = β 1 y n + α 1 y n − 1 + α 1 x n − k , y n + 1 = α 2 x n + β 2 x n − 1 + β 2 y n − k ,

where

α 1 = − b 1 y ¯ [ 1 + ( c 1 + x ¯ + y ¯ ) ln λ ] ( c 1 + x ¯ + y ¯ ) 2 λ x ¯ + y ¯ , α 2 = b 2 λ − ( x ¯ + y ¯ ) c 2 + x ¯ + y ¯ , β 1 = b 1 λ − ( x ¯ + y ¯ ) c 1 + x ¯ + y ¯ , β 2 = − b 2 x ¯ [ 1 + ( c 2 + x ¯ + y ¯ ) ln λ ] ( c 2 + x ¯ + y ¯ ) 2 λ x ¯ + y ¯ .

Note that system (55) is equivalent to W n + 1 = C W n , where

C = 0 β 1 0 α 1 … 0 0 α 1 0 α 2 0 β 2 0 … 0 0 0 β 2 1 0 0 0 … 0 0 0 0 0 1 0 0 … 0 0 0 0 … … … 0 0 … … … 1 0 0 0 0 0 0 0 … 0 1 0 0 , W n = x n y n x n − 1 y n − 1 ⋮ x n − k y n − k .

Then for k > 1 , the characteristic equation of matrix C is

(56) μ 2 k + 2 − α 2 β 1 μ 2 k − ( α 1 α 2 + β 1 β 2 ) μ 2 k − 1 − α 1 β 2 μ 2 k − 2 − ( α 1 + β 2 ) μ k + 1 + α 1 β 2 = 0 .

Also,

(57) x ¯ = a 1 + b 1 y ¯ λ − ( x ¯ + y ¯ ) c 1 + x ¯ + y ¯ and y ¯ = a 2 + b 2 x ¯ λ − ( x ¯ + y ¯ ) c 2 + x ¯ + y ¯ .

Then,

x ¯ ≤ a 1 + a 2 p 1 + p 1 p 2 x ¯ = P + p 1 p 2 x ¯ .

Thus,

(58) a 1 ≤ x ¯ ≤ P 1 − p 1 p 2 .

Similarly from (57) we have

(59) a 2 ≤ y ¯ ≤ Q 1 − p 1 p 2 ,

b 1 ≤ λ 2 ( a 1 + a 2 ) ( c 1 + a 1 + a 2 ) ( c 2 + a 1 + a 2 ) 2 b 2 .

From (53) we obtain,

(60) 0 < p 1 p 2 ≤ 1 2 ⇒ 1 − p 1 p 2 ≥ 1 2 ⇒ 1 1 − p 1 p 2 ≤ 2 and p 1 p 2 1 − p 1 p 2 ≤ 1 .

Then by using (54), (58), (59), and (60) we obtain

∣ α 2 β 1 ∣ + ∣ α 1 α 2 ∣ + ∣ β 1 β 2 ∣ + 2 ∣ α 1 β 2 ∣ + ∣ α 1 ∣ + ∣ β 2 ∣ ≤ p 1 p 2 + Q p 1 p 2 ( c 1 ln λ + 1 ) c 1 ( 1 − p 1 p 2 ) + P p 1 p 2 ( c 2 ln λ + 1 ) c 2 ( 1 − p 1 p 2 ) + 2 P Q p 1 p 2 ( c 1 ln λ + 1 ) ( c 2 ln λ + 1 ) c 1 c 2 ( 1 − p 1 p 2 ) 2 + Q p 1 ( c 1 ln λ + 1 ) c 1 ( 1 − p 1 p 2 ) + P p 2 ( c 2 ln λ + 1 ) c 2 ( 1 − p 1 p 2 ) ≤ [ a 1 + a 2 + a 1 p 2 ] ln λ + 1 c 1 + [ a 1 + a 2 + 2 a 1 p 2 ] ln λ + 1 c 2 + p 1 p 2 + a 2 ( ln λ + 1 c 2 ) + 4 a 2 2 ln λ + 1 c 1 ln λ + 1 c 2 + 2 a 2 ln λ + 1 c 1 + ( 6 a 1 a 2 + 4 p 2 a 1 2 ) ln λ + 1 c 1 ln λ + 1 c 2 < 1 .

By Theorem 1.6 of [15], all the roots of (56) satisfy ∣ μ ∣ < 1 . Hence by the linearized stability theorem and Lemma 3.3, the result follows.□

To illustrate the global asymptotic stability stated in Theorem 3.3, we take λ = 15 , a 1 = 0.0002 , a 2 = 0.0004 , b 2 = 3000.6 , c 2 = 10.1 and we animate bifurcation diagrams of (2) for various values of c between 25 and 80 with bifurcation parameter b 1 . Black and red vertical lines in the animation correspond to the points b 1 = λ 2 ( a 1 + a 2 ) ( c 1 + a 1 + a 2 ) ( c 2 + a 1 + a 2 ) 2 b 2 and

b 1 = ( c 1 + a 1 + a 2 ) λ ( a 1 + a 2 ) p 2 + a 2 ln λ + 1 c 2 + 4 a 2 2 ln λ + 1 c 1 ln λ + 1 c 2 + 2 a 2 ln λ + 1 c 1 × 1 − ( a 1 + a 2 + a 1 p 2 ) ln λ + 1 c 1 − ( a 1 + a 2 + 2 a 1 p 2 ) ln λ + 1 c 2 − ( 6 a 1 a 2 + 4 p 2 a 1 2 ) ln λ + 1 c 1 ln λ + 1 c 2 ,

respectively. We can see in Animation 4, that for values of b 1 between 0 and the vertical line (black or red) near 0, (2) is globally asymptotically stable.

4 Numerical illustrations

Here we give few numerical illustrations and figures to verify our results.

Example 4.1

Let a = 0.03 , b = 2 , c = 2 , λ = 3 , and k = 5 . For initial x values 5.0 , 3.0 , 3.3 , 2.2 , 2.3 we see that after few iterations x n converges to 0.10586093950821493, as shown in Figure 1(a).

When c is changed to 5, we can see that after few iterations x n converges to 0.046480852055143476, as shown in Figure 1(b).

$Figure 1 Plots of x n {x}_{n} of (1) showing the convergence of positive solutions to the equilibrium where a = 0.03 , b = 2 , λ = 3 a=0.03,b=2,\lambda =3 , and k = 5 . k=5. (a) c = 2 c=2 and (b) c = 5 c=5 .$

Figure 1

Plots of x n of (1) showing the convergence of positive solutions to the equilibrium where a = 0.03 , b = 2 , λ = 3 , and k = 5 . (a) c = 2 and (b) c = 5 .

Example 4.2

Let a = 0.03 , b = 9 , c = 7 , λ = 3 , and k = 6 . For initial x values 62 , 11 , 94 , 2.7 , 18.3 , 19.2 we see that after few iterations x n converges to 0.1766657815923124, as shown in Figure 2(a).

When b is changed to 2.5 we can see that after few iterations x n converges to 0.04412460390399768, as shown in Figure 2(b).

$Figure 2 Plots of x n {x}_{n} of (1) showing the convergence of positive solutions to the equilibrium where a = 0.03 , c = 7 , λ = 3 a=0.03,c=7,\lambda =3 , and k = 6 . k=6. (a) b = 9 b=9 and (b) b = 2.5 b=2.5 .$

Figure 2

Plots of x n of (1) showing the convergence of positive solutions to the equilibrium where a = 0.03 , c = 7 , λ = 3 , and k = 6 . (a) b = 9 and (b) b = 2.5 .

Example 4.3

Let a 1 = 0.0002 , a 2 = 0.0004 , b 1 = 0.03 , b 2 = 3000.6 , c 1 = 25 , c 2 = 10.10 , λ = 15 , and k = 5 . For initial x values 6.2 , 10 , 74 , 2.7 , 28.3 and y values 8 . 2 , 11 , 94 , 6.7 , 18.3 we see that after few iterations x n converges to 0.0002660289968728449 and y n converges to 0.06601746237737958, as shown in Figure 3(a).

When c 1 is changed to 35, after few iterations x n converges to 0.0002445102434840107 and y n converges to 0.06148541027270541, as shown in Figure 3(b).

$Figure 3 Plots of x n , y n {x}_{n},{y}_{n} of (2) showing the convergence of positive solutions to the equilibrium where a 1 = 0.0002 , a 2 = 0.0004 , b 1 = 0.03 , b 2 = 3000.6 , c 2 = 10.10 , λ = 15 {a}_{1}=0.0002,{a}_{2}=0.0004,{b}_{1}=0.03,{b}_{2}=3000.6,{c}_{2}=10.10,\lambda =15 , and k = 5 . k=5. (a) c 1 = 25 {c}_{1}=25 and (b) c 1 = 35 {c}_{1}=35 .$

Figure 3

Plots of x n , y n of (2) showing the convergence of positive solutions to the equilibrium where a 1 = 0.0002 , a 2 = 0.0004 , b 1 = 0.03 , b 2 = 3000.6 , c 2 = 10.10 , λ = 15 , and k = 5 . (a) c 1 = 25 and (b) c 1 = 35 .

Example 4.4

Let a 1 = 0.0002 , a 2 = 0.0004 , b 1 = 0.25 , b 2 = 3000.6 , c 1 = 53 , c 2 = 10.10 , λ = 15 , and k = 5 . For initial x values 16 , 31.3 , 17 , 2.07 , 2.8 and y values 1.8 , 1.01 , 4.9 , 6.1 , 1.2 , we see that after few iterations x n converges to 0.0006294612248094049 and y n converges to 0.13000822040427737, as shown in Figure 4(a).

When b 1 is changed to 0.07, after few iterations x n converges to 0.0002743178960513349 and y n converges to 0.06773503675038935, as shown in Figure 4(b).

$Figure 4 Plots of x n , y n {x}_{n},{y}_{n} of (2) showing the convergence of positive solutions to the equilibrium where a 1 = 0.0002 , a 2 = 0.0004 , b 2 = 3000.6 , c 1 = 53 , c 2 = 10.10 , λ = 15 {a}_{1}=0.0002,{a}_{2}=0.0004,{b}_{2}=3000.6,{c}_{1}=53,{c}_{2}=10.10,\lambda =15 , and k = 5 . k=5. (a) b 1 = 0.25 {b}_{1}=0.25 and (b) b 1 = 0.07 {b}_{1}=0.07 .$

Figure 4

Plots of x n , y n of (2) showing the convergence of positive solutions to the equilibrium where a 1 = 0.0002 , a 2 = 0.0004 , b 2 = 3000.6 , c 1 = 53 , c 2 = 10.10 , λ = 15 , and k = 5 . (a) b 1 = 0.25 and (b) b 1 = 0.07 .

5 Conclusion

We derived the conditions for the solutions of (1) and (2) to be bounded, persist, invariant, convergent, and asymptotic in order to compare the results between the equation and system. We found that the lower bound of the invariant intervals (1) and (2) remains constant, while the upper bound varies based on the parameters of the equation and system, respectively. We also determined the conditions for global asymptotic stability for equation (1) and system (2) using the convergence conditions derived for each separately. To further visualize the results, we animated both the plots and bifurcation diagrams of (1) and (2) and demonstrated the boundedness and stability. As we move forward, further exploration and refinement of these models can contribute to a more comprehensive understanding of the dynamics of ecological systems. The insights gained from this research can potentially guide ecological management practices and aid in predicting the behavior of populations in real-world scenarios with mutualistic interactions.

Author contributions: Sibi C. Babu contributed to the main results, S. M. Mathew did the analysis part, and Dr. D. S. Dilip supervised them.
Conflict of interest: The authors state no conflict of interest.
Ethical approval: The conducted research is not related to either human or animals use.
Informed consent: Informed consent has been obtained from all individuals included in this study.
Data availability statement: Data availability is not applicable to this article as no new data were created or analyzed in this study.

Appendix

The animations mentioned in this article can be accessed from the following links:

Animation 1: https://youtu.be/fPPsg6kR67U
Animation 2: https://youtu.be/6YPTKcr1NCg
Animation 3: https://youtube.com/shorts/zH4s7vgsKZ4
Animation 4: https://youtu.be/J3QdDSsfpkc

References

[1] Beverton R. J. H., & Holt S. J. (1957). On the dynamics of exploited fish population, Fishery Investigations (vol. 2, Issue 19), H.M. Stationery Off, London. Search in Google Scholar

[2] Carrillo C., & Fife P. (2005). Spatial effects in discrete generation population models, Journal of Mathematical Biology, 50, 161–188. 10.1007/s00285-004-0284-4Search in Google Scholar PubMed

[3] da Silveira Costa M. I., & dos Anjos L. (2019). Order of events: Optimal harvest fraction in a discrete time model of a spatially structured single population protected by a marine reserve can be overestimated due to an imprecise modeling of harvest timing. Ecological Modeling, 411, 108799. 10.1016/j.ecolmodel.2019.108799Search in Google Scholar

[4] Dilip D. S., Kilicman A., & Babu S. C. (2019). Asymptotic and boundedness behavior of a rational difference equation, Journal of Difference Equations and Applications, 25(3), 305–312. 10.1080/10236198.2019.1568424Search in Google Scholar

[5] Dilip D. S., & Mathew S. M. (2021). Stability analysis of a time varying population model without migration. Journal of Difference Equations and Applications, 27(10), 1525–1536. 10.1080/10236198.2021.1993839Search in Google Scholar

[6] Dilip D. S., & Philip T. (2020). Global stability, periodicity and boundedness behavior of a difference equation, Southeast Asian Bulletin of Mathematics, 44(3), 315–323. Search in Google Scholar

[7] Dilip D. S., & Philip T. (2022). Stability and bifurcation analysis of two spatial population dynamics models, International Journal of Applied and Computational Mathematics, 8, 108. 10.1007/s40819-022-01307-3Search in Google Scholar

[8] Dilip D. S., & Babu S. C. (2021). Rectangular p-periodicity and pq-periodicity of a rational difference system, Examples and counterexamples, 1, 100004. 10.1016/j.exco.2021.100004Search in Google Scholar

[9] Dilip D. S., & Mathew S. M. (2021). Dynamics of a second order nonlinear difference system with exponents. Journal of the Egyptian Mathematical Society, 29, 10. 10.1186/s42787-021-00119-6Search in Google Scholar

[10] Dilip D. S., Mathew S. M., & Elsayed E. M. (2020). Asymptotic and boundedness behavior of a second order difference equation. Journal of Computational Mathematica, 4(2), 68–77. 10.26524/cm82Search in Google Scholar

[11] Elaydi S. N. (2021). Discrete Chaos with Applications in Science and Engineering, New York: Chapman & Hall/CRC Press. Search in Google Scholar

[12] El-Metwally H., Grove E. A., Ladas G., Levins R., & Radin M. (2001). On the difference equation xn+1=α+βxn−1e−xn. Nonlinear Analysis, 47(7) , 4623–4634. 10.1016/S0362-546X(01)00575-2Search in Google Scholar

[13] Gouhier T. C., Guichard F., & Menge B. A. (2010). Ecological processes can synchronize marine population dynamics over continental scales. PNAS, 107(18), 8281–8286. 10.1073/pnas.0914588107Search in Google Scholar PubMed PubMed Central

[14] Khan A. Q., & Qureshi M. N. (2014). Behaviour of an exponential system of difference equations. Discrete Dynamics in Nature and Society, 10.1155/2014/607281. Search in Google Scholar

[15] Grove E. A., & Ladas G. (2005). Periodicities in Nonlinear Difference Equations, Chapman & Hall/CRC, Boco Raton, Florida. 10.1201/9781420037722Search in Google Scholar

[16] Leo Amalraj J., Maria Susai Manuel M., Baleanu D., & Dilip D. S. (2023). Global stability, periodicity, and bifurcation analysis of a difference equation, AIP Advances, 13, 015116. 10.1063/5.0106829Search in Google Scholar

[17] Liebhold A., Koenig W. D., & Bjornstad O. N. (2004). Spatial synchrony in population dynamics. Annual Review of Ecology, Evolution, and Systematics, 35(1), 467–490. 10.1146/annurev.ecolsys.34.011802.132516Search in Google Scholar

[18] Lutscher F. (2019). Integrodifference equations in spatial ecology, Interdisciplinary Applied Mathematics, Springer, Switzerland. 10.1007/978-3-030-29294-2Search in Google Scholar

[19] Mathew S. M., & Dilip D. S. (2007). Dynamics of interspecific k species competition model. Journal of Interdisciplinary Mathematics, 2022, 10.1080/09720502.2021.2012891. Search in Google Scholar

[20] Michael A. Radin, Periodic Character and Patterns of Recursive Sequences. Switzerland: Springer, 2018. 10.1007/978-3-030-01780-4Search in Google Scholar

[21] Ozturk I., Bozkurt F., & Ozen S. (2006). On the Difference Equation yn+1=α+βe−ynγ+yn−1. Applied Mathematics and Computations, 181, 1387–1393. 10.1016/j.amc.2006.03.007Search in Google Scholar

[22] Ozturk I, Bozkurt F., Ozen S. (2009). Global asymptotic behavior of the difference equation yn+1=αe−(nyn+(n−k)yn−k)β+nyn+(n−k)yn−k. Applied Mathematics Letters, 22, 595–599. 10.1016/j.aml.2008.06.037Search in Google Scholar

[23] Papaschinopoulos, Schinas C. J., & Stefanidou G. (2007). On a k-Order System of Lyness Type Difference Equations, Advances in Difference Equations, 31272. 10.1155/2007/31272Search in Google Scholar

[24] Papaschinopoulos G., Ellina G., & Papadopoulos K. B. (2014). Asymptotic behavior of the positive solutions of an exponential type system of difference equations. Applied Mathematics and Computation, 245, 181–190. 10.1016/j.amc.2014.07.074Search in Google Scholar

[25] Papaschinopoulos G., Fotiades N., & Schinas C. J. (2014). On a System of Difference Equations Including Negative Exponential Terms. Journal of Difference Equations and Applications, 20(5–6), 717–732. 10.1080/10236198.2013.814647Search in Google Scholar

[26] Papaschinopoulos G., Radin M. A., & Schinas C. J. (2011). On the system of two difference equations of exponential form. Mathematical and Computer Modelling, 54, 2969–2977. 10.1016/j.mcm.2011.07.019Search in Google Scholar

[27] Papaschinopoulos G., & Schinas C. J. (2012). On the Dynamics of two Exponential Type Systems of Difference Equations. Computers and Mathematics with Applications, 64, 2326–2334. 10.1016/j.camwa.2012.04.002Search in Google Scholar

[28] Papaschinopoulos G., Schinas C. J., & Ellina G. (2014). On the Dynamics of the Solutions of a Biological Model. Journal of Difference Equations and Applications, 20(5–6), 694–705. 10.1080/10236198.2013.806493Search in Google Scholar

[29] Psarros N., & Papaschinopoulos G. (2017). Long-term Behavior of Positive Solutions of an Exponentially Self-regulating System of Difference Equations, International Journal of Biomathematics, 10(3), 1750045. 10.1142/S1793524517500450Search in Google Scholar

[30] Qureshi M. N. M., Qadeer Khan A., Din Q. (2014). Asymptotic behavior of a Nicholson-Bailey model. Advances in Difference Equations, 2014, 62. 10.1186/1687-1847-2014-62Search in Google Scholar

[31] Ricker W. E. (1954). Stock and recruitment. Journal of Fisheries Research Board of Canada, 11, 559–623. 10.1139/f54-039Search in Google Scholar

[32] Robertson S., Cushing J., & Costantino R. (2012). Life stages: interactions and spatial patterns, Bulletin of Mathematical Biology, 74, 491–508. 10.1007/s11538-011-9705-xSearch in Google Scholar PubMed

[33] Tilman D., & Wedin D. (1991). Oscillations and chaos in the dynamics of a perennial grass. Letters to Nature, 353, 653–655. 10.1038/353653a0Search in Google Scholar

[34] Ying Y., Chen Y., Lin L., Gao T., & Quinn T. (2011). Risks of ignoring fish population spatial structure in fisheries management. Canadian Journal of Fisheries and Aquatic Sciences, 68, 2101–2120. 10.1139/f2011-116Search in Google Scholar

[35] Zhao X.-Q. (2017) Dynamical Systems in Population Biology, Second Edition, Springer, Switzerland. Search in Google Scholar

Received: 2023-08-04

Revised: 2024-01-06

Accepted: 2024-01-24

Published Online: 2024-02-29

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

Articles in the same Issue

https://doi.org/10.1515/cmb-2023-0121

Keywords for this article

boundedness; difference equations; persistence; asymptotic behavior

Creative Commons

BY 4.0