Scaling Laws in the Fine-Scale Structure of Range Margins

Abstract

Margins of the geographic distributions of species are important regions in terms of ecological and evolutionary processes, including the species’ response to climate change. This paper reviews some spatially explicit metapopulation models of range margins across environmental gradients (e.g., across latitudes or altitudes). These models share some robust results, which allow for generalizations within a broad variety of species and environments: (1) sharp edges can emerge even across relatively smooth environmental gradients; (2) intraspecific competition combined with dispersal limitation is a sufficient condition for the sharpening; (3) at the margin, the “mainland” of continuous occurrence splits into “islands”. Computer simulations pointed out some characteristic scaling laws in the size distribution of the islands, and in the structure of the hull of the mainland. The hull is a fractal with a dimension 7/4. Its width and length scale with the gradient according to characteristic scaling laws (with exponents 3/7 and 4/7, respectively). These general features follow from a second-order phase transition from a connected to a fragmented state. The results contribute to understanding the origin of vegetation zones and the spatial pattern of ecotones.

1. Introduction

A fundamental goal of ecological research is to understand the distribution of organisms in space and time [1,2]. A lot of effort was devoted to the study of populations at the edges of spatial distributions, where the conditions gradually worsen and, thus, the abundance vanishes. The focus on this region is for many reasons. Firstly, the transition from presence to absence can clearly show the species’ ecological requirements in situ (e.g., its tolerance to extremes in the temperature) and, thus, may help predict the species’ response to climate changes [3,4,5,6,7,8,9,10]. Secondly, the range margins are often special in terms of genetic composition. For example, the range margins, compared to the core, may contain a higher proportion of far-dispersing individuals (e.g., with differences in the seeds, as in Abronia umbellata [5]; see also References [9,11]). Thirdly, the manifestation of phenotypic plasticity may also be stronger at the margins [12,13,14]. A typical example of large phenotypic differences between the core and the margin is the occurrence of a krummholz zone at alpine treelines [12,13].

As the population’s density decreases along the environmental gradient, it is frequently observable that an originally continuous spatial distribution gets fragmented into “islands”. These scattered occurrences are mentioned as spots of potentially rapid evolution for two main reasons: (1) local adaptation; and (2) genetic drift in the relatively small local populations (c.f. [5,9,15,16,17,18]). Both depend on the degree of isolation, i.e., on the gene flow between the islands and the mainland [4,5,16,19,20]. To understand these processes, not only should the dispersal or migration of individuals be taken into consideration, but also the merging and splitting of whole islands with the mainland. This necessitates the study of the spatial structure and dynamics of (meta)populations on a broad scale.

Theoretical modeling contributed significantly to understanding the characteristic patterns and processes at range margins. The aim of this paper was to review some insights from strategic models (sensu [21]), which can be generalized for various species and environments. In terms of modeling techniques, these are spatially explicit metapopulation models with environmental gradients. The usage of the word “gradient” is different in ecology vs. mathematics; I apply the ecological term. An “environmental gradient” is mentioned when an ecological factor changes monotonously along a spatial axis through the studied area (see, for example, the fundamental “gradient analysis” papers by Whittaker [22]). In its simplest usage, the “gradient” means the physical object (e.g., the hillside). In a more abstract usage, the “gradient” is the spatial axis along which the environment changes. Some of the studies yielded general scaling laws to characterize the spatial structure. These shall be discussed in the context of critical phenomena. Finally, I discuss the importance of these findings for the protection of wildlife.

2. Metapopulations along Environmental Gradients

The occupancy of space is generally determined by the balance between two competing processes: colonization and (local) extinction. Metapopulation models are well-developed tools for the study of these events [23]. The first metapopulation model was introduced by Levins [24]. The model describes the occupancy of discrete habitat sites within a geographic region. It assumes spatial homogeneity (no gradient). The proportion of occupied sites ($0\le n\le 1$) changes over time according to the following equation:

$$\frac{dn(t)}{dt}=cn(t)\left[1-n(t)\right]-en(t),$$

(1)

where c and e are the rates of colonization and extinction, respectively. The Levins model inspired a lot of further research on the dynamics of colonization and extinction. An important advance was to represent the spatial locations of the sites (e.g., geographic coordinates [23]). This spatially explicit approach enabled the investigation of limited migration or dispersal distance, which is generally a significant factor in the spreading of species [25]. The spatial positions of the sites can be incorporated in various ways [26]. A typical method is to subdivide the total area into a regular lattice. Distance-dependence can be expressed, for example, by stating that an unoccupied cell can only be colonized from the neighboring occupied cells (e.g., from the four adjacent cells in a square lattice). This “contact” nature of colonization is directly represented in the so-called contact process [27] (CP; see a review about its ecological applications in Reference [28]). It is relatively easy to incorporate an environmental gradient into the CP. The only additional feature is that the colonization rate and/or extinction rate changes along a spatial axis (c(x) and/or e(x)). The gradient contact process (gCP) can be considered as a spatially explicit minimal model for the study of range margins [29].

The earliest gradient metapopulation model I am aware of was published by Noble in 1993 [30]. It was a “contact model”, but was more complex than the gCP, as it included jump dispersal and an Allee effect. Many researchers created and analyzed a large variety of gradient metapopulation models since. Some of them directly applied the gCP [29,31], or a model that could be transformed into a gCP [32,33]; others used more complex approaches (e.g., References [6,9,11,16,34,35,36,37,38,39,40]). The present review does not describe the details of these models; instead, it only highlights some general conclusions.

Due to the diversity of the models and of the fields of sciences involved (from ecology through population genetics to evolutionary biology), many terms are used for naming the same, or nearly the same, geometric objects. It would be a matter of a longer review to clarify whether the terms “range margin”, “range edge”, “range boundary”, and “species border” are synonyms, or whether their usage differs for different authors. In this paper, I consistently use “range margin” when I write about a two-dimensional pattern, and “range edge” when I mean a lower-dimensional object (a straight line or a fractal with <2 dimensions). Finally, I present a remark about the metapopulation’s definition. The metapopulation, in general, is a set of local populations connected by the motion of individuals (including the dispersal of offspring). In many studies, the definition implies that the local populations live in discrete habitat patches (e.g., Reference [41]). In other studies, including the present review, the term is used in a broader sense (e.g., References [6,9,33]), whereby the study area is subdivided into discrete sites (either naturally or artificially), and the organism can (a) move or disperse its offspring between these sites, (b) colonize the empty sites, and (c) get locally extinct from the occupied sites. Therefore, the models I present are essentially spatially discrete, site-based models of population dynamics.

3. Abrupt Changes along Smooth Environmental Gradients

The range margin is usually shaped by many ecological factors (temperature, moisture, availability of resources, etc.). Metapopulation models boil down the plentitude of interactions between the organisms and the environment into the consequent demographic rates (c and e in the CP). This permits comparisons between different species in various habitats (c.f. [38]). The homogeneous (non-gradient) CP can be simplified even further. Taking $1/e$ as a time unit, the system has a single control parameter, the spreading rate $=\frac{c}{e}$.

Let us observe $\widehat{n}(\lambda )$, which is the equilibrium density of occupied sites over the total area in the homogeneous CP. Investigations in statistical physics revealed a second-order (critical) phase transition from the living ($\widehat{n}(\lambda )>0$) to the globally extinct ($\widehat{n}(\lambda )=0$) state at a critical threshold ${\lambda }_c$ [42,43,44]. The value of ${\lambda }_c$ in a square lattice with a four-cell neighbourhood was estimated as ≈1.6488 by Monte Carlo simulations. As the metapopulation approaches the threshold ($\lambda \to {\lambda }_c$), the equilibrium density declines rapidly according to a scaling law,

$$\widehat{n}\propto {(\lambda -{\lambda }_c)}^{\eta },$$

(2)

where ∞ denotes “proportional to”. In the same region, the correlation length in space (${\xi }_s$) and in time (${\xi }_t$, called the relaxation time) diverge. The divergence can also be described by scaling laws,

$${\xi }_s\propto {(\lambda -{\lambda }_c)}^{{}^{-{\nu }_s}}\mathrm{and}{\xi }_t\propto {(\lambda -{\lambda }_c)}^{{}^{-{\nu }_t}}.$$

(3)

In finite samples of size A, the variance to mean ratio (V/M) of the density also diverges,

$$\frac{VA}{M}\propto {(\lambda -{\lambda }_c)}^{-\mu }.$$

(4)

The values of the scaling exponents were estimated numerically: $\eta =0.583$, ${\nu }_s=0.733$, $v_t=1.2956$ [43,44], and $\mu =0.934$ (calculated from data in References [27,42,45], see Reference [28]). Importantly, the scaling exponents are “universal” in the sense that their values depend only on the dimensionality of the system (which is two in the present case), and are independent of the geometry of the lattice (square, hexagonal) of the neighborhood size (4, 8, etc. cells), and other local details. The scaling relationships expressed in Equations (2) and (3) apply to a broader set of processes beyond the CP, within the directed percolation universality class [42,43,44] (see a description for ecologists in Reference [46]). Essentially, the behavior of the system is dominated by long-range correlations in the vicinity of the threshold and, thus, the local details become negligible.

In the gCP, the demographic rates depend on the spatial location: c(x) and/or e(x). It is a fundamental question whether $\lambda (x)=\frac{c(x)}{e(x)}$ is still sufficient for the prediction of the equilibrium density $\widehat{n}(\mathrm{x})$. A complicating factor is that the parent site, from which the colonization occurs, can differ from the offspring’s site in quality (i.e., can be at a higher or lower position along the gradient). Therefore, it is necessary to specify whether the success of colonization is parent-dependent or offspring-dependent. Our Monte Carlo simulations [29] showed that, in the case of offspring dependence, $\lambda (x)$ is always an accurate predictor of $\widehat{n}(\mathrm{x})$. That is, various combinations of c(x) and e(x) with the same $\lambda (x)=\frac{c(x)}{e(x)}$ produce the same $\widehat{n}(\mathrm{x})$. In the case of parent dependence, the whole profile of c(x) and e(x) along the gradient has to be taken into consideration. Nevertheless, $\lambda (x)$ provides a good approximation when the parent-dependent component of colonization does not change abruptly across space.

Another considerable effect is the flow of individuals from the more to the less densely populated region. Due to this flow, the species can be present permanently beyond the extinction threshold (at x, where $\lambda (x)<{\lambda }_c$) [29]. Due to this “rescue effect”, the shape of $\widehat{n}(\lambda )$ in the homogeneous CP is fundamentally different from the shape of $\widehat{n}(\mathrm{x})$ in the gradient CP, even if $\lambda (x)$ is linear ($\lambda (x)=gx$).

The “downslope” flow of individuals is not a specificity of the gCP, but is generally observable in gradient metapopulation models (e.g., References [7,9,11,15,18,32,47]). This phenomenon and, thus, the rescue effect becomes negligible only in the limit of infinite system size when $g\to 0$. However, all real-life systems are finite. Most of them are relatively small compared to the systems used in the Monte Carlo simulations. Since the CP describes spreading by spatial contacts, the most reasonable lattice constant is equal to the maximum dispersal distance in the actual species [28]. Therefore, the finite-size effect is probably considerable in most real-life situations. The sharp extinction threshold at ${\lambda }_c$, suggested by the theory, is presumably blurred by this finite-size effect. Nevertheless, the fact that simple colonization–extinction dynamics produce an abrupt change in the density of occupancy even along smooth environmental gradients is an important result for ecologists.

The abrupt change is essentially caused by dispersal limitation. A species which could disperse its seeds (or other kind of offspring) all over the area, without any spatial limitation, would perfectly “mirror” the spatial pattern of λ in the equilibrium. In contrast, when the distance of dispersal is limited, the spatial distribution is clumped to some extent: the local density of occupied sites around each occupied site is typically higher than the mean (global) density over the area. Therefore, density regulation fails. The failure exacerbates as the discrepancy between the local and global density increases. A diagnostic sign is the increase of fluctuations (Equation (4)). Approaching the threshold can be viewed as getting into a “vertex of extinction”: the lower the density, the more inefficient its regulation. This is the primary reason for the abrupt transition from the living to the extinct state (See more details in Reference [46]).

It is important to note that this kind of abrupt change results from the spatial self-organization of the metapopulation. It is entirely different from those cases in which the range ends because of an external factor, e.g., a geographic barrier. The phase transition, presented here, is produced by the populations’ colonization and extinction events, and it depends on the emergent spatial pattern of occupancy.

The mechanism appears to be fairly general. Colonization is ubiquitous; local extinction also occurs in most populations, at least on fine spatial scales. Dispersal limitation occurs in most species, although its significance depends on the spatial scale of the study. Practically, those cases match into this theoretical framework in which the dispersal distance of the species is much smaller than any linear size of the studied terrain. I think it would be an important subject of future research to delineate the set of population dynamic processes which belong to the directed percolation universality class and, thus, share the characteristics of the CP. The behavior of these should also be tested along environmental gradients, in comparison with the gCP.

The main qualitative result of the gCP is the occurrence of abrupt changes along smooth environmental gradients [29]. This is in agreement with many other gradient metapopulation models. The sudden decline of density is a general observation [4,9,16,18,32,33,34,35], and the increase of random fluctuations toward the edge was also observed [33,35,48].

In concert with the models, empirical studies on various plant and animal species also yielded results about the abrupt decline of density across relatively smooth gradients (e.g., References [49,50,51]). Some ecotones provide well-observable examples. Primarily, those cases suit the metapopulation approach in which the margin is formed by a single, dominant species. In many examples, the dominant species is a tree, and the environment at the margin is too harsh for supporting comparably strong competitor species. Examples include high mountain treelines [13], the transition from taiga to tundra vegetation [38], and the ecotone between woodland and grassland [34,52,53]. In those cases in which the ecotone is an interface zone between two or more similarly strong competitors, a more complex model is needed, e.g., a multi-species gCP. An important, common message from many of the existing gradient metapopulation models (e.g., References [32,33,35,38]) is that dispersal limitation in a single species is sufficient for forming a sharp edge, even in the absence of external factors (e.g., the presence of a barrier or any influence from other species).

4. The Spatial Pattern of Occupied Sites

A common observation in almost all gradient metapopulation models is a transition from continuous to scattered occurrence of the species, i.e., a “mainland” and “archipelago” structure [29,32,33,34,35,36,38,39,54,55] (see also References [47,56] for reviews about empirical studies). Fragmentation is a natural consequence of the decrease of the density of occupied sites. I discuss two questions: (1) Where does the transition occur? (2) How does the size distribution of islands change along the gradient?

Percolation theory provides a firm theoretical basis for the study of the transition from the continuous to the scattered states. In its original, simple version, the theory was formulated for homogeneous (non-gradient) situations [57]. A typical site percolation problem can be described as follows: consider a lattice which consists of N sites; Np sites are occupied and N (1 − p) are empty. The empty and occupied sites are distributed randomly within. We study the connectivity of the occupied sites in the following sense: let us define two sites connected if, and only if, they are neighbors and are both occupied. The neighborhood size Z is arbitrary. A common choice in square lattices is to take the four adjacent cells as neighbors (Z = 4). In spite of the simplicity of this model, there are many exciting questions that can be asked about the overall connectivity of the sites. For example, what is the probability of the occurrence of a connected path from one side of the lattice to the opposite side?

Percolation theory, with more advanced versions of the original model, is broadly applied in natural sciences and engineering [57,58]. A typical usage is the study of the flow of material in porous media. It was also applied in ecology, for example, in the study of landscape connectivity (see Reference [28] for a review). I present here some basic results from the theory, in order to discuss the connectivity of occupied sites at a given time in a metapopulation. A simple example for the practical usage of a percolation model in a population is to investigate the opportunities of spreading (e.g., of a pest) between trees in a forest. However, I am going to use the theory without such interpretation, only for analyzing the spatial structure (thus, the spreading agent can simply be considered as a “dye”).

One of the basic characteristics in the lattice is the probability of percolation. Consider an infinitely large lattice in which a fraction p of the sites is occupied. Starting from a randomly chosen occupied site, what is the probability that the spreading agent can get infinitely far from the point of origin? The percolation probability R was found to undergo a second-order phase transition at a threshold p_c. In the limit of an infinitely large lattice, R is 0 at $p<p_c$, and is positive at $p>p_c$. The threshold value depends on the details of the model. In a square lattice with a four-cell neighborhood, numerical simulations estimated p_c = 0.592746 [57]. In the surviving phase ($p>p_c$) close to the threshold, R can be described by the following scaling law:

$$R\propto {(p-p_c)}^{\beta },$$

(5)

where $\beta =5/36$ [57]; β is a universal scaling exponent (i.e., independent of the neighborhood size and other local details), and is characteristic for the isotropic percolation universality class (which is different from the directed percolation class, mentioned at the CP; see a comparison between these classes in Reference [28]).

It is interesting to observe the change of the sizes of the percolation clusters with the increase of p. Each percolation cluster consists of a set of connected sites, surrounded by non-connected ones. The average cluster size S is defined in the following way: we select an occupied site, and record the size of the percolation cluster to which the site belongs. We calculate the average of these cluster sizes; S is finite below the percolation threshold (at $p<p_c$), and diverges according to a scaling law as the threshold is approached.

$$S\propto {(p_c-p)}^{-\gamma },$$

(6)

where $\gamma =43/18$ [57]. Above the percolation threshold ($p>p_c$), the system contains an infinite cluster by definition (often called the giant component); thus, S is infinite. Disregarding the infinite cluster, the remaining finite ones obey a similar scaling law,

$$S\propto {(p-p_c)}^{-\gamma },$$

(7)

with the same universal exponent $\gamma =43/18$, but with a smaller, non-universal multiplying factor.

In finite lattices, even the largest cluster must be finite and, in general, the scaling laws in Equations (5)–(7) are modified by finite-size effects [57]. Nevertheless, percolation theory has some important messages for finite systems as well, including gradient metapopulations.

In the percolation model presented above, the value of p was the same in every lattice cell. A next step toward gradient metapopulations is to assume that p changes linearly along a spatial axis, $p(x)=ax$. This system is a gradient random map (gRM) [31]. In a further step, we can study the spatial pattern of occupied sites in the gCP. Since the pattern is time-dependent, it is reasonable to study snapshots in the (dynamic) equilibrium. We investigated this system assuming a linear change $\lambda (x)=gx$ as well as other $\lambda (x)$ functions [31,59].

According to the results, the gRM and the gCP share some features. There is a characteristic mainland-with-islands structure along the gradient, where the mainland corresponds to the giant cluster, and the islands to the smaller clusters (Figure 1). Figure 1 shows snapshots from simulations in relatively small lattices (N = 6800). The results presented in References [31,59] were obtained in larger systems (N > 16 million). In the simulations, I used a four-cell neighborhood. Colonization was offspring-dependent; thus, $\lambda (x)$ was sufficient for characterizing the environment. The lower panels demonstrate the abrupt decline of $\widehat{n}(\mathrm{x})$ (green curve) in spite of the smooth change of $\lambda (x)$ (orange line).

The upper panels show the spontaneous fragmentation at the range margin, depending on the slope (compare Figure 1a,b). These observations suggest that the occurrence of population fragments (islands) at a margin does not necessarily indicate any heterogeneity in the background (e.g., patches with particularly favorable microclimate or resource-richness [29,32]). Islands may emerge simply due to the internal dynamics. The two cases, externally vs. internally generated islands, can be distinguished by observing the positions of the islands over time. In the externally generated case, the islands are static (anchored to the favorable habitat patches), while, in the internally generated case, they are mobile.

In concert with the percolation model, many other models of gradient metapopulations produced mainland-and-archipelago structures, both within the gCP modeling framework [32,33], and outside [34,35,36,38,39,54,55]. I am not aware of any empirical study in which the external vs. internal reasons for island-formation were separated. I think this would be an important matter for future research, for example, for better understanding gene flow between the islands and the mainland, which hinges on the frequency of merging and splitting (see Section 1).

The sizes of the islands (population fragments) are also important. For example, genetic drift is more likely in smaller fragments, and local adaptation may be more difficult due to a lower local genetic diversity (c.f. [18]). On the other hand, these islands, as “outposts”, may become sources of population expansion in case of a climate change [51]. Percolation theory has some definite predictions about the sizes of the islands (Equations (6) and (7); see more details in Reference [31])). The sizes are predicted to decrease abruptly with the increase of distance from the mainland, depending on the slope. Observations at an alpine treeline in the Glacier National Park in North America [39] are in agreement with this prediction. The cumulative size distribution of the islands was found to fit a scaling law,

$$F(S\ge z)\propto z^{-\phi },$$

(8)

where the left side is the frequency of the patches that are larger than z. The value of the exponent was estimated to be ϕ = 0.83. The generality of this result in other range margins is yet to be explored. Considering the universality of the scaling laws that characterize the sizes of the islands in the gCP (Equations (6) and (7)), it would be interesting to study the cumulative size distribution of the islands in the gCP, and in other simulated and real range margins.

There are some mechanisms which decrease the likelihood of island-formation, for example, a strong Allee effect [9,35,36,38,40,60] or other local positive feedback mechanisms (self-facilitation [61]). In these cases, the gCP can be considered as a reference to which the other cases, with additional effects, can be compared.

It should be noted that the sizes of the gaps, which are left unoccupied by the species within the mainland, are also important in an ecological and evolutionary sense. The gaps were studied in non-gradient cases both in lattice simulations and on real-life data (in a multi-species situation, in a rainforest in Panama [62,63]). It would be interesting to extend these studies to range margins as well. Equations (2) and (3) predict an abrupt loss of gaps from the edge toward the inside of the mainland, and Equation (7) suggests that the number and area of islands within the mainland also decreases rapidly as the distance from the percolation threshold increases. Gaps in the canopy of the dominant species are important for the survival and reproduction of subordinate species, and they may provide special microclimatic conditions which increase biodiversity. Furthermore, they may be important in the evolutionary adaptation of the dominant species to the environment. Turner and Wong [54] applied stochastic simulations to investigate the problem that a mutation that would be advantageous within a certain range along an environmental gradient may arise further away from that range. The persistence of the mutant genotype hinges on its success of reaching the suitable sites, by percolating through the area covered by the resident genotype(s) (see also References [16,64,65]). This strongly depends on the density of the resident along the environmental gradient. The availability of empty gaps may be decisive in the percolation process.

5. The Hull Edge

The previous chapters presented the abrupt nature of the transition from the living to the extinct state along the gradient (as shown in the lower panels in Figure 1). The transition zone is relatively narrow, but has some interesting properties in magnification. Section 4 already presented some features of the islands. This chapter focuses on the edge of the mainland.

One of the frequently studied geometric structures is the hull. It consists of those sites which belong to the mainland, and are adjacent to sites that belong to the giant unoccupied site (the “ocean”); thus, the hull can be considered as the “coast” (marked by red in Figure 1). The neighborhood size (Z) at determining the clusters is arbitrary again, and is not necessarily the same as the neighborhood size used when producing the pattern in the gCP. The concept of the hull is broadly applied in percolation theory [66], beyond the gCP.

Figure 1 illustrates the winding of the hull in the gCP. The winding is constrained in the x direction by the gradient. The gradient imposes an inevitable finite-size effect, as it occurs even when the total area is infinite. At coarse spatial scales, the hull can be considered as a one-dimensional object (a straight line). At fine spatial scales, however, it has an intriguing fractal structure. The fractal dimension was numerically estimated to be D = 7/4 [31].

Coarsening the scale of observation, there is a crossover from D = 7/4 to 1 (Figure 4 in Reference [31]). Within the fractal region, D = 7/4 is independent of the neighborhood size Z, of the slope g (in $\lambda (x)=gx$), and of some local details of the metapopulation dynamics; both a constant c with linear e(x), and a constant 1/e with linear c(x) yield the same fractal dimension. The independence is observable only within reasonable limits because of the aforementioned finite-size effect. Z or g is obviously too high when islands do not occur beyond the tips of the “peninsulas” of the hull, i.e., it can be suspected that the hull is distorted.

The (relative) independence of D on Z suggests that there is a zone (a range within x) in which the fractal properties are detectable. This may explain the scaling law observed by Zeng and Malanson [39] concerning the size distribution of the islands in a broad range of x (Equation (8)). This is an interesting matter for future research. The fractal nature of range margins is certainly an important finding, especially when the organism under study provides habitat for other organisms (for example, trees are inhabited by many organisms). Several papers in landscape ecology emphasized the importance of flows across boundaries (e.g., the diffusion of nutrients and of living organisms [28,67]). Milne et al. [34] discussed the importance of the hull in this regard. To my knowledge, we were the first to estimate the fractal dimension of the hull [31]. Such a high-dimensional fractal, D = 7/4, may provide a considerable surface for the flows.

Winding hull edges were observed, for example, at the margin of woodland patches in New Mexico [31,34]. The field data were fitted to the predicted D = 7/4 in Reference [31]. The fit was remarkably accurate at three different values of Z, which indicates the relevance of the theoretical results under natural conditions, even in relatively small systems.

The width (u) and length (w) of the hull are also important observables. The width was characterized by the standard deviation of the hull along the x axis. The length was defined as the number of lattice cells that belonged to the hull. Both the width and the length depend on g,

$$u(g)\propto g^{-\zeta };$$

(9)

$$w(g)\propto g^{-\omega }.$$

(10)

The scaling exponents are ζ = 3/7 and ω = 4/7. These values, as well as the fractal dimension D can be directly calculated from the scaling exponents β and γ in ordinary, non-gradient percolation (Equations (5) and (6); see pages 31 and 34 in Reference [31]). This is surprising, since the CP is a rather complex pattern-generating mechanism compared to the conventional generation of a percolation map (Section 4). The CP produces a clumped spatial distribution of occupied sites. The result can probably be explained by the abrupt decrease of spatial correlations away from the extinction threshold (Equation (3)). When g is not very large, i.e., the slope is not very steep, then the hull is at a considerable distance from the extinction threshold (see Figure 1b) and, thus, the spatial correlations may become negligible. Essentially, the system becomes similar to a gradient random map. This hypothesis was confirmed by analyzing the hull in the gRM. Exactly the same values of the scaling exponents (D = 7/4, ζ = 3/7, and ω = 4/7) were found. This suggests that the gRM, in spite of its simplicity, may be a suitable model for the study of the hull across a gradient, even in those cases in which the underlying population dynamic process is more complex. It would be an interesting task for future research to check the validity of this statement in gradient metapopulation models that are more complex than the gCP. The fractal dimension of the hull could certainly be estimated in real range margins as well. Obtaining data for u(g) and w(g) is more difficult.

Finally, the mean position of the hull along the x axis is also an important observable (let us denote it by x_m). This can be considered as the main place of the transition from the connected to the fragmented state. Accordingly, Milne et al. [34] raised the hypothesis that the density of occupied sites at this place is $\widehat{n}(x_m)=p_c$. In square lattices, considering only the most adjacent cells as neighbors (Z = 4), this value is approximately p_c ≈ 0.59 (see Section 4 and Reference [57]). The hypothesis was confirmed by gRM simulations, and similar values were also found in the studied New Mexican woodlands [31,34]. It is remarkable that, in forestry, 60% canopy cover is often used as a limit above which the vegetation can be considered a closed forest [34].

Although gradient metapopulations are fairly complex (see the previous sections), the hull appears to be a simple structure. It is hopeful to understand its main features solely on the basis of classical percolation theory. We investigated the validity of this statement by simulating percolation maps with various p(x) functions, in addition to the previously studied linear $p(x)=ax$. The simulations confirmed the broad validity of D = 7/4. The investigations also provided an asymptotic expression for the probability that a site at a given distance from x_m is on the hull [59]. Within the complex pattern of the range margin (Figure 1), x_m is a well-defined position. It could be used for delineating the range edge unequivocally, and monitoring the shift of the range in case of a climate change [31,34]. I think it is an important matter for future research to explore this opportunity in various species, and to work out a monitoring method that could be used in practice.

6. Summary and Outlook

Understanding the limits of species’ distributions is one of the fundamental goals in ecology. The task is rather difficult, because the range edges are often self-organized, i.e., broader geographic patterns emerge from the local processes in a complex manner [4,9,29,31,32,33,34,38,39,47,52,53,56,68,69]. Comparisons between different species and geographic regions necessitate a common theoretical framework, within which this complexity could be handled. Several papers suggest that the theory of phase transitions provides important help in this attempt [31,34,39,52,53,70,71]. It suggests some broadly generalizable laws, which concern the spatial pattern of occupied sites along environmental gradients.

The theory also suggests some direct simplifications, which reduce complexity. One of these is the use of the hull for delineating the range edge. The hull has some broadly applicable scaling properties, which were proven to be independent of the details of the models, as they could be deduced from ordinary (non-gradient) percolation. These theoretical considerations suggest that the comparison of range margins (between species, places, or years) could be simplified by comparing the hulls. For example, a range shift could be detected sensitively by monitoring the mean position of the hull.

Another opportunity for simplification is suggested by the gCP, and concerns the equilibrium density profile along the gradient, $\widehat{n}(\mathrm{x})$. Computer simulations raised the hypothesis that the spreading rate, $\lambda (x)=\frac{c(x)}{e(x)}$, is sufficient for predicting the densities in a broad range of realistic cases; thus, various combinations of c(x) and e(x) with the same $\lambda (x)$ produce the same equilibrium density profile. The possibility of this kind of simplification, by rescaling time, could be investigated in other gradient metapopulation models as well, beyond the gCP.

I presented two phase transitions: (1) from the living to the extinct state, and (2) from the connected to the fragmented one. Each transition was presented firstly in homogeneous (non-gradient) systems. Then, an environmental gradient was introduced, i.e., it was assumed that the control parameter changed along a spatial axis. The gradient resulted in some inevitable finite-size effects which blurred the original threshold-like transition. Nevertheless, the changes remained abrupt even along smooth gradients. This “sharpening” is an important phenomenon, as it can influence the population’s adaptability to the environment, and the capacity of movement in case of climate change.

This review was based on the gradient contact process (gCP), which is a basic model of colonization and extinction across an environmental gradient. Its main findings are in agreement with many other, more complex gradient metapopulation models. It would be interesting to outline the set of models to which the gCP’s universal features are directly applicable, and to understand those factors that call for fundamental changes in the model’s settings, thus making the gCP’s scaling laws inapplicable. For example, some papers suggest that a strong Allee effect can destroy the CP’s phase transition [40] and can fundamentally change the metapopulations spatial structure [38,40,60]. Positive feedback mechanisms between the organism and the environment (e.g., Reference [61]) require special attention, because they can generate spatial patterns that are fundamentally different from those that are generated by the gCP [72]. The gCP’s transition from the living to the extinct state hinges on the emergent spatial pattern. Background heterogeneity in addition to the gradient (e.g., the occurrence of non-habitable or less habitable patches) can also disturb the phase transition [73], depending on the quality, amount, size, shape, and persistence of the patches (see some examples in References [28,74]). These effects also require further investigation.

Understanding the fine-scale structure and dynamics of range margins is indispensable for predicting the ecological effects of climate changes. Due to the rapidity of anthropogenic changes, we need predictions even on the fine scale of years or decades. Considering the usually limited capacity for data collection in the field, good models are pivotal for predictions and monitoring. I think that linking ecological knowledge with the theory of phase transitions can advance this research significantly.