Dynamical analysis of a discrete two-patch model with the Allee effect and nonlinear dispersal

1.   Introduction

The Allee effect ^[1] refers to a phenomenon in population biology where the fitness and survival rates of individuals decrease when the population size becomes smaller. It suggests that certain species require a minimal population size to effectively find mates, protect against predation, or efficiently gather resources such as food. When the density of a population falls below this critical threshold, reproductive success and overall survival can be hindered (see ^[2,3,4] and the references cited therein).

In recent years, there has been increasing interest in studying the dynamics of patchy populations, where habitats are separated by unsuitable areas. The movement of individuals between patches, known as dispersal, is an important factor in determining the persistence and stability of such populations. Researchers have begun to study the dynamics of two-patch models with strong Allee effect ^[5]. These studies have revealed some interesting behaviors that were not previously understood ^[6,7,8]. Recently, Kang et al. ^[9] studied the following two-patch model with strong Allee effect:

where $D\in[0, 1]$ is the dispersal parameter, representing the fraction of population migration from one patch to the other per unit of time.

In ^[10], the authors proposed the following one species model with additive Allee effect and dispersal:

The results indicate that when the Allee effect constant $a$ increases or $m$ decreases, the total population abundance increases. In addition, when the dispersal rate $D_1$ increases or $D_2$ decreases, the total population density increases. An additive Allee effect can deduce complex dynamics such as saddle-nodes and transcritical bifurcations.

In ^[11], the author proposed the following model with linear dispersal:

where $D$ is the dispersal constant. Later, Allen ^[11] pointed out that the diffusion rate may be influenced by population density and proposed a patchy model with biased (nonlinear) diffusion as follows:

As is also shown in ^[12,13], populations might disperse non-linearly in complex real-life environments. Recently, Xia et al. ^[14] studied the following two-patch model with Allee effect and nonlinear dispersal:

where $u$, $v$ are the densities of the population in the first patch and the second patch, respectively. $A$ is the Allee effect constant. $r$, $d$, and $b$ are the birth rate, natural mortality, and death rate due to intra-prey competition in the first patch, respectively. $a$ and $c$ are the intrinsic growth rate and the death rate due to intra-prey competition of population in the second patch, respectively. $D$ is the dispersal coefficient.

It is well known that discrete models defined by difference equations are preferable to continuous-time ones when a species has non-overlapping generations or a limited population size. Recently, in ^[15,16], the authors studied discrete-time systems with Allee effects. In ^[17,18], the authors investigated the effect of dispersal on asymptotic total population size in the discrete two-patch model. However, to the best of the authors' knowledge, up to now a discrete two-patch model incorporating the Allee effect and nonlinear dispersal has never been put forth. Motivated by the above and with the assistance of the piecewise constant parameter method introduced by Jiang and Rogers ^[19], we convert model (1.1) into a discrete one as follows:

where $0\leq n$ $\leq t < n+1$ and $[t]$ is the greatest integer less than or equal to $t$. Since the right hand side of system (1.2) is constant over the interval $[n, n + 1$), integrating over $[n, t)$ and letting $t\rightarrow n+1$, it yields

Denoting $u(n)$ by $u_{n}$ and $v(n)$ by $v_{n}$, we thus obtain the following discrete-time model:

Different from the continuous two-patch model (1.1) in Xia et al. ^[14], this is the first time that the discrete two-patch model with nonlinear dispersal and the Allee effect has been proposed. The corresponding dynamic behaviors will be investigated in detail. We will also conclude that immediate nonlinear dispersal other than large nonlinear dispersal in ^[14] will be more conducive to the survival of the species.

The rest of the paper is organized as follows. In Section 2, we discuss the existence and stability of fixed points of system (1.4). In Section 3, we present the complete analysis of bifurcation. The influence of the Allee effect and nonlinear dispersal is presented in Section 4. A brief summary and discussion is in Section 5.

2.   Existence and stability of fixed points of (1.4)

To obtain the fixed point of (1.4), we need to solve the following equation:

Clearly, (1.4) always has the boundary fixed points $E_0(0, 0), E_{01}(0, \frac{a}{c+D})$. Moreover, the fixed point $E(u^*, 0)$ on the $u$ coordinate axis exists where $u^*$ satisfies the following equation:

If $d+bA+DA-r\ge 0$, then system (1.4) has no other equilibrium on the coordinate axis. In the following we investigate the case $d+bA+DA-r < 0$, that is, $r > d+bA+DA$. Notice that the discriminant of $F(u)$ is $\Delta = (b+D)^2A^2-2(d+r)(b+D)A+(d-r)^2$. If $\Delta > 0$, then system (1.4) has two boundary fixed points $E_{10}(u^*_{1}, 0)$ and $E_{20}(u^*_{2}, 0)$, where $u^*_{1}$ and $u^*_{2}$ are the positive roots of the equation $F(u) = 0$. If $\Delta = 0$, then system (1.4) has a boundary fixed point $E_{30}(u_3^*, 0)$. In order to simplify the analysis, let

Also, we have

Theorem 2.1 System (1.4) always has two boundary fixed points, i.e., $E_0(0, 0)$ and $E_{01}(0, \frac{a}{c+D}).$ Moreover,

(1) if $r\leq d+bA+DA$, then system (1.4) has no other boundary fixed point.

(2) if $r > d+bA+DA$, then

(i) system (1.4) also has two boundary fixed points $E_{10}(u_{1}^*, 0)$ and $E_{20}(u^*_{2}, 0)$ if $0 < A < A_{1}$, where $u^*_{1} = \frac{-(b+D)A-d+r-\sqrt{\Delta}}{2(b+D)}, u^*_{2} = \frac{-(b+D)A-d+r+\sqrt{\Delta}}{2(b+D)}$;

(ii) system (1.4) also has a boundary fixed point $E_{30}(u_{3}^*, 0)$ if $A = A_{1}$, where $u_{3}^* = \frac{-d+\sqrt{dr}}{b+D}$;

(iii) system (1.4) has no other boundary fixed point if $A > A_{1}$.

A positive fixed point $(u, v)$ of system (1.4) satisfies

that is,

To simplify the analysis, let

Equation (2.3) becomes

Notice that the discriminant of (2.4) is $\Delta_1(m) = (Am+n-r)^2-4mnA = (Am-n-r)^2-4nr$.

If $d = d^{**}$, Equation (2.4) becomes

Therefore, if $r\leq mA$, there is no positive equilibrium; if $r > mA$, Equation (2.5) has a unique positive real root.

If $0 < d < d^{**}$, Equation (2.4) has a unique positive real root.

If $d\geq d^{**}+r$, i.e., $n\geq r$, Equation (2.4) has no positive real root.

In the following, we investigate the case $d^{**} < d < d^{**}+r$. If $m\geq\frac{r-n}{A}$, Equation (2.4) has no positive real root. Next, we consider the case $m < \frac{r-n}{A}$. The discriminant of $\Delta_1(m)$ is $\Delta_2 = 16nrA^2 > 0$. Thus, $\Delta_1(m) = 0$ has two positive real roots, i.e.,

We get $m_1 < \frac{r-n}{A} < \frac{r+n}{A} < m_2$. Therefore, if $m = m_1$, Equation (2.4) has a unique positive real root. If $m < m_1$, Equation (2.4) has two positive real roots. If $m > m_1$, it follows that $\Delta_1(m) < 0$ and then Eq (2.4) has no positive real root. Thus, we have the existence of a positive fixed point as follows:

Theorem 2.2

(1) If $0 < d < d^{**}$, then system (1.4) has a positive fixed point $E_1^*(u_1, v_1)$.

(2) If $d = d^{**}$, then

(i) system (1.4) has a positive fixed point $E_1^*(u_{1}, v_1)$ when $r > mA$;

(ii) system (1.4) has no positive fixed point when $r\leq mA$.

(3) If $d^{**} < d < d^{**}+r$, then

(i) system (1.4) has two positive fixed points $E_1^*(u_{1}, v_1)$ and $E_2^*(u_{2}, v_2)$ when $m < m_1$;

(ii) system (1.4) has a positive fixed point $E_3^*(u_{3}, v_3)$ when $m = m_1$;

(iii) system (1.4) has no positive fixed point when $m > m_1$.

(4) If $d\geq d^{**}+r$, then system (1.4) has no positive fixed point.

And,

Next, we will consider the local stability of the fixed point. The Jacobian matrix of system (1.4) at the equilibrium $E(u, v)$ is

where $N = \exp\big(\frac {ru}{A + u}-d-bu + D(v - u)\big)$ and $G = \exp\big(a-cv+ D(u - v)\big)$. Let $\lambda_1$ and $\lambda_2$ be the two eigenvalues of $J(E).$ To study the local stability of these fixed points, we will use the classification definition of fixed points in ^[20] and obtain the following results:

Theorem 2.3 $E_0(0, 0)$ is always a saddle.

Proof. At the trivial fixed point $E_0(0, 0)$, the Jacobian matrix is

with eigenvalues $\lambda _1 = e^{-d}\in(0, 1)$ and $\lambda _2 = e^a > 1$. Hence, $E_0(0, 0)$ is always a saddle.

Theorem 2.4 For $E_{01}(0, \frac{a}{c+D})$,

(1) it is a sink if and only if $0 < a < \min\big\{2, \frac{d(c+D)}{D}\big\}$;

(2) it is a source if and only if $a > \max\big\{2, \frac{d(c+D)}{D}\big\}$;

(3) it is non-hyperbolic if either $a = 2$ or $a = \frac{d(c+D)}{D}$;

(4) it is a saddle, except for in cases (1)–(3).

Proof. At the boundary fixed point $E_{01}(0, \frac{a}{c+D})$, the Jacobian matrix is

with eigenvalues $\lambda _1 = e^{\frac{aD}{c+D}-d}$ and $\lambda _2 = 1-a$. It is easy to see that

and

Hence, the result is proved.

Theorem 2.5 The boundary fixed point $E_{30}(u_3^*, 0)$ is a non-hyperbolic.

Proof. At the boundary fixed point $E_{30}(u_3^*, 0)$, the Jacobian matrix is

with eigenvalues $\lambda _1 = 1$ and $\lambda _2 = e^{a+\frac{(\sqrt{dr}-d)D}{b+D}} > 1$. Hence, the proof is complete.

Theorem 2.6 For the equilibrium $E_{i0}(u_{i}^*, 0) (i = 1, 2)$,

(1) it is a source if and only if $\frac {rA}{(A +u_{i}^*)^2} > b+D$;

(2) it is non-hyperbolic if and only if $\frac {rA}{(A +u_{i}^*)^2} = b+D$;

(3) it is a saddle if and only if $b+D-2 < \frac {rA}{(A +u_{i}^*)^2} < b+D$.

Proof. At the boundary fixed point $E_{i0}(u_{i}^*, 0), (i = 1, 2)$, the Jacobian matrix is

The corresponding eigenvalues are $\lambda _1 = 1+u_{i}^*(\frac {r}{A + u_{i}^*}-\frac {ru_{i}^*}{(A +u_{i}^*)^2}-b-D)$ and $\lambda _2 = e^{a+u_{i}^*D} > 1$. When $b+D > 2$, it is easy to see that

Hence, the proof is complete.

Theorem 2.7 If $\frac{r}{A} < b$, then the positive fixed point $E^*(u_{*}, v_*)$ is

(1) a sink if and only if $P < 1+Q$ and $Q < 1$;

(2) a source if and only if $P < 1+Q$ and $Q > 1$;

(3) non-hyperbolic if and only if $P = 1+Q$;

(4) a saddle if and only if $P > 1+Q$, where

Proof. At the positive fixed point $E^*(u_{*}, v_*)$, the Jacobian matrix is

The characteristic equation for $J(E^*)$ is $F(\lambda) = \lambda ^2+P\lambda +Q = 0$.

Thus,

If $\frac{r}{A} < b$ holds, then $\dfrac{rA}{(A +u_{*})^2} < \frac{r}{A} < b$ and $\frac{v_*D^2}{a+Du_*} = \frac{D^2}{c+D} < D$, implying that $F(1) > 0$. If the conditions $P < 1+Q$ and $Q < 1$ hold, then $F(-1) > 0, Q < 1$. Hence, according to Lemmas 1 and 2 in ^[21], we obtain that $E^*(u_{*}, v_*)$ is a sink, which is stable. Cases (2)–(4) can be proved in the same way. Hence, Theorem 2.7 is obtained.

The feasibility and local stability criteria of the fixed points of system (1.4) are given in Table 1.

3.   Bifurcation analysis

In the following, we will use the central manifold and bifurcation theories ^[22,23] to discuss bifurcation in system (1.4). In detail, we will analyze the existence of a flip bifurcation at the fixed point $E_{01}(0, \frac{a}{c+D})$ and a fold bifurcation at the fixed points $E_{30}(u_3^*, 0)$ and $E_{3}^*(u_3, v_3)$.

Theorem 3.1 System (1.4) undergoes a flip bifurcation at $E_{01}(0, \frac{a}{c+D})$ when the parameters are varied in a small range of $F_A = \Big\{(a, r, d, D, c):a = 2, a\neq \frac{d(c+D)}{D}, r > 0, d > 0, D > 0, c > 0\Big\}$.

Proof. Theorem 2.4(3) shows that if $a = 2$ and $a\neq \frac{d(c+D)}{D}$, then one of the eigenvalues of the fixed point $E_{01}(0, \frac{a}{c+D})$ is $-1$, and the other eigenvalue is neither 1 nor $-1$. Regarding $a$ as the bifurcation parameter and perturbing $a$ by $\zeta$, system (1.4) can be seen as the two-dimensional map below:

where $|\zeta|\ll1$. By allowing $x_1 = u, y_1 = v-\frac{a}{c+D}$, we then transform the fixed point $E_{01}(0, \frac{a}{c+D})$ of the system (1.4) into the origin and shift model (3.1) into

Here,

and

Next, we use the following transformation:

and then get the normal form of (3.1) as follows:

where

and

By the center manifold theory, the stability of $(x_2, y_2) = (0, 0)$ near $\tilde{\zeta} = 0$ can be determined by studying a one-parameter family of reduced equations on a center manifold, which can be represented as follows:

Here, $y_2$ and $\tilde{\zeta}$ are sufficiently small. We assume that

Then, $h(\tilde{\zeta, }y_2)$ satisfies

Substituting (3.5) into (3.6), we obtain

Thus, the map on the center manifold is

A straightforward calculation gives

Therefore, by ^[24], we can get a flip bifurcation at $E_{01}$.

In Figure 1(a), we can see that the fixed point $E_{01}(0, \frac{a}{c+D})$ is stable if $0 < a < 2$, and it is unstable if $a > 2$. The maximum Lyapunov exponent is shown in Figure 1(b).

Theorem 3.2 System (1.4) undergoes a fold bifurcation at $E_{30}(u_3^*, 0)$ if the parameters vary in the small neighborhood of $F_B = \Big\{(r, d, D, A):r > d > 0, D > 0, A = A_1\Big\}$.

Proof. From Theorem 2.5, one can have that the eigenvalues of $J(E_{30})$ are $\lambda _1 = 1$ and $\lambda _2\neq1, -1$. We choose $A$ as the bifurcation parameter to study fold bifurcation. Let $A = A_1+\xi$, where $\xi$ is a small perturbation and will be treated as a new variable. System (1.4) is then rewritten as

We first use the change of variables $x_3 = u-u_3^*, y_3 = v$ to translate the fixed point $E_{30}(u_3^*, 0)$ to the origin, and then transform system (3.7) into

Here,

and

Then, with the transformation,

system (3.8) becomes

where

and

According to the central manifold theorem, suppose that an approximate representation of the central manifold $W^c_2(0, 0, 0)$ is as follows:

where $x_4$ and $\tilde{\xi}$ are sufficiently small.

By a simple comparison, $k_1 = k_2 = k_3 = 0$ can be obtained. Therefore, the following expression can be easily evaluated:

A straightforward calculation gives

which leads to the existence of a fold bifurcation.

In Figure 2, the fixed points $E_{10}(u_1^{*}, 0)$ and $E_{20}(u_2^{*}, 0)$ bifurcate from $E_{30}(u_3^{*}, 0)$ when $0 < A < A_1$, coalesce at $E_{30}(u_3^{*}, 0)$ when $A = A_1$, and disappear when $A > A_1$.

In addition, the fixed points $E_{1}^*(u_3, v_3)$ and $E_{2}^*(u_3, v_3)$ bifurcate from $E_{3}^*(u_3, v_3)$ when $m < m_1$, merge at $E_{3}^*(u_3, v_3)$ when $m = m_1$, and vanish when $m > m_1$. Thus, we can have the existence of a fold bifurcation as follows:

Theorem 3.3 When $d^{**} < d < d^{**}+r$ and $m = m_1$, system (1.4) undergoes a fold bifurcation at $E_{3}^*(u_3, v_3)$.

Proof. Denote

Since $m = b+\frac{cD}{c+D}$, we choose $b$ as a bifurcation parameter for studying the fold bifurcation of $E_{3}^*(u_3, v_3)$. Let $b = b^*+\delta$, where $b^* = m-\frac{cD}{c+D}.$ $\delta$ is a small perturbation and will be treated as a new variable. System (1.4) is then rewritten as

The eigenvalues of $J_{E_{3}^*}$ are $\lambda_1 = 1, \lambda_2\neq1$ by Theorem 2.7. We first use the change of variables $x_5 = u-u_3, y_5 = v-v_3$ to translate the positive fixed point $E_{3}^*(u_3, v_3)$ to the origin and then transform system (1.4) into

Here,

and

We construct an invertible matrix

and use the translation

and then the map (3.12) becomes

where

and

Next, an approximate representation of the central manifold $W^c_3(0, 0, 0)$ is supposed as follows:

where

Therefore, we consider the map restricted to the center manifold $W^c_3(0, 0, 0)$:

Here,

A straightforward calculation gives

which leads to the existence of a fold bifurcation.

In Figure 3, the parameter values are $d = 0.7, A = 1, D = 0.6, r = 1.2, a = 0.8, c = 0.2$, and the initial value is $(u_0, v_0) = (0.7, 0.6)$. We have $b = 0.456$ and calculate $u_3 = 0.436$ from Theorem 2.2. In other words, the positive fixed points $E_1^*(u_1, v_1)$ and $E_2^*(u_2, v_2)$ bifurcate from $E_3^*(u_3, v_3)$ when $0.38 < b < 0.456$, merge at $E_3^*(u_3, v_3)$ when $b = 0.456$, and vanish when $b > 0.456$.

4.   Effect of the Allee effect and nonlinear dispersal

In this section, we will explore the impact of the Allee effect and nonlinear dispersal on population survival through numerical simulations. As is well known, when population density is low, the Allee effect can lead to a decrease in cooperation and mutual assistance among individuals, thereby reducing the growth rate of the population. We choose

and the initial values of the system (1.4) are $(0.7, 0.2)$ and $(0.5, 0.4)$. As is shown in Figure 4, when the Allee constant $A = 0$, the population density $u$ is relatively high and oscillates greatly. When the Allee constant $A = 0.5$ or $A = 1.5$, the population density decreases and oscillations become less frequent. Increasing Allee effect can effectively reduce the oscillation amplitude of the solution. In this sense, the Allee effect has a stabilizing ability.

In Figure 5, we take $(a, r, c, A, b, d) = (1.5, 1.2, 0.2, 3, 0.3, 0.4)$ and initial values $(0.7, 0.2)$ and $(0.5, 0.4)$. It follows that when the dispersal coefficient $D = 0$, population $u$ goes extinct while population $v$ is permanent. When $D = 0.1$, both populations coexist. As the dispersal coefficient keeps increasing, the fluctuations of both populations become significant. Figure 5 shows that when there is no nonlinear dispersal, the population will tend to be extinct due to the Allee effect. Moreover, large nonlinear dispersal will also lead to fluctuations of the population. That is to say, proper proliferation is conducive to population survival.

5.   Summary

In this paper, we have considered a discrete two-patch model with Allee effect and nonlinear dispersal based on Xia et al. ^[14]. The dynamic behaviors have been analyzed in detail. We discuss the existence of fixed points, the flip bifurcation at the boundary fixed point $E_{01}(0, \frac{a}{c+D})$, the fold bifurcation at the boundary fixed point $E_{30}(u_3^*, 0)$, and the positive fixed point $E_3^*(u_3, v_3)$, respectively. The impact of the Allee effect and nonlinear dispersal on population survival is also presented through numerical simulation.

The topological types of the fixed points are substantially different from those in Xia et al. ^[14] and Grumbach et al. ^[17]. In detail, the model in this manuscript has a maximum of two positive fixed points, whereas the similar discrete model in ^[17] has a unique positive fixed point. Moreover, the discrete model can experience a flip bifurcation at the boundary fixed point $E_{01}$ when the conditions $a = 2, a\neq \frac{d(c+D)}{D}$ hold, while the corresponding boundary equilibrium $E_{v}$ in the continuous model ^[14] is globally asymptotically stable under certain conditions. In addition, the results in Xia et al. ^[14] show that for the continuous case, large nonlinear dispersal may prevent the extinction of species due to the Allee effect. But in this paper, we obtain that small or high nonlinear dispersal will lead to extinction or fluctuations, and only moderate and appropriate nonlinear dispersal can make the species become permanent. In other words, suitable and modest dispersal is advantageous to the survival of the species. The aforementioned point is due to the fact that modest dispersal minimizes intense rivalry with others for resources like food, habitat, and breeding space. But, when the dispersal coefficient is large, the species will go extinct because dispersal increases the complexity of the system, and makes it more susceptible to external influence. Therefore, the conclusion of this article is a complement and promotion of ^[14] and ^[17].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant (11601085) and the Natural Science Foundation of Fujian Province (2021J01614, 2021J01613).

Conflict of interest

The authors declare there is no conflict of interest.