Connecting metapopulation heterogeneity to aggregated lifetime statistics

TL;DR

This paper investigates how heterogeneity in patch properties influences the aggregated lifetime statistics in metapopulations, revealing that variability alone can induce critical-like scaling laws, thus affecting interpretation of such data.

Contribution

It introduces a weakly-coupled metapopulation model linking patch heterogeneity to aggregated lifetime distributions, highlighting the role of variability in critical phenomena.

Findings

Patch heterogeneity can produce critical-like scaling laws in aggregate data.

A connection between patch property variability and lifetime distribution is established.

Heterogeneous systems can be equivalent to homogeneous critical systems under certain conditions.

Abstract

Aggregated metapopulation lifetime statistics has been used to access stylized facts that might help identify the underlying patch-level dynamics. For instance, the emergence of scaling laws in the aggregated probability distribution of patch lifetimes can be associated to critical phenomena, in which the correlation length among system units tends to diverge. Nevertheless, an aggregated approach is biased by patch-level variability, a fact that can blur the interpretation of the data. Here, I propose a weakly-coupled metapopulation model to show how patch variability can solely trigger qualitatively different lifetime probability distribution at the aggregated level. In a generalized approach, I obtain a two-way connection between the variability of a certain patch property (e.g. carrying capacity, environment condition or connectivity) and the aggregated lifetime probability distribution. Furthermore, for a particular case, assuming that scaling laws are observed at the aggregated-level, I speculate the heterogeneity that could be behind it, relating the qualitative features the variability (mean, variance and concentration) to the scaling exponents. In this perspective, the application points to the possibility of equivalence between heterogeneous weakly-coupled metapopulations and homogeneous ones that exhibit critical behavior.

Equations

$${\mathcal{E}}_i=\{T_{i,0}^{†},T_{i,0},T_{i,1}^{†},T_{i,1},\dots \},$$

$$n_i\sim \frac{T_o}{⟨T_i⟩+⟨T_i^{†}⟩}.$$

$$n(m)\sim \frac{N(m)T_o}{\tau (m)+{\tau }^{†}(m)}.$$

$$w(m)=\frac{\mathcal{N}\rho (m)}{\tau (m)+{\tau }^{†}(m)},$$

$$\mathcal{P}(T)={\int }_{\Omega }w(m)p(T|m)dm.$$

$$\mathcal{P}(T)=\mathcal{N}{\int }_{\Omega }\frac{\rho (m)}{\tau (m)+{\tau }^{†}(m)}p(T|m)dm.$$

$$\mathcal{P}(T)=\mathcal{N}{\int }_{\Omega }\frac{\rho (m)}{\tau (m)[\tau (m)+{\tau }^{†}(m)]}{e}^{-T/\tau (m)}dm.$$

$$\mathcal{P}(T)=\mathcal{N}{\int }_{s(0)}^{s(\infty )}\frac{\rho (m(s))m^{\prime }(s)s^2}{1+s{\tau }^{†}(m(s))}{e}^{-Ts}ds,$$

$$\mathcal{P}(T)=-\mathcal{N}{\int }_0^{\infty }\frac{\rho (m(s))m^{\prime }(s)s^2}{1+s{\tau }^{†}(m(s))}{e}^{-Ts}ds.$$

$$\mathcal{P}(T)=-\mathcal{N}\mathcal{L}\left[\frac{\rho (m(s))m^{\prime }(s)s^2}{1+s{\tau }^{†}(m(s))}\right].$$

$$\rho (m)=-\mathcal{N}h(m){\mathcal{L}}^{-1}[\mathcal{P}(T){]}_{s=1/\tau (m)},$$

$$h(m)\equiv \frac{d\tau (m)}{dm}\left(1+\frac{{\tau }^{†}(m)}{\tau (m)}\right),$$

$$\tau (m)=e^{am}-1.$$

$$\mathcal{P}(T)\propto \frac{e^{-T/{\tau }_c}}{({\tau }_0+T{)}^{\beta }},$$

$$\rho (m)\propto h(m)[S(m){]}^{\beta -1}{e}^{-S(m){\tau }_0},\text{for }m\in [0,m^{\star }]$$

$$m^{\star }=\frac{1}{a}\ln ({\tau }_c+1).$$

$$h(m)=a{e}^{am}\left(1+\frac{1}{m(e^{am}-1)}\right)$$

$$\displaystyle\rho(0)\propto\frac{e^{-\frac{\tau_{0}}{am}}}{m^{\beta+1}}\biggr{|}_{m\to 0}\,.$$

$$\rho(m^{\star})\propto(m^{\star}-m)^{\beta-1}\biggr{|}_{m\to m^{\star}}\,.$$

Full text

I Introduction

A metapopulation is composed of subpopulations that live within disconnected habitats (patches) but can interact through individual dispersal. Locally, at each patch, due to demographic stochasticity, populations are vulnerable to extinction. However, the metapopulation as a whole tends to be stable as individual dispersal can balance fluctuation and recolonize patches Hanski (1999). This interplay between local and nonlocal dynamics generates, for each patch, a collection of time intervals that correspond either to an empty or occupied state. With this dataset statistical analysis can be performed to extract relevant information about the metapopulation dynamics Hanski and Ovaskainen (2000a); Hanski (1999); Tilman et al. (1994); Kuussaari et al. (2009).

Aggregated approaches (i.e., grouping data from all patches and species) has been used to access stylized facts, providing a deepen view of inter-patch interaction. Previously, the existence of power-law scaling in the aggregated distribution of species abundance and lifetimes (the occupied time interval) has been associated with critical phenomena that generates interaction between system units across a broad range of scales Allen et al. (2001); Keitt and Stanley (1998); Stanley (1971); Solé et al. (1999). Nevertheless, it has been also pointed out that heterogeneity could mimic similar statistical features, a fact that can compromise the interpretation of the underlying dynamics Allen et al. (2001).

This work deepens previous discussions on the role of heterogeneity in aggregated datasets promoted in Ref. Allen et al. (2001). Focusing on the aggregated lifetime statistics, the following results provide an explicit analytical relation between patch-level variability and the aggregated lifetime probability distribution. For that, I propose a weakly-coupled metapopulation model, in which dispersal timescales are longer than population dynamics correlation times. This way, each patch dynamics is constructed by a local part and an effective term to account for the flux of individuals from other patches (see Fig. 1). In this framework, in Sec. II.1, I start by characterizing the structure of the lifetime dataset. In Sec. II.2, the aggregated lifetime probability distribution is written as a superstatistics of the local dynamics Beck and Cohen (2003), i.e. by mixing each patch exponential distribution with different characteristic timescales Allen et al. (2001). Finally, in Sec. II.3, I obtain the reverse transformation that can reveal patch heterogeneity from the aggregated lifetime probability distribution. Overall, even though, locally, lifetimes are exponentially distributed, due to the variability in the patches’ properties (such as its carrying capacity, environmental condition and connectivity), the aggregated lifetime probability distribution could exhibit forms that deviate from the exponential form Beck and Cohen (2003); Allen et al. (2001).

In Sec. III, by applying the model, I speculate the patch variability behind a generalized power-law form for the aggregated lifetime distribution. This particular form is inspired by the North American Breeding Bird Survey dataset Sauer et al. (1966); Keitt and Stanley (1998) for which was invoked the presence of critical phenomena to justify the scaling laws Allen et al. (2001). Thus, these applied results aim to show the possible equivalence between heterogeneous weakly-coupled metapopulations and homogeneous ones that exhibit critical behavior at the aggregated-level. The relation between the scaling exponents with the qualitative features of the patch-level variability (mean, variance and concentration) is explored in detail.

The idea of connecting patch and global levels has been a central issue in metapopulation dynamics studies Hanski (1999); Hanski and Ovaskainen (2000b). However, a rigorous mathematical treatment for stochastic spatial explicit models faces challenges in providing expressions for macroscopic quantities. Despite recent developments on approximative methods, such as moment closure or perturbative expansions around the mean-field solution, only superficial aspects of the system are straightforward to access Ovaskainen et al. (2014); Meerson and Ovaskainen (2013). The presented weakly-coupled model bypass these challenges, allowing to demonstrate the emergence of scaling purely due to heterogeneity. This is a necessary reference frame to understand the ecosystem across scales Butler et al. (2017); Lafuerza and Toral (2013) and to ensure a proper identification of emergent behavior caused by ecosystem variability (its landscape structure and species), which can play a role in its resilience, stability and functions Loreau et al. (2001, 2003); Solé et al. (1999); Stanley (1971).

II Framework

In the following sections, I start by characterizing the lifetime dataset structure and establishing the direct transformation from patch to the ensemble (aggregated) level, obtaining the lifetime probability distribution considering data from all patches. Then, I find the inverse transformation that reveals patch variability from the aggregated perspective. Later in Sec. III, I apply the approach assuming specific population dynamics models and discuss the appearance of power-law statistics.

II.1 Dataset structure

At patch-level, stochastic fluctuations (demographic and environmental) eventually drive the population to the null (extinction) state, which is locally-absorbing Ovaskainen and Meerson . However, since each patch’s boundary is open, new individuals are able to reestablish the population. The combination of these two forces places the patches out of equilibrium in a succession of colonization and extinction events.

Experimental records for a metapopulation constituted by $N$ patches can be structured as

[TABLE]

where $i\in(1,2,\ldots,N)$, $T_{i}$ is the time interval in which the patch $i$ was occupied, its lifetime, and $T^{\dagger}_{i}$ represent the time interval in which the patch was empty Hanski et al. (1994).

For each patch, the number of records $n_{i}\equiv|\mathcal{E}_{i}|$ is given by

[TABLE]

where $T_{o}$ is the total observation time.

Since patches with same features will exhibit same statistical trends, we can rewrite Eq. (2), defining the mean lifetime as a function of the property $m$,

[TABLE]

$n(m)$ gives the total number of events for all patches with property $m$ within the observation time, being $N(m)$ the number of patches with property $m_{i}\in[m_{i}-\Delta m,m_{i}+\Delta m]$ with $\Delta m/m\ll 1$. In Fig. 2, for a particular model that we will discuss later (see also supplementary material) a typical time series is depicted, highlighting patches’ occupied and empty time intervals.

From the aggregated perspective, we can see that patches with small lifetimes have a higher importance in the statistics, due to large event counting. Besides that, these weights are modulated by the likelihood of having a particular patch property in the metapopulation. Thus, in the large metapopulation limit, $N\gg 1$, the weight $w(m)$ of a particular property $m$ is given by

[TABLE]

where $\rho(m)=\frac{N(m)}{N}$ is the probability density function for observing a certain property $m$ in the metapopulation and $\mathcal{N}^{-1}=\int_{\Omega}\frac{\rho(m)}{\tau(m)+\tau^{\dagger}(m)}dm$ is a normalization factor.

II.2 Superstatistics

In order to construct the aggregated lifetime probability distribution $\mathcal{P}(T)$, it is necessary to mix those of the patches’ accounting for the contributions of different local properties. Namely, in a given domain $\Omega$,

[TABLE]

where $p(T|m)$ is the lifetime probability distribution for a patch having characteristic $m$ and $w(m)$ is the weight of the characteristic in the aggregated set. Substituting the weight given by Eq. (4) in (5),

[TABLE]

It has been shown that lifetimes are, rather generally, exponentially distributed independently of the model specificities Grimm and Wissel . In the framework of Markov processes, being $q(n,t)$ the probability of having $n$ individuals at instant $t$, its temporal evolution is given by $\dot{q}=M(q)q$, where matrix $M$ sets the transition rates for increasing ($n\to n+1$) and decreasing ($n\to n-1$) the number of individuals. Typically, it is found that the leading eigenvalue $\Lambda_{1}$ of $M$ is much smaller than the others, even when including density dependency, environment noise and other features Grimm and Wissel . As a consequence, after a fast transient $q_{n>0}(t)\propto e^{-\Lambda_{1}t}$ the extinction probability becomes $q_{0}(t)\propto 1-Ae^{-\Lambda_{1}t}$. Hence, the lifetime probability distribution is given by $p\propto e^{-\Lambda_{1}t}$ with mean lifetime $\tau\simeq 1/\Lambda_{1}$ Grimm and Wissel . In the supplementary material, I investigate in detail the lifetime distribution for the canonical model Ovaskainen and Meerson . There, exponential laws were also found.

Motivated by this feature, setting $p(\tau|m)=e^{-T/\tau(m)}/\tau(m)$ in Eq. (5), we obtain

[TABLE]

Eq. (7) defines the aggregated lifetime probability $\mathcal{P}$ as a mixing of the subpopulation dynamics in the patch property space, constituting a direct transformation that gives us the aggregated statistics from the patch-level variability.

II.3 Deconvolution of the aggregated data

In this section, I present the final step of the general approach, which involves the inverse transformation that can recover the patch variability from the aggregated statistics. In order to avoid unnecessary complications, I will focus on a single patch property that belongs to a positive interval, $m\in\Omega\equiv[0,\infty)$, and positively influences population persistence. Such property could be resource availability, patch quality or dispersal rate, for instance.

Performing a change of variable in Eq. (7), setting $\tau=1/s$ and writing $m(s)$, we obtain the integral in $s$

[TABLE]

where prime notation means $m^{\prime}(s)=dm(s)/ds$. At this point, note that in order to write $m(s)$, the function set by the dynamics $s(m)$ needs to be invertible.

Since here it is assumed that the mean lifetime increases with the habitat property $m$, we can set the intervals to $s(\infty)=0$ and $s(0)=\infty$. Thus,

[TABLE]

Identifying that Eq. (9) resembles the Laplace transform structure $\mathcal{L}[f](T)=\int_{0}^{\infty}f(s)e^{-Ts}ds$, it is straightforward to obtain

[TABLE]

Finally, within the Laplace transform formalism, it is possible to apply the inverse transformation, recovering the distribution $\rho$ of property $m$ from the aggregated distribution $\mathcal{P}$,

[TABLE]

where

[TABLE]

in which it was used that $s(m)=1/\tau(m)$ and $-s^{2}m^{\prime}(s)=\frac{dm}{d\tau}$.

III Application

Eqs. (10) and (11) provide us a two-way connection between patch-level variability ($\rho$) and the aggregated lifetime probability distribution ($P$). In order to show an application of this result, specific choices for how the mean extinction time $\tau$ depends on the patch property $m$ need to be made. I choose to set that $\tau$ grows exponentially with $m$,

[TABLE]

For this choice, if $m$ is infinitely large, the population persists indefinitely. Conversely, when $m$ is negligible, the mean lifetime goes to zero.

The proposed exponential form in Eq. (13) is motivated by the results found for when $m$ is the carrying capacity in standard population dynamics models, either by analytical or computational means Ovaskainen and Meerson . This form is also observed when investigating the role of migration rates and patch connectivity, as reported in the supplementary material. For an effective single patch model, I show that its lifetime is exponentially distributed with mean lifetime that grows exponentially with the migration rate. For a spatial explicit formulation (under weak-coupling), considering a random patch landscape, analogous results were also observed, where the mean lifetime increases exponentially with connectivity. Thus, for the following results, the property $m$ can be interpreted as the patch carrying capacity, inward flux of individuals or connectivity, since all of them influence the mean lifetime in a similar way.

To specify the function $h$ in Eq. (12) it is also necessary to set how the empty time $\tau^{\dagger}$ depends on $m$. Particularly, I will set $\tau^{\dagger}(m)=1/m$. This extra assumption regarding $\tau^{\dagger}$ applies for the scenarios in which $m$ refers to the migration rate (determining the rate of arrival of new individuals); or represents the carrying capacity when flux is biased towards high-quality patches (informed dispersal) Colombo and Anteneodo (2015); Bengfort et al. (2016). For the case in which $m$ do not couple to the migration process (uninformed dispersal), $\tau^{\dagger}$ will be constant, meaning homogeneous migration rates. For the following results, considering either $\tau^{\dagger}=1/m$ or constant does not produce qualitative changes.

Once the model mediation is set by electing $\tau(m)$ in Eq. (13), given the patch variability $\rho(m)$, it is possible to obtain the aggregated lifetime distribution $\mathcal{P}$ from Eq. (10). More interestingly, it is possible to perform the reverse transformation (shown in Eq. (11)), revealing the patch variability in the aggregated statistics. So, I propose here a generalized power-law form for the aggregated lifetime probability distribution, namely

[TABLE]

where the parameters $\tau_{c}>0$ controls the exponential cut-off, and $\tau_{0}>0$ sets the probability of vanishing lifetimes and ensures normalization for the whole range $\beta>0$, assuming $T\in(0,\infty)$. This choice is motivated by the North American Breeding Bird Survey dataset, for which was found $\beta=1.6$ and $\tau_{c}=14\text{ (years)}$ Keitt and Stanley (1998), and were associated to critical phenomena Keitt and Stanley (1998). The following results aim to obtain the equivalent patch-level variability, that can generate same statistics at the aggregated level.

Substituting Eq. (14) in (11), and noting that Eq. (14) has a common Laplace transform (i.e., the integrate unit step with frequency and time shifts), we finally obtain that

[TABLE]

with $S(m)=1/\tau(m)-1/\tau_{c}$, being

[TABLE]

The function $h$ is obtained from Eq. (12), substituting Eq. (13),

[TABLE]

where we additionally consider that $\tau^{\dagger}(m)=1/m$, as discussed previously.

In Fig. 3, it is shown how the probability distributions for patch variability $\rho(m)$ changes with the scaling exponent $\beta$ of Eq. (14). Different values of the exponent reflects in qualitative differences in patch-level variability. For large $\beta$ values, the probability of observing patches with characteristic $m$ decays with the values of $m$, such that patches with high resource availability, migration rate or connectivity ($m\sim m^{\star}$) are rare in the system. In contrast, for small values of $\beta$ this behavior changes and the system exhibits a maximum at $m^{\star}$, being the wealthy patches (the ones with high connectivity and resource availability) numerous in the metapopulation.

To identify the qualitative different distribution in Fig 3, it is helpful to analyze the behavior of $\rho$ in Eq. (15) under the limits of vanishing and maximum value of $m$. For $m\to 0$, the mean lifetime becomes too small, being experimentally challenging to access this region. In Fig. 3, as a reference, the gray highlighted region corresponds to mean lifetime $\tau<1$, where the behavior of $\rho$ would not be captured by experimental data sampled in intervals with size $\Delta T\geq 1$ (every year, for example). Despite that, analytically, from Eq. (15),

[TABLE]

Since, in Eq. (14) our proposal sets $\tau_{0}>0$ to ensure normalization, $\rho(m\to 0)$ always vanishes, as seen in Fig. 3.

In the limit of maximum $m$ (high-quality limit), taking $m\to m^{\star}$ in Eq. (15),

[TABLE]

Consequently, for $\beta>1$ ($\beta<1$), $\rho(m^{\star})$ goes to zero (diverges), indicating the different regimes of metapopulation heterogeneity (as shown in Fig. 3). Only for $\beta=1$, $\rho(m^{\star})$ assumes a positive finite value. Note that Eq. (19) does not depend on other parameters aside from $\beta$.

More detailed information about the distributions shown in Fig. 3 can be extracted obtaining the moments of the patch property distribution. In Fig. 4, from the probability distribution given by Eq (15), I numerically compute the mean and the standard deviation of patch properties in the metapopulation for different values of $\beta$. As $\beta$ increases, the mean value of $m$ shifts towards small values, while the standard deviation passes through a maximum at $\beta\simeq 1.7$ (as indicated by the vertical dashed line). The particular value of the exponent that indicates maximum variability is sensitive to the parameters choice. Curiously, this optimal value of $\beta$ is close to the one observed in Ref. Keitt and Stanley (1998).

IV Final remarks

The origins of the power-laws observed in the aggregated lifetime probability distribution have been a matter of debate Keitt and Stanley (1998); Allen et al. (2001). On one hand, power-laws can be a consequence of critical phenomenon Stanley (1971), generating interconnected scaling laws for the system macroscopic state variable Stanley (1971); Taubert et al. (2018). On the other hand, however, we see that patch-level variability can also control these power-laws. During the previous sections, I studied the connection between the patch-level variability and the aggregated lifetime probability distribution. Considering a weakly-coupled metapopulation, the results demonstrate explicitly how patch variability can solely trigger different outcomes at the aggregated level. Understanding the intrinsic features of the dataset, I first wrote the aggregated lifetime probability distribution as a superstatistics Beck and Cohen (2003) of the patches dynamics. Following that, an inverse transformation was obtained, revealing patch property variability from the aggregated perspective.

Deepening previous debates Allen et al. (2001); Keitt and Stanley (1998), I applied these results to find the specific class of patch-level variability that generates scaling laws in the metapopulation aggregated lifetime probability distribution, mimicking the outcomes that can also be associated to critical phenomena Allen et al. (2001). Particularly, I investigated a generalized version of the power-law form observed for North American Breeding Bird Survey dataset Keitt and Stanley (1998); Sauer et al. (1966) which was the source of the debate in Ref. Allen et al. (2001). It was then presented that the observed scaling exponent $\beta$ is related to different types of heterogeneity. In the case of $\beta>1$, high-quality patches (well connected or with high resources availability) are rare in the metapopulation. While for $\beta<1$, this feature is inverted. Moreover, for a nontrivial value of $\beta$ ($\beta\simeq 1.7$, in the case depicted in Fig. 4), patch variability achieves a maximum. The results have shown the connection between a macroscopic characteristic of the aggregated set ($\beta$) to patch-level variability ($\rho(m)$).

In sum, the application of the proposed framework indicates that, at the aggregated level, it is possible to have an equivalence between a heterogeneous weakly-coupled system and a homogeneous one exhibiting critical phenomena Allen et al. (2001); Lafuerza and Toral (2013). This stresses the role of patch-level variability in the aggregated data, which cannot be overlooked when aiming to interpret the underlying dynamics Allen et al. (2001).

In order to achieve our results, some specific choices were made. Firstly, I assumed that the metapopulation is weakly coupled. On a first look, this approximation might be seen as a weakness of the results, however, it reveals its strength by demonstrating the emergence of power-laws in the aggregated lifetime statistics in the absence of strong correlations, solely triggered by patch heterogeneity Allen et al. (2001); Keitt and Stanley (1998). Secondly, the developments presented assumed that patch lifetimes are exponentially distributed. It was argued that this aspect is rather general, being present independently of the local population dynamics model. Nevertheless, it is worth to remark that in some special cases this feature might be violated at intermediate and short scales Artime et al. (2018); Ovaskainen and Meerson ; Pigolotti et al. (2005). Then, for these cases, the use of the Laplace transform formalism would be compromised.

It would be relevant to improve certain points of the results for practical situation. For instance, the calculations could be extended to the cases in which variability is present in several patch proprieties. However, to accomplish that, it is necessary to ensure the consistency of the steps done in the Sec. II.3. Also, fortunately, it was straightforward to apply the results to the aggregated probability distribution in the form of Eq. (14), since it has a common Laplace inverse. For special cases, one should develop a generalized transformation such as Mellin transformations Whittaker and Watson (1996).

At last, the results here developed can also be posed in a broader sense, as long as the dataset involves a collection of weakly correlated stochastic processes at different scales. For instance, in the context of animal movement, the aggregated probability distribution of the jump sizes is obtained by mixing data from random walkers that explore space at different spatial and temporal scales Viswanathan et al. (2011); Rodríguez et al. (2017). In this case, power-laws in movement statistics could emerge due to heterogeneity in animal behavior Keitt and Stanley (1998); Viswanathan et al. (2011); Rodríguez et al. (2017).

Supplementary material

S1 Appendix.

Relation between lifetimes and patch properties. Classical scenarios are investigated showing how mean population lifetime depends on patch properties.

Acknowledgments

I would like to thank Profs. Emilio Hernández-Gárcia, Celia Anteneodo and Raúl Toral for the support, comments and suggestions on the manuscript. I also acknowledge funds from the Spanish Research Agency through grant ESOTECOS FIS2015-63628-C2-1-R (AEI/FEDER, EU) and grant MDM-2017-0711 from the Maria de Maeztu Program for units of Excellence in R&D.