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][3][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][10][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][21][22][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][28][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.

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][36][37].

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 i represents the probability that patch i is occupied, and h i the probability that patch i suitable. The notation e p i is the average probability that neighboring patch of patch i is occu- where p j 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 i and h i 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 le p þ m (habitat restoration), where e p 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 ce pe p=ðe p þ aÞ (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].
Note: Here 0, 1, and 2 denotes that a patch is unsuitable, suitable but empty, or occupied, respectively. e p represents the fraction of occupied patches in neighborhood. Other parameters are the same as in the text. https://doi.org/10.1371/journal.pone.0174141.t001

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).
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.

Discussion
Organism-environment positive feedback is ubiquitous in real world, especially in severe environment such as dry land [42][43][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][47][48][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][57][58][59][60][61][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][64][65], biodiversity persistence [66,67], epidemic spread [68][69][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]. 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). https://doi.org/10.1371/journal.pone.0174141.g004

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. 1 p dp dt ) 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 E 0 (0, μ/(μ + d)), and the eigenvalues of the Jacobian matrix at E 0 are − (e + d) and − (μ + d). Thus, the boundary equilibrium E 0 (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 C 1 and C 2 .
If Δ > 0, then the equation f 0 (p) = 0 has two roots, It is obvious that f(p) is strictly decreasing in (p 1 , p 2 ) and increasing in (−1, p 1 ) and (p 2 , +1).  In addition, we get Then, we can easily obtain the following results: 1. If A > 0, B > 0, and A 2 > 3B, or A 2 3B, then system (2) has no positive equilibrium.