## Abstract

A spatial metapopulation is a mosaic of interconnected patch populations. The complex routes of colonization between the patches are governed by the metapopulation's dispersal network. Over the past two decades, there has been considerable interest in uncovering the effects of dispersal network topology and its symmetry on metapopulation persistence. While most studies find that the level of symmetry in dispersal pattern enhances persistence, some have reached the conclusion that symmetry has at most a minor effect. In this work, we present a new perspective on the debate. We study properties of the in- and out-degree distribution of patches in the metapopulation which define the number of dispersal routes into and out of a particular patch, respectively. By analysing the spectral radius of the dispersal matrices, we confirm that a higher level of symmetry has only a marginal impact on persistence. We continue to analyse different properties of the in–out degree distribution, namely the ‘in–out degree correlation’ (IODC) and degree heterogeneity, and find their relationship to metapopulation persistence. Our analysis shows that, in contrast to symmetry, the in–out degree distribution and particularly, the IODC are dominant factors controlling persistence.

## 1. Introduction

The field of metapopulation dynamics was pioneered by Levins [1] more than four decades ago, and it has remained a highly active ecological research agenda ever since [2–5]. Metapopulations are composed of a set of patch populations which are often sustained by colonization between patches. A key area of research requires understanding the factors which determine the fate of a metapopulation, namely whether it will thrive and persist or whether it will become extinct as it evolves in time. This requires studying the structure, landscape and dispersal patterns of the metapopulation, as well as species interactions, and demographic variables (birth and death rates), which may all have a dramatic impact on viability [6–14]. More recently, special attention has been given to the specific issue of connectivity and dispersal pattern and their effects on metapopulation persistence [3,15–18]. In this respect, the network connectivity is determined by the dispersal pathways between patches, so that the ‘dispersal matrix’ plays an important role.

As noted in [16], traditionally, most studies assumed a symmetric dispersal pattern, meaning an equal probability of dispersal from site A to site B as from B to A. However, the assumption of symmetric dispersal is now considered to be too simplistic. In nature, asymmetries in dispersal can occur as a result of ocean currents, elevation gradients, wind and numerous other mechanisms. In the past decade, several studies have investigated the effect of dispersal network symmetry on metapopulation persistence [15–18]. Vuilleumier & Possingham [15] simulated random dispersal patterns and analysed their influence on metapopulation survival. The general conclusion of that study was that an asymmetric dispersal pattern may lead to a decrease in metapopulation viability. Bode *et al.* [16], by analysing correlations between metapopulation viability and different dispersal pattern measures, also found that asymmetric dispersal leads to a higher extinction risk. Vuilleumier *et al.* [17] gave analytical and numerical results for different models of dispersal that supported these previous studies.

More recently, Kleinhans & Jonsson [18] attempted to isolate the effect of asymmetry on metapopulation viability. They formulated a symmetry index that allows a rigorous categorization of the level of dispersal pattern symmetry. By simulating the metapopulation dynamics on networks with different levels of symmetry, they found that symmetry has no effect or only a negligible effect on metapopulation persistence. This finding contradicts many earlier works and has thus led to controversy regarding the effect of symmetry on metapopulation viability. The controversy is further complicated as different researchers use different models and methods in an attempt to generalize the results.

Here, we present a different perspective on the issue which can potentially resolve some of the disagreement and shed new light on the controversy. In order to proceed, we need to characterize the dispersal network in more detail. We suppose that each node in the network represents a patch in the metapopulation, whereas directed links between the nodes represent colonization routes between the patches. A key concept is the notion of a patch's ‘in- and out-degree’. These are defined as the number of colonization routes entering and leaving the patch, respectively. Thus, in the network illustrated in figure 1*a*, the red (lighter) node has an in-degree of *d*_{in} = 3 indicating that it may be colonized from three other patches. The same red (lighter) node has an out-degree of *d*_{out} = 2, indicating that individuals on this patch may colonize to two other patches. The degree heterogeneity (from here on simply ‘heterogeneity’) of the dispersal network is characterized by the variability in the distribution of *d*_{in} and of *d*_{out}. In a heterogeneous network, the in- and out-degrees of the patches can vary considerably from patch to patch. By contrast, in a so-called regular network, all patches have exactly the same in- and out-degrees, implying *d*_{in} = *d*_{out} for all nodes, as illustrated in figure 1*b*. Note that a patch with low in-degree, but high out-degree might be considered a ‘source’, whereas the converse would be considered a ‘sink’, concepts we elaborate on in §5. We also point out that degree heterogeneity is only one of several ways of incorporating heterogeneity into the model. At the end of §2, we shortly discuss a different kind of heterogeneity, namely the possibility of varying extinction probabilities which may also depend on connectivity.

An important finding in the study of metapopulations is that persistence of a metapopulation is determined by the spectral radius of the dispersal matrix [11,19], i.e. the absolute value of the eigenvalue of greatest magnitude. However, with some exceptions, e.g. [3,11], most of these studies of the spectral radius have dealt with undirected dispersal networks in which a colonization route from patch-*i* to *j* guarantees colonization from patch *j* to *i*. A more realistic analysis should cover the important case of directed dispersal networks. Nevertheless, this has been avoided because of the difficulties in characterizing the spectral radius of directed networks. Here, we demonstrate that the spectral radius of the dispersal matrix is determined by the in- and out-degree distributions of the nodes in the network. More specifically, the spectral radius can be well approximated using the so-called in–out degree correlation (IODC) of the dispersal network [20]. In the simplest interpretation, a high IODC means that patches with high in-degree tend to have high out-degree and vice versa, whereas a low IODC means that patches having high in-degree tend to have low out-degree and vice versa as is illustrated in figure 2.

We explore the IODC in depth shortly, and show that metapopulation viability strongly depends on it. As we explain, this insight can potentially disentangle the contradicting results regarding the effect of symmetry on persistence.

## 2. The stochastic metapopulation model and spectral radius criterion

To analyse metapopulation dynamics, we use the framework of stochastic patch occupancy models (SPOMs) from previous studies [15–18]. In these models, one assumes a network of patches where each patch is either occupied (1) or extinct (0). In each time step of the model, each patch can go extinct with probability *e*. We also assume that dispersal events between connected patches, from colonized patches to extinct patches, occur with probability *c*. For a metapopulation consisting of *N* patches, the dispersal network, *M*, is an *N* × *N* matrix in which *m _{ij}* = 1 if patch-

*j*is connected to patch-

*i,*whereas

*m*= 0 otherwise. Note that the network is directed in the sense that a connection from patch-

_{ij}*i*to

*j*does not necessarily imply a connection from patch-

*j*to

*i*(i.e.

*m*does not necessarily equal

_{ij}*m*). An illustration of a random directed network is given in figure 1

_{ji}*a*.

We first provide a formulation of the metapopulation model and then demonstrate that the spectral radius of the dispersal network *M* governs metapopulation persistence. Ovaskainen *et al.* [21] present a related model and eigenvalue result, but only provide a brief sketch of the details which we find instructive and useful for understanding the key ideas of this paper.

We begin by defining *p _{i}*(

*t*) to be the probability that patch

*i*is occupied at time step

*t*. In a mean-field approximation, we assume that the probabilities

*p*(

_{i}*t*+ 1) depend only on the probabilities

*p*(

_{i}*t*). In the case of random networks, it was shown that this approximation is quite accurate [22,23]. Therefore, we can write

*p*(

_{i}*t*+ 1) as a sum of two contributions: either patch

*i*was colonized at time

*t*, and then, we multiply by the probability of it not becoming extinct during one time step, or it was extinct at time

*t*and was recolonized in the last time step 2.1

By the definitions above
2.2a
2.2b
2.2cEquation (2.2*b*) is simply the complement of the probability of patch *i not* being recolonized. Note that equation (2.2*c*) is a first-order approximation in *p _{i}*(

*t*) which is accurate when . We obtain the following expression for

*p*(

_{i}*t*+ 1): 2.3

Now, we are interested in whether a metapopulation will persist or not, and we therefore check the linear stability of the extinction state (*p _{i}*(

*t*) = 0, ∀

*i*).

If the extinction state is stable, then the metapopulation is prone to extinction, as this state will be an attractor for the dynamics of *p _{i}*(

*t*). On the other hand, its instability implies the existence of a meta-stable state in which the metapopulation will persist for a long time. In this case, the first-order approximation in

*p*(

_{i}*t*) will no longer hold, and we can speculate that the extinction time will be exponential in the number of occupied patches in the equilibrium. To complete the stability analysis, we retain only linear terms in

*p*(

_{i}*t*) (i.e. close to the extinction state), which gives 2.4where Δ

*p*(

_{i}*t*) =

*p*(

_{i}*t*+ 1) −

*p*(

_{i}*t*). In matrix form, the last equation becomes 2.5where

**M**is the dispersal network matrix,

**I**is the identity matrix and

**p**and

**Δp**are the vectors of

*p*(

_{i}*t*) and

*p*(

_{i}*t*+ 1) −

*p*(

_{i}*t*), respectively. The last expression implies that the extinction state is unstable and the metapopulation will persist if 2.6Here,

*λ*is the spectral radius of

_{m}*M*, which according to the Perron–Frobenius theorem [24] is necessarily real for all irreducible dispersal matrices, symmetric or asymmetric. By irreducible it is meant that any patch population in the metapopulation is connected to, and may be reached from, any other patch population. When the eigenvalue criterion is met, the metapopulation will persist for long timescales, but will eventually become extinct owing to demographic fluctuations. Otherwise, it will go extinct with the number of occupied patches declining as an exponential function of time,

*N*

_{occupied}∼ exp(−

*βt*), (implied by equation (2.5)), i.e. exponentially fast.

Combined with equation (2.6), this shows that a larger spectral radius of the dispersal matrix leads to persistence for a larger range of the ratio between the extinction and colonization parameters. In the case of a time varying environment, in which *e* and *c* are subject to noise, we expect that a higher spectral radius can lead to more robustness and resilience of the metapopulation. Therefore, for the purpose of assessing the long-time behaviour of a metapopulation, it suffices to know its spectral radius.

As mentioned, we assumed a constant extinction probability, *e*. In a more realistic model, the extinction probability can vary between patches, depending on the conditions on a specific patch, and might also be correlated to other connectivity properties such as in or out degrees. However, as we wish to focus solely on connectivity effects, we assume a constant extinction probability for all patches, a practice that is not uncommon in similar metapopulation modelling approaches [15,16,18].

## 3. To what degree does symmetry impact metapopulation persistence?

Having described the metapopulation model and assumptions and derived the spectral radius criterion, here we intend to examine the impact of network symmetry on metapopulation persistence. Symmetry is the tendency of patches to be connected in both directions, i.e. if there is a dispersal route from patch *i* to patch *j*, then there will be a dispersal route from *j* to *i*. We first rigorously define an index for the level of symmetry and afterwards relate it to metapopulation persistence through a numerical calculation of dispersal matrices' spectral radii and through simulations.

We follow [18] and define the degree of symmetry of a dispersal matrix, *γ*, as the proportion of symmetric links out of all links. That is, *γ* is equal to the number of links where patch *i* is connected to *j* and *j* to *i*, divided by the total number of links. Thus, *γ* = 0 represents a fully asymmetric matrix and *γ* = 1 corresponds to a completely symmetric pattern. More rigorously
3.1In equation (3.1), the numerator is the number of symmetric links, and the denominator is the total number of links. In order to check the effect of symmetry on metapopulation persistence, we relate the spectral radius of the dispersal matrix, *λ _{m}*, to the degree of symmetry,

*γ*. The algorithm in [18] was used to generate 100 matrices of dimension

*N*×

*N*= 1000 × 1000 for each value of

*γ*. Figure 3

*a,b*shows the mean spectral radius (solid line), obtained by averaging over 100 matrices, as a function of symmetry,

*γ*, for matrices with mean degree

*d*= 5 and

*d*= 10, respectively. Also plotted are 95% CIs (dotted-dashed line; obtained by removing the upper and lower extreme 2.5% of the same 100 matrices). The linear increase of the mean spectral radius with

*γ*is evident, however, it is rather minor in extent given there is only a 10–20% change exhibited over the full range of symmetry.

To confirm the mild effect of symmetry on the persistence expected by the relatively small increase in spectral radii, stochastic simulations were also performed by using an SPOM. Again, for varying values of *γ*, 50 dispersal patterns with *N* = 1000 were generated, and the dynamics was simulated for different values of *c*, the colonization parameter.

Simulations were initiated by choosing 10 random patches. In each time step, a colonized patch was randomly chosen, and all possible future events involving it were enumerated and listed, namely its extinction or colonization of any neighbouring patch. Then, one of these possible events was randomly drawn and executed with probability *e* for extinction and *c* for colonization. This simulation protocol mimics the more accurate Gillespie stochastic simulations [25] yet uses only a fraction of the computing resources, because it is not necessary to enumerate all possible events of the full metapopulation and calculate their probabilities in each time step. The simulation run ends in one of three conditions: (i) all patches are extinct; (ii) the number of occupied patches reaches 10% of the total (100 patches in the case of *N* = 1000); and (iii) a threshold time of 100 *N* (10^{5} steps for *N* = 1000) is reached without all nodes becoming extinct. In the second and third cases, the run is deemed persisting, whereas in the first, it is extinct. The rationale is that cases (ii) and (iii) indicate that the extinction state is unstable and that there exists a persisting non-zero equilibrium as expected by the mean-field derivation. Note that there is no need to discuss here the unusual and unlikely case of marginal stability, in which one eigenvalue of the system is zero. Suffice to say, that in such a situation, the solution will slowly reach the extinction state by the stochastic fluctuations.

Figure 4 shows the results for mean degree *d* = 5 and *d* = 10, respectively. Figure 4 shows the fraction of persisting runs versus the colonization *c* for different values of *γ*. A slight trend can be seen with more symmetric matrices having a slightly lower threshold for persistence. While notable by eye, the results in figure 4 agree with [18] in that raising the symmetry does not result in a substantial increase in persistence.

## 4. The effects of the in–out degree distribution

In §3, we showed that the effect of symmetry on metapopulation persistence is negligible. Here, to obtain better ecological insight, we explore other properties of the dispersal network and their effect on persistence. Recall that for undirected networks the degree of a node is represented by a single number and therefore is determined by a single distribution, *p _{k}*, the probability of a random node to have degree

*k*. However, for directed networks, one needs to consider a two-dimensional distribution,

*p*the probability of a node having in-degree

_{jk}*j*and out-degree

*k*. The in–out degree distribution can be characterized by several attributes: 4.1 4.2 4.3 4.4The first and coarsest property of the in–out degree distribution is its mean (equation (4.1)). Because every link goes out from a patch to a different patch, the mean in-degree must equal to the mean out-degree. Network properties more subtle than the mean take into account higher moments. Second moment properties are in or out coefficients of variation ( or , equations (4.2) and (4.3)) and in–out degree correlation (IODC , equation (4.4)). Strictly speaking, equation (4.4) does not define a true Pearson correlation (e.g. it does not vary solely in the interval [−1,1]) or even a covariance (the mean is not subtracted from

*d*

_{in}and

*d*

_{out}). We use the term ‘correlation’ loosely to convey the main concept. Namely, when the positive quantity is large, the variables

*d*

_{in}and

*d*

_{out}are highly correlated, and when small, the variables are anti-correlated.

The coefficients of variation determine the heterogeneity of the in- and out-degree distributions. For example, a regular network, i.e. a network with all nodes having equal in- and out-degrees will have an in- and out-CV^{2} of zero. If the in- or out-degree varies a lot between different patches, then the network will have a high in or out CV^{2}, respectively. The IODC is a measure of the extent to which a patch has its in-degree close in number to its out-degree. Figure 2 shows example components of networks with a high and low IODC.

There are also other properties or measures of between-patch connectivity, for example, assortativity which is the tendency of high (in or out) degree nodes to be connected to other high (in or out) degree nodes. Such properties of the dispersal pattern are beyond the scope of this work. In §4a,b, we focus on the specific effects of IODC and in/out-degree heterogeneity.

### (a) The effect of the in–out degree correlation

Here, we examine the IODC of a dispersal network and its effect on the spectral radius of the network and also perform explicit simulations to corroborate our predictions. In [20], the following approximation for the spectral radius of a matrix of a large random directed network is derived:
4.5Here, *d*_{in}/*d*_{out} are the in/out degrees of a node. The brackets denote averaging over all nodes and is the mean degree of the dispersal matrix. The numerator on the right-hand side of equation (4.5) is the IODC of the network which is the tendency of nodes to have similar in- and out-degrees. This approximation implies that for networks with equal mean degree, the IODC is perhaps the most dominant factor which determines the spectral radius of the dispersal matrix and therefore metapopulation persistence. In practice, the mean degree and IODC of a given dispersal network are approximated by
4.6

As with the symmetry, to isolate the strong relationship between the spectral radius and the IODC, it is necessary to generate matrices with a constant mean degree and varying IODC. The task of generating large random matrices adhering to the aforementioned constraints is quite challenging. To accomplish it, we used the algorithm for generating a network given an arbitrary joint in–out degree distribution [26] combined with the following original functional form for *p _{jk}*:
4.7The parameters

*a*,

*b*and

*c*make it possible, in practice, to tune IODC while keeping an almost constant mean degree.

*D*is a normalization factor.

The main parameter for tuning the IODC is *c* as the last term in the multiplication represents a ‘repulsion’ or ‘attraction’ between *d*_{in} and *d*_{out} (*j* and *k* in the formula) according to the sign of *c*, thereby controlling the tendency of equality of in- and out-degrees. By taking different values for the parameters, we managed to obtain an IODC in the range between approximately 50 and 120 for matrices with *N* = 1000 and *d* ∼ 10. Figure 5*a,b* illustrates *p _{jk}* in the extreme cases of , 120 with mean degree

*d*∼ 10. In figure 5

*a*, preference is given to pairs with low in-degree and high out-degree or vice versa, whereas in figure 5

*b,*the converse is true, a preference is given to pairs in which the in-degree is equal to the out-degree. The above method generates suitable degree distributions

*p*, and with these established, the next step is to generate random connectivity matrices satisfying these degree distributions. As mentioned the algorithm in [26] was used for this purpose.

_{jk}Figure 6 shows the results of a numerical calculation of the spectral radius for the matrices of varying IODC. As expected, when the IODC increases 2.5-fold, the spectral radius also increases almost threefold. From this, it is clear that the IODC is crucial for determining the fate of a metapopulation to persistence or extinction.

Finally, to corroborate our predictions we also performed simulations for varying , *N* = 1000 and a close to constant degree of *d* ∼ 10. The simulations were performed as described in §3, and the results can be seen in figure 7. The dramatic expected increase in persistence, or decrease in the colonization threshold, with higher is evident.

### (b) The effect of in/out degree heterogeneity

In §4a, we saw the dramatic effect the IODC can have on metapopulation persistence for networks with the same mean degree. In this part, we examine how degree heterogeneity affects these previous results. As with the symmetry and the IODC, to single out the effect of heterogeneity, we must generate dispersal patterns with constant mean in/out degrees and varying CV^{2}. To this end, we used the gamma distribution, which makes it possible to specify the mean as well as the variance by tuning its scale and shape parameters. For the sake of simplicity, we restricted ourselves to the case where the in-degree distribution is identical to the out-degree distribution. Using the same in- and out-degree distribution leaves us with the question of how to pair the in- and out-degrees in each patch. We answered this question by examining three extreme cases: (i) the in-degree is exactly equal to the out-degree; (ii) match descending ordered in-degrees with ascending order out-degrees, so that a node with the highest in-degree has the lowest out-degree and vice versa; and (iii) random matching-choose random pairs, so that the probability of a node having out degree *j* is independent of the probability of a node having in degree *k* (*p _{jk}* =

*p*). In figure 8, we plotted the spectral radius (solid lines) and spectral radius approximation (dotted-dashed; equation (4.5)) of dispersal patterns with mean degree and CV

_{j}p_{k}^{2}varying between 0.1 and 4. The blue (uppermost) lines correspond to case (i) and show an approximately fourfold increase in spectral radius, indicating that dispersal matrices with high heterogeneity and high IODC will be the most robust. The red (lowermost) lines show case (ii), namely anti-correlated in- and out-degrees. In this case, it is notable that the point with highest spectral radius is also the point with minimal heterogeneity. This implies that in anti-correlated matching, heterogeneity in the degree distribution is detrimental to metapopulation survival. The green (middle) lines show case (iii), i.e. the case of random pairing of in- and out-degrees. In this case, spectral radius remains constant and approximately equals as can be expected by the eigenvalue approximation equation (4.5) along with equation (4.4) when

*p*=

_{jk}*p*.

_{j}p_{k}These extreme cases show that increasing CV^{2} is accompanied by an increase in the possible range of the spectral radius and IODC. Thus, a complex picture emerges, one cannot claim that heterogeneity on its own leads to better or worse persistence. Rather, heterogeneity carries potential for better persistence or for worse, where the direction to which this potential will be exercised is determined by the matching of the in- and out-degrees. In the case that in- and out-degrees are anti-correlated the network will become, so to speak, stuck. In the case that in- and out-degrees are highly correlated, the metapopulation will benefit a high spectral radius of the dispersal pattern and therefore increased resilience. Finally, the random matching case remains at a constant fixed value.

## 5. Summary and discussion

Here, we summarize our main findings, interpret our results, examine a limiting case of regular networks and give a heuristic explanation concerning the impact of IODC. In the above sections, we showed that symmetry of a dispersal pattern has only a very limited impact on metapopulation persistence. We also examined in detail the effect of IODC and in–out degree heterogeneity and specifically showed the dramatic effect they may have on persistence.

A limiting case which is useful to examine is that of regular dispersal networks. By definition, in a regular network of degree *d*, every single node has the same in- and out-degree *d*, and as such the spectral radius is always *λ _{m}* =

*d*. Note that there are many types of regular networks of degree

*d*, from random to lattice-like structures, and they may take on many values of

*γ*. However, they will all have a single value for the spectral radius and the IODC. This example shows that the spectral radius of a regular network, and hence, the persistence of a metapopulation, cannot significantly depend on the symmetry index,

*γ*. Indeed, in [18], regular matrices are explicitly examined by simulations, and no effect of symmetry is found; a result which is in good agreement with our findings.

To deepen our understanding of the results, we ask the question of why should the IODC affect metapopulation persistence? An ecological perspective on metapopulation dynamics gives interesting insights. Focus first on an individual patch, which has just been recolonized. Because the patch is always prone to extinction, it will have a finite lifetime (typically approximately 1/*e*). Within its lifetime, it can recolonize neighbouring patches. If it has a large number of incoming links and a small number of outgoing links, then it will not be able to recolonize many other patches and in a sense the flow is bottlenecked or stopped. Such a node will act as a sink to the dynamics. If the opposite happens, , then the focal node can recolonize many patches though once it is extinct, it will be hard for it to be recolonized by another patch. In this case, the node acts as a source, yet effectively also bottlenecks the flow. Mathematically, is maximal when the in- and out-degrees of nodes are equal, but maximally varied between nodes, as we have shown in §4b. Therefore, for this case, the spectral radius will be maximal. This heuristic interpretation matches the analysis in [6]. Symmetry of the dispersal matrix is a specific case where the in- and out-degrees tend to be equal though it can happen also for fully asymmetric matrices. This explains why symmetry slightly enhances persistence and also shows why it is not the intrinsic property of symmetry in itself, but actually the IODCs combined with high degree heterogeneity.

Besides gaining a better theoretical understanding of the factors contributing to metapopulation persistence, these results may have practical application for conservation action and for terrestrial and marine protected areas' design. For example, designing a network of marine protected areas, so that the in- and out-degree correlations are higher, may lead to better persistence of metapopulations of different marine species.

To summarize, we have presented a hitherto neglected aspect of the question of the effect of dispersal pattern connectivity on metapopulation persistence. By analysing the spectral radius of dispersal matrices and through simulations, we came to the conclusion that symmetry is not a significant factor in determining metapopulation persistence. Rather, the IODC and degree heterogeneity were found to be very important factors in sentencing a metapopulation to sustained persistence or extinction. We also provided a heuristic explanation of why the IODCs have such an influence. The results may have practical significance for environmental policy making and conservation action.

## Funding statement

This research was supported by the Porter School of Environmental Studies at Tel Aviv University. The research was also supported under Australian Research Council's Discovery Projects funding scheme (project number DP150102472).

## Authors' contributions

E.S. and L.S. participated in the design and theoretical foundation of the study, and drafted the manuscript. E.S. additionally carried out analytical mathematical derivations, implemented simulations and numerical procedures.

- Received January 29, 2015.
- Accepted March 10, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.