The influence of density in population dynamics with strong and weak Allee effect

Abstract

In this paper, we consider a reaction–diffusion model in population dynamics and study the impact of different types of Allee effects with logistic growth in the heterogeneous closed region. For strong Allee effects, usually, species unconditionally die out and an extinction-survival situation occurs when the effect is weak according to the resource and sparse functions. In particular, we study the impact of the multiplicative Allee effect in classical diffusion when the sparsity is either positive or negative. Negative sparsity implies a weak Allee effect, and the population survives in some domain and diverges otherwise. Positive sparsity gives a strong Allee effect, and the population extinct without any condition. The influence of Allee effects on the existence and persistence of positive steady states as well as global bifurcation diagrams is presented. The method of sub-super solutions is used for analyzing equations. The stability conditions and the region of positive solutions (multiple solutions may exist) are presented. When the diffusion is absent, we consider the model with and without harvesting, which are initial value problems (IVPs) and study the local stability analysis and present bifurcation analysis. We present a number of numerical examples to verify analytical results.

Introduction

In population biology, a small or sparse population for some species is difficult to study because many factors such as mating may be difficult to analyze. A small or sparse population is one most important reasons for the extinction of the population and is known as the Allee effect. Allee effects are broadly defined as a decline in individual fitness at low population size or low density, which can result in critical population thresholds below which populations crash to extinction. There are many reasons by which an Allee effect can arise including

• an Allee effect begins with reproductive mechanisms, including fertilization efficiency in broadcast spawners. For example, eggs of aquatic animals (such as fishes or oysters) that lay many small eggs;

• pollen limitation and mate finding;

• sperm limitation, reproductive facilitation by conspecifics, and female choice;

• mechanisms related to survival: environmental conditioning and particularly predation such as flocking, coloniality, and group vigilance;

• finally, Allee effects are in social and cooperative species, where group size is important for both reproduction and survival.

The American ecologist, Warder Clyde Allee [1], was the first researcher who presented experimental studies on Alle effects in 1931. Since then many researchers studied Allee effects in different frameworks. For two reasons, many of the early results on Allee effects were theoretical. First, demonstrating an Allee effect in the wild usually requires a long and complete data set of population dynamics; Second, population dynamicists were not interested in studying under-crowding, which was regarded as an interesting but anecdotal process. The studies during the early nineties (before the real start of the upsurge), most of the researchers were concerned about marine invertebrates, plants and insects [2,3,4]. Only very few studied about vertebrates, except for results on fisheries [5]. During this time, a few theoretical studies were also published [6, 7].

The Allee effect is an interesting and challenging topic in population dynamics and many researchers studied Allee effects in many different techniques, both theoretically and experimentally. After the dominance of experimental and theoretical studies on Allee effects, researchers now focusing their attention on two new directions: (i) theoretical studies via mathematical models of ever increasing details; and (ii) empirical studies with the unveiling of Allee effects in natural populations. Much progress has been made by researches under these new ideas. Researchers now have a better understanding of Allee effects and they can draw better distinctions. Note that individual fitness is assumed higher at low densities because of lower intra-specific competition.

There are several kinds of Allee effects such as component Allee effect, demographic Allee effect, dormant Allee effect, etc. and there are different types of relationships between these Allee effects. A dormant Allee effect is a component Allee effect that either does not result in a demographic Allee effect or results in a weak Allee effect. If interacting with a strong Allee effect, causes the overall Allee threshold to be higher than the Allee threshold due to the strong Allee effect alone. Also, a strong Allee effect is a demographic Allee effect with an Allee threshold whereas a weak Allee effect is a demographic Allee effect without an Allee threshold. For details see Fig. 1.

Fig. 1

Classical negative density dependence (solid) compared to strong (dashed) and weak (dotted) Allee effects [8]

A demographic Allee effect is weak if at low density the per capita population growth rate is lower than at higher densities, but remains positive (dotted curve on Fig. 1). On the contrary, it can become so low as to become negative below a certain value, called the Allee threshold (dashed curve on Fig. 1): it is then a strong Allee effect ( [8], chapter 3). If a population subject to a strong Allee effect drops below that threshold, the population growth rate becomes negative and the population will get smaller at a hasten rate until it reaches zero or extinction. Very recently a group of ecologists experimented demographic and component Allee effect in group-level species and they found that at all levels, component Allee effect may generate a demographic Allee effect. A group-level demographic Allee effect always implies at least one group or subgroup-level component Allee effect [9].

Experimentally, it is easier to demonstrate that an Allee effect is strong than that an Allee effect is weak. With adequate time, we can easily trace the density for which the per capita population growth rate becomes negative. Estimates of the Allee threshold are usually approximate because low abundance will yield high observation error, fluctuations due to demographic stochasticity, and a significant proportion of counts at zero [10]. Although having an impressive dataset, the study on the island fox could not identify the value of the Allee threshold even though the Allee effect was strong in some populations because these populations were already declining [11]. A strong Allee effect has also been shown in the gypsy moth (Lymantria dispar), an invasive pest spreading across eastern North America. Also see the recent studies for competition model with the Allee effect [12,13,14,15].

Despite the history of the Allee effect in mathematical biology, the Allee effect is broadly discussed in [16,17,18,19,20,21]. Variation of the Allee effect as well as growing function exhibits ecological phenomena more precisely, especially when a plant or, animal, or organism is about to extinct. In 2016, Rocha et al. discussed the bifurcation and extinction cases of species with a weak Allee effect, where the growth follows Richards’ growth law (details in [22]). Also see [23, 24] for fishery and genetic model with Allee effect and history within. Populations differ in the diffusion strategies as they employ as well as in their environmental intensities [25,26,27]. Further, the competition between populations also has influences on diffusion [28,29,30].

In this paper, our main goal is to evoke the importance of Allee effect in reaction-diffusion equation and show how it effects populations in ecology. We study the impact of multiplicative Allee effect in classical diffusion when the sparsity is either positive or negative. We show that negative sparsity implies weak Allee effect and the population survives in some domain and diverges otherwise. Also, positive sparsity gives strong Allee effect and the population extinct without any condition. The stability conditions and the region of positive solution are presented. We present a number of numerical examples.

This paper is organized as follows. In “Mathematical model” section, we present mathematical models and the fundamental assumption on per capita growth function, f(t, x, u). Allee effects without diffusion are presented in “Allee effects without diffusion” section, where in “Logistic model with Allee effects and harvesting” section we have discussed the harvesting model without diffusion. Bifurcation point and steady state for the diffusion coefficient, $d \rightarrow 0$ with and without harvesting are presented in this portion. In “Allee effects with diffusion” section we present main results of this paper. We present the existence of positive solution and the region where the solution exists for model (11). In “Examples and applications”section, we present some numerical examples where we present bifurcation points and regions of positive solution for particular functions. Finally, in “conclusion” section we present the conclusion and some open problems for further research.

Mathematical model

Describing the Allee effect in a mathematical model is more challenging. We will use the form (1) for the logistic equation in cases where the coefficients are constant, but we want to consider the situations when the intrinsic growth rate r changes sign (concerning time, or in spatial models for space).

$$\begin{aligned} \frac{\displaystyle du}{\displaystyle dt} =ru \left( 1-\frac{\displaystyle u}{\displaystyle K} \right) . \end{aligned}$$

(1)

The assumptions that increase population density lead to decreases in the birth rate and increase the death rate may not always be valid.

Allee (1931) observed that many animals engage in social behavior such as cooperative hunting or group defense that can cause their birth rate to increase or their death rate to decrease with population density [1]. Also, the rate of predation may decrease with prey density in some cases, as discussed by Ludwig in [31] (see recent studies on predation under the Allee effect in [42] and references therein). In the presence of such effects, which are typically known as Allee effects, the model (1) will take a more general pattern

$$\begin{aligned} \frac{\displaystyle du}{\displaystyle d t}=uf(t,x,u), \end{aligned}$$

(2)

where f(t, x, u) may be increasing for some values of u and decreasing for others. A simple case of a model with an Allee effect is

$$\begin{aligned} \frac{\displaystyle du}{\displaystyle d t}=ru(K-u)(u-M), \end{aligned}$$

(3)

where, $r > 0$ and $0< M < K$. The model (3) implies that the density of u will decrease if $0< u < M$ or $u>K$ but increase if $M< u < K$. Introducing diffusion in this model is more challenging and effective in spatial heterogeneity.

To demonstrate the model with diffusion and the Allee effect, it is essential to introduce the following short description. Mathematically the growth rate per capita function f(t, x, u) will not have maximum value at $u=0$. If f(t, x, u) is negative when u is small, we define such a growth pattern has a strong Allee effect or mandatory extinction. If f(t, x, u) is decreasing but still positive for the very low density of u, then this growth has a weak Allee effect or extinction-survival situation. Considering diffusion in ecology model has open a new era.

Let us now consider the general reaction-diffusion mathematical model with homogeneous Neumann boundary conditions

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{\displaystyle \partial u}{\displaystyle \partial t}=d \Delta u+u(t,x)f(t,x,u),&{}t>0,~ x\in \Omega ,\\ \nabla u\cdot n =0,&{}t>0, ~x\in {\partial \Omega },\\ u(0,x)=u_0(x),&{}x\in {\Omega }. \end{array}\right. } \end{aligned}$$

(4)

Here d is the diffusion rate. The habitat $\Omega$ is a bounded region in ${\mathbb {R}}^n$ with ${\partial }{\Omega }\in C^{2+\beta }$, $\beta >0$ and n is the outward normal vector. The functions u(t, x) represent the population density of the species and it’s migration rate is positive, $d > 0$. We assume that all the functions are in the class of $C^{1+\beta } ({\overline{\Omega }})$, $\beta >0$ for any $x\in {\overline{\Omega }}$.

The similar type model (4) was studied by Shi and Shivaji in [16] for weak Allee effects with global bifurcation analysis. They considered the following functions as a reaction term:

$$\begin{aligned} f(t,x,u)=r(x)\left( 1-\frac{u(t,x)}{K(x)}\right) - \frac{A}{1+Bu(t,x)},\;\;\text {or}\;\;f(t,x,u)=m(x)u(t,x)-b(x) u^2(t,x), \end{aligned}$$

where the Allee effect is due to a type-II functional response.

It is important to note that the function f(t, x, u) in (4) will be negative or positive such that the Allee effect will be seen depending on the choice of species mating difficulty or, low growth rate or high death rate as well as competition, etc. There are several kinds of growth function with the Allee Effect has been discovered since 1931. There are numerous example of f(t, x, u) with Allee effects, for examples:

In 1988, Leah Edelstein-Keshet use a quadratic polynomial of the form [32]

$$\begin{aligned} f(t,x,u)=r(x)-b(u(t,x)-a)^2. \end{aligned}$$

(5)

The negative quadratic Eq. (5) describes the per capita growth rate under the Allee effect. Here the function r(x) is the intrinsic growth rate of species and bounded. Also $a,b>0$, and depending on the choice of a there could be unconditional stability; extinction-stability situation; and unavoidable extinction which we defined as no Allee effect (similar to logistic growth); weak Allee effect; and strong Allee effect. We will now define an important function for further studies.

Definition 1

A sparse function is a function $f: \Omega \rightarrow {\mathbb {R}}$ such that $f(\Omega )<<|\Omega |$, where $|\Omega |$ is the magnitude of a set (i.e., the number of elements in the set). The function f is also in the class of $C^{1+\alpha }(\Omega ),\;0<\alpha <1$.

Biologically sparse means thinly scattered (species are set or planted here and there or, not being dense or close together), and any species of the population having this property is known as a sparse population. This idea leads to a sparsity function, which gives biological meaning to a mathematical expression.

The next two growth function with Allee effects are a variation of (5) and are obtained by modifying the Verhulst logistic model

$$\begin{aligned} f(t,x,u)= & {} r(x)\left( 1-\frac{u(t,x)}{K(x)}\right) \left( \frac{u(t,x)-M(x)}{K(x)}\right) . \end{aligned}$$

(6)

$$\begin{aligned} f(t,x,u)= & {} r(x)\left( 1-\frac{u(t,x)}{K(x)}\right) \left( \frac{u(t,x)}{M(x)}-1 \right) . \end{aligned}$$

(7)

The function (6) was broadly discussed in [33, 34]. Here M(x) is a sparse function, defined as $M(x)<K(x)$, with K(x) being the carrying capacity. And by choice of different M(x), different type of Allee effects occurs. Similarly (7) was discussed in [35].

Another choice of f(t, x, u) shown in (8), which is more appealing choice for Allee effect than the previous Eqs. (6) and (7) and is broadly discussed in [36, 37].

$$\begin{aligned} f(t,x,u)=r(x)\left( 1-\frac{u(t,x)}{K(x)}\right) \left( 1- \frac{M(x)+C}{u(t,x)+C}\right) . \end{aligned}$$

(8)

Here M(x) has the same interpretation as Eq. (6) and (7). Also by choice of constant C, (8) is a captivating choice, since the dynamics approaches logistic growth for large u. Model (8) has already been used with $C\equiv 0$ in competition, predation and meta-population models ([36, 38] and references therein). Therefore, considering model (8) with $C>0$ yields a more flexible model with variable extinction rates.

The problem (8) with $C=\sup \limits _{\Omega } \left[ \frac{\displaystyle r(x)}{\displaystyle K(x)}\right]$ and re-parameterization gives

$$\begin{aligned} f(t,x,u)&=r(x)\left( 1-\frac{u(t,x)}{K(x)}\right) \left( 1- \frac{M(x)+C}{u(t,x)+C}\right) \\&=r(x)\frac{K(x)-u(t,x)}{K(x)} \cdot \frac{u(t,x)+C-M(x)-C}{u(t,x)+C}\\&=\frac{r(x)}{K(x)}\cdot \frac{(K(x)-u(t,x))(u(t,x)-M(x))}{u(t,x)+C}\\&=C \cdot \frac{u(t,x)(K(x)+M(x)-u(t,x))-K(x)M(x)}{u(t,x)+C}\\&=\frac{-K(x)M(x)}{1+\frac{1}{C}u(t,x)}+u(t,x)\frac{(K(x)+M(x)-u(t,x))}{1+\frac{1}{C}u(t,x)}. \end{aligned}$$

Now choosing $b=-\frac{K(x)M(x)}{1+\frac{1}{C}u(t,x)}$, $c=\frac{1}{C}$ and $a=K(x)+M(x)$, we get the following equation

$$\begin{aligned} f(t,x,u)=b+u(t,x)\frac{a-u(t,x)}{1+cu(t,x)}. \end{aligned}$$

(9)

In 1996, Takeuchi proposed and discussed this form (9) of Allee effect with $a, b, C > 0$.

The per capita growth function (10) proposed by Jacobs in 1984 combines features of the heuristic approach with biological reasoning; (see [40])

$$\begin{aligned} f(t,x,u)=r_0+ \alpha \beta \frac{u^{\omega }}{u^{\omega }+ \beta }-cu^{\gamma }, \end{aligned}$$

(10)

where $\alpha , \beta , \omega , \gamma >0$. The function defined in (10) treats the Allee effect and negative density dependence separately and may combine them in a wide range of biologically plausible relationships. On the other hand, the model is quite complicated, and consequently, much work has not been studied.

Since each functions (5)–(10) are biologically meaningfull and generalized, see [41] for more details, here we mainly consider the function (7) with our diffusion model (4). Using Eq. (7) is more appropriate to this paper since the main concept of this paper is to show both strong and weak Allee effect in a heterogeneous environment with and without diffusion. The sparsity function here gives the demographic Allee effect and made extinction of species within a short period in the presence of diffusion. The following sections will cover the analytic discussions and numerical approximations for further possibilities of this kind of Allee effect. By considering suitable sparse function and parameters the function of the other could give similar results and also broadly discussed and established in the references provided. Numerically we also consider the cases where functions (6) and (8) appears.

Let us now consider the following mathematical model of (4) with f(t, x, u) as defined in (7). Then we get the following regular (classical) diffusion model with Allee effect:

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{\displaystyle \partial u}{\displaystyle \partial t}=d \Delta u(t,x) +r(x)u(t,x)\left( 1-\frac{\displaystyle u(t,x)}{\displaystyle K(x)}\right) \left( \frac{\displaystyle u(t,x)}{\displaystyle M(x)}-1\right) ,&{}t>0,~ x\in \Omega ,\\ \nabla u\cdot n =0,&{}t>0 ,~x\in {\partial }{\Omega },\\ u(0,x)=u_0(x),&{}x\in {\Omega }. \end{array}\right. } \end{aligned}$$

(11)

For spatial positive functions M(x) and K(x), it is assumed that $0<M(x)<K(x)$ for any $x\in {\Omega }$. Here, M(x) is the sparsity function and $r(x)>0$ is the intrinsic growth rate. Because of the choice of sparsity function, a species could have an Allee effect or not. Therefore, the maximum population density is K(x) and is known as the environment carrying capacity. Throughout the paper, we have assume that the growth rate per capita (7) [as well as (5)–(10)] satisfies:

• (h1) For any $u \ge 0$, $f(t,\cdot ,u) \in { C^{1 +\beta } ({\overline{\Omega }})}$ for $\beta \in {(0,1)}$, and for any $x \in {{\overline{\Omega }}}$, $f(t,x,\cdot ) \in {C^1({\mathbb {R}}^+)}$;

• (h2) For any $x \in {{\overline{\Omega }}}$, there exists $u_2(x) \ge 0$ such that $f(t,x,u) \le 0$ for $u(t,x) > u_2(x)$, and $u_2(x) \le K(x)$ for all $x \in {{\overline{\Omega }}}$; and

• (h3) For any $x \in {{\overline{\Omega }}}$, there exists $u_1(x) \ge 0$ such that $f(t,x,\cdot )$ is increasing in $[0,u_1(x)]$, $f(t,x,\cdot )$ is decreasing in $[u_1(x), \infty )$, and there exists $N>0$ such that $N \ge f(t,x,u_1)$ for all $x \in {{\overline{\Omega }}}$.

From the assumptions, the function $u_2(x)$ indicates the crowding effects of the population at x, which may vary by location but it has a uniform upper bound, K(x). The function $u_1(x)$ is where f(t, x, u) attains the maximum value. Here we still allow logistic growth in cases where $u_1(x) = 0$. The constant N is the uniform upper bound of the growth rate per capita.

Fig. 2

The graphs on top row are growth rate uf(t, x, u), and the ones on lower row are growth rate per capita f(t, x, u), where a logistic for $M \rightarrow 0$; b weak Allee effect for $M=1.1+ \cos (\pi x)$; c strong Allee effect for $M=-(1.1+ \cos (\pi x))$; with $K(x)=2.0+ \cos (\pi x)$, $d=1.0$, $r(x)=1.2$ over $\Omega =[0,1]$. In all figure f(t, x, u) follows (7)

The function f(t, x, u) could be logistic or having Allee effect based on the following assumptions (see also Fig. 2):

• Logistic type when $f(t,x,0)>0$, $u_1(x)=0$ and $f(t,x,\cdot )$ is decreasing in $[0,u_2(x)]$ (Fig. 2a).

• Weak Allee effect type when there exists $x_s>0$ such that $f(t,x,0)\ge 0$ for $x>x_s$ and $f(t,x,0)\le 0$ for $x<x_s$ in a non-empty open domain $\Omega _s\subset \Omega$. Also $u_1(x)>0$, $f(t,x,\cdot )$ is increasing in $[0,u_1(x)]$, and $f(t,x,\cdot )$ is decreasing in $[u_1(x),u_2(x)]$ (Fig. 2b).

• Strong Allee effect type when $f(t,x,0)< 0$, $u_1(x)>0$, $f(t,x,u_1)>0$, $f(t,x,\cdot )$ is increasing in $[0, u_1(x)]$, $f(t,x,\cdot )$ is decreasing in $[u_1(x),u_2(x)]$ (Fig. 2c).

Throughout the paper, we will consider f(t, x, u) as defined in (7), which satisfies (h1)–(h3). To estimate the long-time dynamics of the species whether they will persist or become extinct; steady state solution and eigenvalue-eigenfunction method are used to determine the dynamical behavior.

Allee effects without diffusion

Logistic model with Allee effects

Consider the case when diffusion is ignorable; that is $d \rightarrow 0$ (for historical interest, see [42]), then the species birth (or death) is only dependent on growth rate (or death rate). If we allow Allee effects in this type of model, we get the following reduction model of (11) with $d=0$:

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{\displaystyle d u}{\displaystyle d t}= ru(t)\left( 1-\frac{\displaystyle u(t)}{\displaystyle K}\right) \left( \frac{\displaystyle u(t)}{\displaystyle M}-1 \right) ,&{} t>0,\\ u(0)=u_0,&{} \end{array}\right. } \end{aligned}$$

(12)

where $M < K$. The parametric values, r, is defined as the intrinsic growth rate, K is the carrying capacity, and M is the Allee threshold.

Let $g(u)=r u(t)\left( 1-\displaystyle \frac{\displaystyle u(t)}{\displaystyle K}\right) \left( \displaystyle \frac{ u(t)}{\displaystyle M}-1 \right)$. To find the stationary solutions, $u^*$ of (12), directly we set $\frac{\displaystyle du^*}{\displaystyle dt}=0$ such that

$$\begin{aligned}&\therefore ru^* \left( 1-\frac{\displaystyle u^*}{\displaystyle K}\right) \left( \frac{\displaystyle u^*}{\displaystyle M}-1 \right) =0\\&\Rightarrow u^*= \left\{ 0,\; M,\; K \right\} . \end{aligned}$$

Also we get

$$\begin{aligned} \frac{d g}{d u}=r\left\{ \left( 1-\frac{\displaystyle u}{\displaystyle K}\right) \left( \frac{\displaystyle u}{\displaystyle M}-1 \right) + \frac{u}{M} \left( 1-\frac{\displaystyle u}{\displaystyle K}\right) - \frac{u}{K} \left( \frac{\displaystyle u}{\displaystyle M}-1 \right) \right\} . \end{aligned}$$

(13)

Since parameters r and K, except M, are strictly positive, and M is the Allee threshold (positive or negative threshold) and $M<K$, then depending on the choice of M, we have the following two cases:

Case 1 When $0<M<K$ (weak Allee effect):

Using $u^*= \left\{ 0,\; M,\; K \right\}$ in (13), we get,

$$\begin{aligned} \frac{d g}{d u}(0)= & {} -r<0;\\ \frac{d g}{d u}(M)= & {} r\left( 1-\frac{M}{K} \right) >0;\;\text {and}\\ \frac{d g}{d u}(K)= & {} -r\left( \frac{K}{M}-1 \right) <0. \end{aligned}$$

Therefore, $u^*=\left\{ 0,\; K \right\}$ are stable solution of (12) and $u^*=M$ is the unstable state of (12) when the species are facing weak Allee effects.

Case 2 When $M<0<K$ (strong Allee effect):

Putting $u^*= \left\{ 0,\; M,\; K \right\}$ in (13), where $M=-|M|$, we get that the equilibrium 0 does not change its stability. But when, $M=-|M|$,

$$\begin{aligned} \frac{d g}{d u}(M)= & {} \frac{d g}{d u}(-|M|)=r\left( 1+\frac{|M|}{K} \right)>0;\;\text {and}\\ \frac{d g}{d u}(K)= & {} r\left( 1+\frac{K}{|M|} \right) >0. \end{aligned}$$

Therefore, $u^*=0$ is the only stable solution of (12) and $u^*=M, \; K$ are unstable state of (12), when the species are facing strong Allee effects.

For a single species population in ecology, the biological significance of these stable-unstable strategies is that either the species repelled or converge to the carrying capacity. The present study employs that the stable situation is either $u=0$ or $u(t)=K$, i.e., the population will die out or survive, respectively, and the maximum population density is K. A bifurcation point or threshold level can be found in an unstable solution, $u(t)=M$.

For the sake of diversity and to cover the wide range of studies, in the following section, we introduce the harvesting term in absence of diffusion.

Logistic model with Allee effects and harvesting

If harvesting is involved then the model can be formulated in the following way from (12):

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{\displaystyle d u}{\displaystyle d t}= ru(t)\left( 1-\frac{\displaystyle u(t)}{\displaystyle K}\right) \left( \frac{\displaystyle u(t)}{\displaystyle M}-1 \right) -Hu(t),&{} t>0,\\ u(0)=u_0&{}. \end{array}\right. } \end{aligned}$$

(14)

Here the harvesting function H is constant, $H \ge 0$ and $H < r$, and in case of Allee effects, $M < K$. Let us define the following function

$$\begin{aligned} g^*(u)= u(t)\left[ r \left( 1-\frac{ u(t)}{ K}\right) \left( \frac{ u(t)}{ M}-1 \right) -H\right] , \end{aligned}$$

and

$$\begin{aligned} \frac{d g^*}{d u}=r\left\{ \left( 1-\frac{\displaystyle u}{\displaystyle K}\right) \left( \frac{\displaystyle u}{\displaystyle M}-1 \right) + \frac{u}{M} \left( 1-\frac{\displaystyle u}{\displaystyle K}\right) - \frac{u}{K} \left( \frac{\displaystyle u}{\displaystyle M}-1 \right) \right\} -H. \end{aligned}$$

(15)

Now letting $\frac{\displaystyle d u}{\displaystyle d t}=0$, we get

$$\begin{aligned}&u\left[ r \left( 1-\frac{ u}{ K}\right) \left( \frac{ u}{ M}-1 \right) -H\right] =0\\&\quad \Rightarrow ru^2-r(K+M)u+KM(r+H)=0; \; \text {and} \; u=0\\&\quad \Rightarrow u=\frac{\displaystyle r(K+M) \pm \sqrt{r^2(K-M)^2-4rKMH}}{\displaystyle 2r}; \; \text {and} \; u=0. \end{aligned}$$

Therefore, the three equilibrium states of (14) are

$$\begin{aligned} u_1&=0,\\ u_2&=\frac{\displaystyle 1}{\displaystyle 2r}\left[ r(K+M) + \sqrt{r^2(K-M)^2-4rKMH}\right] ,\;\text {and} \\ u_3&=\frac{\displaystyle 1}{\displaystyle 2r}\left[ r(K+M) - \sqrt{r^2(K-M)^2-4rKMH}\right] . \end{aligned}$$

The equilibrium states $u_2$ and $u_3$ exist as long as

$$\begin{aligned} r^2(K-M)^2-4rKMH>0. \end{aligned}$$

Let us define $D:= \sqrt{r^2(K-M)^2-4rKMH}$ which modify $u_2$ and $u_3$ in the following mathematical structures

$$\begin{aligned} u_2= & {} \frac{\displaystyle 1}{\displaystyle 2}\left[ (K+M)+\frac{\displaystyle D}{\displaystyle r}\right] ,\\ u_3= & {} \frac{\displaystyle 1}{\displaystyle 2}\left[ (K+M)-\frac{\displaystyle D}{\displaystyle r}\right] . \end{aligned}$$

Now from Eq. (15), we get

$$\begin{aligned} \frac{\displaystyle d g^*}{\displaystyle d u}(u_1=0)=-(r+H) <0, \;\;\text {since}\;r>0,\;H\ge 0. \end{aligned}$$

Thus we conclude that $u_1=0$ is a stable solution either there is a strong Allee effect or weak Allee effect. The rest two equilibria states are $u_2$ and $u_3$, and in terms of H and M, the following two situations arise:

When $0<M<K$, i.e., weak Allee effect, then we get the following three cases in terms of H:

$$\begin{aligned} H<\frac{\displaystyle r(K-M)^2}{\displaystyle 4KM},\;\; H=\frac{\displaystyle r(K-M)^2}{\displaystyle 4KM},\;\text {and}\;\; H>\frac{\displaystyle r(K-M)^2}{\displaystyle 4KM}. \end{aligned}$$

Case I When $H>\displaystyle \frac{\displaystyle r(K-M)^2}{\displaystyle 4KM}$, that is if $D^2<0$ then $u_2, u_3 \in {\mathbb {C}}\setminus {\mathbb {R}}$ hence their stability is irrelevant and we have no conclusion.

Case II If $H=\displaystyle \frac{\displaystyle r(K-M)^2}{\displaystyle 4KM}$ then $D=0$ and the system (14) has one identical solution, $u_2=u_3=\frac{\displaystyle K+M}{\displaystyle 2}$ and from (15) it is seen that

$$\begin{aligned} \frac{\displaystyle d g^*}{\displaystyle d u}\left( \frac{\displaystyle K+M}{\displaystyle 2}\right) =0. \end{aligned}$$

Therefore, in this situation, the system has a unique positive equilibrium, which is semi-unstable.

Case III If $H<\displaystyle \frac{\displaystyle r(K-M)^2}{\displaystyle 4KM}$ then $D>0$ and the problem (14) has two distinct solutions $u_2 \; \text {and} \; u_3$, where

$$\begin{aligned} u_2=\frac{\displaystyle 1}{\displaystyle 2}\left[ (K+M)+\frac{\displaystyle D}{\displaystyle r}\right] \; \text {and} \; u_3=\frac{\displaystyle 1}{\displaystyle 2}\left[ (K+M)-\frac{\displaystyle D}{\displaystyle r}\right] . \end{aligned}$$

Since $D>0$, $u_2$ is always positive. We recall (15) and get the simplified form as follows:

$$\begin{aligned} \frac{d g^*}{d u}&=r\left[ -\frac{\displaystyle 3}{\displaystyle KM} u^2+ \frac{\displaystyle 2(K+M)}{\displaystyle KM}u -1 \right] -H, \end{aligned}$$

such that

$$\begin{aligned} \frac{d g^*}{d u}(u_2)&=r\left[ -\frac{\displaystyle 3}{\displaystyle KM} u_2^2+ \frac{\displaystyle 2(K+M)}{\displaystyle KM}u_2 -1 \right] -H. \end{aligned}$$

Using $u_2=\frac{\displaystyle 1}{\displaystyle 2}\left[ (K+M)+\frac{\displaystyle D}{\displaystyle r}\right]$, we get

$$\begin{aligned} \frac{d g^*}{d u}(u_2)=r\left[ \frac{\displaystyle (K+M)^2}{\displaystyle 4KM} -1 \right] -H -\left[ \frac{K+M}{2KM}D+\frac{3}{4rKM}D^2 \right] . \end{aligned}$$

(16)

It is noted that $\frac{\displaystyle d g^*}{\displaystyle d u}\left( \frac{\displaystyle K+M}{\displaystyle 2} \right) =r\left[ \frac{\displaystyle (K+M)^2}{\displaystyle 4KM} -1 \right] -H=0$ when $H=\frac{\displaystyle r(K-M)^2}{\displaystyle 4KM}$. Hence there exits at least one set of positive parameters, $r,\;K,\;M,\;\text {and}\;H$ with $M<K$ such that we assume $H=\frac{\displaystyle r(K-M)^2}{\displaystyle 4KM}-\epsilon ,\;\epsilon >0$ for $H<\frac{\displaystyle r(K-M)^2}{\displaystyle 4KM}$, where $\epsilon$ is small enough. Using the value of H, it is found that $D=2\sqrt{r\epsilon KM}>0$. Substituting H and D in (16) yields

$$\begin{aligned} \frac{d g^*}{d u}(u_2)&=-r+\epsilon -\left[ \frac{K+M}{2KM}2\sqrt{r\epsilon KM}+\frac{3}{4rKM} 4r\epsilon KM \right] \\&=-r+\epsilon - \frac{K+M}{KM}\sqrt{r\epsilon KM}- 3\epsilon \\&=-r-2 \epsilon - \frac{(K+M)\sqrt{r\epsilon KM}}{KM}<0. \end{aligned}$$

In case of equilibrium point $u_3$, it is seen that $u_3$ is positive only for $r(K+M)>D$. Then using the fact that, $\frac{\displaystyle dg^*}{\displaystyle du}\left( \frac{\displaystyle K+M}{\displaystyle 2} \right) =0$ and $H<\frac{\displaystyle r(K-M)^2}{\displaystyle 4KM}$ and from the above constructive manner, we found that $\frac{d g^*}{d u}(u_3)>0$.

All the results from case III are summarized as follows:

• 1. For $D>0$, it is obvious that $u_2 >0$. Using the fact that $\displaystyle \frac{\displaystyle d g^*}{\displaystyle d u}\left( \frac{\displaystyle K+M}{\displaystyle 2}\right) =0$ and $H<\frac{\displaystyle r(K-M)^2}{\displaystyle 4KM}$, we have $\displaystyle \frac{\displaystyle d g^*}{\displaystyle d u}(u_2) <0.$

• 2. When $r(K+M)>D$ then $u_3 >0$, and we have $\displaystyle \frac{\displaystyle d g^*}{\displaystyle d u}(u_3) >0.$ If $r(K+M)<D$, the solution $u_3$ is strictly negative and the solution is not acceptable.

Fig. 3

Numerical solutions of (14) for $r = 0.8147,\; K = 1,\; u_0=0.5$ with a $0<M = 0.30<K,\; H = 0.2661$, and a $M = -0.9058<0<K,\; H = 0.1270$

The stable and convergent solutions are presented in Fig. 3 both for weak and strong Allee effects. For large time, t, the solution $u(t)\rightarrow u_2$ when M is positive (Fig. 3a). Similarly, in case of strong Allee effect, the solution u(t) is converging to $u_1=0$ (Fig. 3b).

In conclusion, for weak Allee effect, we have two stable equilibria $u_1=\left\{ 0, \; u_2\right\}$, and one unstable equilibrium state $u_3$. Hence the solution of the model (14) is decreasing for $0<u(t)<u_3$, and increasing for $u_3<u(t)<u_2$.

In case of strong Allee effects, that is when $M<0<K$, we put $-|M|$ instead of M, then we get the following equilibria

$$\begin{aligned} u_1= & {} 0,\\ u_2= & {} \frac{\displaystyle 1}{\displaystyle 2} \left[ (K-|M|) + \frac{\displaystyle D^*}{\displaystyle r}\right] ,\;\text {and}\\ u_3= & {} \frac{\displaystyle 1}{\displaystyle 2}\left[ (K-|M|)-\frac{\displaystyle D^*}{\displaystyle r}\right] . \end{aligned}$$

Where $D^*:= \sqrt{r^2(K+|M|)^2+4rK|M|H}$. Here we also get three cases in terms of H:

$$\begin{aligned} H<\frac{\displaystyle -r(K+|M|)^2}{\displaystyle 4K|M|},\;\; H=\frac{\displaystyle -r(K+|M|)^2}{\displaystyle 4K|M|},\;\text {and}\;\; H>\frac{\displaystyle -r(K+|M|)^2}{\displaystyle 4K|M|}. \end{aligned}$$

Since $D^*$ is not valid for this $H<\displaystyle \frac{\displaystyle -r(K+|M|)^2}{\displaystyle 4K|M|}$, so there is no equilibrium point, that is $u_2, u_3 \in {\mathbb {C}}\setminus {\mathbb {R}}$.

On contrary, the value $-\frac{\displaystyle r(K+|M|)^2}{\displaystyle 4K|M|}$ is strictly negative since $r, \; K$ and |M| are strictly positive. Hence H is negative for either, $H=\displaystyle \frac{\displaystyle -r(K+|M|)^2}{\displaystyle 4K|M|}$ or, $H>\displaystyle \frac{\displaystyle -r(K+|M|)^2}{\displaystyle 4K|M|}$, which is a contradiction of our earliest assumption that harvesting can not be negative. Considering this cases, we can conclude that there is no condition of $D^*$ found, such that a positive equilibrium exists when the harvest is positive.

Therefore, if we allow harvesting in a strong Allee effect, the resulting behavior is diverging, and does not exist a stable point except possibly, $u_1 = 0$.

Let us now introduce the diffusion term in our considered model for further study.

Allee effects with diffusion

In the dynamical behavior of (11) without an Allee effect, the results are well established. When the species has logistic growth, there is a critical value $d_1>0$ such that $0<d<d_1$, i.e. for slow diffusion, there is a unique positive steady state $u^*(x)$, which is the asymptotic limit for any non-negative initial distribution except $u_0 \equiv 0$, thus the persistence of the population is achieved. If $d>d_1$, the diffusion is fast and the only non-negative steady state solution is $u^*=0$, thus the extinction is inevitable; the readers can check the literature ( [16], Section 1) for more details about the inevitable extinction. If u(x) is any stationary solution of (11), consider the following eigenvalue problem to analyze the linear stability of the system,

$$\begin{aligned} \Delta \phi + \sigma \; r(x) \phi [f+uf_u]=0,~ x\in \Omega ,\;\nabla \phi =0,~x\in {\partial }{\Omega }, \end{aligned}$$

where $f_u$ is the derivative of f with respect to u, and

$$\begin{aligned} f(t,x,u)=r(x) \left( 1-\frac{\displaystyle u^*(x)}{\displaystyle K(x)}\right) \left( \frac{\displaystyle u^*(x)}{\displaystyle M(x)}-1\right) . \end{aligned}$$

The function $\phi (x)$ is the eigenfunction and $\sigma$ is the corresponding eigenvalue. Here the constant $d_1=1/\sigma _1(f,\Omega )$, and $\sigma _1(f,\Omega )$ is the principle eigenvalue of the linearized system of

$$\begin{aligned} \Delta \phi + \sigma \; r(x) \phi \left( 1-\frac{\displaystyle u^*(x)}{\displaystyle K(x)}\right) \left( \frac{\displaystyle u^*(x)}{\displaystyle M(x)}-1\right) =0,~ x\in \Omega ,\;\nabla \phi =0,~x\in {\partial }{\Omega }. \end{aligned}$$

(17)

Bifurcation analysis

Let $\sigma = d^{-1}$ and $u^*(x)$ be the stationary solution of (11) then

$$\begin{aligned} {\left\{ \begin{array}{ll} \Delta u^*(x) + \sigma \; r(x)u^*(x)\displaystyle \left( 1-\frac{\displaystyle u^*(x)}{\displaystyle K(x)}\right) \left( \frac{\displaystyle u^*(x)}{\displaystyle M(x)}-1\right) =0,&{} x\in \Omega ,\\ \nabla u^*\cdot n =0,&{}x\in {\partial }{\Omega },\\ u(0,x)=u_0(x),&{}x\in {\Omega }. \end{array}\right. } \end{aligned}$$

(18)

When a stable solution of (18) is viewing as an equilibrium solution to (11) can be determined by the eigenvalue problem:

$$\begin{aligned}&\Delta \phi + \sigma \; r(x) \phi \left( 1- \frac{\displaystyle u^*(x)}{\displaystyle K(x)}\right) \left( \frac{\displaystyle u^*(x)}{\displaystyle M(x)}-1\right) =0 \\&\quad \Rightarrow - \int \limits _{\Omega }|\nabla \phi |^2 \; dx+\sigma \int \limits _{\Omega } r(x) \phi ^2 \left( 1-\frac{\displaystyle u^*(x)}{\displaystyle K(x)}\right) \left( \frac{\displaystyle u^*(x)}{\displaystyle M(x)}-1\right) \; dx=0 \\&\quad \Rightarrow \sigma = \frac{\int \limits _{\Omega }|\nabla \phi |^2 \; dx}{\int \limits _{\Omega } r(x) \phi ^2 \left( 1-\frac{\displaystyle u^*(x)}{\displaystyle K(x)}\right) \left( \frac{\displaystyle u^*(x)}{\displaystyle M(x)}-1\right) \; dx}. \end{aligned}$$

Hence the principle eigenvalue is defined by

$$\begin{aligned} \sigma _1 = \sup \limits _{u \in {H_0^1(\Omega )}} \frac{\int \limits _{\Omega }|\nabla u^*|^2 \; dx}{\int \limits _{\Omega } r(x) u^{*2}(x) \left( 1-\frac{\displaystyle u^*(x)}{\displaystyle K(x)}\right) \left( \frac{\displaystyle u^*(x)}{\displaystyle M(x)}-1\right) \; dx}. \end{aligned}$$

Thus the bifurcation point for positive solutions is defined by

$$\begin{aligned} d_1=\frac{1}{\sigma _1(f,\Omega )}=\sup \limits _{u \in {H_0^1(\Omega )}}\left\{ \int \limits _{\Omega } f(t,x,u^*) u^{*2}(x) \;dx : \int \limits _{\Omega }|\nabla u^*|^2 \; dx =1\right\} . \end{aligned}$$

(19)

Here $f(t,x,u^*)$ follows assumptions (h1)–(h3) and $\sigma _1 \equiv \sigma _1(f,\Omega )$ is a bifurcation point, where nontrivial solution of (18) bifurcate from trivial solution. The classical type of bifurcation diagram at $\sigma _1$ is shown in Fig. 4. The more details about this type of bifurcation can be found in ([16], Section 2).

Fig. 4

Classical bifurcation diagram from left to right showing for logistic; strong Allee effect and weak Allee effect, respectively. Here $\sigma =d^{-1}$, where d is the diffusion constant

Existence of a positive stationary solution

Consider the steady state of (11), if $u^*(x)$ being the positive solution then we get

$$\begin{aligned} 0=d \Delta u^*(x) + r(x)u^*(x) \left( 1-\frac{\displaystyle u^*(x)}{\displaystyle K(x)} \right) \left( \frac{\displaystyle u^*(x)}{\displaystyle M(x)}-1\right) ,~ x\in \Omega , \end{aligned}$$

which implies

$$\begin{aligned} {\left\{ \begin{array}{ll} - \Delta u^* = \displaystyle \frac{\displaystyle 1}{\displaystyle d} ru^*\left( 1-\frac{\displaystyle u^*}{\displaystyle K} \right) \left( \frac{\displaystyle u^*}{\displaystyle M} -1 \right) , &{} x\in \Omega ,\\ \nabla u^*\cdot n =0,&{}x\in {\partial }{\Omega },\\ u^*(0,x)=u_0(x),&{}x\in {\Omega }. \end{array}\right. } \end{aligned}$$

(20)

Now the following situations occur for a steady-state solution.

Theorem 1

Let $\sigma _1 \ge r_1 \displaystyle \frac{\displaystyle (K_1-M_1)^2}{\displaystyle 4K_1M_1}$. Then (20) has no positive solution.

Proof

Let $\sigma _1$ be the principle eigenvalue of the operator $-\Delta$ with Neumann boundary conditions and $\phi >0$ in $\Omega$ be a corresponding eigenvector. Assume that u be a positive solution of (20) then by the Green’s identity, we have

$$\begin{aligned}&\int \limits _{\Omega } (u \Delta \phi - \phi \Delta u) \; dx =0 \nonumber \\&\quad \Rightarrow \int \limits _{\Omega } \left[ \frac{\displaystyle 1}{\displaystyle d} r(x)u(t,x) \left( 1-\frac{\displaystyle u(t,x)}{\displaystyle K(x)} \right) \left( \frac{\displaystyle u(t,x)}{\displaystyle M(x)}-1 \right) \phi -\sigma _1 \phi u\right] \; dx=0 \nonumber \\&\quad \Rightarrow \int \limits _{\Omega } \phi u(t,x)\left[ \frac{\displaystyle 1}{\displaystyle d} r(x) \left( 1-\frac{\displaystyle u(t,x)}{\displaystyle K(x)} \right) \left( \frac{\displaystyle u(t,x)}{\displaystyle M(x)}-1 \right) -\sigma _1 \right] \; dx=0. \end{aligned}$$

(21)

Since f(t, x, u) follows assumption (h1)–(h3), let ${\tilde{f}}(u)=r_1\displaystyle \left( 1-\frac{\displaystyle u}{\displaystyle K_1} \right) \left( \frac{\displaystyle u}{\displaystyle M_1}-1 \right)$, where $r_1, \; K_1$ are maximum value of the functions $r(x),\; K(x)$ and $M_1$ is the minimum value of M(x), respectively in the domain $\Omega$. Then

$$\begin{aligned} \frac{\displaystyle \partial {\tilde{f}}(u)}{\displaystyle \partial u}=r_1\left[ \frac{\displaystyle K_1+M_1}{\displaystyle K_1M_1}-\frac{\displaystyle 2}{\displaystyle K_1M_1}u \right] . \end{aligned}$$

For critical value we know $\displaystyle \frac{\displaystyle \partial {\tilde{f}}(u)}{\displaystyle \partial u} \equiv 0$, which gives $u = \frac{\displaystyle K_1+M_1}{\displaystyle 2}$.

Since $\displaystyle \frac{\displaystyle \partial ^2 {\tilde{f}}(u)}{\displaystyle \partial u^2}=-\frac{\displaystyle 2r_1}{\displaystyle K_1M_1} <0$, the function ${\tilde{f}}(u)$ has a maximum value when $u = \displaystyle \frac{\displaystyle K_1+M_1}{ \displaystyle 2}$. Thus

$$\begin{aligned} {\tilde{f}}_0&:= \sup \limits _{\Omega } {\tilde{f}}(u)\\&= r_1\left( 1-\frac{\displaystyle K_1+M_1}{\displaystyle 2K_1}\right) \left( \frac{\displaystyle K_1+M_1}{\displaystyle 2M_1}-1\right) \\&=r_1 \frac{\displaystyle K_1-M_1}{\displaystyle 2}\cdot \frac{\displaystyle K_1-M_1}{\displaystyle 2} \\&= r_1\frac{\displaystyle (K_1-M_1)^2}{\displaystyle 4K_1M_1}, \end{aligned}$$

and $u(t,x)f(t,x,u) \le u(t,x){\tilde{f}}(t,x,u)$ for $u(t,x) \ge 0$. Thus Eq. (21) implies

$$\begin{aligned}&0=\int \limits _{\Omega } (u \Delta \phi - \phi \Delta u) \; dx \le \int \limits _{\Omega } \phi u \; \left[ \frac{\displaystyle 1}{\displaystyle d}{\tilde{f}}(t,x,u) - \sigma _1\right] \; dx \nonumber \\&\therefore \int \limits _{\Omega } \phi u \; \left[ \frac{\displaystyle 1}{\displaystyle d} {\tilde{f}}(t,x,u) - \sigma _1\right] \; dx \ge 0. \end{aligned}$$

(22)

But if $\sigma _1 \ge r_1\displaystyle \frac{\displaystyle (K_1-M_1)^2}{\displaystyle 4K_1M_1}$ then

$$\begin{aligned} \int \limits _{\Omega } \phi u \; \left[ \frac{\displaystyle 1}{\displaystyle d}{\tilde{f}}(t,x,u) - \sigma _1\right] \; dx < \int \limits _{\Omega } \phi u \; \left( \frac{\displaystyle 1}{\displaystyle d}{\tilde{f}}_0 - \sigma _1\right) \; dx \le 0, \end{aligned}$$

which is a contradiction of (22). Hence (20) has no positive solution for $\sigma _1 \ge r_1 \displaystyle \frac{\displaystyle (K_1-M_1)^2}{\displaystyle 4K_1M_1}$. $\square$

Theorem 2

There exists positive constants $M_0, \; r_1$, $\tilde{K_0}$ and $\mu \in {(0,1]}$ such that for $K_0>\tilde{K_0}$, (20) has a positive solution and the positive solution occurs in the region $\displaystyle \left[ \frac{ K_0}{ 2} \mu ^2 , \frac{ K_0}{ 2}\right]$.

Proof

At first we will prove an existence result for (20). Let $\sigma$ be the principal eigenvalue of the operator $-\Delta$ with Neumann boundary conditions and $\phi >0$ is the corresponding eigenfunction such that $||\phi ||_{\infty }=1$. Hence there exists $\delta >0, \; \mu \in {(0,1]}$ and $m>0$ such that

$$\begin{aligned}&|\nabla \phi |^2 -\sigma \phi ^2 \ge m \; \text {for} \; x\in {\Omega _{\delta }}, \end{aligned}$$

(23)

$$\begin{aligned}&\nabla \cdot \phi \ge \mu \; \text {for} \; x\in {\partial \Omega _{\delta }}, \end{aligned}$$

(24)

where $\Omega _{\delta } :=\left\{ x\in {\Omega }\; |\; d(x, \partial \Omega ) < \delta \right\}$. From (20), we have the reaction function as defined in a simplest form such that

$$\begin{aligned} g(u)=r_1u \left( 1-\frac{\displaystyle u}{\displaystyle K_0} \right) \left( \frac{\displaystyle u}{\displaystyle M_0}-1\right) , \end{aligned}$$

where $K_0 \le K_1$, $M_0 \le M_1$ and $K_1, M_1 \; \text {and} \; r_1$ are the maximum value of $K(x), M(x) \; \text {and} \; r(x)$, respectively as defined in Theorem 1.

It is easily observed that the zeros of g are $0, \; M_0$ and $K_0$ where $K_0>M_0$ and hence

$$\begin{aligned} g(u):=-\frac{\displaystyle r_1}{\displaystyle K_0M_0}u(u-K_0)(u-M_0). \end{aligned}$$

Let $u^*$ be the first positive zero of $\displaystyle \frac{\displaystyle \partial g}{\displaystyle \partial u}$ in which

$$\begin{aligned} u^* = \frac{\displaystyle 1}{\displaystyle 3} \left[ (K_0+M_0) - \sqrt{K_0^2+M_0^2-K_0M_0} \right] < \frac{\displaystyle K_0+M_0}{\displaystyle 3}. \end{aligned}$$

But g is convex on $\displaystyle \left( 0, \frac{\displaystyle K_0+M_0}{\displaystyle 3} \right)$ and hence

$$\begin{aligned} \lambda := - \inf \limits _{u \in {[0, K_0]}} g(u) < \frac{\displaystyle r_1}{\displaystyle 3} \left[ (K_0+M_0) - \sqrt{K_0^2+M_0^2-K_0M_0} \right] =r_1u^*. \end{aligned}$$

If we recall the second set of Fig. 2, which describes the behavior of $u \; \text {and} \;u^*$ (Fig. 5). We first claim that

$$\begin{aligned} \frac{\displaystyle \lambda }{\displaystyle K_0}&<\frac{\displaystyle 1}{\displaystyle K_0} \left[ \frac{\displaystyle r_1}{\displaystyle 3}\left\{ (K_0+M_0)-\sqrt{K_0^2+M_0^2-K_0M_0} \right\} \right] \\&=\frac{\displaystyle r_1M_0}{\displaystyle (K_0+M_0)+\sqrt{K_0^2+M_0^2-K_0M_0}}. \end{aligned}$$

Thus $\displaystyle \frac{\displaystyle \lambda }{\displaystyle K_0} \rightarrow 0$ as $K_0$ tends to infinity. Hence there exists $K_0^{(1)}$ such that for every $K_0 >K_0^{(1)}$, we have

$$\begin{aligned} m > \frac{\displaystyle \lambda }{\displaystyle K_0}. \end{aligned}$$

(25)

Next, we also claim that $\displaystyle \frac{\displaystyle K_0}{\displaystyle M_0} \rightarrow \infty \; \text {as} \; K_0 \rightarrow \infty$. Hence there exists $K_0^{(2)}$ such that for every $K_0 >K_0^{(2)}$, we have

$$\begin{aligned} \left[ \frac{\displaystyle K_0}{\displaystyle 2} \mu ^2 , \frac{\displaystyle K_0}{\displaystyle 2}\right] \subset (M_0, K_0), \end{aligned}$$

(26)

and $K_{\mu } := \inf \limits _{u \in {\left[ \frac{K_0}{2} \mu ^2 , \frac{K_0}{2}\right] }} g(u) >0$. Finally

$$\begin{aligned} \frac{\displaystyle K_{\mu }}{\displaystyle K_0}&=\frac{\displaystyle 1}{K_0}\left[ \min \left\{ g\left( \frac{\displaystyle K_0}{\displaystyle 2} \mu ^2\right) ,g\left( \frac{\displaystyle K_0}{\displaystyle 2}\right) \right\} \right] \\&=\frac{\displaystyle 1}{\displaystyle K_0}\left[ \min \left\{ -\frac{\displaystyle r_1 \mu ^2}{\displaystyle 4M_0}K_0(\mu ^2-2) \left( \frac{\displaystyle K_0 \mu ^2}{\displaystyle 2}-M_0 \right) , \frac{\displaystyle r_1}{\displaystyle 4M_0}K_0 \left( \frac{\displaystyle K_0}{\displaystyle 2}-M_0 \right) \right\} \right] \\&= \min \left\{ -\frac{\displaystyle r_1 \mu ^2}{\displaystyle 4M_0}(\mu ^2-2)\left( \frac{\displaystyle K_0 \mu ^2}{\displaystyle 2}-M_0\right) , \frac{\displaystyle r_1}{\displaystyle 4M_0}\left( \frac{\displaystyle K_0}{\displaystyle 2}-M_0\right) \right\} , \end{aligned}$$

which tends to infinity as $K_0 \rightarrow \infty$. Thus there exists $K_0^{(3)}>K_0^{(2)}$ such that for every $K_0 >K_0^{(3)}$, we have

$$\begin{aligned} \sigma < \frac{\displaystyle K_{\mu }}{\displaystyle K_0}. \end{aligned}$$

(27)

Define $\tilde{K_0} := \max \left\{ K_0^{(1)}, K_0^{(3)}\right\}$ and for some $M_0>0, \; r_1>0$ and $K_0>\tilde{K_0}$, (20) has a positive solution. At this stage, we will use the upper and lower solution method to show the existence of a positive integral which lies in the interval $\left[ \displaystyle \frac{\displaystyle K_0}{\displaystyle 2} \mu ^2 , \displaystyle \frac{\displaystyle K_0}{\displaystyle 2}\right]$; $\mu \in {(0,1]}$.

Fig. 5

Growth rate function uf(t, x, u) with weak Allee effect for $M=1.1+ \cos (\pi x)$, $K(x)=2.0+ \cos (\pi x)$, $d=1.0$ and $r(x)=1.2$ over the domain $\Omega =[0,1]$

First, we prove that $\psi = \displaystyle \frac{\displaystyle K_0}{\displaystyle 2} \phi ^2$ is a sub-solution of (20). For $x\in {\Omega _{\delta }}$,

$$\begin{aligned} -\Delta \psi&=K_0 (\sigma \phi ^2 -|\nabla \phi |^2)\\&\le -K_0 m; \;\text {using}\; (23) \\&\le -\lambda ; \; \text {using} \; (25)\\&\le \inf \limits _{u \in {[0, K_0]}} g(u)\\&\le r_1u \left( 1-\frac{\displaystyle u}{\displaystyle K_0} \right) \left( \frac{\displaystyle u}{\displaystyle M_0}-1\right) . \end{aligned}$$

Also for the boundary:

$$\begin{aligned} -\Delta \psi&=K_0 (\sigma \phi ^2 -|\nabla \phi |^2)\\&\le K_0 \sigma \\&\le K_{\mu }= \inf \limits _{u \in {\left[ \frac{K_0}{2} \mu ^2 , \frac{K_0}{2}\right] }} g(u); \; \text {using} \; (27)\\&\le r_1u \left( 1-\frac{\displaystyle u}{\displaystyle K_0} \right) \left( \frac{\displaystyle u}{\displaystyle M_0}-1\right) . \end{aligned}$$

Hence $\psi = \displaystyle \frac{\displaystyle K_0}{\displaystyle 2} \phi ^2$ is a sub solution of (20). We also note that $\phi =K_0$ is a super solution for (20) since K(x) is the carrying capacity of the environment, therefore the solution would not exceed it. However, if the solution exceeds the carrying capacity for sometimes, the solution will not be stable, and within a short period, it will again converge to K(x). Thus there exists a positive solution of (20) in $\left[ \displaystyle \frac{K_0}{2} \mu ^2 , \frac{K_0}{2}\right]$. $\square$

Examples and applications

To learn the importance of application of Allee effects in biological invasions, we can consider two situations: In one case, consider a species, which is in red list, that is the population density is very low, so there is not enough mating and the environment is not favorable, so ultimately in a short time, the species will go extinct. But if we can improve the sparse function by giving their favorable sparse (in this case, the sparse function will be very low), then the species could survive.

In other cases, if a species is invaded by another species (e.g., bacteria or virus), then the sparse function for invader species will play an important role to immune the invaded species. In this situation, we have to higher the parametric value of sparse function so the invader species could not find favorable space to spread in invaded species. Consequently, the invader species will go extinct and the invaded body will be immune.

Considering all types of similar cases and form the numerical assumptions, the following numerical behaviors are simulated.

Example 1

Consider the carrying capacity $K=2.5$ with an initial density of any species $u_0=0.5$ and intrinsic growth rate of the species, $r=1.0$. Then the numerical behavior of solutions of (12) shows that if the species has a strong Allee effect and the sparse function M is negative, then the solution will always be repelled to zero, which is following in Fig. 6a.

Fig. 6

Solutions of (12) with respect to time t, for $K=2.5$, $u_0=0.5$ and $r=1.0$ where in (a) $M=-1.5, -1.0, -0.5,-0.2, -0.1$ and in (b) $M=1.5, 1.0, 0.5,0.2, 0.1$

If the species are following a weak Allee effect and the sparse function M is positive but $M<K$, then we see that for any $M\in {(0, u_0)}$, the solution of (12) converges to its maximum yields, towards carrying capacity. But if $u_0< M< K$, then the solution is repelling and converging to zero. For $M=u_0$, there exists a constant solution. See Fig. 6b which follows our assumption (h1)–(h3) and the stable conditions of (12).

Example 2

When diffusion is ignorable, that is in the model (12) we have two stable solutions, one is repelling and the other one is asymptotically stable. In numerical behavior, we see similar results where either the species is extinct or survives. In numerical behavior at first we vary initial density, and see that the solution, u(t) of (12) for $K=2.5$, $r=1.0$ and $M=1.15$ either tends to zero or, converge to K. So the solutions of (12) are locally stable. Also there is a constant solution if $u_0 \equiv M$ as seen in Fig. 7a.

Fig. 7

Solutions of (12) with respect to time t, where $r=1.0$ and in (a) $K=2.5$, $M= 1.15$ and $u_0=0.5, 0.9, 1.15, 1.5$ and in (b) $K=2.0$, $u_0=0.9$ and $M=0.6, 0.9, 1.2, 1.8$

In Fig. 7b, we exhibit that if sparse function M is having lower significance and the initial population is in a better position then the species will exist and will converge to K. Contrariwise if the sparse function is more significant in the domain and initial densities are very low, then the species will die out in a short time. In Fig. 7b we vary sparse function as $M=\left\{ 0.6,\; 0.9, \; 1.2, \; 1.8 \right\}$ with $u_0=0.9$, $K=2$ and growth rate is same as above. Here we also see one stable solution and one repelled solution.

Example 3

Consider the Allee effect model with harvesting as in (14). Here we see that for a small harvesting rate and sparse function, positive and lower than initial density gives stability at its utmost point. And on the contrary situation, the solutions converge to zero.

Fig. 8

Solutions of (14) with respect to time t for $K=4.0$, $u_0=1.0$ and $r=1.0$ with various harvesting rate as $H=0.1,\; 0.3,\;0.7,\; \text {and} \; 1.0$; where a $M= 0.6$ (weak Allee effect), and b $M=- 0.6$ (strong Allee effect)

Again we consider various harvesting rates; lower or higher than initial density for strong Allee effects, that is the sparse function is being negative. In this case, zero is the only converging solution. The following Fig. 8 is showing the significance of harvesting functions while we choose $K=4.0$, $u_0=1.0$, $r=1.0$ and for $M= 0.6$ (Fig. 8a) and for $M=- 0.6$ (Fig. 8b). This is also a case of local stability and both figure shows the conditions of Allee effects as described for model (14).

Example 4

Assume now that diffusion is not ignorable, consider in (11), $K(x)=2.5+\cos (\pi x)$ and $M(x)=1.2+\cos (\pi x)$ with migration rate $d=0.1$, growth rate $r=1.0$ and initial density $u_0=0.5$. Then the model (11) has positive solution in the domain $\Omega =(0,1)$, which is shown in Fig. 9.

Fig. 9

a Solutions of (11) with respect to space at $t=50$ and b average solutions of (11) with respect to time t for $K(x)=2.5+ \cos (\pi x)$, $M(x)= 1.2+ \cos (\pi x)$ and $r(x)=1.0$ on $\Omega =(0,1)$, with $u_0=0.5$, $d=0.1$

Example 5

While considering diffusion is not zero and not close to zero if we choose some particular function then we see that there exists one $d^*$ for which, when $d\ge d^*$ there exist no positive solution, positive solutions occur only on $0<d<d^*$. Particularly in (11) let $K(x)=2.5+ \cos (\pi x)$, $M(x)= 1.2+ \cos (\pi x)$, $r(x)=1.0$ and $u_0=0.5$. Then numerically we get a bifurcation point $d^*=0.154$.

Fig. 10

Dynamical behavior of solutions of (11) for various diffusion coefficients at $t=50$ in (a) and average solutions of (11) in (b) for $K(x)=2.5+ \cos (\pi x)$, $M(x)= 1.2+ \cos (\pi x)$ and $r(x)=1.0$ on $\Omega =(0,1)$, with initial density $u_0=0.5$

Any value greater than this $d^*$ is always tending to zero, and the species extinct after a short time. The following Fig. 10 is showing the dynamical behavior of $d^*$, that is when diffusion or migration is very slow, the species will persist. This numerical behavior is analytically proven in Theorem 1.

Example 6

For model (11), let the distribution functions are $K(x)=2.5+ \cos (\pi x)$, $M(x)= 1.2+ \cos (\pi x)$, growth rate is $r(x)=1.0$ on the domain $\Omega =(0,1)$ and initially crowding density is $u_0=0.5$.

Fig. 11

Region of positive solutions of (11) at $t=50$ for diffusion coefficients a $d=0.05$, b $d=0.1$, and c $d=0.153$, where $K(x)=2.5+ \cos (\pi x)$, $M(x)= 1.2+ \cos (\pi x)$ and $r(x)=1.0$ on $\Omega =(0,1)$, with initial density $u_0=0.5$

Then the positive solutions of (11) lies on the region $[0,\frac{K_0}{2}]\subset \left[ \frac{ K_0}{ 2} \mu ^2 , \frac{ K_0}{ 2}\right]$; $\mu \in {(0,1]}$, which is analytically asserted in Theorem 2. In numerical assertion, from Fig. 11 when $d=0.05$ (Fig. 11a), region of positive solution is $[0.4,0.9] \subset [0,\frac{K_0}{2}]$; when $d=0.1$ (Fig. 11b), region of positive solution is $[0.4,1.2] \subset [0,\frac{K_0}{2}]$ and when $d=0.154$ (Fig. 11c), region of positive solution is $[0,0.4] \subset [0,\frac{K_0}{2}]$. Here $d=0.154=d^*$ is a bifurcation point as described in Fig. 10.

Fig. 12

Dynamical behavior of solutions of (11) for various diffusion coefficients at $t=200$ in (a) solutions of (4) for $u_0=0.5$ (b) solutions of (4) for $u_0=0.9$ with $K(x)=2.5+ \cos (\pi x)$, $M(x)= 1.2+ \cos (\pi x)$ and $r(x)=1.0$ on $\Omega =(0,1)$

Example 7

Considering diffusion model (4) with growth rate function with Allee effect as in (6), then the behavior of solution diverges for certain $d>d^*$. If we increase the initial value then $d^*$ changes for the same sparse function and carrying capacity. For an initial value $d^*$ is not unique and each cases solution exists for $0<d \le d^*$. In (4) with (6), let $K(x)=2.6+ \cos (\pi x)$, $M(x)= 1.1+ \cos (\pi x)$, $r(x)=1.0$. Then numerically we get a bifurcation point $d^*=0.021$ for $u_0=0.5$, see Fig. 12a. For $u_0=0.9$, we get bifurcation point $d^*=0.038$ as in Fig. 12b. The numerical illustrations also indicate that there is an impact of initial density when the Allee effect occurs. The simulating result ensures the theoretical outcome as shown in Theorem 1.

Conclusion

In this paper, we studied a generalized reaction-diffusion model for single species having an Allee effect, where the species has carrying capacity, K, growth rate function, r and sparse function, $M, \; 0<M<K$. For this generalized model, we discussed various Allee effects, and regarding these Allee effects, we discussed different types of growth function and their behavior. While considering a particular growth with the Allee effect (either strong or, weak, or both), we first considered diffusion is ignorable and found the bifurcation point and stable steady state when $d \rightarrow 0$. We also find a stable state involving harvesting. When $d \not \equiv 0$ and not close to zero, we have found a bifurcation point $d^*$ for which the steady states changes to unsteady. For some particular function, we see that there exists one $d^*$ for which when $d\ge d^*$ there exists no positive solution, positive solutions occur only on $0<d \le d^*$. For any greater value of $d^*$, the solution tends to zero and the species is in extinction at a short time. The Fig. 10 is showing the dynamical behavior of $d^*$ as shown analytically in Theorem 1. The existence of a solution (might have multiple solutions) is proved and the region of the positive solution is also found. Some numerical examples were presented to justify the analytic study in a non-empty open domain.

Finally, let us outline some open problems and topics for research and discussion

• 1. Consider (11) with the Dirichlet and the Robin boundary conditions.

• 2. For variable periodic parameters, does a positive periodic solution depend on the diffusion coefficient?

• 3. Originally some of the models described in terms of ordinary differential equations incorporate some kind of delay, for example, the Nicholson’s blowflies equation. Extend the results of the present paper to delay models, explore the influence of delays and compare results with classical and directed diffusion strategies.

• 4. It is known that, for example, for the logistic growth with a periodic carrying capacity, the average of the positive periodic solution over the time cycle is less than the average resource function. Is this true for the periodic model with the directed diffusion?

• 5. Consider the modified model of (11) as

  $$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{\displaystyle \partial u}{\displaystyle \partial t}=d \Delta \left( \frac{u(t,x)}{P(x)}\right) +r(x)u(t,x)\left( 1-\frac{\displaystyle u(t,x)}{\displaystyle K(x)}\right) \left( \frac{\displaystyle u(t,x)}{\displaystyle M(x)}-1\right) ,&{}t>0,~ x\in \Omega ,\\ \nabla (u/P)\cdot n =0,~t>0, &{}x\in {\partial \Omega },\\ u(0,x)=u_0(x),&{}x\in {\Omega }. \end{array}\right. } \end{aligned}$$

  (6.1)

  where P(x) is the resource distribution function of species u(t, x), and other functions have the same interpretation as noted for (11).