Next Article in Journal
H and Passive Fuzzy Control for Non-Linear Descriptor Systems with Time-Varying Delay and Sensor Faults
Previous Article in Journal
On an Anti-Torqued Vector Field on Riemannian Manifolds
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Metapopulation Persistence and Extinction in a Fragmented Random Habitat: A Simulation Study

by
Hashem Althagafi
1,2,* and
Sergei Petrovskii
1,3,*
1
School of Computing and Mathematical Sciences, University of Leicester, Leicester LE1 7RH, UK
2
Department of Mathematics Science, Umm Al-Qura University, Mecca 24381, Saudi Arabia
3
Peoples Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya St, Moscow 117198, Russia
*
Authors to whom correspondence should be addressed.
Mathematics 2021, 9(18), 2202; https://doi.org/10.3390/math9182202
Submission received: 20 June 2021 / Revised: 1 September 2021 / Accepted: 3 September 2021 / Published: 8 September 2021
(This article belongs to the Section Mathematical Biology)

Abstract

:
Habitat fragmentation is recognized as the most serious threat to biodiversity worldwide and has been the focus of intensive research for a few decades. Due to the complexity of the problem, however, there are still many issues that remain poorly understood. In particular, it remains unclear how species extinction or persistence in a fragmented habitat consisting of sites with randomly varying properties can be affected by the strength of inter-site coupling (e.g., due to migration between sites). In this paper, we address this problem by means of numerical simulations using a conceptual single-species spatially-discrete system. We show how an increase in the inter-site coupling changes the population distribution, leading to the formation of persistence domains separated by extinction domains. Having analysed the simulation results, we suggest a simple heuristic criterion that allows one to distinguish between different spatial domains where the species either persists or goes extinct.

1. Introduction

Factors and processes that either enhance the persistence of natural populations or, on the contrary, make them prone to extinction have been a major focus of attention over the last few decades [1]. This interest has been fuelled by the ongoing global environmental change, e.g., the global warming and its impact on populations an communities [2]. Among other specific effects, the global warming is known to alter species ranges and/or lead to habitat fragmentation. Habitat fragmentation can also result from a direct anthropogenic activity such as forest logging, building new roads, etc. In its turn, habitat fragmentation as well as habitat loss more generally often has distinctly negative effect on corresponding populations and may result in species extinctions [3,4]; in fact, it is regarded as the most serious threat to biodiversity worldwide [5,6].
Understanding the properties of population dynamics in a complex or fragmented habitat is therefore a problem of high practical importance and there is indeed a large number empirical and theoretical studies addressing this issue [7,8,9,10,11,12,13,14,15,16]. A commonly accepted theoretical framework to study population dynamics on a fragmented habitat is the concept of metapopulation [4,17,18,19,20,21,22]. The fragmented habitat is considered as an ensemble of separate sites; the subpopulations of a given species inhabiting the sites may be coupled either by the inter-site dispersal or by a common spatially-correlated external factor such as weather fluctuations, the latter being known as the Moran effect [23,24]). Metapopulation models can be either spatially implicit where the state of metapopulation is described by a single ‘global’ variable such as the fraction of occupied sites [4,17,18,19] or spatially explicit where each site is described by its own ‘local’ variables, e.g., the probability of a given patch being inhabited (e.g., see [25]) or by the size of a given subpopulation (e.g., [10,16,26]. In the latter case, especially where the relative location of sites is described explicitly, the metapopulation model becomes similar to lattice models [27,28] and network models [29,30].
There is a considerable number of studies where the problem of metapopulation persistence/extinction is addressed subject to habitat geometry [25,31,32] and/or environmental and demographic stochasticity [33,34]. There is, however, much less studies where the issue of persistence/extinction is considered in connection with the strength of coupling between different sites (albeit it may sometimes be implicitly included into the habitat geometry, as the strength would normally decay with an increase in the inter-site distance). Meanwhile, this issue is important, as there is some evidence that the strength of inter-site coupling may change as a result of the climate change [35,36].
It must be mentioned here that the possibility of extinction depends on the type of density-dependence in the local population growth. In particular, in deterministic models, in a closed domain (i.e., without outwards migration), the population with the logistic growth cannot get extinct as the extinction state is unstable [37,38,39] and hence is, in the dynamical systems sense, a repeller. However, we mention here that natural populations rarely exhibit logistic growth: much more common is the growth rate with the Allee effect [40,41,42,43], also see [44] and references there. The presence of a strong Allee effect makes the population dynamics significantly different [40,42,43,45,46,47,48]. In particular, the extinction state becomes stable; hence, the extinction of a closed population becomes possible.
In this paper, we address the issue of metapopulation persistence/extinction subject to a change in the strength of the inter-site coupling by considering the dynamics of a single-species system in a 1D chain of sites connected by isotropic population dispersal. Since natural environment is rarely uniform, we consider sites with random properties, i.e., with the carrying capacity being a random variable drawn from a certain distribution. In order to make the dynamics nontrivial and also more biologically realistic (see the previous paragraph), we assume that the local growth is impeded by a strong Allee effect. Our approach is to a large extent inspired by an earlier study, see [49], where it was shown for a conceptual two-site system that an increase in the coupling strength may lead, under some constraints, to a population outbreak when the system ‘jumps’ from the low-density steady state to the high-density one. In mathematical terms, such a transition follows a saddle-node bifurcation when the low-density steady state disappears as a result of the increase in coupling. Although in this paper we are more interested in the extinction rather than in an outbreak, our working hypothesis is that the extinction may follow a sufficiently large increase in the coupling strength due to essentially the same mechanism as in [49].

2. Two-Patch System

We begin with a generic system consisting of two coupled sites or patches, say, Patch 1 and Patch 2. We consider an unstructured population (thus neglecting its age or stage structure) which is fully described by the population size in each patch, i.e., u 1 ( t ) and u 2 ( t ) , respectively, at time t. We consider a species with overlapping generations, so u 1 ( t ) and u 2 ( t ) are continuous functions of time. We neglect demographic stochasticity and hence consider a fully deterministic model where the local population growth is described by an ordinary differential equation [37,38,39]. We assume that the individuals moves around and hence can migrate between the patches. We further assume that the individual movement is random, i.e., diffusion-like. The dynamics of the two-patch system is then described by the following equations:
d u 1 d t = f 1 ( u 1 ) + q 12 u 2 u 1 , d u 2 d t = f 2 ( u 2 ) + q 21 u 1 u 2 .
Here f 1 , 2 describe the local population growth. Since both subpopulations are of the same species, functions f 1 ( u 1 ) and f 2 ( u 2 ) are of the same type, e.g., logistic of with the Allee effect, but the parameters (such the per capita growth rate, the carrying capacity, and the Allee threshold where relevant) can be different if the patches are of different quality. The linear terms in the right-hand side of Equation (1) describe the diffusive coupling where coefficients q 12 and q 21 quantify the strength of the inter-patch coupling, e.g., see [26,50]. In the discussion below, we further assume that the individual movement is isotropic and does not have any preferred direction; hence, q 12 = q 21 = q .
The properties of system (1) are determined by its steady states; in particular, a long-term persistence of the two subpopulations is only possible if there exists a stable ‘coexistence’ steady state, i.e., a positive solution of the following system:
f 1 ( u 1 ) + q u 2 u 1 = 0 , f 2 ( u 2 ) + q u 1 u 2 = 0 .
Technically, the existence of steady states can be established from considering the isoclines of the system, which are described by the following equations:
u 1 = u 2 + 1 q f 1 ( u 1 ) , u 2 = u 1 + 1 q f 2 ( u 2 ) ,
Equations (3) are readily obtained from Equations (2). Obviously, the graph of each isocline in (3) is a deviation from the bisector line u 2 = u 1 , the deviation being given by the corresponding growth rate function (i.e., f 1 or f 2 ) scaled by the coupling strength q.
Any intersection of the two curves defined by Equations (3) is a steady state of the system. The existence of a positive steady state depends both on the type of the growth rate and on the coupling strength. Let us first have a brief look at the simpler case of logistic growth; see Figure 1. It is readily seen that, since the isoclines are convex curves, for any value of q the two curves will always have an intersection point in the interior of the first quarter of the phase plane ( u 1 , u 2 ) . Therefore, for any biologically meaningful parameter values, there exists a unique positive steady state (shown by the black dot in Figure 1), its stability can be shown straightforwardly by considering the direction of phase flow. Correspondingly, starting from any nontrivial initial condition, the subpopulation sizes converge to their steady state values.
The properties of system (1) become different in the presence of a strong Allee effect [40,41,42,44,45,48], i.e., when there exists an intermediate steady state—say, β (often called the Allee threshold)—so that the growth rate becomes negative if the population density falls below the threshold, f ( u ) < 0 for u < β . A commonly used parameterization of the growth rate is a cubic polynomial, so that the generic two-patch model (1) takes the following specific form:
d u 1 d t = α 1 u 1 u 1 β 1 1 u 1 k 1 + q u 2 u 1 , d u 2 d t = α 2 u 2 u 2 β 2 1 u 2 k 2 + q u 1 u 2 .
In this case, since the isoclines have sigmoidal shape, an increase in inter-patch coupling can result in a disappearance of the positive steady state; see Figure 2. While for a sufficiently small q the two isoclines (3) will always have a positive intersection point (Figure 2a), for a sufficiently large q the intersection (and hence the positive steady state) may disappear through a fold bifurcation (Figure 2b). Therefore, the temporal dynamics of the system will be qualitatively different for weak and strong coupling, with the subpopulation sizes converging to a positive steady state for small q but decaying to zero, thus signifying extinction, for a large q.
A crucial issue is to identify the properties of the growth functions f 1 ( u 1 ) and f 2 ( u 2 ) that may lead to the extinction of the two-patch metapopulation. The precise answer to this depends not only on the system’s parameters but also on the initial conditions. Let us consider the case where at the beginning one patch is inhabited and the other one is empty, so that for one subpopulation its initial size is at the corresponding carrying capacity, say u 1 ( 0 ) = k 1 , and for the other it is at zero, u 2 ( 0 ) = 0 . For an uncoupled system ( q = 0 ), since the initial values correspond to stables steady states, there will be no dynamics, so that the subpopulation sizes remain at their initial values forever. For a small positive q, the steady state values are changing as follows: for the second subpopulation, it becomes positive but remains small, 0 < u ¯ 2 β 1 , and for the first subpopulation, it decreases only slightly, u ¯ 1 k 1 , u ¯ 1 < k 1 .
A further increase in q tends to homogenize the system: indeed, it is readily seen from Equations (2) that in the limit q , ( u 1 u 2 ) 0 . Therefore, with an increase in q, u ¯ 1 is decreasing and u ¯ 2 is increasing. However, we further notice that, for any q, the steady state values, if they exist, must satisfy the following equation:
f 1 ( u ¯ 1 ) + f 2 ( u ¯ 2 ) = 0 ,
(which is obtained by adding together Equations (2)). From (5), we readily conclude that u ¯ 1 , although being a decreasing function of q, cannot approach zero if the following condition holds:
max f 1 ( u ¯ 1 ) + min f 2 ( u ¯ 2 ) > 0 or max f 1 ( u ¯ 1 ) > | min f 2 ( u ¯ 2 ) | .
As long as the properties of the growth rate are such that condition (6) holds, the two-patch metapopulation cannot get extinct.
Condition (6) is general, i.e., does not depend on the specific form of functions f 1 , 2 . In case of a given parameterization, the terms in (6) can be further specified. In particular, in case the growth rate is described by a cubic polynomial, condition (6) can, in principle, be expressed via the roots of the polynomials, i.e., via β 1 , 2 and k 1 , 2 . For a general relation between the roots, the condition is given by a bulky algebraic formula which is not instructive; we do not show it here for the sake of brevity. However, in a special case of equal carrying capacities, k 1 = k 2 = k , condition (6) takes a simple, elegant form: the first subpopulation cannot get extinct as along as
k β 1 > β 2 .
By introducing the average value of Allee threshold as β = 0.5 ( β 1 + β 2 ) , relation (7) can be written in a more generic form as
k 2 β > 0 .
Note that relation (8) does not depend on q. Therefore, it gives a sufficient condition of the two-patch metapopulation persistence, but not a necessary one.

3. Extinction and Persistence in a Multi-Patch System

Now, we extend the baseline two-patch model onto a more general multi-patch case. The generalization is straightforward: for an N-patch system, we consider the following system of N differential equations:
d u i d t = f i u i + j = 1 N q i j u j u i , i = 2 , , N 1 ,
where, similarly to the above, u i is the size of the ith subpopulation (i.e., the population inside the ith patch), Obviously, system (9) describes a 1D string of patches where the subpopulation in each patch is coupled (because of diffusive movement of its individuals) to its immediate neighbours. The equations at the end of the chain, i.e., for i = 1 and i = N can be slightly different depending on the boundary conditions. Here we assume that the system as a whole is closed and there are no flows outwards. Therefore, the two boundary equations—the two equations to complement system (9)—are
d u 1 d t = α 1 u 1 ( u 1 β 1 ) ( 1 u 1 / k 1 ) + q ( u 2 u 1 ) ,
d u N d t = α N u N ( u N β N ) ( 1 u N / k 100 ) + q ( u N 1 u N ) .
With regard to the growth rates f i , for the reasons discussed above we focus on the case of the strong Allee effect. Correspondingly, the local population growth is parameterized as follows:
f i ( u i ) = α i u i u i β i 1 u i k i ,
where, generally speaking, the parameters are space-dependent, i.e., different in different patches. Correspondingly, the system as a whole is described by the following three arrays:
{ α i , β i , k i , i = 1 , , N } ,
We consider patches to be random, in the sense that the parameters of the growth rate are random functions of space. However, having all three arrays as actually random is technically challenging as the system will have too many parameters, their value being largely hypothetical. Moreover, it may be excessive as different factors may have a similar effect on the population dynamics, so that, for instance, the dynamics depends on their linear combination. e.g., see relation (7). Correspondingly, below We consider only { k i } as a random array and assume that the per capita growth rate and the Allee threshold are fixed parameters, α 1 = = α N α and β 1 = = β N β . For { k i } , we consider it to be uniformly distributed in [ k m i n , k m a x ] , so that it can be simulated as follows:
k i = k m i n + ( k m a x k m i n ) ψ , i = 1 , , N ,
where ψ is a random number uniformly distributed in [ 0 , 1 ] and k m i n , k m a x are parameters.

Simulation Results

System of Equations (9)–(11) was solved numerically using Runge-Kutta method (using function ode45 in Matlab). We considered a system of one hundred patched, i.e., N = 100 . Each patch is quantified by its carrying capacity and the Allee threshold. as is described above. Simulations were run for the strength of coupling in the range 0 q 20 varying q with the step Δ q = 0.1 . We performed simulations for two different types of initial conditions. For the first type, for each new value of q, the initial subpopulation size in all patches is considered to be the same, that is
u i ( 0 ) = u 0 , i = 1 , , N .
For the second type, the initial condition is defined by ‘continuation’. This is particularly relevant when the system properties are investigated subject to a change in coupling. In this case, simulations are first run for a given value of q for a sufficiently long time (up to t = 200 ) to make sure that the system reaches its large-time asymptotic equilibrium state:
lim t u i ( q ; t ) = u ¯ ( q ) i , i = 1 , , N .
For the next value of q q + Δ q , the initial conditions are taken as the ‘final’ value reached in the previous simulation, that is
u i ( q + Δ q , 0 ) = u ¯ ( q ) i , i = 1 , , N .
For instance, starting from q = 0 , the above approach provides the initial conditions for each considered value of q > 0 .
We begin with the first type initial conditions, cf. (15). In order to provide the first insight into the system properties, We first run a number of simulations for different combinations of q and β . We mostly focus on the asymptotical steady-state behaviour of the system and hence show the metapopulation distributions obtained for a large time, hence ignoring the transient behaviour; it was checked in simulations that t = 200 was sufficiently large to ensure the system’s convergence to its steady state distribution.
Figure 3 shows several typical examples of metapopulation distribution over space for different coupling intensities. We readily observe that, for a small q, the subpopulation sizes fluctuate erratically between the neighbouring sites, as depicted in Figure 3 (i,iv). However, an increase in the coupling strength apparently brings a correlation between neighbouring sites, resulting in the formation of clusters: while the subpopulations can disappear in some subdomains, it will persist in others, as depicted in Figure 3 (ii,v). With a further increase in q, the metapopulation may go to extinction over the whole domain, see Figure 3 (v,vi).
The observation that an increase in q may lead to metapopulation collapse deserves a closer look. Figure 4 shows the results obtained for q = 20 and two other values of β . We observe that an increase in the coupling strength does not necessarily lead to the global metapulation extinction but can lead the formation of persistence clusters instead, provided β is sufficiently large.
We now recall that the carrying capacities are random. Since k is an array made of random numbers, it can be regarded as a realization of a certain random process. We then ask a question: whether a different realization of that hypothetical process, i.e., k where its elements are different random numbers (albeit generated with the same probability distribution (14)), may change the system properties? In order to make an insight into this issue, we generate two more random arrays, say k ( m ) where m = 2 , 3 , and perform simulations respectively. Results are shown in Figure 5. We readily observe that, for those two realizations of the carrying capacity, there is no metapopulation extinction. From the results of Figure 3, Figure 4 and Figure 5, we therefore provisionally conclude that a strong coupling is more likely to lead to a pattern formation (persistence clusters) rather than to a global metapopulation collapse.
While the spatial metapopulation distribution shown at different moments of time and/or for different parameter values provide an exhaustive information about the system state, they make it difficult to follow a change in the system properties resulting from a change in the controlling parameter—in our case, q. We therefore complement the results shown in Figure 4 with another figure (see Figure 6) where the state of the metapopulation is quantified by a single global variable, i.e., the average subpopulation size < u i > . Figure 6 (obtained for the same parameters as in Figure 4) shows how the average subpopulation size gradually decreases along with an increase in q, eventually going to extinction.
A monotonous decrease in the average density is, however, only observed for a smaller value of the Allee threshold ( β = 6 in Figure 6). For a larger β , a change in < u i > is not necessarily monotonous, and in fact not necessarily lead to species extinction. This is illustrated in Figure 7 where the dependence of average subpopulation size on the coupling strength q is shown for a few different values of the Allee threshold ( β = 7 , 8 , 9 ). We readily observe that the dependence < u i > ( q ) reveals a variety of scenarios. For β = 7 (red line), it shows a strongly non-monotonous behaviour first decreasing steadily to a very small value but increasing back to a ‘safe’ (not small) value with a further increase in q. For the critical threshold β = 8 (black line), < u i > ( q ) shows only relatively small fluctuations around a certain intermediate value of < u i > 3 . For β = 9 (blue line), < u i > ( q ) shows a tendency to increase slightly before stabilizing at a higher value < u i > 5 .
In Figure 8, we now focus on the second scenario whereby the subpopulations go extinct only in certain subdomains (but remain persistent in other subdomains) as we vary q and β .
A important point is to understand why in some subdomains the subpopulations are at high risk of extinction while in others they remain persistent. Being inspired by the results obtained in the baseline case of two-patch metapopulation, We hypothesize that there may be a certain relation between the parameters such as the Allee threshold β , the coupling strength q and the average value of the carrying capacity k. In order to verify this, we first compute the average carrying capacity < k ( m ) > in the domains where the subpopulations persist or, on the contrary, go extinct.
Table 1 shows the average subpopulation size < k > focusing on the case where subpopulations go extinct in some subdomains (numbered for convenience). Additionally, Table 2 presents the average subpopulation size < k > for the subdomains where local populations are persistent.
We now recall that, for the baseline case of two-patch system, we obtained a condition that makes impossible the metapopulation persistence; see (8). A question arises as to whether relation (8) may apply to a more complicated multi-patch system, perhaps in a somewhat modified form. In order to check this hypothesis, we consider condition (8) in the following slightly generalized form:
p 2 β k < p * ,
where
k = 1 N N k i , β = 1 N N β i ,
Here p * is a certain ‘critical’ value of quantity p. Note that parameter q is absent from (18). Similarly to the two-patch metapopulation, it indicates that condition is sufficient but not necessary and persistence may become impossible before condition (18) is met (which is indeed the case, see Figure 9 below).
The last column in Table 1 provides the values of p calculated for the spatial domains (shown in Figure 8) where the subpopulations go extinct. Similarly, the last column in Table 2 provides the values of p calculated for some of the persistence clusters. These results are brought together in Figure 9. We readily observe that, indeed, there is a critical value of p: a lower bound of the persistence interval. The critical value can be estimated from the simulation results as p * 8.4 ; it remains unclear how it can be estimated analytically.
Now we are going to make a more systematic look at the effect of coupling strength on the emerging spatial structure of the metapopulation. For this purpose, we run simulations for q being varied in the range 0 q 20 with a small step Δ q = 0.1 and combine the results into a single ‘blue-yellow’ plot; see Figure 10, Figure 11 and Figure 12. For all values of q, the initial condition is used as u i ( 0 ) = c o n s t = 5 . We observe that an increase in coupling tends to decrease the total metapopulation size (cf. Figure 10 middle, Figure 11 middle and Figure 12 right). However, this is not always the case (see Figure 11 left and Figure 12 left): domains where subpopulations went extinct for a certain value of q may become persistence clusters for a larger q, e.g., see Figure 11, left. We also notice that an increase in β leads to an increase in the number of the persistence clusters as well as their size. We further notice that the effect of randomness in the carrying capacity, although altering the specific distributions, does not bring any qualitative changes, as the tendency in all three cases (i.e., left, middle and right) is roughly the same in all three cases shown in Figure 10, Figure 11 and Figure 12.
An interesting question is whether the stationary large-time spatial distribution may depend on the initial conditions. In order to check it, we consider a different initial condition: simulations are now run with the second type condition; see Equations (16) and (17). Results are shown in Figure 13, Figure 14 and Figure 15. We readily observe that in this case the dependence on q is strictly monotonous, which makes the emerging distributions quite different from the case of the constant initial distribution. We therefore conclude that the system exhibits bistability: there exist two (at least) stationary asymptotic distributions, so that it depends on the initial conditions which one actually emerges in the course of dynamics. We also notice that the tendency of having more persistence clusters over the spatial domain, and of a larger size, occurs as well in the case of second type initial condition. This seem to suggest its generality.

4. Discussion and Concluding Remarks

Understanding population dynamics in complex and/or fragmented environments is a major challenge in ecology [11,12,15,22,32,46,51]. In particular, conditions resulting in population collapse and species extinction in a fragmented habitat have long been a focus of the metapopulation theory; previous research suggested specific habitat geometry and demographic/environment stochasticity as factors that can, under some conditions, lead to metapopulation collapse [31,32,33,34] This paper intends to contribute to this discussion by revealing an overlooked factor potentially resulting in metapopulation extinction. We investigates the dynamics of a single-species metapopulation in a one-dimensional random environment; our main focus is the possibility of global extinction subject to the strength of inter-habitat coupling (e.g., as a result of intensified dispersal/migration between the patches) and related factors that can make the extinction more likely.
We first considered a baseline two-patch metapopulation. For this model, we have obtained a sufficient condition preventing metapopulation extinction, see Equation (7). We also showed that, in case the local population growth is hampered by the strong Allee effect, an increase in the strength of inter-patch coupling (due to the inter-patch dispersal) may lead to metapopulation collapse and species extinction, see Figure 2.
We then considered a 1D random metapopulation: a string of patches coupled by a short-distance dispersal (i.e., where each patch is coupled to its immediate neighbours) where the carrying capacity of the local population growth is a random function of space. For this system, using numerical simulations, we showed that an increase in coupling may either lead to metapopulation collapse and global species extinction or to the formation of ‘persistence clusters’ (groups of patches where the subpopulations persist) separated by large stretches of empty space where the subpopulations go extinct. Such persistence clusters usually consist of 10–25 individual patches, so that their spatial size corresponds to an intermediate scale: they are much larger than a single patch but considerably less than the whole system. We emphasize that the persistence clusters are completely self-organized, as our model does not include any long-distance correlations. We also suggested a semi-empirical criterion to identify the persistence clusters based on the relation between an average carrying capacity and the Allee threshold (see Equation (7)) and, using the simulation results, showed that it worked rather well, although it remains unclear how the critical value of quantity p can be calculated analytically.
Thus, in this paper we have examined a new mechanism that may lead to, on one hand, metapopulation extinction or, on the other hand, pattern formation through creating persistence clusters. These results allow for multiple ecological interpretations as well as questions. An increase in individual mobility can lead to an increase in the diffusive interplay between the individual habitats (patches). However, to demonstrate this, we used a very simple, conceptual single-species model. In a real ecosystem, due to the complex nature of inter-specific interactions not only there can be multiple mechanisms that can lead to population extinction, there are also numerous conditions that determine which mechanism takes place. A single-species model is typically only applicable on certain timescales [52]. A question therefore arises whether a more complex model, e.g., explicitly including multiple interacting species, would exhibit similar feedbacks. Despite useful insights from previous work [15,16,19], this issue remains controversial. The complexity even in the dynamics of a few-species community can create periodical and chaotic regimes. Coupling different habitats may result in synchronization/desynchronization of the subpopulations’ dynamics. This issue is of high ecological importance and should be a focus of future research.

Author Contributions

S.P. conceived the ideas and designed methodology; H.A. performed numerical simulations (as part of his PhD Programme under the supervision of S.P.); H.A. and S.P. analysed the data; H.A. prepared the first draft of the manuscript; S.P. finalized the manuscript. Both authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowledgments

H.A. expresses his gratitude to the University of Leicester for a comprehensive research support to his work. He also thanks the Saudi Government and Umm Al-Qura University for both financial and academical support. S.P. was supported by the RUDN University Strategic Academic Leadership Program.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ehrlich, P.; Ehrlich, A. Extinction; Ballantine Books: New York, NY, USA, 1981. [Google Scholar]
  2. Keith, D.A.; Akçakaya, H.R.; Thuiller, W.; Midgley, G.F.; Pearson, R.G.; Phillips, S.J.; Regan, H.M.; Araújo, M.B.; Rebelo, T.G. Predicting extinction risks under climate change: Coupling stochastic population models with dynamic bioclimatic habitat models. Biol. Lett. 2008, 4, 560–563. [Google Scholar] [CrossRef] [PubMed]
  3. Brooks, T.M.; Mittermeier, R.A.; Mittermeier, C.G.; Fonseca, G.A.B.D.; Rylands, A.B.; Konstant, W.R.; Flick, P.; Pilgrim, J.; Oldfield, S.; Magin, G.; et al. Habitat loss and extinctions in the hotspots of biodiversity. Conserv. Biol. 2002, 16, 909–923. [Google Scholar] [CrossRef] [Green Version]
  4. Tilman, D.; May, R.; Lehman, C.L.; Nowak, M.A. Habitat destruction and the extinction debt. Nature 1994, 371, 65–66. [Google Scholar] [CrossRef]
  5. Fahrig, L. Relative Effects of Habitat Loss and Fragmentation on Population Extinction. J. Wildl. Manag. 1997, 61, 603–610. [Google Scholar] [CrossRef]
  6. Wilcox, B.A.; Murphy, D.D. Conservation Strategy: The Effects of Fragmentation on Extinction. Am. Nat. 1985, 125, 879–887. [Google Scholar] [CrossRef]
  7. Kimura, M. “Stepping stone” model of population. Ann. Rep. Nat. Inst. Genet. Jpn. 1953, 3, 62–63. [Google Scholar]
  8. Renshaw, E. A survey of stepping-stone models in population-dynamics. Adv. Appl. Probab. 1986, 18, 581–627. [Google Scholar] [CrossRef]
  9. Cox, J.T.; Durrett, R. The stepping stone model: New formulas expose old myths. Ann. Appl. Probab. 2002, 12, 1348–1377. [Google Scholar] [CrossRef]
  10. Jansen, V.A.A.; Lloyd, A.L. Local stability analysis of spatially homogeneous solutions of multi-patch systems. J. Math. Biol. 2000, 41, 232–252. [Google Scholar] [CrossRef]
  11. DeAngelis, D.L.; Zhang, B.; Ni, W.-M.; Wang, Y. Carrying Capacity of a Population Diffusing in a Heterogeneous Environment. Mathematics 2020, 8, 49. [Google Scholar] [CrossRef] [Green Version]
  12. Kareiva, P. Population dynamics in spatially complex environments: Theory and data. Philosophical Transactions of the Royal Society of London. Ser. B Biol. Sci. 1990, 330, 175–190. [Google Scholar]
  13. Nisbet, R.M.; Briggs, C.J.; Gurney, W.S.; Murdoch, W.W.; Stewart-Oaten, A. Two-patch metapopulation dynamics. In Patch Dynamics; Springer: Berlin/Heidelberg, Germany, 1993; pp. 125–135. [Google Scholar]
  14. Blasius, B.; Huppert, A.; Stone, L. Complex dynamics and phase synchronization in spatially extended ecological systems. Nature 1999, 399, 354–359. [Google Scholar] [CrossRef] [PubMed]
  15. Solé, R.V. Chaos, dispersal and extinction in coupled ecosystems. J. Theor. Biol. 1998, 193, 539–541. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. McCann, K.; Hastings, A.; Harrison, S.; Wilson, W. Population outbreaks in a discrete world. Theor. Popul. Biol. 2000, 57, 97–108. [Google Scholar] [CrossRef] [PubMed]
  17. Amarasekare, P. Allee effects in metapopulation dynamics. Am. Nat. 1998, 152, 298–302. [Google Scholar] [CrossRef] [PubMed]
  18. Levins, R. Some demographic and genetic consequences of environmental heterogeneity for biological control. Bull. Entomol. Soc. Am. 1969, 15, 237–240. [Google Scholar] [CrossRef]
  19. Levins, R. Extinction. In Some Mathematical Problems in Biology; Gersten-haber, M., Ed.; American Mathematical Society: Providence, RI, USA, 1970; pp. 77–107. [Google Scholar]
  20. Pires, M.A.; Queiros, S.M.D. Optimal dispersal in ecological dynamics with Allee effect in metapopulations. PLoS ONE 2019, 14, e0218087. [Google Scholar] [CrossRef] [Green Version]
  21. Wang, W.D. Population dispersal and Allee effect. Ric. Mat. 2016, 65, 535–548. [Google Scholar] [CrossRef]
  22. Hanski, I. Metapopulation Ecology; Oxford University Press: Oxford, UK, 1999. [Google Scholar]
  23. Moran, P.A.P. The statistical analysis of the Canadian lynx cycle II: Synchronization and meteorology. Aust. J. Zool. 1953, 1, 291–298. [Google Scholar] [CrossRef]
  24. Royama, T. Analytical Population Dynamics; Chapman & Hall: London, UK, 1992. [Google Scholar]
  25. Hanski, I.; Ovaskainen, O. The metapopulation capacity of a fragmented landscape. Nature 2000, 404, 755–758. [Google Scholar] [CrossRef]
  26. Namba, T.; Umemoto, A.; Minami, E. The effects of habitat fragmentation on persistence of source–sink metapopulations in systems with predators and prey or apparent competitors. Theor. Popul. Biol. 1999, 56, 123–137. [Google Scholar] [CrossRef]
  27. Allen, J.C.; Schaffer, W.M.; Rosko, D. Chaos reduces species extinction by amplifying local population noise. Nature 1993, 364, 229–232. [Google Scholar] [CrossRef] [PubMed]
  28. Roughgarten, J.; Iwasa, Y. Dynamics of a metapopulation with space-limited subpopulations. Theor. Pop. Biol. 1986, 29, 235–261. [Google Scholar] [CrossRef]
  29. Bodin, O.; Saura, S. Ranking individual habitat patches as connectivity providers: Integrating network analysis and patch removal experiments. Ecol. Mod. 2010, 221, 2393–2405. [Google Scholar] [CrossRef]
  30. Urban, D.; Keitt, T. Landscape Connectivity: A Graph-Theoretic Perspective. Ecology 2001, 82, 1205–1218. [Google Scholar] [CrossRef]
  31. Kininmonth, S.; Drechsler, M.; Johst, K.; Possingham, H.P. Metapopulation mean life time within complex networks. Mar. Ecol. Prog. Ser. 2010, 417, 139–149. [Google Scholar] [CrossRef]
  32. With, K.A.; King, A.W. Extinction thresholds for species in fractal landscapes. Conserv. Biol. 1999, 13, 314–326. [Google Scholar] [CrossRef]
  33. Harrison, S.; Quinn, J.F. Correlated Environments and the Persistence of Metapopulations. Oikos 1989, 56, 293–298. [Google Scholar] [CrossRef]
  34. Legendre, S.; Schoener, T.W.; Clobert, J.; Spiller, D.A. How Is Extinction Risk Related to Population-Size Variability over Time? A Family of Models for Species with Repeated Extinction and Immigration. Am. Nat. 2008, 172, 282–298. [Google Scholar] [CrossRef]
  35. Croteau, E.K. Causes and Consequences of Dispersal in Plants and Animals. Nat. Educ. Knowl. 2010, 3, 12. [Google Scholar]
  36. Travis, J.M.J.; Delgado, M.; Bocedi, G.; Baguette, M.; Bartoń, K.; Bonte, D.; Boulangeat, I.; Hodgson, J.A.; Kubisch, A.; Penteriani, V.; et al. Dispersal and species’ responses to climate change. Oikos 2013, 122, 1532–1540. [Google Scholar] [CrossRef] [Green Version]
  37. Edelstein-Keshet, L. Mathematical Models in Biology; McGraw-Hill: New York, NY, USA, 1988. [Google Scholar]
  38. Murray, J.D. Mathematical Biology; Springer: Berlin/Heidelberg, Germany, 1989. [Google Scholar]
  39. Kot, M. Elements of Mathematical Ecology; Cambridge University Press: Cambridge, UK, 2001. [Google Scholar]
  40. Dennis, B. Allee effects: Population growth, critical density, and the chance of extinction. Nat. Resour. Model. 1989, 3, 481–538. [Google Scholar] [CrossRef]
  41. Stephens, P.A.; Sutherl, W.J. Consequences of the Allee effect for behaviour, ecology and conservation. Trends Ecol. Evol. 1999, 14, 401–405. [Google Scholar] [CrossRef]
  42. Lidicker, W.Z., Jr. The Allee effect: Its history and future importance. Open Ecol. J. 2010, 3, 1. [Google Scholar] [CrossRef]
  43. Berec, L. Allee effects under climate change. Oikos 2019, 128, 972–983. [Google Scholar] [CrossRef]
  44. Courchamp, F.; Berek, L.; Gascoigne, J. Allee Effects in Ecology and Conservation; Oxford University Press: Oxford, MA, USA, 2008. [Google Scholar]
  45. Lewis, M.A.; Kareiva, P. Allee dynamics and the spread of invading organisms. Theor. Popul. Biol. 1993, 43, 141–158. [Google Scholar] [CrossRef]
  46. Keitt, T.H.; Lewis, M.A.; Holt, R.D. Allee effects, invasion pinning, and species’ borders. Am. Nat. 2001, 157, 203–216. [Google Scholar] [CrossRef]
  47. Boukal, D.S.; Berec, L. Single-species models of the Allee effect: Extinction boundaries, sex ratios and mate encounters. J. Theor. Biol. 2002, 218, 375–394. [Google Scholar] [CrossRef]
  48. Sun, G.Q. Mathematical modeling of population dynamics with Allee effect. Nonlinear Dyn. 2016, 85, 1–12. [Google Scholar] [CrossRef]
  49. Petrovskii, S.; Li, B.L. Increased coupling between subpopulations in a spatially structured environment can lead to population outbreaks. J. Theor. Biol. 2001, 212, 549–562. [Google Scholar] [CrossRef]
  50. Amarasekare, P. Interactions between local dynamics and dispersal: Insights from single species models. Theor. Popul. Biol. 1998, 53, 44–59. [Google Scholar] [CrossRef]
  51. Seno, H. Effect of a singular patch on population persistence in a multi-patch system. Ecol. Model. 1988, 43, 271–286. [Google Scholar] [CrossRef]
  52. Ludwig, D.; Jones, D.D.; Holling, C.S. Qualitative analysis of insect outbreak systems: The spruce budworm and forest. J. Anim. Ecol. 1978, 47, 315–332. [Google Scholar] [CrossRef] [Green Version]
Figure 1. The phase plane of system (1) in case the growth rate is logistic, parameterized here as a square polynomial f i ( u i ) = r i u i ( 1 u i / k i ) , i = 1 , 2 . Red and blue curves show the isoclines (3) obtained for hypothetical parameter values r 1 = r 2 = 1 , k 1 = 0.5 , k 2 = 2 and q = 0.5 . The black dot shows the positive steady state, the black line is the bisector line.
Figure 1. The phase plane of system (1) in case the growth rate is logistic, parameterized here as a square polynomial f i ( u i ) = r i u i ( 1 u i / k i ) , i = 1 , 2 . Red and blue curves show the isoclines (3) obtained for hypothetical parameter values r 1 = r 2 = 1 , k 1 = 0.5 , k 2 = 2 and q = 0.5 . The black dot shows the positive steady state, the black line is the bisector line.
Mathematics 09 02202 g001
Figure 2. (a,b) The phase plane of system (4) shown for different values of the coupling strength: (a) q = 0.5 and (b) q = 1 . The black dot shows the positive steady state, the black line is the bisector line. The positive steady state exists in (a) but disappears in (b). Other parameters are k 1 = 1 , k 2 = 3 , α 1 , 2 = 1 , β 1 = 0.2 and β 2 = 1.2 . (c,d) The temporal dynamics of the system obtained for cases (a,b), respectively, for the initial conditions u 1 ( 0 ) = 1.1 and u 2 ( 0 ) = 3.1 .
Figure 2. (a,b) The phase plane of system (4) shown for different values of the coupling strength: (a) q = 0.5 and (b) q = 1 . The black dot shows the positive steady state, the black line is the bisector line. The positive steady state exists in (a) but disappears in (b). Other parameters are k 1 = 1 , k 2 = 3 , α 1 , 2 = 1 , β 1 = 0.2 and β 2 = 1.2 . (c,d) The temporal dynamics of the system obtained for cases (a,b), respectively, for the initial conditions u 1 ( 0 ) = 1.1 and u 2 ( 0 ) = 3.1 .
Mathematics 09 02202 g002
Figure 3. The subpopulation size (the black circles) vs space, i.e., at different patches, obtained for the initial condition u i ( 0 ) = 5 (the dashed green line), for different coupling ( q = 0.1 in (i,iv), q = 2 in (ii,v) and q = 20 in (iii,vi)) and different Allee threshold ( β = 7 in (iiii) and β = 6 in (ivvi), as shown by the orange line). The thin blue broken line show the random carrying capacity.
Figure 3. The subpopulation size (the black circles) vs space, i.e., at different patches, obtained for the initial condition u i ( 0 ) = 5 (the dashed green line), for different coupling ( q = 0.1 in (i,iv), q = 2 in (ii,v) and q = 20 in (iii,vi)) and different Allee threshold ( β = 7 in (iiii) and β = 6 in (ivvi), as shown by the orange line). The thin blue broken line show the random carrying capacity.
Mathematics 09 02202 g003
Figure 4. Variation of subpopulation size over space with a large coupling q = 20 obtained for β = 8 in (i) and β = 9 in (ii).
Figure 4. Variation of subpopulation size over space with a large coupling q = 20 obtained for β = 8 in (i) and β = 9 in (ii).
Mathematics 09 02202 g004
Figure 5. Variation of subpopulation size over space in case of a large coupling q = 20 and β = 7 obtained for two other realizations of the random carrying capacity (generated according to (14)), k ( 2 ) and k ( 3 ) (shown by the thin blue broken line) for (i,ii), respectively.
Figure 5. Variation of subpopulation size over space in case of a large coupling q = 20 and β = 7 obtained for two other realizations of the random carrying capacity (generated according to (14)), k ( 2 ) and k ( 3 ) (shown by the thin blue broken line) for (i,ii), respectively.
Mathematics 09 02202 g005
Figure 6. Dependence of the average subpopulation size < u i > on the coupling strength, q, for different random carrying capacities, k ( 1 ) , k ( 2 ) , k ( 3 ) , and constant β = 6 .
Figure 6. Dependence of the average subpopulation size < u i > on the coupling strength, q, for different random carrying capacities, k ( 1 ) , k ( 2 ) , k ( 3 ) , and constant β = 6 .
Mathematics 09 02202 g006
Figure 7. Dependence of the average population size < u i > on the coupling strength q obtained for different values of the Allee threshold as shown by different colours: β = 7 , β = 8 and β = 9 for red, black and blue, respectively.
Figure 7. Dependence of the average population size < u i > on the coupling strength q obtained for different values of the Allee threshold as shown by different colours: β = 7 , β = 8 and β = 9 for red, black and blue, respectively.
Mathematics 09 02202 g007
Figure 8. Subpopulation size vs space for different coupling strength q, different values of β and different random carrying capacities.
Figure 8. Subpopulation size vs space for different coupling strength q, different values of β and different random carrying capacities.
Mathematics 09 02202 g008
Figure 9. Values of p = 2 β < k > calculated according to (18): red arrows mark the values of p i corresponding to domains with subpopulation extinction (see Table 1), blue arrow mark the values of p i corresponding to persistence clusters (see Table 2); the position of the arrows is approximate. Numbers correspond to case numbering in the tables. The corresponding critical value (shown by the double red-blue arrow) is estimated as p * 8.4 .
Figure 9. Values of p = 2 β < k > calculated according to (18): red arrows mark the values of p i corresponding to domains with subpopulation extinction (see Table 1), blue arrow mark the values of p i corresponding to persistence clusters (see Table 2); the position of the arrows is approximate. Numbers correspond to case numbering in the tables. The corresponding critical value (shown by the double red-blue arrow) is estimated as p * 8.4 .
Mathematics 09 02202 g009
Figure 10. Stationary metapopulation distributions over space as a function of q obtained for the first type initial conditions (15), β = 6 , and three random realizations of the carrying capacity (left to right): k ( 1 ) , k ( 2 ) and k ( 3 ) , respectively.
Figure 10. Stationary metapopulation distributions over space as a function of q obtained for the first type initial conditions (15), β = 6 , and three random realizations of the carrying capacity (left to right): k ( 1 ) , k ( 2 ) and k ( 3 ) , respectively.
Mathematics 09 02202 g010
Figure 11. Same as in Figure 10 but for β = 7 .
Figure 11. Same as in Figure 10 but for β = 7 .
Mathematics 09 02202 g011
Figure 12. Same as in Figure 10 but for β i = 8 .
Figure 12. Same as in Figure 10 but for β i = 8 .
Mathematics 09 02202 g012
Figure 13. Stationary metapopulation distributions over space as a function of q obtained for the second type initial conditions (16) and (17), β = 6 , and three random realizations of the carrying capacity (left to right): k ( 1 ) , k ( 2 ) and k ( 3 ) , respectively.
Figure 13. Stationary metapopulation distributions over space as a function of q obtained for the second type initial conditions (16) and (17), β = 6 , and three random realizations of the carrying capacity (left to right): k ( 1 ) , k ( 2 ) and k ( 3 ) , respectively.
Mathematics 09 02202 g013
Figure 14. Same as in Figure 13 but for β = 7 .
Figure 14. Same as in Figure 13 but for β = 7 .
Mathematics 09 02202 g014
Figure 15. Same as in Figure 13 but for β = 8 .
Figure 15. Same as in Figure 13 but for β = 8 .
Mathematics 09 02202 g015
Table 1. Average carrying capacity < k > and quantity p (see relation (18)) calculated for the extinction domains shown in Figure 8.
Table 1. Average carrying capacity < k > and quantity p (see relation (18)) calculated for the extinction domains shown in Figure 8.
Case No.Domain β qWhich k ( m ) <k>p
1i a β = 6 q = 1 k ( 1 ) 6.34355.6565
2i b6.40225.5978
3 i i a β = 7 5.68418.3159
4 i i b6.20107.7990
5 i i c6.84557.1545
6 i i i a β = 8 6.34639.6537
7 i i i b7.64848.3516
8 i v a β = 9 6.346311.6537
9 i v b7.648410.3516
10v a β = 7 q = 20 6.00957.9905
11 v i a β = 8 6.53319.4669
12 v i i a β = 9 6.459011.5410
13 v i i i a β = 6 q = 1 k ( 2 ) 6.42365.5764
14 v i i i b5.86916.1309
15 v i i i c5.63496.3651
16 i x a β = 7 6.88907.1110
17 i x b6.17727.8228
Table 2. Average carrying capacity < k > and quantity p (see relation (18)) calculated for the persistence clusters shown in Figure 8.
Table 2. Average carrying capacity < k > and quantity p (see relation (18)) calculated for the persistence clusters shown in Figure 8.
Case No.Domain β qWhich k ( m ) <k>p
18v x 1 β = 7 q = 20 k ( 1 ) 5.57098.4291
19 v i x 1 β = 8 5.269710.7303
20 v i x 2 5.298310.7017
21 v i i x 1 β = 9 5.362112.6379
22 v i i x 2 5.289512.7105
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Althagafi, H.; Petrovskii, S. Metapopulation Persistence and Extinction in a Fragmented Random Habitat: A Simulation Study. Mathematics 2021, 9, 2202. https://doi.org/10.3390/math9182202

AMA Style

Althagafi H, Petrovskii S. Metapopulation Persistence and Extinction in a Fragmented Random Habitat: A Simulation Study. Mathematics. 2021; 9(18):2202. https://doi.org/10.3390/math9182202

Chicago/Turabian Style

Althagafi, Hashem, and Sergei Petrovskii. 2021. "Metapopulation Persistence and Extinction in a Fragmented Random Habitat: A Simulation Study" Mathematics 9, no. 18: 2202. https://doi.org/10.3390/math9182202

APA Style

Althagafi, H., & Petrovskii, S. (2021). Metapopulation Persistence and Extinction in a Fragmented Random Habitat: A Simulation Study. Mathematics, 9(18), 2202. https://doi.org/10.3390/math9182202

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop