The effect of colonization dynamics in competition for space in metacommunities

One of the main factors that determines habitat suitability for sessile and territorial organisms is the presence or absence of another competing individual in that habitat. This type of competition arises in populations occupying patches in a metacommunity. Previous studies have looked at this process using a continuous-time modeling framework, where colonizations and extinctions occur simultaneously. However, different colonization processes may be performed by different species, which may affect the metacommunity dynamics. We address this issue by developing a discrete-time framework that describes these kinds of metacommunity interactions, and we consider different colonization dynamics. To understand potential dynamics, we consider specific functional forms that characterize the colonization and extinction processes of metapopulations competing for space as their limiting factor. We then provide a mathematical analysis of the models generated by this framework, and we compare these results to what is seen in nature and in previous models.


Introduction
Metacommunities are a theoretical concept expressed as a set of local communities linked by dispersal of multiple species ). Using metacommunities' integration of local and regional processes, we can explore community processes that occur at a more regional scale, such as species distribution in the landscape or patch occupancy. This concept has seen an increase in interest in recent years (Grainger and Gilbert 2016;Leibold and Chase 2017;Guichard 2017), including possible applications in conservation biology ).
An important example of this framework is competitive metacommunities, which are used to explain coexistence mechanisms among competing populations ). These competitive metacommunities allow us to explain species richness in several communities, as they explain how local dynamics can affect regional presence of species (Harrison and Cornell 2008). Some examples where this concept has been used to explain species richness are plant communities (Kuglerova et al. 2015) or crustacean communities (Cottenie and Meester 2003).
These competitive interactions are caused by overlaps in the niches of the populations. Habitat suitability is a key component of a population's niche, which determines its ability to survive in a given environment (Kearney 2006). This suitability is defined by environmental conditions and the population's life and evolutionary history. One of the fundamental components of habitat suitability is the presence or absence of another individual in a given region of space. This is more critical in sessile organisms such as plants, where the presence of an individual in a specific point of space solely limits the possibility of other organisms to exist in that same point (Aarssen et al. 2006).
If one focuses only on the presence or absence of a population in space, the environment can be understood as an arrangement of patches distributed in space, where a patch is not suitable if it's occupied, and suitable otherwise. This patch occupancy idea was originally studied for metapopulations by Levins (1969). These metapopulations consist of populations of a single species distributed through a landscape of patches, which are then connected by dispersal. Later, Hastings (1980) and Tilman (1994) expanded these ideas to competitive metacommunities, where multiple metapopulations occupy patches in space and can be excluded from those patches, leading to local extinctions. These local extinctions can be caused by dispersal failure (Mouquet and Loreau 2003), high competition strength , or disturbance (Sousa 1984). The latter is particularly useful when studying how disturbance affects hyperdiverse ecosystems such as tropical forests, where single population dynamics are not useful to understand coexistence of many species (ter Steege et al. 2019).
One of the key processes a metapopulation goes through is the colonization of unoccupied patches, which has an impact on the fate of the metacommunity structure (Lobel et al. 2009;Howeth and Leibold 2010;Wisnoski et al. 2019). These colonization processes might be impacted by the presence of competitors and can lead to different dynamics of patch occupancy and overall diversity at a landscape level. Although empirical evidence of the importance of the different dynamics has appeared, the theoretical ideas behind those results have been less explored.
These processes can be easily described in a discrete-time modeling framework. In discrete-time models, the state of the system is measured at discrete, regular time intervals, compared to continuous-time models, where the state of the system is continuously measured. This difference allows us to model underlying processes of the dynamics as separate substeps and focus our analysis on the step of the colonization dynamics. Discrete-time models have been applied to metacommunities before, putting more emphasis on the local competition dynamics (Wilson 1992), coexistence mechanisms (Shoemaker and Melbourne 2016), or realism with a simulation based analysis (Thompson et al. 2020). In this work we put more emphasis on the colonization dynamics and how different colonization processes affect the metacommunity structure. This approach has been recently considered to study metapopulations (Marculis and Hastings 2021), and here we expand this idea to metacommunities.
In this paper we construct and analyze a simple discretetime metacommunity model using difference equations that captures the colonization and local extinctions separately and use it to analyze how different colonization, competition, and disturbance effects impact the persistence of several metapopulations in an environment. In the next section, we describe the framework that considers colonization and extinction processes separately in time, and we use functional forms that describe the patch occupancy in metacommunities where space availability is the limiting factor. In Section 3 we present a mathematical analysis of the model with the given functional forms. Finally, in Section 4 we discuss this analysis, focusing on how it relates to the continuous-time models, and what other dynamics can be explored using this framework.

Model and methods
To describe competition for space, we expand on the metapopulation framework presented by Marculis and Hastings (2021). Suppose we have a system with n species, where the i-th species at time t occupies a proportion p (i) t of the available patches in a given region. Then, at the next time step, the proportion of occupied patches will change following two processes in sequence: Colonization of unoccupied patches and local extinctions. Suppose species i colonizes unoccupied patches following a function g i ( t ) , where t is a vector where the j-th entry corresponds to p (j) t . After that, species i faces local extinctions following a function f i (g 1 ( t ), … , g n ( t )) . Then, if and are vectors where the i-th entry corresponds to g i and f i respectively, then the proportions of the species at time t + 1 follow the difference equation: In order to use this framework in a competition context, we will arrange our metacommunity in the same way as Hastings (1980) and Tilman (1994). Suppose species are ordered by competitive capacity, i.e., species i is outcompeted by species j for j < i . In addition, assume that the growth and extinction processes of patches by species i are determined by the proportion of patches occupied by species j ≤ i . Then, Equation 1 simplifies to the system of equations: Next, we focus our attention on reasonable functions that describe the colonization process. To do this, we follow the ideas of Hastings (1980) and Marculis and Hastings (2021) of using discrete-time competition models to describe the proportion of patches colonized.

Leslie-Gower model
The Leslie-Gower model, further analyzed in Cushing et al. (2004), is an extension to the single-population Beverton-Holt to include the effect of competing species over the carrying capacity of the population of interest. Similar to the model presented by Hastings (1980), which extends logistic growth to competitive metacommunities in continuous time, we propose the Leslie-Gower model as an extension of the Beverton-Holt to competitive metacommunities in discrete time. In this study, we frame this model as follows: Suppose species i has a colonization (1) t+1 = t . (2) rate d i > 1 , then the proportion of patches colonized by species i at time t is given by

Ricker competition model
In some cases, colonization processes may not follow a monotonic growth, but rather present periods of overcompensation due to interference between colonizers (Goubault et al. 2005) or due to boom and bust dynamics (Simberloff and Gibbons 2004). This formulation expresses the idea that when too many individuals simultaneously attempt to colonize an empty patch, none may survive. So if the number of occupied sites is too high, it is possible that the probability of colonization of empty sites decreases. In this case, modeling the process as a Ricker model makes biological sense. Similar to the Leslie-Gower model, the Ricker model can be expanded to include interspecific competition (Cushing et al. 2004). In addition, we normalize the Ricker model to make sure that our values remain in [0, 1], as they represent proportions. Suppose that species i has a colonization rate d i > 0 . Then, the proportion of patches colonized by species i at time t is given by

Leslie-Gower with allee effect
Alternatively to the aspects described by the Ricker description, too few colonizers may make colonization rates per occupied patches lower, which would resemble an Allee effect in population dynamics. We can describe this process by expanding the Leslie-Gower model to include an Allee effect term (Chow and Jang 2014). Suppose species i has a colonization rate d i > 1 and a probability 1∕m i of finding a new patch to colonize. Then, the proportion of patches colonized by species i at time t is given by

The extinction function
Next, we describe the extinction process as follows: Local extinctions of populations will occur due to disturbance or competitive exclusion (Hastings 1980). Suppose the environment presents a homogeneous disturbance that causes a proportion e of patches to become extinct. In addition, suppose the presence of better competitions causes extinctions with a strength c. Then, if x i = g i p (1) t , … , p (i) t , the extinction process is described by the function

Model analysis
Using the framework presented in System 2, we fix the extinction function to follow Equation 6 and analyze the asymptotic dynamics of the framework when the colonization function follows Equations 3, 4, and 5. The analyses of said dynamics can be achieved by using a standard linear analysis of the fixed points. We consider then how the different parameters d i , e, m i , c affect the stability of the fixed points through bifurcation analysis, and if nonlinear attractors can arise from the given colonization functions. In the following section, we present the results of this analysis and their proofs are found in Appendix A.

Leslie-Gower colonization
The Leslie-Gower colonization process provides a Lotka-Volterra-type competition. This allows us to compare it directly to the model presented in Hastings (1980). We first check its convergence to a single configuration.
Theorem 1 For any given parameters e, c, d i for i = 1, … , n , the System 2 with colonization function given by Equation 3 and extinction function given by Equation 6 has a globally asymptotically stable equilibrium in the hypercube [0, 1] n . The metapopulations with positive equilibria p(i) also satisfy that Theorem 1 shows that, using a Leslie-Gower type of colonization process, we recover the original results as Hastings (1980). The last inequality of Theorem 1 implies that there must be a trade-off between competitive advantage and colonization ability, regardless of the amount of disturbance in the environment. In general, in this type of dynamics one would expect that the worse competitors are better colonizers.
We can see how this behavior looks like starting from an almost empty landscape in Fig. 1. Notice that the dynamics towards the equilibrium are not monotonic, and the better competitors that are not the best (in Fig. 1, i = 2 ) will even invade a larger proportion of the patches than the top competitor before being outcompeted by it and start decreasing to their equilibrium value.
Compare this time series with a similar one produced with the continuous model of Hastings (1980). This model is rescaled to follow the following equation: and the produced time series is presented in Fig. 2. By comparing this time series with that presented in Fig. 1, we can see that they have a similar qualitative behavior, especially in the fact that each of the metapopulations converge to their equilibrium value.
Our model also recovers the results by Tilman (1994), where in the presence of disturbance (i.e., e > 0 ), there is Time series for a system with 3 species using a colonization process given by Equation 3 and extinction process given by Equation 6. In this simulation we use e = 0.2, c = 0.01, d i = 1.5 i . Notice how for all species, the time series converges to a positive equilibrium value, following Theorem 1

Fig. 2
Time series for a system with 3 species using the model presented in Equation 7.
In this simulation we use e = 0.2, c = 0.01, d i = 1.5 i . If we compare it to the time series presented in Fig. 1, we get that they get qualitatively similar results, and they both converge asymptotically to their equilibrium values always a proportion of patches at equilibrium that will not contain any population. This result can be written as follows.
Proposition 1 In the System 2 with colonization function given by Equation 3 and extinction function given by Equation 6, regional coexistence of more than one metapopulation is only possible in the presence of disturbance.
We can also see how total proportion of patches occupied grows as we introduce more metapopulations into our system in Fig. 3. Notice that although numerically our choice of parameters allow for many species to persist (84 in total), most of the available patches are already occupied by the first top competitors. In a real scenario where the number of available patches is finite, the lesser competitors might not be able to persist at all, significantly reducing the effective number of coexisting species.

Ricker colonization
In the case of the colonization function given by Equation 4, the persistence equilibrium for a given population might not be stable, which could lead us to a cyclic or chaotic behavior similar to that of the classical Ricker model (Ricker 1954). To further study this behavior and keep our results analytically tractable, we limit the following analysis to the simplest case of two competing species (i.e., n = 2 ). This gives us the following result. We can see this effect when comparing the bifurcation diagrams for species 1 with different intensities of disturbance in Fig. 4. Notice that as disturbance increases, the period doubling bifurcation occurs at a later point and allows the proportion of occupied patches to stay stable for higher colonization rates. Notice as well that before the bifurcation point, we see a decrease in the total number of patches occupied at equilibrium. This would suggest a mechanism of interference between the different populations trying to colonize the same patches.
If we vary d 1 and e, we get bifurcation diagrams for i = 2 as shown in Fig. 5. Notice that when the equilibrium proportion for species 1 is stable, species 2 seems to behave stable as well for small colonization rates. However, there exists certain threshold (presumably the period doubling bifurcation) under which species 2 collapses, leading to the regional extinction as the only possibility. In the case species 1 shows chaotic behavior, this behavior passes on to species 2, which presents high Bifurcations diagrams for p(1) using a Ricker colonization process, with e = 0.2 (left), and e = 0.5 (right). This figure shows how disturbance produces a stabilizing force into species 1, as the period doubling bifurcation occurs at a higher colonization rate d 1 compared to the case with no disturbance shows that disturbance tends to stabilize p(2) . Notice that there is a limit on how well of a colonizer can species 2 be, as for d 2 bigger than some value around 7.5 will collapse into the regional extinction equilibrium only variation in its proportions, before collapsing for high enough colonization rates. As expected, at low competition parameters, disturbance is a destabilizing factor for species 2, as can be seen in Fig. 5c. Notice as well that if species 1 is a better colonizer (when d 1 = 3 ), the chaotic attractor changes less abruptly, which follows the proof of Theorem 2, which shows that the colonization rate of species 1 is a stabilizing force for species 2. Also notice that we restrict our choice of c to a value far away from the threshold we found in Theorem 2 as bigger values seem to not provide a positive equilibrium value, suggesting that in reality the only biologically feasible condition for whether disturbance has a destabilizing effect for species 2 is if In the case of Figs. 5a and 5c, d 1 does not satisfy this inequality, and thus, an increase in disturbance will cause more unstable dynamics. On the other hand, in Figs. 5b and 5d, d 1 does satisfy this inequality, and an increase in disturbance leads to more stable dynamics.
We show how the dynamics between the two metapopulations unfold in the time series in Figs. 6 and 7. Figure 6 shows that when the bottom competitor is too strong of a colonizer, it will present overshoots of the patches occupied (bottom). Notice that species 2 presents an overshooting process, where they invade a higher proportion of patches at the beginning, followed with a sudden collapse of the metapopulation. These time series also show that an increase in disturbance prolongs the collapse of the metapopulation at the beginning of the time series, followed by a sudden collapse and regional extinction. This is an indicator of boom and bust dynamics. In addition, notice that since d 1 is small, disturbance leads to more unstable dynamics for the bottom competitor. Figure 7 shows the effect that the oscillatory behavior of species 1 has over species 2. In the case where species 1 does not present oscillatory behavior, the proportion of patches occupied by species 2 at equilibrium is low. On the other hand, when species 1 has an oscillatory behavior, it allows species 2 to access more patches. An interesting observation is that this phenomenon occurs in synchrony, i.e., the bottom competitor will occupy more patches when the top competitor does. This could be explained by the effect of interference the top competitor has at high colonization intensities.

Allee effect
A more general version of the model provided by Equation 5 for two competing populations was studied in Chow and Jang (2014). Based on the framework provided by System 2, we can get the following result for any number of species.
Theorem 3 For any given parameters e, c, d i , m i for i = 1, … , n , there exists a value i (that depends on p (1) 0 , … , p (i) 0 ) such that, for the System 2 with colonization function given by Equation 5 and extinction function given by Equation 6 we get that where p(i) + is the biggest positive fixed point whenever it exists. In addition, a necessary condition for p(i) + to exist is that We can see in Fig. 8 how these dynamics play out. Notice that with similar parameters, the total amount of metapopulations that the environment can sustain are less than when using Equation 3. In this case, we still have the same competition/colonization trade-off as in the Leslie-Gower model.
However, being a better colonizer also consists of being better at finding suitable patches, unlike in the previous cases, where being a better colonizer consisted only of being faster at colonizing unoccupied patches. This can be seen in the case of species 3, which although a worse competitor than species 2, can persist due to a better ability of finding a patch to occupy (by having a probability of finding a patch one order of magnitude higher). This figure also shows that this capability is also dependent on the initial conditions, as in the bottom figure a lower initial occupancy of species 3 will eventually cause a regional extinction.

Discussion
In this paper we develop a simple yet general model for patch occupancy of metacommunities, which is then used to study competition for space in several scenarios. In the simplest case, we use Leslie-Gower dynamics as the colonization function, which has been shown to be a discrete-time analogous model to the classic Lotka-Volterra competition model (Liu and Elaydi 2001). We find that using this function as our colonization function generates the same results as those presented in Hastings (1980) and Tilman (1994). Hence, for one description of competition and colonization (Equation 3), our model is a discretetime analogue of the model used in these previous papers.
However, the generality of the framework allows us to consider alternative colonization functions that capture different types of dynamics. The probability of colonization of empty patches may not depend linearly on the number of occupied patches if successful colonization is greatly influenced by the number of individuals arriving. If there is interference, then too many occupied locations will reduce the colonization probability, leading to overcompensation which we described using the Ricker model. Too few occupied patches may also reduce the likelihood of colonization, which we described using a model incorporating an Allee effect. These processes can be naturally included in a discrete-time model.
When we use the Ricker function as our colonization function, the model dynamics are drastically different when the colonization rates of the metapopulations are high. In our runs, we see a boom and bust of the lower competitor followed by a regional collapse of the metapopulation, or a cyclic behavior in its patch occupancy. The boom and bust behavior in the lower competitor when colonization rates are higher can be explained either due to resource depletion at high densities, or due to a high number of invasions that generate patches with smaller densities and thus a lower competitive capability (Simberloff and Gibbons 2004;Sagata and Lester 2009).
These biological explanations make sense in the context of the model. Since the top competitor is not affected by the dynamics of the bottom competitor, these patches initially Fig. 7 Time series for p t using a Ricker colonization process, with c = 0.01 , and e = 0.2, d 1 = 5, d 2 = 2 (top), and e = 0.5, d 1 = 5, d 2 = 2 (bottom). When the top competitor metapopulation is oscillating, this allow the bottom competitor to colonize a bigger proportion of patches in synchrony with the top competitor. These time series show that an increase in disturbance allows for more stable dynamics to occur, thanks to the higher value of d 1 ◂ occupied by the bottom competitor will soon be occupied by a better competitor, which will cause a collapse of the metapopulation. As seen in Fig. 5, if the colonization rate is too high, this "bust" process will cause a regional collapse of the metapopulation, without any possibility for persistence. In other scenarios, persistence of the bottom competitor is enhanced thanks to oscillatory behavior of the top competitor metapopulation. This type of persistence induced by nonlinear attractors has been previously studied at the community level (Huisman and Weissing 1999) and metacommunity level (Koelle and Vandermeer 2005) using continuous time models. Here we have been able to find it at the metacommunity level using a discrete-time model and possibly caused by interference of the top competitor with itself. In the context of competition for space, this cyclic behavior shows up in competition caused by allelopathy (Nakamaru and Iwasa 2000;Lenski and Riley 2002). This behavior also resembles the effects of competition between two populations over a single resource, which is possible when the top competitor has oscillatory behavior (Armstrong and McGehee 1980;Adler 1990). This suggests that, when considering space as a resource, the metapopulations behave in a similar manner to competing populations for a single resource.
In the case of the colonization process with an Allee effect, we see the same colonization-competition trade-off as in the Leslie-Gower colonization process. However, being a good colonizer in this context not only consists of having high reproductive or dispersal rates (having a higher d i ), but also being good at finding a suitable patch to colonize (having a smaller m i ) which mitigates the negative effects that occur at low patch occupancy. This scenario is important in highly fragmented habitats, where search efficiency is as important as dispersal capabilities (Niebuhr et al. 2015). In this similar context, the dependency on the initial conditions of the system can lead to alternate stable states, where the amount of species in the environment depends on this probability of finding a good patch, which is linked to the fragmentation of the habitat (Buenau et al. 2007).
Using these functional forms may seem an artificial way to make metacommunity dynamics resemble single population dynamics in discrete time. However, evidence of such behaviors occurring in natural systems exists. In the case of Leslie-Gower, we find that the dynamics are equivalent to those observed in continuous time models (Hastings 1980;Tilman 1994), which represent an important element of literature in ecological theory. In the case of the Ricker colonization, evidence of interference in colonization processes has been observed in several species of insects, especially when intraspecific competition between larval stages is significant (Dye 1982;Balmer et al. 2009). In the case of the Leslie-Gower with Allee effect, habitat degradation can lead to regional extinction of metapopulations, which we find can be further intensified when metapopulations are competing for space Bascompte and Sole (1996). In addition, more degraded habitats have been observed to have less diversity at the metacommunity level than better suited habitats (Halme et al. 2013).
Previous works have shown that dispersal between patches in a population leads to more stable dynamics when this population follows a Ricker growth function (Doebeli 1995). In the metacommunity context, we observe that when our metapopulations follow a Ricker-like colonization function, increasing dispersal can have varying effects, depending on the amount of disturbance that leads to local extinctions. In our analysis the interaction effects of dispersal and disturbance over the top competitor follow similar dynamics to those found in Marculis and Hastings (2021), where high dispersal intensities lead to instabilities in the system. However, when a metapopulation is outclassed in competition for space, these interactions lead to more complex dynamics, and dispersal may either stabilize the metapopulation or destabilize it, following certain thresholds which depend on the dispersal of the better competitors, and the disturbance and competition intensities.
With regard to the effect of competition, we see that in the three models analyzed, competition directly hinders the amount of space available for each species at equilibrium. One important aspect of how we modeled competition is that competition only occurs after populations have settled in their respective patches; an important feature of territorial animals (Stams and Krishnan 2001). However, this behavior can be also seen in plant communities, where competition for space occurs at the root level (Gersani et al. 2001). We suspect that the dynamics will differ if competitive exclusion occurs before the colonization, especially in colonization models with alternative stable states such as the Leslie-Gower with Allee effect.
Notice that in this work we have considered the events of colonization and extinction to occur in a different order than the events in Marculis and Hastings (2021). Since we only have two events occurring (colonization and extinction), Fig. 8 Time series for a system with 3 species using a colonization process given by Equation 5 and extinction process given by Equation 6. In this simulation we use e = 0.2, c = 0.1, d i = 10 i , m i = 1 10 i and start with (top) p (1) Notice how a lower initial density for species 3 causes it to go extinct. In both simulations species 2 was not able to persist in the landscape, which could mean that their threshold go to their positive equilibrium needed to be even higherTime series for a system with 3 species using a colonization process given by Equation 5 and extinction process given by Equation 6. In this simulation we use e = 0.2, c = 0.1, d i = 10 i , m i = 1 10 i and start with (top) p (1) Notice how a lower initial density for species 3 causes it to go extinct. In both simulations species 2 was not able to persist in the landscape, which could mean that their threshold go to their positive equilibrium needed to be even higher ◂ the order at which these events happen does not affect the qualitative dynamics of the models. This is because one can think that which event happens first will be determined by when will the dynamics be observed. If one observes the metacommunity before the colonization process occurs, one would get the model seen in this work, whereas if one observes the metacommunity after the colonization process, a model that resembles more that of Marculis and Hastings (2021) would show up. In both cases, the observed dynamics would be the same, just measured at different moments in time.
This analysis shows that by using a discrete-time framework for a patch occupancy model, we are able to recover a variety of dynamics that are not possible in a continuoustime framework using the simple models based on the Levins metapopulation model. The generality of this framework allows us to study other types of metacommunity interactions, such as spatial heterogeneity over a landscape in systems with alternative stable states  or metacommunities with multiple trophic levels (Guzman et al. 2019). The simplicity of this model allows a deeper analysis of these dynamics, which will then be useful to unveil further properties of these metacommunities.

Proof of theorem 1
First we will prove the following Lemma. If we consider the colonization function to satisfy Equation 3, for each i, up to two nonnegative fixed points p(i) show up, p(i) = 0 and, if p(j) are known for j < i , let (k) , then: To check for stability, System 2 has a triangular Jacobian matrix, which implies that the eigenvalues of the system are given by: If p(i) = 0 , then our eigenvalue simplifies to which is less than 1 provided (1 − e)d i < This implies that the regional extinction equilibrium will be unstable as long as the persistence equilibrium exists. In the case of the persistence equilibrium, notice that the expression further simplifies to In this case, since s i > s i−1 , this eigenvalue i is smaller than 1. Following the proof of this Lemma, we can prove the theorem. We will prove this Theorem by induction over the number of species n. Notice that if n = 1 , then the model simplifies to a Beverton-Holt-type system of the form Since this model is analogous to the Beverton-Holt model, if (1 − e)d 1 > 1 , the system has a single globally asymptotically stable equilibrium Otherwise, the system only has one equilibrium in [0, 1] which corresponds to regional extinction, p = 0 . In both cases, the system has a globally asymptotically stable equilibrium. Suppose now that the theorem holds for some number of species n and consider a system with n + 1 species. This implies that the system with the first n species has a globally asymptotically stable equilibrium in the hypercube [0, 1] n .
Based on this observation, now notice that for big t, the metapopulation p (n+1) t will behave similarly to a single population model. Suppose without loss of generality that p (i) t =p (i) for i ≤ n , where p(i) is the stable equilibrium for species i, which corresponds to either the regional extinction or the value given by Equation 8, depending on the values of e, c, d i for i = 1, … , n + 1 . It is enough to show that p (n+1) t →p (n+1) as t → ∞ . This is the case because, given our assumptions, the equation for species n + 1 simplifies to the Beverton-Holt-type equation given by , then the system converges to the single stable equilibrium 0 < � p (n+1) ≤ 1 given by Equation 8. Otherwise, it converges to the regional extinction p(n+1) = 0.

Proof of proposition 1
This can be easily noticed from Equation 8, which gives us an expression for s i : In this expression, when e > 0 we get that the numerator becomes smaller than the denominator for any values of d i , c , which implies that s i < 1 for all i. Even further, in the absence of disturbance, the only possible scenario is that the top competitor will take all of the patches (i.e., p(1) = 1 ). This proves the proposition.

Proof of theorem 2
If our colonization process follows Equation 4, we have up to two nonnegative fixed points p(i) , the regional extinction p(i) = 0 and the persistence equilibrium given by The eigenvalues of the system are given by: For the regional extinction equilibrium p(i) = 0 , this eigenvalue becomes which is less than 1 provided (1 − e)d i < Similar to the Leslie-Gower colonization, this implies that the regional extinction will be unstable as long as the persistence equilibrium exists. In the case of the . (10) persistence equilibrium, the expression simplifies even further, as the expression = 1 . This implies that the eigenvalue is given by For i = 1 , the equilibrium has the form: The eigenvalue associated with this equilibrium is given by This equilibrium will be linearly stable provided that d 1 < exp (1) 1−e . As a function of e, this expression is increasing in [0,1]. This shows that for species 1, disturbance acts as a stabilizing force, by increasing the upper gap under which the population growth does not produce cyclic behavior. Now we focus our attention to the behavior of species 2. For i = 2 , we first find x 1 , which is given by the expression: which gives us an expression for p(2) using Equation 11: This has the eigenvalue Then, this equilibrium will be stable provided that d 2 satisfies Notice that for species 2, the effect of disturbance over its stability depends on the effect of competition of species 1. The right-hand side of the inequality will be decreasing in e if the inequality is satisfied. If d 1 > exp − 1 2 1 1−e , then for any value of c > 0 , this expression cannot be satisfied. Otherwise, if p (1) = 1 + ln((1 − e)d 1 ) d 1 .
c 2 log(d 1 (1 − e)) + 1 + d 1 (1 − e) < 0 then the inequality is satisfied. If this inequality is satisfied, then as disturbance increases, p(2) can be stable at higher values of d 2 . Thus, at high dispersal intensities of the top competitor, increasing disturbance will tend to a less stable equilibrium for the bottom competitor, whereas at low dispersal intensities, for high enough competition values, decreasing disturbance will tend to a more stable equilibrium for the bottom competitor.

Proof of theorem 3
In order to simplify our analysis, we rescale our variables in Equation 5 to have the following functional form for the colonization function: In this case, the Allee effect colonization function is provided by Equation 13; solving for fixed points p(i) ≠ 0 we get that those fixed points are the real roots to the parabola in p where Notice that the terms A, C are both positive. This implies that the real roots of this parabola, if any, share the same sign. In particular, there exist two positive fixed points for species i if B > 0 and B 2 > 4AC . We can explicitly write these fixed points as We can also find an explicit expression for the i-th eigenvalue as c > −d 1 (1 − e) 2 ln(d 1 (1 − e)) + 1 For any i, in the case of the regional extinction p(i) = 0 we have that i = 0 < 1 . Now, let r be the root term in Equation 15. Then, if we let p(i) =p (i) + , then we can write i as follows: The first term can be verified to be smaller than 1/2, whereas the second term is clearly smaller than 2. Therefore, i < 1 . A similar procedure, using the fact that the sign of r changes throughout the calculation, shows that if p(i) =p (i) − , then i > 1 . This shows that the Allee effect behavior is preserved.
The case for n = 1 was previously studied in Chow and Jang (2014). Using those results and following a similar argument to Theorem 1, we prove the theorem.
need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.