Interaction between Allee effects caused by organism-environment feedback and by other ecological mechanisms.

Understanding Allee effect has crucial importance for ecological conservation and management because it is strongly related to population extinction. Due to various ecological mechanisms accounting for Allee effect, it is necessary to study the influence of multiple Allee effects on the dynamics and persistence of population. We here focus on organism-environment feedback which can incur strong, weak, and fatal Allee effect (AE-by-OEF), and further examine their interaction with the Allee effects caused by other ecological mechanisms (AE-by-OM). The results show that multiple Allee effects largely increase the extinction risk of population either due to the enlargement of Allee threshold or the change of inherent characteristic of Allee effect, and such an increase will be enhanced dramatically with increasing the strength of individual Allee effects. Our simulations explicitly considering spatial structure also demonstrate that local interaction among habitat patches can greatly mitigate such superimposed Allee effects as well as individual Allee effect. This implies that spatially structurized habitat could play an important role in ecological conservation and management.

Introduction

Allee effect, the positive relationship between per capita growth rate (individual’s fitness) and population density at low population density [1], has received considerable attention in ecology and conservation because it is directly related to population extinction [2-4] and can incur complicated spatial pattern [5, 6]. Many mechanisms may account for Allee effect such as organism-environment feedback [7, 8], mate-finding limitations [9-11], cooperative defense [12, 13], social dysfunction [1, 14, 15], inbreeding depression [16, 17], and predator avoidance or defense [1, 18, 19]. Allee effect also occurs at a metapopualtion level similar to local populations, which may be because of dispersal cost or colonization difficulty [15, 20–23]. Due to various ecological mechanisms leading to Allee effect, it is necessary to study the influence of multiple Allee effects on population dynamics and persistence. Berec et al. [24] has suggested that two or more Allee effects could affect the dynamics of a population, and more empirical evidence for multiple Allee effects have been provided [10]. We here studied organism-environment positive feedback when Allee effect caused by other ecological mechanism is involved in also.

Organism-environment positive feedback (i.e. habitat modification caused by organism) is an important ecological mechanisms leading to Allee effect [7, 25]. Growth, survival and metabolism of organisms are limited by the physical environments (e.g. temperature, wind, and turbulence) and chemical environments (e.g. oxygen, toxins, and hormones) [26]. Individuals in larger groups will live better if organisms can improve their environment in some ways [10, 27–29]. Such an organism-environment positive feedback is a crucial ecological process in ecosystems [27, 30], which could lead to Allee effect because the fitness of a population may indirectly depend on population density through environmental intermediary [7, 8, 31, 32]. However, the dynamics and consequence of the organism-environment feedback is still unclear when the population is also subjected to Allee effect by other ecological mechanism.

In this paper, we therefore investigate multiple Allee effects respectively caused by the organism-environment feedback (AE-by-OEF) and other ecological mechanisms (AE-by-OM), and specifically analyzed how the superimposed and interacted Allee effects affect the dynamics and persistence of population. Through simulations considering explicitly spatial structure, we also demonstrate the effects of spatial structure on Allee effects.

Models

Mean-field model

With mean-field assumption, the organism-environment positive feedback in patchy habitat can be described by following ordinary equations [8] The variable p is the fraction of patches occupied by the species, and h is the fraction of suitable patches (thus 1 − h indicates the fraction of unsuitable patches, or the fraction of habitat loss), where obviously p ≤ h. Parameter c and e are the colonization and extinction rate, respectively; d is the habitat destruction rate due to human activities or natural causes; λ and μ are the internal (organismic) and external habitat restoration rate, respectively. Parameter c, e, d, λ and μ are all non-negative constants. The system (Eq 1) returns to the classical Levins model when λ = 0 [33]. The organism-environment positive feedback, reflected by parameter λ > 0 in Eq 1, can generate Allee effect (i.e. AE-by-OEF, [8, 27, 34]).

To investigate how the AE-by-OEF functions when the population is also subjected to Allee effect by other mechanism (AE-by-OM), we follow classical method to incorporate AE-by-OM into Eq 1 [11, 35–37]. The term p / (p + a) indicate AE-by-OM, where a ≥ 0 denotes the intensity of the Allee effect [23, 35–38]. In particular a = 0 excludes AE-by-OM (Eq 2 returns to Eq 1), and λ = 0 means AE-by-OEF disappears (see Results). Our model (Eq 2) thus potentially includes the two types of Allee effects when a > 0 and λ > 0.

Probability transition model

In order to reveal the effect of spatial structure on Allee effect, we constructed probability transition model (PTM) based on Markov process, which captures the dynamical change of state probability of each patch over time depending on the state of local neighboring patches. Because there are all three possible states for each patch in our modeling system (i.e. occupied, suitable but empty, and unsuitable), we can only focus on two independent probabilities that a patch is occupied or suitable. Specifically, the dynamical changes of the two probabilities for patch i were given by [39] where p represents the probability that patch i is occupied, and h the probability that patch i suitable. The notation is the average probability that neighboring patch of patch i is occupied, i.e. , where p is the probability that a neighboring patch j of patch i is occupied and z neighborhood size (z = 4 called von Neumann neighborhood and z = 8 Moor neighborhood). In particular Eq 3 will return to Eq 2 when we average p and h over all patches with infinite neighborhood size. This shows that Eq 3 can not only handle situations with various neighborhood sizes but also includes Eq 2. Because it is very difficult to solve the set of ordinary equations (Eq 3) because of its extremely high dimension, we numerically analyzed it on lattice space with periodic boundary.

Cellular automata model

Using a cellular automata model (CAM) with considering stochasticity [12, 40], we also simulated the system mentioned above on lattice space. Each patch is either occupied (denoted by 2), empty but suitable (denoted by 1), or unsuitable (denoted by 0). We define state transition rule (i.e. update rule) corresponding to above models. An unsuitable patch becomes suitable (0 → 1) at rate (habitat restoration), where represents the fraction of occupied patches in its neighborhood; an empty but suitable patch becomes unsuitable (1 → 0) at rate d (habitat destruction) or occupied (1 → 2) at rate (recolonization); and an occupied patches becomes unsuitable (2 → 0) at rate d (habitat destruction) or empty but still suitable (2 → 1) at rate e (local extinction). The update rules are summarized in Table 1. We used the periodic boundary condition and synchronous update in our simulations [12, 41].

Table 1

Transition rate for the lattice model.

State transition	Rate
0 → 1	λp˜+μ
1 → 0	d
1 → 2	cp˜p˜p˜+a
2 → 0	d
2 → 1	e

Note: Here 0, 1, and 2 denotes that a patch is unsuitable, suitable but empty, or occupied, respectively. represents the fraction of occupied patches in neighborhood. Other parameters are the same as in the text.

Results

Organism-environment positive feedback will individually incur Allee effect (i.e. AE-by-OEF) when λ > (μ + d)2 / d, of which inherent characteristic depends on model parameters (Eqs 1 or 2 with a = 0). For example, the population will be subjected to weak, strong, and fatal Allee effect in turn with increasing habitat destruction rate d (see Table 2, Figs 1 and 2). Individual AE-by-OM in our model is strong when a < (δ–β/(β + 1))2/(4δ) and δ < β/(β + 1) (where δ = (e + d)/c, α = λ/d, β = μ/d), and fatal otherwise (Fig 1). This system becomes more complicated when the two types of Allee effects are superimposed and interacted. Through analyzing per capita growth rate when habitat dynamics is at equilibrium (i.e. dh/dt = 0), we found that the combinations of strong AE-by-OM and either strong or weak AE-by-OEF, or even organism-environment feedback without Allee effect (λ > (μ + d)2/d), are all able to lead to strong Allee effect with higher Allee threshold (Fig 3A, 3C and 3E). These combinations can also incur fatal Allee effect which means population extinction in any case (Fig 3B, 3D and 3F). Phase plane analysis also confirmed these results (Fig 4).

Table 2

Condition for Allee effect caused by organism-environment positive feedback.

Characteristic of Allee effect	Parameter condition
No	λ<(μ+d)^2d,δ<ββ+1
Weak	λ>(μ+d)^2d,δ<ββ+1
Strong	λ>(μ+d)^2d,ββ+1<δ<(α−1)^2+βα
Fatal	λ>(μ+d)^2d,δ>(α−1)^2+βα)

Note: δ = (e + d)/c, α = λ/d, β = μ/d)

Fig 1

Interaction between organism-environment positive feedback and Allee effect caused by other ecological mechanisms.

Organism-environment feedback results in weak, strong, and fatal Allee effect as well as the situation of no Allee effect (denoted by green, red, blue, and cyan area, respectively). Allee effect caused by other ecological mechanisms is strong (left of black line) and fatal (right of black line) Allee effect. The interaction of the two kinds of Allee effects leads to strong (below yellow line in (a), above yellow line in (b) and (c)) and fatal (above yellow line in (a), below yellow line in (b) and (c)) Allee effect. For the parameters, We take: μ = 0.03 and e = 0.1 in (a)-(c); c = 0.5 in (a) and (b); d = 0.05 in (b) and (c); λ = 0.8 in (a); c = 0.3 in (c); Point A-F in (a) correspond to panel A-F in Fig 3.

Fig 2

Definitions of various Allee effects.

Note: function f(p) is per capita growth rate of population, and p* > 0 its maximum point. In schematic diagram, the horizontal axis is population size or density, and vertical axis is per capita growth rate of population.

Fig 3

Relationship between population size and per capita growth rate.

Green and red lines represent organism-environment feedback (a = 0) and Allee effect caused by other mechanism (d = 0), respectively; blue lines represent the combination of the two types of Allee effects. Panel A-F corresponds to point A-F on Fig 1a.

Fig 4

Phase plane analysis.

The curves indicate zero isoclines of corresponding system. Filled and empty circles represent stable and unstable equilibriums, respectively. The left column correspond to Fig 3A and 3B (also point A and B in Fig 1a).

The stability analysis also roots for the above results (Appendix). There is only one boundary equilibrium (0, β / (β + 1)) for both Eqs 1 and 2, which represents population extinction. The equilibrium is stable for Eq 1 if δ > β / (β + 1) and unstable otherwise, while it is always stable for Eq 2 (Appendix). For the both equations, when the boundary equilibrium is stable, the two equations have either two interior equilibriums (the smaller unstable and the larger stable) or no interior equilibrium (Appendix), which means strong or fatal Allee effect. When the boundary equilibrium is unstable for Eq 1, the equation has unique interior equilibrium, meaning weak Allee effect or no Allee effect (Table 2). Totally, multiple Allee effects either enlarge Allee threshold or change the inherent characteristic of Allee effect, and such enlargement and change will be enhanced dramatically with increasing individual Allee effects.

Our simulations based on both PTM and CAM with explicitly spatial structure demonstrated that spatially local interaction between organism and environment can make both individual and superimposed Allee effects mitigate greatly and even disappear (see Fig 5, notably a = 0 and λ = 0 respectively represent individual AE-by-OEF and AE-by-OM). The spatial heterogeneity of population distribution could play important role for the consequence of the local effect. Although low-density local regions are subjected to Allee effect and make sub-populations extinct, high-density local regions save the population globally. In our simulations we obviously observed the aggregated distribution pattern with distinct spatial boundary not only for population but also suitable habitat (Figs 6 and 7). However, because colonization and habitat restoration incurred by organism are limited, local interaction leads to lower global population density than the predictions from mean-field model (Fig 5). This shows that spatially limited interaction is good for preventing population extinction although it decreases population density potentially.

Fig 5

The dependence of population dynamics on its initial state when population dispersal locally.

(A) and (B) based on probability transition model, (C) and (D) based on cellular automata model (neighborhood size n = 4 for A and C, n = 8 for B and D). Red areas indicate the initial values leading to the population extinction, while green areas represent the initial values from which the population persists stably. The simulations ran on 100 × 100 lattices with periodic boundary condition. Parameter λ = 0.4, μ = 0.1, d = 0.1; and other parameters are the same as in Fig 1a.

Fig 6

The distribution of population subjected to Allee effects on two-dimensional space, which are snapshots of 100 × 100 lattice when simulations arrive at 3000 time step.

The color hot indicates the probability of patch occupancy (A) and suitable patch (B). Parameters: d = 0.1 and other parameters are the same as in Fig 1a.

Fig 7

Spatial distribution of population subjected to Allee effect.

Red, green, and blue colors indicate occupied cells, suitable but no occupied cells, and destroyed (unsuitable) cells. The simulations ran on 100 × 100 lattices with periodic boundary condition and 4 neighborhood. Parameters are taken as a = 0, 0.1, and 0.15 from left to right, λ = 0.3, 0.5, and 0.8 from bottom to top, respectively. Other parameters values are the same as in Fig 4.

Discussion

Organism-environment positive feedback is ubiquitous in real world, especially in severe environment such as dry land [42-44] and intertidal zone [27]. Allee effect is a distinct feature caused by the positive feedback. Due to diverse mechanism leading to Allee effect, it is necessary to estimate the effect of Allee effect caused by the organism-environment feedback on population dynamics and persistence when another mechanism generating Allee effect also act on the population [10, 24]. We here study theoretically the interaction of the two type of Allee effect to reveal their effect on population dynamics and persistence. The combination of multiple Allee effects could largely increase the extinction risk of population either due to the enlargement of Allee threshold or the change of inherent characteristic of Allee effect (Fig 3). We further found that the spatial range of the organism-environment feedback also profoundly affect the appearance of Allee effect. Both combined and individual Allee effect greatly mitigates when interaction is limited locally (Fig 5). This implies spatial scale could play an important role in ecological conservation [12, 15, 39, 45].

Mathematical models are useful in helping to understand the dynamic processes of involved populations and in making practical predictions [46-49]. In general, we can classify modeling methodologies for Allee effect into two distinctive categories: phenomenological model and mechanistic model. The former directly incorporates Allee effect into the dynamic models of populations, and the modelers are often interested in the consequence of Allee effect (e.g. [35, 36, 50]); the later models the underlying mechanism that incurs Allee effect, and generally tests how Allee effect is generated (e.g. [8, 51, 52]). In our model, organism-environment feedback which can result in Allee effect was described in details (i.e. mechanistic modeling), while another Allee effect was phenomenologically modeled using a frequently-used classical method (Eqs 1 & 2). Therefore, our models totally include the two kinds of modeling methodologies for the two types of Allee effects.

Much attention has been paid to Allee effect not only because it is an important ecological problem itself, but also it is directly related to population extinction and thus significant to ecological conservation [53, 54]. Recent years, multiple Allee effects, resulting from various mechanisms but acting on the same population, have become a hot point for study, with examples from plants, invertebrates and vertebrates, from natural and exploited populations, and from terrestrial and marine ecosystems [24, 55]. Mathematical models can also increase our understanding of the interaction among multiple Allee effects [10, 24]. We here constructed a set of mathematical models for organism-environment feedback system subjected to Allee effect from other mechanism, without and with spatial structure. Agreeing with the results of Berec et al. [24], our model demonstrated that multiple Allee effects could be stronger than the simple sum of involved individual Allee effects (Fig 3). This implies that the interaction of multiple Allee effects could be nonlinear (i.e. nonadditivity), thus we need to carefully refer to those predictions from linear models in conservation and management.

Reaction-diffusion equations are used traditionally to model the spatial dynamics of populations [5, 6, 56–62]. However, it is difficult to handle various neighborhood sizes (e.g. changeable dispersal kernel) for them. We thus constructed our model system based on lattice space in order to conveniently assess the ecological significance of spatial scale [39]. Spatial scale is an important ecological factor in ecology, which plays a key role in the coexistence of competitive species [39, 63–65], biodiversity persistence [66, 67], epidemic spread [68-70], evolutionary consequence [12, 71], pattern formation [13, 39, 72] and so on. We here found that spatial scale has important significance for the appearance of Allee effect. When organism-environment interaction is limited locally, Allee effect will mitigate although the global size of population decrease (Fig 5). The emergence of spatial pattern resulting from local interaction could account for this mitigation of Allee effect (Figs 6 and 7). Although low-density local regions are subjected to Allee effect and make sub-populations extinct, high-density local regions save the population globally. Pattern formation could prevent population extinction [13, 39, 73].

Appendix: Stability analysis for our system

When a = 0, the Eq 2 returns to Eq 1. When habitat dynamics is at equilibrium (i.e. dh/dt = 0), the per capita growth rate of population (i.e.) is From the condition of Fig 2, we deduce that organism-environment positive feedback will individually incur Allee effect (AE-by-OEF) when λ > (μ + d)2/d. We summarize parameter conditions of various Allee effects (Details are shown in Table 2). When λ = 0, the AE-by-OEF disappears and the population is only subjected to AE-by-OM in Eq 2. Through simple calculation, the interior equilibrium must be the solution of following equation Thus, When λ = 0, Eq 2 has at most two interior equilibriums, and necessary conditions for the existence of two interior equilibriums are a < (δ–β/(β + 1))2/(4δ) and δ < β/(β + 1). When δ > β/(β + 1), there is no interior equilibrium.

When λ > 0 and a > 0, Eq 2 includes the two types of Allee effects. Eq 2 has one boundary equilibrium E0(0, μ/(μ + d)), and the eigenvalues of the Jacobian matrix at E0 are − (e + d) and − (μ + d). Thus, the boundary equilibrium E0(0, μ/(μ + d)) is locally asymptotically stable.

The zero isocline for dp/dt = 0, denoted by where δ = (e + d)/c.

And the zero isocline for dh/dt = 0, denoted by The interior equilibriums of system (2) are the points of intersections of two curves C1 and C2. It is convenient to express this as where A = δ– 1 + (μ +d)/λ, B = δ(μ +d)/λ +a–μ/λ, C = a δ(μ +d)/λ.

In order to determine whether there exists interior equilibriums, let After direct computation, we have f(1) > 0, f(0) > 0, f′(1) > 0, f″(1) > 0, f‴(p) >0.

We can deduce that f″(p), f′(p) and f(p) are monotonically increasing functions for p ∈ (1, +∞). From the monotonicity of f(p), we can see that the Eq (A5) has no positive root when p > 1. That is to say, all positive solutions of Eq (A5) exist in (0, 1).

Note Δ = A2 – 3B, if Δ ≤ 0, f(p) is strictly increasing in (0, +∞) since f′(p) ≥ 0, which yields f(p) ≥ f(0) = C > 0. And Eq (A5) has no positive root.

If Δ > 0, then the equation f′(p) = 0 has two roots, . It is obvious that f(p) is strictly decreasing in (p1, p2) and increasing in (−∞, p1) and (p2, +∞).

In addition, we get

Then, we can easily obtain the following results:

If A > 0, B > 0, and A2 > 3B, or A2 ≤ 3B, then system (2) has no positive equilibrium.

If either of the following inequalities holds: B < 0, f(p2) < 0 or A < 0, B > 0, A2 > 3B, f(p2) < 0, system (2) has two positive equilibriums.

If either of the following inequalities holds: A > 0, B < 0, f(p2) = 0; or A ≤ 0, A2 > 3B, f(p2) = 0, system (2) has an unique positive equilibrium.

Denote the right-hand sides of the two equations of system (2) by f(p, h) and g(p, h) respectively. At any interior (positive) equilibrium (p*, h*), we can write the Jacobian matrix as where subscripts denote partial derivatives. Obviously, three elements of J satisfy f > 0, g > 0 and g < 0.

Implicit differentiation of the equation for the p− null line, f(p, h1(p)) = 0, gives the slope

Similarly, from g(p, h2(p)) = 0, we obtain

We use Eqs (A10) and (A11) to write the determinant as . Then if f < 0 and , therefore the Routh-Hurwitz conditions are satisfied and the equilibrium is stable [74]. If system (2) has two positive equilibriums, we denote them as (p1*,h1*) and (p2*,h2*) (without loss of generality, we suppose that p1* < p2* (see Fig 4). By the properties of h1(p) and h2(p), we can obtain , (p2*)2 > aδ, , thus the equilibrium (p2*,h2*) is asymptotically stable. Whereas and the equilibrium (p1*,h1*) is a saddle point and unstable (see Fig 4).