Modeling the effects of density dependent emigration, weak Allee effects, and matrix hostility on patch-level population persistence

Abbreviations: DDE: Density dependent emigration; DIE: Density independent emigration; +DDE: Positive relationship between density and emigration; −DDE: Negative relationship between density and emigration; UDDE: U-shaped relationship between density and emigration; hDDE: Hump-shaped relationship between density and emigration; AER: Allee effect region

1.   Introduction

1.1.   Background and motivation

Allee effects, the positive effects of increasing density on fitness, were first described in the early 1930s for cooperatively breeding species ^[1,2]. Scarcity of reproductive opportunities at low densities are thought to be a common cause for Allee effects such that there is a population threshold below which the population will go extinct ^[3,4]. However, Allee effects can also arise when increasing density facilitates feeding, habitat modification, defense against predators, or reduces the likelihood of inbreeding depression ^[5,6,7]. Empirical support for Allee effects spans a wide diversity of taxa ^[6,7], although it can be difficult to detect (e.g., ^[8]). The potential for positive density dependence in small populations to increase extinction risk has made Allee effects an important subject of study in conservation biology ^[6], fisheries management ^[9], invasion biology ^[10,11] and insect pest management ^[12,13,14]. Allee effects have also been recognized as important in tumor growth and epidemiology ^{[15,16,17,18,19]}.

Allee effects are particularly important in the context of metapopulation or landscape ecology and there is a growing list of studies that have examined the interplay between Allee effects and dispersal (for recent review, ^[20]). Under some circumstances, dispersal can increase persistence times in metapopulations with Allee effects ^[21]. More often, elevated dispersal rates exacerbate the negative influences of Allee effects on local or regional persistence ^{[22,23,24,25]} or the spatial spread of populations ^[12,13,26].

The relationship between conspecific density and the probability of emigrating from a patch can play an essential role in determining the population-dynamic consequences of an Allee effect. In a population subjected to Allee effects at the metapopulation-level, positive density-dependent emigration (+DDE), whereby individuals remain in patches at low density (Figure 1), will allow for Allee mechanisms to operate ^[22]. In contrast, negative density-dependent emigration (-DDE; Figure 1) would allow for individuals to escape Allee effects ^[27,28]. Although the most widely accepted view of emigration behavior is that species should exhibit +DDE ^[29,30,31], other forms of density-dependent emigration (DDE), including -DDE, exist (Figure 1). In a recent review of the empirical literature in ^[32], the authors found that $35\%$ of the cases exhibited +DDE, $30\%$ were density independent (DIE), $25\%$ were -DDE, $6\%$ were U-shaped (UDDE) and $4\%$ were humped shaped (hDDE). Importantly, recent mathematical models have revealed that DDE forms with a negative slope (-DDE and UDDE) can also induce something similar to an Allee effect at the patch-level (a version of bi-stable population dynamics such that the trivial steady state and a positive steady state are both stable) ^[32,33,34]. Matrix hostility, which relates to the probability that individuals die when exiting the patch and entering the matrix was found to have little effect on the generation of Allee-like effects ^[32], but we do not know how it would affect stability in the presence of an Allee effect. To date, no studies have examined how weak Allee effects (i.e., those that, on their own, do not result in a critical population size below which persistence is not possible ^[35]; Figure 2), the form of DDE and matrix hostility interact to affect population stability.

In this paper, we explore the effects of density dependent emigration on a population that exhibits a demographic weak Allee effect. In particular, we analyze a reaction diffusion model for a population with weak Allee effect per-capita growth inside a patch surrounded by a hostile matrix with boundary conditions that correspond to the situation where the probability of an individual remaining in the patch is a function of conspecific density at the point where the individual reaches the patch boundary. We do not consider a strong Allee effect in our analysis since previous studies have shown that such a per-capita growth rate will generate a patch-level Allee effect (i.e. bi-stable population dynamics) for all patch sizes (see, for example, ^[36]). Figure 3 shows dynamical differences between a population with a patch-level Allee effect generated from either a weak Allee effect per-capita growth rate with DIE or a logistic per-capita growth rate with -DDE and one with no patch-level Allee effect. In the Allee effect case, the population will exhibit a patch-level Allee effect for patch sizes with corresponding $\lambda$-values ($\lambda$ here is a composite parameter proportional to patch size squared) in what we will call the Allee effect range, i.e. $(\lambda_m, E_1)$ (Figure 3). The Allee effect range (AER) can be used to quantify the strength of a patch-level Allee effect in the sense that a large AER will indicate that a population will exhibit a patch-level Allee effect over a wide range of patch sizes. In contrast, a population with a small AER may theoretically exhibit a patch-level Allee effect, but practically will not exhibit such an Allee effect since the proportion of patches with a size in the AER will be almost zero.

We are interested in the interaction/interplay of a weak Allee effect per-capita growth rate (i.e. an initial positive relationship between density and per-capita growth rate), form of the DDE relationship (e.g. +DDE, -DDE, UDDE and hDDE), and matrix hostility. Particularly, 1) will a certain form of DDE enhance or even diminish the patch-level Allee effect, as measured by the AER, 2) will the interaction between the weak Allee effect per-capita growth rate and a particular DDE form create other versions of dynamical bi-stability besides a patch-level Allee effect, and 3) what effect will DDE have on the structure of the steady state density profile?

1.2.   Model formulation

We will present and study a model built upon the reaction diffusion framework which will incorporate a weak Allee effect per capita growth rate and density dependent boundary conditions. Reaction diffusion models have been extensively studied in the literature, see ^{[37,38,39,40,41,42,43]} and references therein for a detailed history of the reaction diffusion framework. In this model, $u(t, x)$ represents the normalized density (i.e. carrying capacity is equal to one) of a population in the patch $\Omega = (0, \ell)$ with patch size $\ell > 0$, surrounded by a hostile matrix, denoted by $\Omega_M = \mathbb{R}\setminus \overline{\Omega}$ (Figure 4). The boundary of $\Omega$ is denoted by $\partial \Omega$. Here, the variable $t$ represents time and $x$ represents spatial location within the patch. The model is then

where $D > 0$ is the patch diffusion rate, $r > 0$ is the patch intrinsic growth rate, $a \in (0, 1)$ is a parameter measuring the strength of the weak Allee effect (in the sense that per-capita growth rate is increasing for $u \in [0, \frac{1 - a}{2})$), $u_0(x)$ is the initial population density distribution in the patch, and $\alpha(u)$ denotes the probability of an individual remaining in the patch upon reaching the boundary. The term $\frac{\partial u}{\partial \eta}$ denotes the outward normal derivative of $u$. Here, the parameter $S^* \geq 0$ is a measure of the hostility of the matrix towards the organism, has units of length by time, and can assume different forms depending upon the patch/matrix interface assumptions ^[44]. If $\alpha(u) \equiv 0$ then the boundary is absorbing, i.e. all individuals that reach the boundary will emigrate, whereas if $\alpha(u) \equiv 1$ then the boundary is reflecting, i.e. the emigration rate is zero. Also, a given relationship between density and emigration can be included in the model by selecting an appropriate $\alpha(u)$ (see, for example, ^{[33,34,44,45,46]}).

We now introduce a standard scaling,

It is easy to see that after applying this scaling and dropping the tilde, (1.1) becomes

with corresponding steady state equation:

where $\lambda = \frac{r \ell^2}{D}$ and $\gamma = \frac{S^*}{\sqrt{r D}}$ are both unitless, $f(u) = \frac{1}{a} u (u + a) (1 - u)$, $g(u) = \frac{1 - \alpha(u)}{\alpha(u)}$, and $\Omega_0 = (0, 1).$ Thus, for fixed $r, D$, the composite parameter $\lambda$ is proportional to patch size squared and $\gamma$ to matrix hostility.

We next choose prototypical functions for the five most common DDE forms reported in the recent literature review in ^[32]. In order to remain consistent in choosing these forms, we employ a single $\alpha(u)$ template and it's mirror image, namely

where $M_1 > 0$ and $m(u) \geq 0$ with $m(0) = 0$ are appropriately chosen to model a given DDE form. Note that the emigration rate at zero will be the same across all forms, i.e. $1 - \alpha_i(0) = 0.5$, $i = 1, 2$. Table 1 lists the exact $m(u)$'s that were chosen to model the five DDE forms (also, Figure 1).

1.3.   Structure of the paper

We will present some preliminary mathematical results in section 2. An evolution of the structure of positive steady states of (1.3) as $\gamma$ is varied is given in section 3, followed by an analysis of the AER in section 4. Section 5 covers the affect of DDE on density profile of the steady states of (1.3). Finally, we discuss some consequences of our results in section 6.

2.   Preliminaries

In this section, we state and prove several mathematical results that will aid in the study of the model (1.3). First, we consider the stability of the trivial steady state, $u(x) \equiv 0$, of (1.3). Let $E_1(\gamma)$ be the principal eigenvalue of the boundary value problem:

We now obtain the following theorem which connects $E_1(\gamma)$ to the stability of $u(x) \equiv 0$.

Theorem 2.1. The trivial solution of (1.4) is asymptotically stable if $\lambda < E_1(\gamma)$, and it is unstable if $\lambda > E_1(\gamma)$.

Before presenting a proof of Theorem 2.1, we recall the following results from ^[34,47].

Lemma 2.1. ^[47] Let $\sigma_1$ be the principal eigenvalue of the linearized equation associated with (1.4), namely

where $u$ is any solution of (1.4). Then the following hold.

$a)$ If $\sigma_1 > 0$, then $u$ is stable. Furthermore, if $u$ is isolated then it is asymptotically stable.

$b)$ If $\sigma_1 < 0$, then $u$ is unstable.

Lemma 2.2. ^[34] Let $u$ be a solution of (1.4) and $\sigma_1^*$ be the principal eigenvalue of the following boundary value problem

Then, $sign(\sigma_1^*) = sign(\sigma_1)$ for $\sigma_1^*, \sigma_1 \neq 0.$

In the light of Lemma 2.2, it suffices to study the relationship between $\sigma_1^*$ and $\lambda$ in order to prove Theorem 2.1.

Proof of Theorem 2.1:

Let $\lambda < E_1(\gamma)$. By Lemma 2.2, we see that the zero solution is asymptotically stable if the principal eigenvalue $\sigma_1^{*}$ of (2.3) with $u \equiv 0$ is positive. Note that, for $\lambda < E_1(\gamma)$, the zero solution is isolated since $\lambda$ is not a bifurcation point on the solution curve $(\mu, 0)$. Let $\mu_1 = \mu_1(\beta)$ be the principal eigenvalue of:

where $\beta \geq 0.$ Then, $\mu_1(\beta)$ is a strictly increasing concave function which passes through the origin, bounded above by $\lambda_D$, the principal eigenvalue of the Dirichlet problem (see ^[44,48]):

Let $\beta = \gamma\sqrt{\lambda}g(0)$. Since $\mu_1(\beta)$ is a strictly increasing concave function of $\beta$ and $\frac{\beta^2}{\gamma^2 g(0)^2}$ is a strictly increasing convex function of $\beta$ which passes through the origin they intersect at exactly two points, namely at $(0, 0)$, and say at $(\beta^*, \mu_1(\beta^*))$ for $\beta^* > 0$ (Figure 5). From (2.1), we can easily see that $\mu_1(\beta^*) = E_1(\gamma)$ and $\beta^* = \gamma \sqrt{E_1(\gamma)}g(0)$. Further, $\lambda+\sigma_1^{*} = \mu_1(\gamma \sqrt{\lambda}g(0))$, where $\sigma_1^{*}$ is the principle eigenvalue of (2.3). Thus, if $\lambda < E_1(\gamma)$ then $\gamma \sqrt{\lambda}g(0) < \beta^*$ and $\mu_1(\gamma \sqrt{\lambda}g(0)) > \lambda$, implying $\sigma_1^{*} > 0$. By Lemma 2.2 the zero solution is asymptotically stable if $\lambda < E_1(\gamma)$.

Next, let $\lambda > E_1(\gamma)$. By Lemma 2.2, the zero solution is unstable if the principle eigenvalue $\sigma_1^{*}$ of (2.3) is negative. But when $\lambda > E_1(\gamma)$, $\gamma \sqrt{\lambda}g(0) > \beta^*$ and $\mu_1(\gamma \sqrt{\lambda}g(0)) < \lambda$ implying $\sigma_1^{*} < 0$ (Figure 5). Hence, Theorem 2.1 is proven.

The second result gives a sufficient condition for the model (1.3) to exhibit a patch-level Allee effect which only requires knowledge of the existence of a positive steady state of (1.3) and not its stability properties.

Lemma 2.3. Let $\gamma > 0$ and $a \in (0, 1)$ be given. If (1.3) has at least one positive steady state for $\lambda < E_1(\gamma)$ then the model (1.3) will exhibit a patch-level Allee effect if the patch size is $\ell = \sqrt{\frac{\lambda D}{r}}$.

Proof of Lemma 2.3: Assume that $\gamma > 0$, $a \in (0, 1)$, and $\lambda < E_1(\gamma)$, are given and $u_1(x)$ is a positive solution of (1.4). By Theorem 2.1, the trivial steady state, $u(x) \equiv 0$, of (1.3) is asymptotically stable. Since $f(s) < 0$ for all $s > 1$, any constant $M \geq 1$ is a supersolution for (1.4) and a strict supersolution if $M > 1$ (see Definition 4.1 of ^[49]). Thus, any positive solution, $u(x)$, of (1.4) must satisfy $0 < u(x) < 1$ for $x \in [0, 1]$. Now, since $u_1(x)$ is a positive solution of (1.4) then it is also a subsolution and satisfies $u_1(x) \leq 1$. For any $u_0(x)$ such that $u_1(x) \leq u_0(x) \leq 1$ for $x \in (0, 1)$, Theorem 6.6 of ^[49] guarantees that the solution of (1.3), $u(t, x)$, with $u(0, x) = u_0(x)$ for $x \in (0, 1)$ must satisfy $0 < u_1(x) < u(t, x) < 1$ for all $x \in [0, 1], t \geq 0$. It is now clear that the model (1.3) will predict extinction for initial population densities, $u_0(x)$, with $\|u_0\|_\infty \approx 0$, whereas the model will predict persistence for $u_0(x)$ satisfying $u_1(x) \leq u_0(x) \leq 1$ for $x \in (0, 1)$. This establishes a patch-level Allee effect proving the lemma.

The final result is referred to as a time map analysis and will allow study of the structure of positive steady state solutions of (1.3) as the composite parameters $\lambda$ and $\gamma$ vary.

Theorem 2.2. A positive solution, $u(x)$, of (1.4) with $\rho = ||u||_\infty$, $n = u(0)$, and $q = u(1)$ exists if and only if $\lambda > 0$, $\rho \in (0, 1)$, and $n, q \in [0, \rho)$ satisfy:

and

where $F(s) = \int_0^s f(t)dt$.

Remark 2.1. For $\rho \in (0, 1)$, since $f(\rho) > 0$, it can be shown that the improper integral in (2.4) is convergent.

Figure 6 illustrates a prototypical positive solution of (1.4). We now provide a proof of Theorem 2.2.

Proof of Theorem 2.2: Assume that $u(x)$ is a positive solution to (1.4) with $\rho: = ||u||_\infty$, $n: = u(0)$, $q = u(1)$. Since (1.4) is an autonomous differential equation, if there exists an $x_0 \in (0, 1)$ such that $u'(x_0) = 0$ then $v(x) : = u(x_0 + x)$ and $w(x) : = u(x_0 - x)$ will both satisfy the initial value problem,

for all $x \in [0, d)$ with $d = \min\{x_0, 1 - x_0\}$. Picard's Existence and Uniqueness Theorem asserts that $u(x_0 + x) \equiv u(x_0 - x)$. Hence, $u(x)$ must be symmetric about $x_0$, $u'(x) \geq 0;\ [0, x_0]$, and $u'(x) \leq 0;\ [x_0, 1]$.

Multiplying both sides of (1.4) by $u'$ we obtain

and integrating both sides gives

Substituting $x = x_0$, $x = 0$, and $x = 1$ into (2.8) gives

Combining (2.9) with (2.10) and (2.11) we have,

Now substitution of (2.9) into (2.8) yields

Solving for $u'(x)$ in (2.14) and using the fact that $u'(x) > 0;\ [0, x_0)$ and $u'(x) < 0;\ (x_0, 1]$ we have

Integration of (2.15) from $0$ to $x$ and (2.16) from $x_0$ to $x$ yields

Through a change of variables and using the fact that $u(0) = q$ and $u(x_0) = \rho$ we have

Substituting $x = x_0$ into (2.19) and $x = 1$ into (2.20) gives

Now subtraction of (2.22) from (2.21) yields,

Next, assume $\lambda > 0, \; \rho \in (0, 1)$, and $n, q \in [0, \rho)$ satisfy (2.4) and (2.5). Define $u(x):[0, 1] \rightarrow \mathbb{R}$ by

We will now show that $u(x)$ is a positive solution to (1.4). It is easy to see that the turning point given by $x_0 = \frac{1}{\sqrt{2\lambda}}\int_n^{\rho}\frac{dt}{\sqrt{F(\rho) - F(t)}}$ is unique for fixed $\lambda$-, $n$-, and $\rho$-values. The function,

is a differentiable function of $u$ which is strictly increasing from $0$ to $x_0$ as $u$ increases from $n$ to $\rho$. Thus, for each $x \in [0, x_0]$, there is a unique $u(x)$ such that

Moreover, by the Implicit Function theorem, $u(x)$ is differentiable with respect to $x$. Differentiating (2.26) gives,

Through a similar argument, $u(x)$ is a differentiable, decreasing function of $x$ for $x \in (x_0, 1)$ with

This implies that we have,

Differentiating again, we have,

Thus, $u(x)$ satisfies the differential equation in (1.4). It only remains to be seen that $u(x)$ satisfies the boundary condition in (1.4). However, from (2.24) and (2.25) it is clear that $u(0) = n$ and $u(1) = q$. Since $n$ is a solution of (2.12), we have

Substituting $x = 0$ into (2.27) gives,

Combining (2.29) and (2.30) we have,

A similar argument shows that

Hence, $u(x)$ satisfies (1.4) and the proof is complete.

3.   Structure of positive steady states of (1.3) as $\gamma$ varies

In this section, we will present numerically generated bifurcation curves of positive steady states of the model (1.3) and some corresponding biological implications. We employ (2.4) and (2.5) from Theorem 2.2 to numerically obtain bifurcation curves via Mathematica (Wolfram Inc., ver. 12.0). First, we briefly discuss our methodology. Fixing $a \in (0, 1)$ and $\gamma > 0$ we let $x_i = \frac{i}{n+1}; \; i = 1, ..., n$ for some $n\geq 1$. Setting $\rho = x_1$, we solve the Eq (2.5) for $n$ and $q$ using the $\texttt{FindRoot}$ command in Mathematica. The values of $n, q$ and $\rho$ are then substituted into (2.4) to find the corresponding value of $\lambda$. Repeating this procedure for $\rho = x_i$, $i = 2, ......, n$, we obtain $(\lambda, \rho)$ points generating a bifurcation diagram of $\lambda$ vs. $\rho = \|u\|_\infty$ for positive solutions of (1.4).

Recall from Lemma 2.3 that if there exists a range of $\lambda < E_1(\gamma)$ for which a positive solution of (1.4) exists then the model will predict an Allee effect at the patch-level for patch sizes corresponding to these $\lambda-$values. In this Allee effect case, the population density must surpass a certain threshold in order for persistence to be predicted. Since our growth rate $f(u)$ is taken to be of a weak Allee effect form, we would expect model predictions of an Allee effect at the patch-level in the case of a DIE. We are particularly interested in model predictions of bi-stability scenarios other than a patch-level Allee effect in the case of density dependent emigration. We will present an evolution of the bifurcation curves for all five DDE forms as $\gamma$ increases for the cases: 1) where the forms of DDE are relatively weak and parameter values are $a = 0.5$, $M_1 M_2 = 1.1$ and $M_3 = 0.5$ (Figure 7), and 2) where the forms of DDE are relatively strong and pronounced with parameter values $a = 0.5$, $M_1 M_2 = 0.08$ and $M_3 = 0.25$ (Figure 8). In both cases, an $a$-value of $0.5$ gives a substantial weak Allee effect, i.e. the per-capita growth rate will increase for $u$-values in $[0, 0.25)$. Note that presentation of an exploration of the entire parameter space would be quite challenging and is outside of the scope of this work.

As shown in both Figures 7 and 8, the bifurcation curves' starting value, $E_1(\gamma)$, satisfies $E_1(0) = 0$, $E_1(\gamma)$ is strictly increasing in $\gamma$, and $E_1(\gamma) \rightarrow \pi^2$ as $\gamma \rightarrow \infty$ (see ^[48] or ^[50], for example). The positive relationship between density and emigration probability in +DDE and initially in hDDE cause the maximum steady state values of these two forms to be much less than the DIE case, whereas the negative relationship in $-$DDE and initially in UDDE cause an increase in maximum steady state values as compared with the DIE form. The difference in maximum steady state values appears to be greatest for intermediate values of $\gamma$ and the least when $\gamma$ is large. Notice that as $\gamma \rightarrow \infty$, i.e. when the matrix is completely hostile, the +DDE, $-$DDE, UDDE, and hDDE curves all converge to the DIE form as illustrated in Figures 7c and 8c. A patch-level Allee effect is present in all values of $\gamma$ for Figure 7, but the initial positive relationship between density and emigration probability of hDDE is able to completely counteract the patch-level Allee effect in Figure 8a, even though the +DDE case does not. This discrepency is due to the positive relationship being much stronger (at least initially) in the hDDE case versus the +DDE case, as shown in Figure 8d. In Figure 7, the only bi-stability of steady states predicted by the model is the aforementioned patch-level Allee effect. In contrast, Figure 8 shows examples of other types of bi-stability in the case of hDDE in (a) and -DDE in (b). Though not shown here, a similar non-Allee effect bi-stability also appeared in the UDDE case for the same parameter values in Figure 8. In fact, any S-shaped bifurcation curve (or even a more complicated shape) occurring for $\lambda > E_1(\gamma)$ will not qualify as an Allee effect since by Theorem 2.1, the trivial steady state, $u(x) \equiv 0$, is unstable. In both cases, model predictions of persistence vary over a wide range as the matrix hostility, as measured in the composite parameter $\gamma$, varied.

4.   Allee effect region length

In this section, we explore the relationship between DDE form and the strength of the patch-level Allee effect predicted by the model (1.3). In order to accomplish this, we study the length of the AER, defined as $E_1(\gamma) - \lambda_m(\gamma)$, for fixed values of $M_1, M_2, M_3,$ and $a$ (Figure 3). We calculate $\lambda_m(\gamma)$ by employing Theorem 2.2 and Mathematica (Wolfram Inc., ver. 12.0) to numerically generate the bifurcation curve of positive solutions of (1.4) (see section 3) for a fixed set of parameters. The smallest $\lambda$-value on the curve is then $\lambda_m$. Using the Mathematica command $\texttt{NDEigensystem}$, we numerically estimate the value of $E_1(\gamma)$. If $\lambda_m(\gamma) < E_1(\gamma)$ then for $\lambda \in (\lambda_m(\gamma), E_1(\gamma))$, there is at least one positive solution, and by Lemma 2.3 the model predicts a patch-level Allee effect. However, if $\lambda_m(\gamma) \geq E_1(\gamma)$ then no such patch-level Allee effect can exist, since by Theorem 2.1, the trivial steady state is unstable for $\lambda \geq E_1(\gamma)$. In what follows, we will first compare the length of the AER for all the DDE forms given in Table 1 and then explore the possibility of the +DDE form counteracting a patch-level Allee effect.

4.1.   Qualitative connection between AER length and DDE form

Choosing $M_3 = 0.25$ and $a = 0.5$, we computed the AER length for different $\gamma$-values for each of the five DDE forms. This choice of $a$ will ensure a substantial weak Allee effect, i.e. the per-capita growth rate will increase for $u$-values in $[0, 0.25)$, whereas, $M_3 = 0.25$ will cause the minimum and maximum emigration probabilities to occur at $u = 0.25$ for UDDE and hDDE, respectively. We evaluated many other parameter values for $M_3$ and $a$ but obtained similar results. Although a full exploration of the entire parameter space is outside the scope of this work, we aim to provide a qualitative connection between the form of DDE and length of AER as the matrix hostility is varied via $\gamma$. Figures 9–11 illustrate this connection for $M_1 M_2 = 0.1, 0.5,$ and $1$. These $M_1 M_2$-values produce DDE forms that are somewhat different from DIE when $M_1 M_2 \approx 0$ to almost identical to DIE when $M_1 M_2$ is large.

In all three cases of $M_1 M_2$-values, the model always exhibited a patch-level Allee effect in the DIE, +DDE, $-$DDE, and UDDE cases. Also, when $\gamma$ is large, the length of the AER is virtually identical to DIE across all DDE forms. The AER length approached zero in all DDE forms and in all parameter choices as $\gamma$ approached zero. The +DDE form partially counteracted the patch-level Allee effect by slightly lowering the AER length for all $\gamma$-values, though for these parameter choices, the +DDE relationship was not strong enough to fully counteract the Allee effect. In contrast, the hDDE form, which is initially a positive relationship between density and emigration rate, was able to completely counteract the patch-level Allee effect for $\gamma$ approximately in $[1.5, 2.5]$ in Figure 9 and in $(0, 0.5]$ in Figure 10. This discrepancy between the +DDE and hDDE forms is due to the positive relationship in the hDDE being clearly stronger than the one in +DDE in both Figures 9 and 10. Due to switching from a positive relationship to a strong negative one, a patch-level Allee effect reappeared for $\gamma < 1.5$ for hDDE in Figure 9. However, this switch in the relationship in hDDE was not sufficient to allow the Allee effect to reappear in Figure 10.

In all three Figures, both $-$DDE and UDDE forms caused an increase in length of the AER as compared to the DIE case. In fact, in Figures 9 and 10, the AER length initially increased as $\gamma$ decreased but then began to decrease as $\gamma$ became small for both $-$DDE and UDDE, even boasting a peak value of almost four-times the DIE AER length in Figure 9. In Figure 11, all DDE forms had strictly decreasing AER length as $\gamma$ decreased. Interestingly, in Figure 9, the AER length for $-$DDE exhibited a steep increase from around one for $\gamma \approx 4$ to around four for $\gamma \approx 2.5$. A positive relationship between density and emigration rate at least partially counteracted the patch-level Allee effect for +DDE and hDDE forms, whereas the negative relationship enhanced the Allee effect for $-$DDE and UDDE forms. Also, this counteraction and enhancement of the patch-level Allee effect is dependent upon the hostility of the surrounding patch matrix, as measured by the parameter $\gamma$.

4.2.   Counteracting a patch-level Allee effect with +DDE

Our analysis of the structure of positive steady states for the model indicates that DDE forms containing a negative slope can increase the strength of the patch-level Allee effect as measured by the AER length, whereas, a positive slope can counteract the Allee effect. Even though both +DDE and hDDE have the potential to completely counteract a patch-level Allee effect for small patch sizes, the hDDE form's negative slope for $u > M_3$ will allow the Allee effect to reappear as the patch size approaches zero (Figure 9). Thus, we chose to focus on +DDE in an attempt to quantify when a patch-level Allee effect will be completely counteracted by a DDE relationship containing a positive slope. To that end, we again employed Theorem 2.2 and Mathematica (Wolfram Inc., ver. 12.0) to numerically generate bifurcation curves of positive solutions for (1.4) for fixed sets of parameter values. To establish the existence of a patch-level Allee effect in the +DDE case, it suffices to show that the slope of the bifurcation curve is negative for $\rho \approx 0$, i.e. we consider $\lambda = \lambda(\rho)$ ($\rho$ denotes the maximum steady state value) and numerically calculate $\lambda'(0)$. Figure 12 illustrates the parts of the parameter space for which a patch-level Allee effect is predicted by the model, i.e $\lambda'(0) < 0$, (Region Ⅰ) and parts where an Allee effect is not predicted, $\lambda'(0) > 0$, (Region Ⅱ) for the case of $a = 9.$ Notice that the boundary between Region Ⅰ and Ⅱ is comprised of the $M_1 M_2$- and $\gamma$-values such that $\lambda'(0) = 0$.

There is clearly a maximal $M_1 M_2$-value, such that for $M_1 M_2$ larger than this value the model will predict a patch-level Allee effect for all $\gamma > 0$. In contrast, it appears that for any $\gamma > 0$, there is always a small range of $M_1 M_2$-values such that no patch-level Allee effect is present.

Figure 13 compares the boundary curve separating parameter space into Region Ⅰ and Ⅱ for $a = \; 0.5, \; 0.75\;$ and $0.9$. Recall that $a \in (0, 1)$ measures the strength of the demographic weak Allee effect in the model via the per-capita growth rate. Thus, the demographic Allee effect varies from almost not present for $a \approx 1$ to substantial for $a \approx 0$ (Figure 13 (left)). Figure 13 shows that for smaller $a$-values, the +DDE response must become correspondingly stronger as indicated in the smaller $M_1 M_2$-values. Figure 14 illustrates this point for fixed $\gamma = 0.59275$ and $a = 0.75$, in which we compare the +DDE forms from Region Ⅰ and Ⅱ. Notice that for $M_1 M_2$-values that are sufficiently small (corresponding to solid curves in Figure 14) the positive relationship between density and emigration probability is strong enough to completely counteract the demographic Allee effect in the per-capita growth rate to produce no patch-level Allee effect. In contrast, the remaining dashed curves in Figure 14 represent +DDE forms that only partially counteract the patch-level Allee effect. Figure 15 further illustrates this point by comparing the actual bifurcation curves for +DDE forms belonging to Region Ⅰ (dashed) and Region Ⅱ (solid). Notice that, initially, the Region Ⅰ +DDE form bifurcation curves all decrease in $\lambda$ (i.e. $\lambda'(0) < 0$), while the Region Ⅱ +DDE form bifurcation curves increase in $\lambda$ (i.e. $\lambda'(0) > 0$), both as the maximum steady state value increases.

5.   Steady state density profile

In this section, we explore the effects of DDE on the density profile of positive steady states of (1.3). Figure 16 illustrates the two possible density profiles for a positive steady state of the model, namely asymmetric (whenever $n = u(0) \neq u(1) = q$) and symmetric (whenever $n = q$). Typically in a model with DIE, positive steady states would intuitively be symmetric since organisms would interact with both boundary points in a similar way and experience emigration into the hostile matrix. However, when DDE is present it is not clear when or even why asymmetric steady states could arise. To help understand the connection between DDE form and density profile of the model steady states, we state and prove sufficient conditions for all positive steady states of the model (1.3) to be symmetric, begining with a useful lemma.

Lemma 5.1. If $h(s) : = g(s)s$ is increasing for all $s > 0$ then every positive solution of (1.4) is symmetric about $x = \frac{1}{2}$.

Proof: Let $u(x)$ be a positive solution of (1.4) such that $n = u(0)$ and $q = u(1)$. From Theorem 2.2, $n, q$ must satisfy

Hence,

or equivalently,

Since $F(s)$ is increasing for $s > 0$, (5.1) can hold only if $n = q$, proving the lemma. We now state and prove the main result of this subsection.

Theorem 5.1. Let $m(s) \geq 0$ for $s \geq 0$.

$(a)$ If $\alpha(u) = \alpha_1(u) = \frac{M_1}{2M_1 + m(u)}$ then

$(i)$ If $m(s) \equiv 0$ (DIE) then all positive solutions of (1.4) are symmetric.

$(ii)$ If $m'(s) \geq 0$ (+DDE) then all positive solutions of (1.4) are symmetric.

$(iii)$ If $m(s) = \frac{s^2-2 M_3 s}{M_2}$ (UDDE) and $M_1 M_2 > \frac{4 M_3^2}{3}$ then all positive solutions of (1.4) are symmetric.

(b) If $\alpha(u) = \alpha_2(u) = \frac{M_1 + m(u)}{2M_1 + m(u)}$ then

$(i)$ If $m(s) = \frac{s^2}{M_2}$ ($-$DDE) and $M_1 M_2 > 1$ then all positive solutions of (1.4) are symmetric.

$(ii)$ If $m(s) = \frac{s^2-2 M_3 s}{M_2}$ (hDDE) and $M_1 M_2 > 1$ then all positive solutions of (1.4) are symmetric.

Proof: To prove (a), we first note that:

Thus, if $m'(s) \geq 0$ then we must have $h'(s) > 0$ for all $s > 0$, and (ⅰ) and (ⅱ) hold by Lemma 5.1. To show (ⅲ), we again calculate $h'(s)$

It is easy to see that if $4 M_3^2 -3 M_1 M_2 < 0$, or equivalently, $M_1 M_2 > \frac{4 M_3^2}{3}$ then we must have that $h'(s) > 0$ for $s > 0$ and hence (ⅲ) holds by Lemma 5.1.

To show (b), we calculate $h'(s)$ when $m(s) = \frac{s^2-2 M_3 s}{M_2}$ for $M_3 \geq 0$:

From the proof of Lemma 2.3, we have that $s < 1$. Hence, if $M_1 M_2 > 1$ then $h'(s) > 0$ for $s > 0$ and hence (b) holds by Lemma 5.1 and the result is proved.

Theorem 5.1 gives a fairly strong result for DIE and +DDE forms in that for any positive (or independent) relationship between density and emigration probability, asymmetric positive steady states of the model are not possible. Notice that the theorem gives results for the remaining DDE forms ($-$DDE, UDDE, and hDDE) that are dependent on the actual shape of the form via $m(s)$. Also, in these three cases where the DDE forms contain a negative relationship between density and emigration the result requires that $M_1 M_2$ not be too small. Note that in each of these three forms, placing a lower bound on $M_1 M_2$ restricts the minimum emigration probability. Since the $-$DDE and hDDE forms achieve their lowest emigration probability when $u = 1$, it is easy to see that $M_1 M_2 > 1$ gives that both forms' minimum emigration probability is bounded from below by $\frac{1}{3}$ for $-$DDE and $\frac{1}{3-2M_3}$ for hDDE. The UDDE form achieves its lowest emigration probability at $u = M_3$, and thus requiring $M_1 M_2 > \frac{4 M_3^2}{3}$ will ensure that the minimum emigration probability is bounded below by $\frac{1}{5}$. This analysis implies that asymmetric positive steady states for the model occur only when a negative relationship between density and emigration probability is present and sufficiently strong.

6.   Discussion

If at low densities a population exhibits a decline in its per-capita growth rate as density decreases then the population is said to exhibit an Allee effect. ^[1,6,51] provided a definition for Allee effects distinguishing between so-called component Allee effects, in which some component of individual fitness is in a positive relationship with either numbers or density of conspecifics, and demographic Allee effects, in which a positive relationship between total fitness and either numbers of density or conspecifics exists. They also mentioned that causes for an Allee effect may be dependent on the spatial scale of consideration for the population. In this article, we have been concerned with demographic Allee effects at the scale of a patch. In the case of a sufficiently strong demographic Allee effect, small populations can be expected to go extinct while larger populations expected to persist. In the context of mathematical models, this phenomenon corresponds to a bi-stability situation where both a positive population density and the zero population density are both stable (see ^[33], for example). We denote the range of patch sizes for which this type of bi-stability is predicted by the model as the Allee effect range and employ AER length as a metric to quantify overall strength of patch-level Allee effect.

In this article, we analyzed a model built on the reaction diffusion framework for a single species in a patch of habitat with a nonlinear boundary condition designed to model density dependent emigration which occurs in various forms in nature ^[32]. This boundary condition accounted for loss due to hostility of the matrix surrounding the patch and an edge-mediated increase or decrease to emigration probability in response to conspecific density. Inside the patch we assumed that population growth was governed by a weak Allee effect with diffusion. Our analysis of the model showed that for the DIE case a patch-level Allee effect occurred for a certain range of patch sizes for any combination of matrix hostility and strength of the weak Allee effect built into the per-capita growth rate. Using this case for comparison, the overall strength of the patch-level Allee effect, measured by the length of the AER, was enhanced by DDE forms with a negative slope ($-$DDE and UDDE). Interestingly, $-$DDE and UDDE can generate a patch-level Allee effect on its own ^[32] and this may explain, in part, why the AER for these two forms of DDE are enhanced. In contrast, when a sufficiently strong positive relationship between density and emigration probability exists (+DDE and hDDE) the patch-level Allee effect can be completely counteracted. For the +DDE case and for several different parameter values, we numerically quantified the strength of positive relationship required to ensure Allee effects were completely counteracted. The hDDE form's negative slope for large population densities prevented a complete counteraction to the patch-level Allee effect for any matrix hostility values. The length of AER was also observed to be dependent on the hostility of the surrounding matrix. For DIE and +DDE forms, the AER length was strictly decreasing as matrix hostility decreased and approached zero when the matrix becomes non-hostile. However, $-$DDE, UDDE, and hDDE showed a more complex relationship with matrix hostility, increasing for certain ranges of matrix hostility and decreasing for others. In this case, AER length strictly decreased as matrix hostility decreased and approached zero when the matrix became non-hostile.

Our analysis of the structure of positive steady states of the model indicated that DDE forms sufficiently similar to DIE only exhibited an Allee effect type of bi-stability. However, if the maximum or minimum emigration probability for a DDE form other than +DDE was sufficiently close to one or zero, respectively, then the model predicted more complicated dynamics with multiple positive steady states for the model in a range of patch sizes where a patch-level Allee effect is not possible, i.e. in a range of patch sizes where the zero population density is unstable. The dynamics of the model were noted to be dependent on the level of matrix hostility. In fact, for the DDE forms, structure of positive steady states of the model was identical to the DIE case when matrix hostility was high, whereas, population dynamics much more complicated than that of the DIE form were exhibited for intermediate to low matrix hostility levels. Our analysis of steady state density profiles of the model under DDE revealed that symmetry of the steady states about the center of the patch is ensured for all parameter values in the case of DIE and +DDE. In this case, density levels at either boundary point are equal. Notwithstanding, we provided sufficient conditions limiting the lower bound of emigration probability of the forms $-$DDE, UDDE and hDDE under which positive steady states of the model were guaranteed to be symmetric about the center of the patch. Our numerical results indicated that violating these sufficient conditions could produce asymmetric positive steady states in which the population density was higher at one boundary point than the other. We also note that this asymmetric steady state phenomenon only occurred when the level of matrix hostility was sufficiently large.

Our results have several consequences for studying populations with a demographic weak Allee effect and DDE. First, seemingly contradictory results were observed when comparing Allee effects at the patch-level versus the metapopulation-level. For example, +DDE will allow for Allee mechanisms to operate when considered at the metapopulation-level ^[22], whereas, +DDE will counteract Allee effects, even alleviate them, if the positive slope is sufficiently strong when considered at the patch-level. This difference rather serves as theoretical evidence to support the caution suggested in ^[51] that Allee effect mechanisms may operate differently depending on the spatial scale. Second, DDE can have profound effects on dynamics of a population by creating, enhancing, and counteracting patch-level Allee effects, causing more complicated dynamics than a patch-level Allee effect via the existence of multiple positive steady states of the model, and creating asymmetric density profiles for steady states of the model. For example, with a moderately hostile matrix, UDDE can cause multiple stable states for a relatively broad range of patch sizes. Bi-stability and multiple stable states, brought about by demographic Allee effects, may contribute to the high variability observed in patch occupancy and density in real metapopulations (e.g., ^[52]). Although a few theoretical studies have explored Allee effects and multiple stable states ^[32,33,34], this subject has scarcely been considered by ecologists but see ^[2]. Our results may also provide some explanation on why it is difficult to empirically verify an Allee effect type of bi-stability at the patch-level. Determining that an organism exhibits a demographic weak Allee effect may not be enough to guarantee existence of a patch-level Allee effect since the presence of density dependent emigration that has a positive slope may alleviate or even completely counteract the Allee effect. Even if the patch-level Allee effect is not completely counteracted by the DDE form, the Allee effect range may become so small that practically finding patch sizes where an Allee effect type of bi-stability occurs is near impossible.

Finally, we note that all of these phenomena are intricately dependent upon the hostility of the matrix surrounding the patch. At low to moderate matrix hostilities, the form of DDE appears to be critically important in determining population dynamics but at high hostility, all forms of DDE have similar dynamic consequences. In human-altered landscapes, where the matrix may be completely inhospitable (e.g., urban areas, clearcuts), the relationship between dispersal and density may matter far less than the magnitude of the dispersal rate in determining metapopulation persistence. Future theoretical and empirical studies on habitat fragmentation should consider, not only the size of remnant patches, but also the deterioration of the patch matrix.

Acknowledgments

This material is based upon work supported by the National Science Foundation under Grant No. DMS-1516519, DMS-1516833, & DMS-1516560.

Conflict of interest

All authors declare no conflicts of interest in this paper.