On an age structured population model with density-dependent dispersals between two patches

1.   Introduction

Among the various features for population dynamics of a single species are the heterogeneity of different habitats for the species and the age structure of species. Due to some natural barriers such as rivers and mountains, man-made barriers such as buildings and highways, and the local habitation of insects and small animals, a patchy environment is often considered in population models. Random migrations betweens patches are often assumed which adopts constant dispersal rates between patches. However, for many species, their intention of dispersal is heavily influenced by the population distribution among the patches and the variance of resources in all patches. This is because individuals typically seek habitats with abundant sources and mates and flee from the threat of predators ^[1]. On the other hand, most species have a clear and well-defined age-structure consisting of immature (juvenile) and mature (adult) stages, and behave differently at different stages. Therefore, it is vital to explore the effects of density dependent dispersal rates on population dynamics that involve both heterogeneity among habitats and the age structure.

In ^[2], Levin proposed a two-patch model to explore the spatial heterogeneity of a single species population. In ^[3], the authors systematically developed age-structured population models to describe the evolution of populations over multiple life stages. Combining the two aspects in ^[2] and ^[3] and assuming constant dispersal rates between patches, So et. al. ^[4] derived a system of delayed differential equations (DDEs) to describe the adult populations of a single species living on two patches; and by analyzing the derived model, they were able to illustrate some effect of the immature death rate on global dynamics. In subsequent studies in ^[5,6,7] on the model assuming identical patches, rich dynamics including Hopf bifurcation, synchronized periodic oscillations and unstable phase-locked oscillations were obtained by applying local bifurcation theory to the DDE model; and global convergence to the trivial state or a positive homogeneous state was established under certain conditions in ^[6]. A recent study ^[8] examined the model with two general patches and three typical birth functions (linear, Ricker and Allee functions). For these types of birth function and in two patches, the author investigated the extinction and persistence of the species under consideration, the boundedness of the solutions to the model, and certain attractors of the model system.

In reality, the mechanisms of dispersion/migration of many species can be very complicated, and constant dispersal rate would be too simple and too ideal. In ^[1] the authors introduced an adaption mechanism of dispersion in a predator-prey model by assuming that according dispersal rates are dependent of the population densities; more specifically, they assumed that a higher density of prey population in a patch results in a lower migration rate of predators in the same patch, and a prey tends to flee from the regions with a high density of predators. They incorporated adaptive migration with multiple timescales, and numerically observed co-dimension two Bautin bifurcation. Recently, the authors in ^[9,10,11] constructed some mathematical models to examine the adaptive dispersal rates of predator-prey populations and compared the adaptive dispersal rate of competitive species between their habitats. The authors in ^[9] proposed a patch payoff concept to describe adaptive dispersal; this concept entails assuming that the net movement between patches is always toward the patch with the higher per capita growth rate. In ^[10,11], the authors described the fitness of individuals according to the local growth rate and assumed that each individual tends to move to a patch with higher fitness. By analyzing the local stability of equilibria, they obtained some results that can hardly occur in patch models with constant dispersal rates, such as an adaption-induced change in persistence patterns ^[10,11], oscillations and Neimark-Sacker bifurcation ^[9]. The terms of adaptive dispersal typically result in complex nonlinearity, and hence, the global dynamics remains an open and challenging problem in the above mentioned three studies.

Because of the limitation in resources, biological species tend to leave a crowded habitat. In this study, we incorporate a density-dependent dispersal strategy into the prototypical age-structured model proposed in ^[3]. Motivated by and following the derivation of the model in ^[4], we derive a DDE model for the adult population in the case of density dependent dispersal. In contrast to the model in ^[4], there are two additional difficulties in theoretically studying this new DDE model: the nonlinear instantaneous dispersal terms and the nonlinear delayed birth terms mediated by the nonlinear dispersals. Allowing a non-monotone birth function such as the Ricker or Allee functions is another common challenge in population models.

Global dynamics of a DDE system is generally very difficult to investigate, even in a scalar DDE (see, e.g., ^[12,13,14]). The authors in ^[12,13] used a domain decomposition method to study mono-stability and bi-stability in scaler DDEs respectively, and obtained some nice results, including a Poincaré-Bendixson type result, heteroclinic orbits, and basins of attraction. This method generally fails for systems of more than one DDE, and hence, the convergent dynamics in a two-variable DDE system is still a hard problem, especially when it involves a non-monotone feedback. Among the known attempts, Xu ^[6] investigated a population model in two identical patches and attacked this problem by seeking an attractor for the system and applying the monotone dynamical system theory; and for a species living in a general patchy environment with heterogeneity between two patches, Terry ^[15] examined the population dynamics of the model with impulsive culling of the adults by the comparison with certain linear DDE system.

The main goal of this study is to explore the impact of density-dependent migration on the dynamics of an age-structured population model in two heterogeneous patches. In Section 2, we derive our model which turns out to be a system of DDEs with nonlinearity in the dispersion terms and spatial non-locality in the delayed nonlinear birth terms. In Section 3, we show well-posedness of the proposed model, and focus on the case with only adult dispersal. In Section 4, by analyzing this model, we obtain some criteria for the existence of positive equilibrium and determine the uniform persistence of populations in the two patches. Making use of these results, we discuss how the density-dependent dispersal affects persistence or extinction of populations. In addition, we investigate the global convergence of the model by dividing the phase space, identifying an attractor for the model system, and applying the theory of monotone dynamical systems. We also identify a positive invariant set within which possible sustained oscillations may be induced by the maturation delay. We conclude the paper by Section 5, in which we present some numerical simulations and offer some discussion of our results.

2.   Model derivation

We follow the procedure in ^[4] to derive our model with density dependent dispersals. Let $u_i(t, a)$ denote the density of a single species at time $t$ ($t\geq 0$) of age $a$ ($a\geq 0$) on patch $i$ ($i = 1, 2$). In the case of density-independent dispersal rates between the patches as considered in ^[4], by the basic equation for age structured populations in ^[3], there hold

Here $d_i(a)$ is the death rate of the individuals of age $a$ in patch $i$, and $D_j(a)$ denotes the dispersal rate of the species of age $a$ from patch $j$ to patch $i$, for $1\leq j\neq i\leq 2$. So et. al. ^[4] considered two age stages: immatures and matures, and assumed that

and

where $r\geq 0$ denotes the maturation age, $I$ and $M$ stand for immature and mature respectively. Let $w_i(t)$ be the total number of adults at time $t$ on patch $i$, that is,

Noting that $u_i(t, \infty) = 0$ and that only adults produce and hence

where $b_i$ is a birth function of the population on the $i$-th patch, So et. al. derived the following system in ^[4] for the adult populations

where $\sigma = e^{-\int^r_0d_I(a)da}$ and $\hat{D}(a) = D_1(a)+D_2(a)$.

Now we modify (2.1) by incorporating density dependent dispersion between patches. In general, the dispersal rate of individuals may depend on the population of all ages. For simplicity of mathematics, in this paper we only consider a scenario that the dispersal rate of individuals depend only on the mature population and this leads to the following revision of (2.1):

where

with the functions $f_i$ satisfying

Such a scenario can be justified as below. Firstly, for many species, each single mature individual reproduces a fixed numbers of offsprings. For examples, a couple of mature shorebirds lay fixed number of eggs in a clutch and then hatch and breed them ^[16,17]; some particular species of geckos and lizards also lay fixed number of eggs to maximize their reproductive success ^[18]. For such a species, it is reasonable to assume that the population at age $a$ is proportional to the mature population. Secondly, the density-dependent dispersal is proposed to investigate the strategy of dispersing toward a relatively abundant resource, and many species compete for resources mainly for their mature individuals because the immature individuals consume relatively little or none resource directly but obtain all nutrition from their raisers. With the above two considerations, the dispersal rate given by (2.7) is reasonable for such species.

The value of $\eta_i\geq 0$ in (2.8) represents the strength of adaptive dispersal in patch $i$. We call $D_i(a)$ and $\tilde{D}_i(a, w_i)$ the intrinsic and adaptive dispersal rates respectively. When $\eta_1 = \eta_2 = 0$, (2.6) reduces to (2.1) reflecting density-independent dispersal rates. In this work, we consider (2.6) with the assumptions (2.2), (2.3) and the typical function

By integrating (2.6) with respect to $a$ from $r$ to $\infty$ and making use of the fact $u_i(t, \infty) = 0$, we obtain

To describe $u_i(t, r)$ in terms of the adult populations, we fix $s$ and consider the function

From (2.6),

for $s\leq t\leq s+r$ and $1\leq i\neq j\leq 2$. Since $d_i(t-s) = d_I(t-s)$ for $i = 1, 2$ and $s\leq t\leq s+r$, we have

and then, by (2.4),

For $s\leq t\leq s+r$, we rewrite (2.10) as

where $K^s(a) = d_I(a)+D_i(a)f_i(w_i(a+s))+D_j(a)f_j(w_j(a+s))$. Solving (2.11), we obtain

Thus,

where

In addition,

Combining this with (2.12), we have

with $\sigma = \exp(-\int_0^r d_I(a) \, da)$ as before. Plugging this into (2.9), we finally obtain the following system

where $w_{it} \in C([-r, 0], \mathbb{R}) \to \mathbb{R}$ is defined by $w_{it}(\theta) = w_i(t+\theta)$ for $\theta \in [-r, 0]$, $i = 1, 2$, and

Here $b_i(\cdot): \mathbb{R}_+ \to \mathbb{R}_+, \, \, i = 1, 2,$ are still the birth functions which are assumed to be continuously differentiable. Depending on the particular species, they may have different forms. Comparing with (2.5), we note that both instantaneous and delayed dispersal terms have been modified in (2.13) by nonlinear functions $f_i(w_i), \, \, \, i = 1, 2$, representing the density dependent adaptive dispersal strategy.

3.   Well-posedness of the model

The model system (2.13) is a system of delay differential equations, which is associated with the phase space $C: = C([-r, 0], \mathbb{R}^2)$. By the fundamental theory of functional differential equations (see, e.g., ^[19]) or directly by an argument of steps, one knows that with the initial condition

(2.13) has a unique solution $w(t) = (w_1(t), w_2(t))$ for $t > 0$. Next, we will show that if the initial functions $\phi_1$ and $\phi_2$ are nonnegative, that is,

then the unique solution also remains nonnegative; moreover, if the birth function is bounded, then the solution is bounded.

Lemma 3.1. If $\phi \in C_+$, then the solution of (2.13) remains non-negative for all $t \ge 0$.

Proof. By the non-negativity of the functions $f_i(\cdot)$ and $b_i(\cdot), \, \, \, i = 1, 2$, [Theorem 2.1, P81] in ^[20] applies, leading to the conclusion of the lemma.

Lemma 3.2. Assume that the birth functions $b_i(\cdot), \, \, i = 1, 2$, are bounded and let $\phi \in C_+$. Then the corresponding solution of (2.13) is bounded.

Proof. Let $B = \sup \{b_i(u); u \in \mathbb{R}_+, \, i = 1, 2.\}$ and $\underline{d} = \min \{d_{1, m}, d_{2, m} \}$. Then the total adult population $W(t) = w_1(t)+w_2(t)$, satisfies

This leads to

implying that $W(t)$ is bounded. This together with the non-negativity of $w_1(t)$ and $w_2(t)$ then further implies the boundedness of $w_1(t)$ and $w_2(t)$, completing the proof.

In the model system (2.13)- (2.14), the general forms of the dispersion rate functions $D_i(a)$ seem to prevent us from going further in analyzing the mode, due to the occurrence of multiple iterated integrals. In the sequel, we only consider the following simple case for the dispersion rate functions:

With these simple dispersion rate functions, the delayed term (2.14) reduces to

Also, the birth functions $b_i(\cdot)$ have many choices depending on the particular species. In the rest of this paper, we will chose the form of the Ricker's birth function for both $b_i(\cdot), \, \, \, i = 1, 2$:

For the above choices for the birth and dispersion rate functions, the model (2.13) becomes

4.   Analysis of the simpler version (3.3)

First, we study the equilibria of the system (3.3) which are governed by the following system of algebraic equations:

For $w\geq 0$ and $i = 1, 2$, denote

Then (4.1) is equivalent to $G_j(w_j) = F_i(w_i)$ for $1\leq i\neq j\leq 2$. Note that $G_i(0) = 0$ and $\lim_{w\rightarrow\infty}G_i(w) = \infty$ for $i = 1, 2$. Obviously, each function $G_i$ has an increasing inverse function $G_i^{-1}$ on $[0, \infty)$ satisfying $G^{-1}(0) = 0$ and $\lim_{w\rightarrow\infty}G_i^{-1}(w) = \infty$. Hence, the curve $L_i$ in the first quadrant is described by

for $w_i\geq 0$ and $F_i(w_i) > 0$. According to the feature of functions $F_i(w)$, we obtain the following properties of the curves $L_i$.

Lemma 4.1. The following hold:

(i) if

then the function $H_i$ is increasing on $[0, \infty)$ and only has the trivial root $w = 0$;

(ii) if

then in addition to the trivial root, $H_i$ also has a positive root $w_i^r > 0$ and is strictly increasing on $[w_i^r, \infty)$.

Proof. First, we explore features of the functions $F_i$. Write it as $F_i(w) = w [h_i(w)-k_i(w)]$ where $h_i(w): = d_{i, m}+D_{i, m}f_i(w)$ is a nondecreasing and $k_i(w): = \sigma b_i(w)/w = \sigma \beta_i e^{-\gamma_i w}$ is a nonincreasing. Thus, if $\sigma\beta_i \leq d_{i, m}+D_{i, m}$, then $h_i(0) \ge k_i(0)$ and $F_i(w) = 0$ has no positive root and $F_i(w)$ is increasing on $[0, \infty)$; and if $d_{i, m}+D_{i, m} < \sigma\beta_i$, then in addition to the trivial root, $F_i(w) = 0$ also has an unique positive root $w_i^r > 0$ and $F_i(w)$ is increasing and positive on $[w_i^r, \infty)$ and negative on $[0, w_i^r)$. Now combining these features of $F_i(w)$ with the properties of $G_j^{-1}(\cdot)$ immediately leads to the assertions on $H_i$, completing the proof.

Lemma 4.2. The graph of $L_i$ in the first quadrant admits an asymptote $w_j = p_iw_i+q_i$, $1\leq i\neq j\leq 2$, for some constant $q_i$ and

Proof. In the first quadrant, a point $(w_i, w_j)$ on the graph of $L_i$ obeys

Note that $\lim_{w_i\rightarrow\infty}F_i(w) = \infty$ and $\lim_{w\rightarrow\infty}G_j(w) = \infty$. In addition, a point $(w_i, w_j)$ on the curve $L_i$ satisfies $w_j\rightarrow\infty$ as $w_i\rightarrow\infty$. The assertion follows from the above and the form of $f_i(w)$, completing the proof.

Based on previous results, we discuss existence of positive equilibria in the following.

Theorem 4.1. The system (3.3) always has the trivial equilibrium $E^0 = (0, 0)$. In addition, it admits an unique positive equilibrium when one of the following conditions holds:

($\mathcal{A}$) $\sigma\beta_i \le d_{i, m}+D_{i, m}$ for $i = 1, 2$, and $\Pi_{k = 1, 2}(d_{k, m}+D_{k, m}-\sigma\beta_k) < D_{1, m}D_{2, m}$;

($\mathcal{B}$) $\sigma\beta_i > d_{i, m}+D_{i, m}$ for $i = 1$ or $2$.

Proof. $E^0 = (0, 0)$ is obviously always an equilibrium of (3.3). An non-trivial equilibrium is governed by

for $w_i > 0$ with $H_i(w_i) > 0$. Note that each $G_i^{-1}$ is increasing on $[0, \infty)$. As explored in Lemma 4.1, $F_i$ is also increasing on $[0, \infty)$ under ($\mathcal{A}$) (respectively on $[w_i^r, \infty)$ under ($\mathcal{B}$)), when $F_i(w) = 0$ has one (respectively, two) nonnegative root. Therefore, the function $H_i$ is increasing on $[0, \infty)$ under ($\mathcal{A}$) (respectively on $[w_i^r, \infty)$ under ($\mathcal{B}$)). Next, we discuss the existence of roots to (4.3).

First, we consider the system (3.3) under the criterion($\mathcal{A}$). From Lemma 4.2, it always holds that

where $H_i'(\infty): = \lim_{w\rightarrow\infty}H_i'(w)$. The condition $\sigma\beta_i < d_{i, m}+D_{i, m}$ for $i = 1, 2$ implies that each $F_i(w) = 0$ only has the trivial root and $F_i(w) > 0$ for all $w > 0$, and then the function $G_i^{-1}$ is defined and increasing on $[0, \infty)$. Therefore, each $H_i$ is defined, continuous and increasing on $[0, \infty)$. The further condition in ($\mathcal{A}$),

is equivalent to

Since $H_1(0) = H_2^{-1}(0)$ and $H_1'(0) < (H_2^{-1})'(0)$, there is a small $\check{w} > 0$ such that $H_1(\check{w}) < H_2^{-1}(\check{w})$. From $H_1'(\infty) > (H_2^{-1})'(\infty)$, there is a large enough $\hat{w}$ such that $H_1(\hat{w}) > H_2^{-1}(\hat{w})$. Since functions $H_1$ and $H_2^{-1}$ are continuous on $[0, \infty)$, there exists at least one point $(w_1^*, w_2^*)$ in the first quadrant such that $H_1(w_1^*) = H_2^{-1}(w_1^*) = w_2^*$, giving a positive root to (4.3). From (4.3) and the monotonicity of each function $H_i$, $i = 1, 2$, the positive equilibrium is unique when it exists. See Figure 1(a) for a demonstration of this case.

Assume that the criterion($\mathcal{B}$) holds only for $i = 1$ then $F_1(w) = 0$ has two nonnegative roots, 0, $w_1^r$, and $F_1(w) > 0$ for all $w > w_1^r$. Hence, the function $H_1$ is defined, continuous and increasing on $[w_1^r, \infty)$. Note that $H_1(w_1^r) = 0\leq (H_2)^{-1}(0) < (H_2)^{-1}(w_1^r)$ and the fact (4.4) implies a large enough $\hat{w}$ such that $H_1(\hat{w}) > H_2^{-1}(\hat{w})$. Similar to the case with ($\mathcal{A}$), but considering $[w_1^r, \infty)$ as the domain of $H_1$, there exists a unique positive root to (4.3). See Figure 1(b) for an illustration of this case. Similarly, the assertion holds when ($\mathcal{B}$) holds only for $i = 2$ or for both $i = 1, 2$. The proof is completed.

Theorem 4.2. Consider the system (3.3) and assume that neither($\mathcal{A}$) nor ($\mathcal{B}$) holds. Then the population in (3.3) goes to extinction in both patches, that is, every solution $(w_1(t), w_2(t))$ of (3.3) with $\phi\in C_+$ satisfies $\lim_{t\rightarrow\infty}w_i(t) = 0, \; {\rm for}\; i = 1, 2$.

Proof. From (3.3), it holds that

where

The auxiliary system

is cooperative and irreducible ^[20] and admits the unique trivial equilibrium when neither ($\mathcal{A}$) nor ($\mathcal{B}$) holds. Similar to the proof of Lemma 3.2, the solutions of (4.5) with initial conditions in $C_+$ are uniformly bounded. From [20,Theorem 5.4.1], each solution of (4.5) converges to the unique trivial equilibrium, i.e. $\lim_{t\rightarrow\infty}x_i(t) = 0, \; {\rm for}\; i = 1, 2$. From the comparison principle [20,Theorem 5.1.1], $w_i(t)\leq x_i(t)$ for $t\geq 0$, $i = 1, 2$, and then it completes the proof.

Next, we consider permanence of the populations described by (3.3).

Theorem 4.3. Consider the system (3.3) and assume that either ($\mathcal{A}$) or ($\mathcal{B}$) holds. Then (3.3) is uniformly persistent in the sense that there is a positive constant $\rho^*$ such that every solution $(w_1(t), w_2(t))$ of (3.3) with $\phi\in C_+\setminus\{\hat{0}\}$, where $\hat{k}$ denotes the constant function on $[-\tau, 0]$ taking value $k$, satisfies

Proof. First, we suppose that ($\mathcal{A}$) holds. Define

Obviously, both $X$ and $X_0$ are positive invariant under the semiflow of (3.3) and $X\setminus X_0$ is relatively close in $X$. From Lemma 3.2, system (3.3) is point dissipative. Set

where $\Phi(t)$ is the semiflow generated by (3.3). Then obviously $M_{\partial} \subset \partial X_0$. We claim that $M_{\partial} = \{ \hat{0} \}$. It is clear that $\{\hat{0}\} \subset M_{\partial}$, so it suffice to show $M_{\partial} \subset \{\hat{0} \}$. Assume the opposite, there is $\phi = (\phi_1, \phi_2) \in M_{\partial}$ with $\phi \neq \hat{0}$. Without loss of generality, we suppose $\phi_1\neq \hat{0}$, $\phi_2 = \hat{0}$. (ⅰ) If $\phi_1(0) > 0$, we have $\frac{dw_2(0)}{dt} > 0$ and then there is a small $t_0 > 0$ such that $w_2(t) > 0$ for $t\in(0, t_0)$. Since $\phi_1(0) > 0$, there is a $t_1\leq t_0$ such that $w_1(t) > 0$ for $t\in(0, t_1)$. Thus, $\Phi(t)\phi\in X_0$ for $t\in(0, t_1)$, a contradiction to the fact that $\phi\in M_{\partial}$. (ⅱ) If $\phi_1(0) = 0$ and $\phi_1(-r) > 0$, we see that

From continuity of solutions to (3.3) and $\phi\in C$, there is a $t_2 > 0$ such that

for $t\in[0, t_2]$. Then $w_1(t) > K_1t$, for $t\in[0, t_2]$. In addition,

where

is finite. By a comparison theorem,

for $t\in(0, t_2]$, also a contradiction to $\phi\in M_{\partial}$. (ⅲ) If $\phi_1(0) = 0$ and $\phi_1(-r) = 0$, we set $r^* = \sup\{ -\theta | \theta\in[-r, 0], \; \phi_1(\theta)\neq 0 \} < r$. Then $\frac{dw_1(t)}{dt} = 0$ and $w_1(t) = 0$ for $t\in[0, r-r^*]$. From the assumption $\phi_2 = \hat{0}$, it also holds that $w_2(t) = 0$ for $t\in[0, r-r^*]$. Since $\phi\in C$, there is a small $\epsilon_1 > 0$ such that $\phi_1(-r^*+\epsilon_1) > 0$. Define

then $\psi_2(\theta)\geq \hat{0}$ for $\theta\in[-r, 0]$ and

By the comparison principle (for $\psi_2(\theta)\geq \tilde{\psi}_2(\theta): = \hat{0}$) and previous result (for $\phi_1 = \psi_1$ and $\phi_2 = \tilde{\psi}_2$), it yields that $\psi\notin M_{\partial}$. Note that $\Phi(t)\psi = \Phi(t+r-r^*+\epsilon_1)\phi$. By the positive invariance of the set $M_{\partial}$, it leads to $\phi\notin M_{\partial}$, a contradiction again. From the contradictions in all three cases (ⅰ)-(ⅲ), we conclude that $M_{\partial} \subset \{ \hat{0} \}$, and hence $M_{\partial} = \{ \hat{0} \}$ indeed, proving the claim.

From ($\mathcal{A}$), there exist constants $0 < \epsilon_2 < 1$ and small $\rho_1 > 0$ such that for $0 < \rho < \rho_1$

For this $\epsilon_2$, there is a small $0 < \rho^* < \rho_1$ such that

Now, we claim that

Suppose, for the sake of contradiction, that there exist an initial condition $\phi\in X_0$ and a $t_3 > 0$ such that $|w_i(t)|\leq \rho^*$, $i = 1, 2$, for $t\geq t_3-r$. From (4.7), for $t\geq t_3$,

We consider an auxiliary equation

and its associated ordinary differential equation

From the comparison theory in [20,Theorem 5.5.1],

where $\tilde{\Phi}(t)$ is the semiflow of (4.9). In the ordinary differential equation (4.10), the characteristic equation of the related Jacobian matrix at the trivial equilibrium is

From (4.6), it admits two roots, $\lambda_1 < 0 < \lambda_2$, i.e. the stability modulus of (4.10) is positive. Also note that (4.9) is a cooperative irreducible system. From Theorem 5.5.1 and Corollary 5.5.2 in ^[20], the system (4.9) also admits a positive stability modulus associated with a positive eigenvector $u$. Note that the semiflow of (3.3) is eventually strong monotone in $[\hat{0}, \hat{\rho}^*]$, see [20,Corollary 5.3.5], and $\hat{0}$ is an equilibrium therein. There exist a $t_4 > t_3$ and a small $\alpha > 0$ such that

Hence, we have

for $t\geq 0$, which is a contradiction to boundedness of the semiflow $\Phi(t)$, and this contradiction proves (4.8).

Obviously, $\{\hat{0}\}$ is an isolated invariant set in $X$, and then the set $M_{\partial}$ consists with an acyclic equilibrium point. From (4.8), $W^s(\hat{0})\cap X_0 = \emptyset$, where $W^s(\hat{0})$ denotes the stable manifold of $\hat{0}$. By the persistence theory in [21,Theorem 4.6], the system (3.3) is uniformly persistent with respect to $(X_0, \partial X_0)$ and the assertion is proved.

Next, we suppose that ($\mathcal{B}$) holds, i.e. $\sigma\beta_i > d_{i, m}+D_{i, m}$ for some $i$. Then there is a positive constant $\rho_2$ such that

and

for $0\leq w_i(t) < \rho_2$. For this $i$, consider an auxiliary equation

From (4.11), it obviously admits an unique positive equilibrium $z^*$. According to $z^*\leq 1/\gamma_i$ or $z^* > 1/\gamma_i$, we refer to Proposition 3.2 or Theorem 3.5 in ^[13] respectively to derive the existence of positive constants $\rho_3 < \rho_2$ and $t_5$ such that $z(t) > \rho_3$ whenever $t > t_5$. The comparison principle in [20,Theorem 5.1.1] implies that $w_i(t)\geq z(t)\geq\rho_3$ whenever $0\leq w_i(t)\leq \rho_2$. Thus $\liminf_{t\rightarrow\infty}w_i(t)\geq\rho_3 = :\rho^*_i$. If, in addition, $\sigma\beta_j > d_{j, m}+D_{j, m}$ for $j\neq i$, then we have also $\liminf_{t\rightarrow\infty}w_j(t)\geq \rho^*_j$ and the assertion holds for $\rho^* = \min\{\rho^*_1, \rho^*_2\}$. If $\sigma\beta_j\leq d_{j, m}+D_{j, m}$ for $j\neq i$, then from the equation of $w_j$ in (3.3),

for $t > t_5$. Let

which is obviously strictly increasing in $w \in [0, \infty)$ and satisfies $\tilde{F}_j(w) = 0$ and $\tilde{F}_j(\infty) = \infty$. Thus, there is a unique $\tilde{\rho}^*_j > 0$ such that $\tilde{F}_j(\tilde{\rho}^*_j) = D_{i, m}f_i(\rho^*_i)\rho^*_i$. Hence, for $t > t_5$,

whenever $0\leq w_j(t) < \tilde{\rho}^*_j$. Therefore, $\liminf_{t\rightarrow\infty}w_j(t)\geq \tilde{\rho}^*_j$, and the assertion holds for $\rho^*: = \min\{\rho^*_i, \tilde{\rho}^*_j\}$.

To demonstrate the convergent dynamics in system (3.3) that involves the non-monotone feedback, we identify an compact set that is invariant and attracting for (3.3) and in which, (3.3) is monotone dynamics and hence the monotone dynamical system theory is applicable therein. For convenience of formulation, we denote $\tilde{w}_i = \frac{1}{\gamma_i}$ for $i = 1, 2$, $\tilde{w} = (\tilde{w}_1, \tilde{w}_2)$ and an order interval

Theorem 4.4. Consider the system (3.3). Suppose that either ($\mathcal{A}$) or ($\mathcal{B}$) holds. Suppose that $(w_1^*, w_2^*)\leq (\tilde{w}_1, \tilde{w}_2)$ and

where $\tilde{F}_i$ is defined in (4.12). Then the positive equilibrium, $(w_1^*, w_2^*)$ attracts all solutions in $X$, i.e. $\lim_{t\rightarrow\infty}w_i(t) = w_i^*$ for $i = 1, 2$.

Proof. Denote, for $l\geq \tilde{w}_1+\tilde{w}_2$,

see Figure 2. We claim that each $\Omega_l$ is positive invariant under (3.3). Except for those on the axes, there are three segments of the boundary of $\Omega_l$,

Denote the right side of (3.3) by $J = (J_1, J_2)^T$, where $T$ indicates the transpose. If $\phi\in\Omega_l$ and $\phi(0)\in B^l_1$, then we have

where the last inequality follows from (4.13). Similarly, the case that $\phi\in\Omega_l$ and $\phi(0)\in B^l_2$ leads to $J_2(\phi)\leq 0$. If $\phi\in\Omega_l$ and $\phi(0)\in B^l_3$, then we have

Combining the result in Lemma 3.1 with [20,Remark 5.2.1], $\Omega_l$ is positive invariant under the semiflow of (3.3).

We write $\tilde{l} = \tilde{w}_1+\tilde{w}_2$ and note that $\lim_{l\rightarrow \tilde{l}}\Omega_l = I$. Next, we prove that $I$ attracts each solution of (3.3). It suffices to show that the omega limit set of each solution is contained in $I$. From boundedness in Lemma 3.2, the solution to (3.3) with initial value $\phi$ has a nonempty, compact and positively invariant omega limit set $\omega(\phi)$. Define

Then, there exists a $\psi = (\psi_1, \psi_2)\in\omega(\phi)$ such that $\psi(0)\in B^{l_0}_1, B^{l_0}_2$ or $B^{l_0}_3$. Also note that $\varphi_i(\theta)\leq l_0-\tilde{w}_j$ for $1\leq i\neq j\leq j$, and $\varphi_1(\theta)+\varphi_2(\theta)\leq l_0$ for each $\varphi = (\varphi_1, \varphi_2)\in\omega(\phi)$. From the invariance of $\omega(\phi)$, there exists a $\xi = (\xi_1, \xi_2)\in\omega(\phi)$ such that $w_r(\xi) = \psi$. Now, suppose that $l_0 > \tilde{l}$. If $\psi(0)\in B^{l_0}_1$, it holds that

where the strict inequality follows from $l_0 > \tilde{l}$. It implies a $t_1 < r$ such that

which contradicts to $w_{t_1}(\xi)\in\omega(\phi)\subset\Omega_{l_0}$. Similarly, the case $\psi(0)\in B^{l_0}_2$ also leads to a contradiction. If $\psi(0)\in B^{l_0}_3$, it holds that

It implies a $t_2 < r$ such that

which contradicts to $w_{t_2}(\xi)\in\omega(\phi)\subset\Omega_{l_0}$. Therefore, the omega limit set $\omega(\phi)\in \Omega_{\tilde{l}} = I$.

From Theorem 4.3, each omega limit set $\omega(\phi)$ is in fact contained in $I^*: = [\hat{\rho}^*, \hat{\tilde{w}}]$ for some $\rho^* > 0$. In addition, we see that $(w_1^*, w_2^*)$ is the unique equilibrium of (3.3) in $I^*$ and this system is cooperative and irreducible over the space $I^*$. By [20,Proposition 5.4.2], the unique equilibrium $(w_1^*, w_2^*)$ attracts all solutions in $I^*$, and therefore, it attracts all solutions in $X$.

When $(w_1^*, w_2^*) \leq (\tilde{w}_1, \tilde{w}_2)$ does not hold, numerical simulations (see Section 5) show that it is possible for (3.3) to have oscillatory solutions. To conclude this section, we establish some estimates for the dynamics of (3.3) for the case $(w_1^*, w_2^*) > (\tilde{w}_1, \tilde{w}_2)$, and identify an invariant set for (3.3) for this case. To this end, we first have the following lemma.

Lemma 4.3. When $(w_1^*, w_2^*) > (\tilde{w}_1, \tilde{w}_2)$, the system

has a unique positive solution for $(p_1, q_1, p_2, q_2)$. Moreover, there holds the dichotomy: either $p_i < w_i^* < q_i$ for $i = 1, 2$, or $p_i = q_i = w_i^* = \tilde{w}_i$ for $i = 1, 2$.

Proof. Consider the equation

which is equivalent to

Substituting $y = (M-d_{1, m}x)/d_{2, m}$ into (4.15) leads to the equation

From the facts that $P(0) = -G_2(\frac{M}{d_{2, m}})-\sigma b_1(\tilde{w}_1) < 0$, $P(\frac{M}{d_{1, m}}) = \tilde{F}_1(\frac{M}{d_{1, m}})-\sigma b_1(\tilde{w}_1) > 0$ and the continuous function $P$ is strictly increasing on $[0, \frac{M}{d_{1, m}}]$, there is a unique positive root satisfying $P(q_1) = 0$, and thus there are unique $q_1$ and $q_2 = (M-d_{1, m} q_i)/d_{2, m}$ satisfying the two left equations in (4.14). A similar argument shows that the two right equations in (4.14) has a unique positive solution $(p_1, p_2)$.

If $w_i^* > q_i$ for $i = 1, 2$, we derive

which is a contradiction. If $w_i^* > q_i$ and $w_j^*\leq q_j$ for some $i$ and $j\neq i$, then strict monotonicity of functions $\tilde{F}_i$ and $G_j$ implies

It contradicts to the fact from (4.14),

Therefore, $w_i^*\leq q_i$ for $i = 1, 2$.

To show $p_i\leq w_i^*$ for $i = 1, 2$, define the strictly increasing functions

Then

The strict monotonicity of $Q_i$ implies $p_i\leq w_i^*$. From previous argument, we see that either $p_i < w_i^* < q_i$ for $i = 1, 2$; or $p_i = q_i = w_i^* = \tilde{w}_i$ for $i = 1, 2$.

We now show that the constants $p_i, q_i$ in (4.14) can help us estimate all omega limit sets and identify a positive invariant set for (3.3).

Theorem 4.5. Assume that either ($\mathcal{A}$) or ($\mathcal{B}$) holds and $(w_1^*, w_2^*) > (\tilde{w}_1, \tilde{w}_2)$. Let $p = (p_1, p_2)$ and $q = (q_1, q_2)$ be as in Lemma 4.3. Then for each $\phi \in C_+$, $\omega(\phi)\leq \hat{q}$; moreover, if $b_i(q_i)\leq b_i(p_i)$ for $i = 1, 2$, then the set $[\hat{p}, \hat{q}]$ is positively invariant under the semiflow of (3.3).

Proof. When $w_i(t) > q_i$ and $w_j(t)\leq q_j$ for $1\leq i\neq j\leq 2$,

When $w_i(t) > q_i$ for $i = 1, 2$,

Therefore, it holds that $\omega(\phi)\leq \hat{q}$ for all $\phi\in C$.

Suppose that $b_i(q_i)\leq b_i(p_i)$ for $i = 1, 2$, let $\phi\in[\hat{p}, \hat{q}]$ and $w(t) = (w_1(t), w_2(t))$ be the corresponding solution with initial function $\phi$. If there is a $t_1 > 0$ such that $w_i(t_1) = p_i$, then

If there is a $t_2 > 0$ such that $w_i(t_2) = q_i$, then

By [20,Theorem 5.2.1 and Remark 5.2.1], $[\hat{p}, \hat{q}]$ is positively invariant under the semiflow of (3.3).

5.   Summary and discussion

In this study, we incorporate adaptive dispersal, in the form of density dependent dispersal, into a two-patch population model with a maturation delay derived in ^[4] which assumes constant dispersion rates between patches, leading to a more realistic model. The improved model turns out to be a system of delay differential equation with spatial non-local birth terms resulted from the dispersals of the immature individual. The density dependent dispersals not only affect the instantaneous migrations of the mature population, they also have an impact on the nonlocal birth terms. For mathematical tractability, we have analyzed the case when the immatures only have constant dispersal rates, adopting the Ricker's birth function and a particular form for the density dependent dispersal rate functions. We have addressed the well-posedness of the model, structure of equilibria, threshold dynamics (in the sense of extinction and persistence) in terms of the conditions ($\mathcal{A}$) and ($\mathcal{B}$) (see Theorems 4.1-4.3). When the population is persistent, under a condition accounting for a monotone scenario (i.e., ($w_1^*, w_2^*) < (\tilde{w}_1, \tilde{w}_2$)), we have shown that the positive equilibrium is globally attractive; while when this condition is not satisfied (i.e., $(w_1^*, w_2^*) > (\tilde{w}_1, \tilde{w}_2)$), we have established an upper bound for omega limit sets of solutions and identified an invariant set within which, periodic oscillations may occur.

We remark that surprisingly the threshold conditions $\mathcal{A}$ and $\mathcal{B}$ for extinction/persistence do not depend on the $\eta_i$'s. However, $\eta_i$'s that reflect the level of density dependent dispersals do affect the value of the positive equilibrium ($w_1^*, w_2^*$), and accordingly, affect whether or not the condition $(w_1^*, w_2^*) \le (\tilde{w}_1, \tilde{w}_2)$ holds. Thus, $\eta_i$'s do have an impact on the long time dynamics of the solutions in the case of persistence. Below we numerically explore such an impact. To this end, we fix $d_{1, m} = 0.4, D_{1, m} = 0.6, d_{2, m} = 0.2, D_{2, m} = 0.6, \beta_1 = 7, \gamma_1 = 2, \beta_2 = 4, \gamma_2 = 2, d_I(a)\equiv 0.05$ and $r = 6$. Then, by choosing various combinations for $\eta_1$ and $\eta_2$, we may observe different persistent dynamics. Figure 3 shows the results for three sets of ($\eta_1, \eta_2$) values, revealing that they may cause stability switch of the positive equilibrium.

Since the model is a DDE system, it is not surprising that the large delay $r$ may destroy the stability of an equilibrium leading to stable periodic solutions with the period depending on the value of $r$. This is demonstrated in Figure 4, in which $d_{1, m} = 0.45, D_{1, m} = 0.65, d_{2, m} = 0.35, D_{2, m} = 0.55, \beta_1 = 8, \gamma_1 = 2, \beta_2 = 7, \gamma_2 = 2, d_I(a)\equiv 0.05, \eta_1 = \eta_2 = 1$ and $r = 4, 10$ respectively. We point out that although the amplitudes of the oscillations change as $r$ varies, these periodic solutions are contained in the invariant set $[p, q] = [p_1, q_1] \times [p_2, q_2]$ (see Section 3) which depends on the value of $r$ since the value of $\sigma$ in (4.14) depends on $r$.

It was discovered both in a two-patch discrete time competitive model ^[22] and a diffusion model for multiple competitive species ^[23], with random dispersal in spatially varying but temporally constant environments, that lower dispersal rates are beneficial to the population persistence. Our adoption of a density-dependent dispersal assumes a relatively higher migration rate than that of the random dispersal. This naturally raises a question as to whether or not density-dependent dispersal is advantageous for population's persistence in a single species model. Model (3.3) contains a non-monotone birth function in (3.2), potentially yielding a result different from that of competition models. As remarked above, according to Theorems 4.2 and 4.3, the levels of density dependence of the dispersal ($\eta_i$'s) do not affect the threshold conditions for extinction/persistence. However, when it persists, the values of $\eta_i$'s do play a significant role in affecting the final size of population, as the value of positive equilibrium $(w_1^*, w_2^*)$ depends on $\eta_i$'s. In Figure 5, we present some numeric results on the impact of $\eta_i$'s on $(w_1^*, w_2^*)$, with choices of balanced adaptive dispersals ($\eta_1 = \eta_2$) and imbalanced adaption ($\eta_1 \neq \eta_2$), revealing diverse consequences. In Figure 5-(a)-(b), we consider two identical patches by assuming $D_{1, m} = D_{2, m} = 0.2, d_{1, m} = d_{2, m} = 0.1, \beta_1 = \beta_2 = 1, \gamma_1 = \gamma_2 = 1$. Take $d_I(a)\equiv 0.1$ and $r = 10 \ln 5$ leading to $\sigma = 0.2$. Then, we can see in this case that balanced adaption does not influence the final sizes of population; and for the distribution among two patches while imbalanced adaption (with adaption in patch-1 and without adaption in patch-2) not only changes the distribution of population but also decreases the final total population ($w_1^*+w_2^*$) as adaption ($\eta_1$) is increased. In Figure 5-(c)-(d), we set the same parameter values except $d_{2, m} = 0.3$ and $\beta_1 = 2$ which means that patch-1 has lower mortality and higher reproduction and is thus a better habitat than patch-2. In Figure 5-(c), even if two groups adopt the same dispersion strategy, the final population decreases in (better) patch-1, so does the total population with respect to $\eta_1$. The imbalanced adaption leads to a similar but worse consequence. In contrast, when patch-1 is a relative poorer habitat (with the same parameters as in (a) except $d_{1, m} = 0.3$ and $\beta_2 = 2$), although Figure 5-(e) shows that the balanced adaption is also disadvantageous to the final population, the imbalanced strategy with only adaption in patch-1, on the other hand, does help boost the species's prosperity (total population) shown in Figure 5-(f).

The intrinsic age-structured dispersal rate, $D_i(a)$, varies among organisms. The general case involving a non-constant dispersal rate at each life stage is required for a population to exhibit its sensitive age structure. Our original model incorporating both nonlinear dispersal terms and highly complex delayed birth terms is a biologically meaningful and a mathematically intractable problem. In addition, the derived model with an Allee birth function, for example $b_i(w) = \delta_i w^2e^{-\varepsilon_i w}$ where $\delta_i$ and $\varepsilon_i$ are constants, can involve the dynamics of bi-stability, and the convergent dynamics in such a non-monotone feedback DDE is also interesting but very challenging. We will address these questions in future research.

Acknowledgments

C. Y. Cheng and S. S. Chen were partially supported by the Ministry of Science and Technology, Taiwan, R.O.C. (Grant No. MOST 106-2115-M-153-003 and 105-2115-M-003-008); X. Zou was partially supported by NSERC of Canada (RGPIN-2016-04665).

Conflict of interest

All authors declare no conflicts of interest in this paper.