Optimal dispersal and diffusion-enhanced robustness in two-patch metapopulations: origin’s saddle-source nature matters

A two-patch logistic metapopulation model is investigated both analytically and numerically focusing on the impact of dispersal on population dynamics. First, the dependence of the global dynamics on the stability type of the full extinction equilibrium point is tackled. Then, the behaviour of the total population with respect to the dispersal is studied analytically. Our findings demonstrate that diffusion plays a crucial role in the preservation of both subpopulations and the full metapopulation under the presence of stochastic perturbations. At low diffusion, the origin is a repulsor, causing the orbits to flow nearly parallel to the axes, risking stochastic extinctions. Higher diffusion turns the repeller into a saddle point. Orbits then quickly converge to the saddle’s unstable manifold, reducing extinction chances. This change in the vector field enhances metapopulation robustness. On the other hand, the well-known fact that asymmetric conditions on the patches is beneficial for the total population is further investigated. This phenomenon has been studied in previous works for large enough or small enough values of the dispersal. In this work, we complete the theory for all values of the dispersal. In particular, we derive analytically a formula for the optimal value of the dispersal that maximizes the total population.


Introduction
Metapopulation theory has provided key results into the dynamics of species inhabiting fragmented populations (patches).Since the seminal work by Levins (1969Levins ( , 1970)), metapopulation models have been widely used to investigate the dynamics and persistence of fragmented populations under different scenarios (Clobert et al. 2012;Abbott 2011).The stability of metapopulations relies on the dynamics of the constituent subpopulations and their synchrony, where dispersal plays a significant role by affecting both subpopulation dynamics and such synchrony (Abbott 2011).Experimental evidence has confirmed the role of dispersal in mediating metapopulation stability (Ellner et al. 2011;Bonsall et al. 2002;Dey and Joshi 2011;Smith 2022).While some studies specifically examined the effects of dispersal rates on the dynamics and synchrony in single-species systems (Dey and Joshi 2011;Smith 2022), others focused on local extinctions with varying degrees of linkage between patches hosting multiple interacting species (Fahrig and Merriam 1985;Ellner et al. 2011;Bonsall et al. 2002;Ruiz-Herrera 2018).
Metapopulations are common and widespread in both terrestrial and marine ecosystems, especially for species relying on dispersal to maintain interconnected populations.In terrestrial systems, some plant species in grasslands rely on seed dispersal to colonize new patches, and their populations go through cycles of local extinctions and recolonizations (Honnay et al. 2002).The bog copper butterfly (Lycaena epixanthe) inhabits wetlands and bogs, which are often isolated from one another, undergoing local extinctions and recolonizations within these isolated habitats (Hanski and Ovaskainen 2003).The Glanville fritillary butterfly (Melitaea cinxia) inhabits a mosaic of meadows and patches of suitable habitat in Europe, and its populations experience regular extinctions and recolonizations (Hanski and Gilpin 1994).The spotted owl (Strix occidentalis) in North America exhibits a metapopulation structure across its range, as it relies on different forest patches for nesting and foraging habitats.These forest patches are often separated by unsuitable habitats, leading to a fragmented distribution of the populations (Franklin et al. 2000).More recently, the transient dynamics of a metapopulation of the colonial coastal bird Audouin's gull have been studied under the framework of nonlinear collective dispersal responses to biotic perturbations (Oro et al. 2023).Metapopulations are also common in marine species such as fishes in estuaries and both rocky and coral reefs, seagrass, intertidal invertebrates and coastal decapodes, among others (Kritzer and Sale 2006).
These previous examples, among many others, suggest that metapopulation theory can play a crucial role in delving into the dynamics of spatially-distributed populations having applications for conservation.Metapopulation dynamics have been studied with discrete- (Allen et al. 1993;Dey et al. 2014) and continuous-time (Levin 1974;Pulliam 1988;Sardanyés and Fontich 2010;Sardanyés et al. 2019) models.For instance, the study of two Ricker maps with symmetric and asymmetric coupling showed that dispersal rates stabilized chaotic behaviour (Dey et al. 2014).In a similar direction, Allen and co-workers identified that chaotic dynamics provided robustness under local and global perturbations in coupled subpopulations modelled with logistic and Ricker maps (Allen et al. 1993).Another model considering twopatch discrete models using coupled logistic maps studied the sensitive dependence on initial conditions for the basin of attraction of the periodic orbits, showing that chaos in one patch can be stabilized by dispersal from the other patch (Hastings 1993).Two-patch time-continuous models have been also thoroughly investigated (see, e.g.Fang et al. (2020) and references therein).For instance, using two coupled logistic systems to study the total population for arbitrarily large dispersal rates and the impact of key parameters such as intrinsic growth rates and carrying capacities (Arditi et al. 2015).The same system of coupled logistic models was later inspected by using the so-called balanced dispersal model instead of linear diffusion (Arditi et al. 2016).The exploration of two-patch models with a generic growth function indicated some conditions of optimality also showing that the total population can be higher than the addition of the carrying capacities of each independent patch (Holt 1985).Two-patch models have been also investigated considering local autocatalytic growth (Sardanyés and Fontich 2010).
Despite the intensive research on two-patch metapopulation models, such works often explored only a limited range of dispersal rates and focused on local dynamics.This approach hindered a comprehensive examination of potential interactions between dispersal rates and asymmetry reflected in different local dynamics within the subpopulations.For instance, some studies focused on homogeneous patches with symmetric dispersal (Gonzalez-Andujar and Perry 1993; Gyllenberg et al. 1993;Hastings 1993;Lloyd 1995), while others allowed variations in the parameter determining dynamics between the subpopulations but maintained symmetric dispersal (Kendall and Fox 1998).Conversely, certain studies examined asymmetric dispersal but limited their analysis to cases with identical qualitative dynamics in both subpopulations (Doebeli 1995;Ylikarjula et al. 2000).Other authors studied the interaction between several species in a rock-paper-scissor interaction (Park 2022;Wang et al. 2011).To gain a deeper understanding in metapopulation dynamics, it is necessary to investigate a broader range of local dynamics and dispersal rates in a systematic manner, considering their potential interactions with spatial heterogeneity and asymmetry in dispersal.Moreover, most of these models lack results on global dynamics and have not considered stochastic perturbations, either intrinsic or extrinsic, in the overall metapopulation dynamics.As far as we know, few works have explored the interplay between dispersal, noise and metapopulations' persistence, mainly providing numerical results in discrete-time models (Allen et al. 1993;Sardanyés et al. 2019).
This paper is organized as follows.In "Mathematical model and dynamical aspects" section, we present the model and its adimensionalization of units, which allows to reduce the system from five to three parameters without loss of generality.In subsequent subsections, we discuss the existence of equilibrium points and their stability, bifurcations and global dynamics.Most of the results regarding local dynamics were already known by the community (this is pointed out along the text).The main take-home message of these sections is that global dynamics are affected remarkably by the stability of the origin.We prove that a certain region of the phase space is foliated by heteroclinic connections between the global extinction and the coexistence equilibria if the origin is a source, but, if the origin is a saddle, there is a single heteroclinic connection connecting both equilibria.This phenomenology is further explored in "On the role of the heteroclinic connection" section.In "Role of diffusion in metapopulation robustness to perturbations" section, we provide numerical evidence on how the results identified in "Local and global dynamics" and "On the role of the heteroclinic connection" sections impact on the persistence of the metapopulation under stochastic population fluctuations.In "Optimal dispersal rate" section, we conduct an analytical study of the dependence of the total population with respect to dispersal.This problem has been addressed previously in the literature in the limit case when the dispersal tends to infinity (see (Holt 1985;Arditi et al. 2015) and references therein), and for a small values of the dispersal (Ruiz-Herrera 2018).In the infinite case, the authors showed that having a certain asymmetric conditions on the patches lead to the total population to overcome the sum of the carrying capacities of both patches.In Ruiz-Herrera (2018), a different and less restrictive conditions were shown to be beneficial for small enough values of the dispersal.We here analyse the general case (any value of the dispersal) integrating the previously studied hypotheses in a unified theory.Moreover, we derive an analytic formula to compute the optimal value of the dispersal that maximizes the total population.Finally, "Conclusions" section is devoted to main conclusions.

Mathematical model and dynamical aspects
We introduce the two-patch metapopulation model given by two logistic models coupled by linear diffusion as the simplest way to model within-patch population dynamics and dispersal of individuals between patches.The model is given by the next couple of autonomous ordinary differential equations (ODEs): State variables x i represent the population numbers at patch i, r i the within-patch intrinsic growth rate, k i being the the carrying capacity at patch i. Parameter D denotes the diffusion (dispersal) among patches, which is assumed to be symmetric and follow Fick's law.Let us consider the unit of the population as the carrying capacity of the first patch (dividing all the population variables by k 1 ) and the unit of time such that the rate of birth in the first patch is equal to 1 (dividing r i by r 1 ).This fixes r 1 = k 1 = 1 without loss of generality considering dimensionless variables.We let the parameters r 2 , k 2 and D free for the present study.The role of the asymmetric configurations can be studied by fixing those free parameters bigger or smaller than 1.
In order to avoid cumbersome notation, we rename r 2 as r and k 2 as k.With these modifications, the system reads as (1) with i, j = 1, 2;i ≠ j.
Notice that the parameters k and r are non-dimensional.As we have mentioned in the Introduction, this model has been largely studied.The motivation of this work is to obtain an analytical estimate of the dispersal rate optimizing the size of the population at the subpopulation level thus maximizing global density.Also, we will provide a detailed investigation of the dynamics close to the origin (involving full extinction) considering random fluctuations and seek how dispersal determines the fate of the population with low initial conditions in both patches and in the metapopulation under random fluctuations.
Let us now investigate the dynamics of Eqs.
(2), focusing on the global aspects of the phase space, namely how trajectories of the system connect the origin and the coexistence equilibrium point.These analyses provide an important theoretical support of "Role of diffusion in metapopulation robustness to perturbations" section.Some of the results presented here are known but they are included here for completeness.For example, the existence of a coexistence equilibrium point was given in Freedman and Waltman (1977) for small values of the parameter D. The existence and local stability for any dispersal rate has been proved in several papers, as well as the global stability of the coexistence equilibrium (Angelis et al. 1979;Holt 1985;Arditi et al. 2015).

Local and global dynamics
We first inspect the number of equilibrium points for Eq.(2).Since the vector field consists of two quadratic polynomials with degree two and no independent term, it is trivial to prove that the origin is always an equilibrium point and, by Bezout's theorem, at most, there exist four equilibria.We now prove this elementary result since it will be useful for further analysis that we will develop below.Proof Let us assume, first, that D ≠ 0 (the case D = 0 can be studied separately).The nonlinear system of equations verified by the equilibrium points is: Here, = r∕k .We isolate x 2 with respect to x 1 using the first equation: (2) Now, we substitute x 2 in the second equation, obtaining: Obviously, x 1 = 0 is always a solution of this quartic.Substi- tuting x 1 = 0 in (3), we get x 2 = 0 and we recover the equi- librium point corresponding to the extinction in both patches (the origin).We can factor the later equation as x 1 p(x 1 ) = 0 , where p is a cubic polynomial (which always has a real solution).Two additional real solutions may appear depending on the parameters r, k, and D.
The system for D = 0 has four equilibrium points given by (0, 0), (1, 0), (0, k) and (1, k).◻ As we said, the existence of at least two equilibrium points was given in Freedman and Waltman (1977); Angelis et al. (1979);Holt (1985); Arditi et al. (2015).One of those points is the origin and the other one is located at the interior of the first quadrant.The typical argument to prove the existence of the later is to show that the two nullclines of the system cross in the interior of the first quadrant (we reproduce the argument in Lemma 2.3).
The behaviour of the coexistence equilibrium point when D → ∞ is a recurrent issue being thoroughly discussed in references (Freedman and Waltman 1977;Angelis et al. 1979;Holt 1985;Arditi et al. 2015).These authors were interested in the limit case as a model for perfectly mixing conditions and, more precisely, in the sum of the coordinates of the coexistence point (namely, the numbers at which the metapopulation is stabilized).This aspect of system (2) will be further explored in "Optimal dispersal rate" section.For the moment being, we let the following remark, stating that the limit of the total population can be obtained from the quartic expression given by (4) recovering the results of previous papers in an alternative way.Before analysing the limit cases, numerical results displaying the interior equilibrium points are shown in Fig. 1 together with some orbits in phase portraits.These orbits actually provide information on how they move away from the origin and from the axes.This behaviour, as we discuss below, plays a crucial role in the fate of the species and the entire metapopulation under stochastic perturbations (see "Role of diffusion in metapopulation robustness to perturbations" section).The equilibrium values of the species at each patch are also displayed for different values of D and r in Fig. 2. (4)

Its non-trivial root is
On the other hand, Eq. ( 3) implies that It holds that, In Angelis et al. (1979), the stability of the origin is explored.There, the authors proved that there always exist a repelling direction and that the origin is never a sink.Thus, extinction can never be achieved under the deterministic setting modeled with Eqs.(2).From Eq. ( 3), it is clear that the multiplicity of the origin as an equilibrium point is the same as the multiplicity of x 1 = 0 as a root of the quartic (4).
Remark 2 (Multiplicity of zero) If r = D∕(1 − D) , zero has multiplicity two as solution of the quartic Eq. ( 4).If, besides, k = 1∕r 2 , it has multiplicity three.Proof The Jacobian matrix of the vector field evaluated at the point (0, 0) is given by: An elementary computation shows that the eigenvalues are: Notice that + is always positive.On the other hand, − is negative for r > r * ∶= D∕(1 − D) and positive for r < r * .Therefore, the origin changes from repelling node to saddle as r crosses the curve r * from above (see Fig. 1).The type of bifurcation is determined by the multiplicity of the origin as equilibrium point (see Remark 2): If the zero has multiplicity two, the bifurcation is transcritical while if it has multiplicity three, it is, generically, a pitchfork.◻ The following result establishes that the nullclines enclose a region of the phase space which is positively invariant by the flow.This is sketched in Fig. 3.

Lemma 2.3 (Intersection of the nullclines) For any positive values of r, k and D, the horizontal and vertical nullclines define a positively invariant region.
Proof The vertical nullcline, d dt x 1 = 0 , is described by the parabola: The horizontal nullcline, d dt x 2 = 0 , is given by: Notice that the vertical nullcline crosses the horizontal axis at N 1 = ((1 − D), 0) while the horizontal nullcline crosses the vertical axis at N 2 = (0, (r − D)∕ ) .It holds that, and The two parabolas intersect in a point contained in the first quadrant.Indeed, if r > D (see Fig. 3, left), at the vertical axis, the horizontal nullcline g h is above the vertical one g v .Moreover, it behaves as √ x 1 when x 1 is large.On the other hand, the vertical nullcline behaves as x 2 1 .A direct application of the Bolzano theorem shows that there exists a crossing point for x 1 > 0 and x 2 > 0 .If r < D (see Fig. 3, right), the horizontal nullcline crosses the vertical axis for x 2 = 0 and for a negative value of x 2 .At the origin, the slope of the horizontal nullcline is positive and the slope of the vertical nullcline is negative.Therefore, near zero, the horizontal nullcline is above the vertical one and the same argument can be used to establish the existence of the crossing point at the first quadrant.
From the fact that the two nullclines always cross, we get a region that separates the phase space.Let us see that this region is positively invariant.If r > D the first compo- nent of the vector field is positive along the vertical axis for values of 0 < x 2 < (r − D)∕ while the second component  2) in the parameter space (r, D) for x 1 (left) and x 2 (right) computed numerically with k = 0.5 and x 1 (0) = 0.45 and x 2 (0) = 0.2 .Overlapped we display the curve separating the source from the saddle of the origin shown in Fig. 1 Fig. 3 Sketch of the positive invariant region (blue) defined by the vertical and horizontal nullclines.Left: r = 1 > D .Right: r = 0.3 < D .In both cases D = 0.4 , k = 0.4 .Observe that the blue region is surrounded by the isoclines and their direction points inwards of the vector field is positive along the horizontal axis for x 1 < (1 − D) .This, together with the fact that the parabo- las are nullclines establishes the positive invariance of the region.If r < D the result follows from the second compo- nent of the vector field being positive along the horizontal axis for x 1 < (1 − D) .
◻ Remark 3 (Positively invariant region) The first quadrant is also positively invariant.The region enclosed by the nullclines is not contained in it (see Fig. 3).For the next result, we shall focus on the intersection between the region enclosed by the isoclines and the first quadrant.This intersection, that we name as I + , is also positively invariant.
Remark 4 (Unstable eigenvector of the origin) The eigenvector associated to the + eigenvalue always points towards I + .Indeed, the eigenvector associated to + is given by: Observe that the first component of v + is always positive, therefore, pointing to the first quadrant.If r > D , the iso- clines cross the axis for positive values.Therefore, always point to I + .If r < D , then the slope of the eigenvector must be between the slopes of g v and g h : Both inequalities are true while Theorem 2.4 (Global dynamics) The following sentences hold for system (2): 1.There are no periodic orbits.
2. It has an equilibrium point ( P 3 ) inside the first quadrant that is globally asymptotically stable (GAS) in the first quadrant.3.If D < r∕(1 + r) , there exist infinitely many heteroclinic connections between the origin ( P 0 ) and P 3 .4. If D > r∕(1 + r) there is a unique heteroclinic connec- tion between P 0 and P 3 .
Proof We shall prove the statements by order.
1. We have shown that the two nullclines cross at a point in the first quadrant.A periodic orbit in the first quadrant should encircle the crossing points of the vertical and horizontal nullclines.But this would contradict the fact that they are positively invariant.
2. We discuss first the case D = 0 .It is elementary to show that it has four equilibrium points, (0, 0) being a repelling node, (1, 0), (0, k) being saddle points; and (1, k) being an attracting node.The lines x 1 = 0 and x 1 = 1 are vertical nullclines and, similarly, x 2 = 0 and x 2 = 1 are horizontal nullclines.Therefore, the rectangle defined by the four equilibrium points is positively invariant.This means that there are no periodic orbits around none of the equilibria and, therefore, the -limit for a full measure set of initial conditions is (1, k).For the case D ≠ 0 we argue that, because of the positively invariant region defined by the nullclines, the coexistence point has to be an attractor.Finally, since there are no periodic orbits (in the first quadrant), it is a global attractor.3.If P 0 is a source, any initial condition on I + with -limit P 0 must have P 3 , which is GAS, as -limit.4. If P 0 is a saddle, the unstable manifold has P 3 as -limit.
Any other initial condition in I + has its -limit outside the first quadrant.

On the role of the heteroclinic connection
In this section, we study the fact that the existence of a unique heteroclinic connection affects the global dynamics.This is done by studying the attracting character of the connection.Lemma 2.2 establishes that the stability type of the origin changes from a repelling node (a source) to a saddle point (as D crosses a certain critical value depending on r).This phenomenon has a significant impact on the dynamics near the origin, and as we will see below, will play a central role in the robustness of the metapopulation under stochastic fluctuations.Indeed, by Hartman's theorem, the linearized system provides an approximation on the trajectories nearby.When the origin is a source, every initial condition in the positively invariant region lies on a heteroclinic connection between the origin and the coexistence equilibrium point.Additionally, two initial conditions that are very close to the origin tend to separate from each other exponentially over time.Conversely, when the origin is a saddle point, there is a single locally attracting heteroclinic connection.In this case, nearby solutions to the origin converge rapidly to its unstable manifold and then follow the heteroclinic connection towards the coexistence equilibrium point.
To gain a better understanding of this phenomenon, we restrict ourselves to the symmetric case with r = 1 and k = 1 , which albeit being the simplest case gathers all the elements playing a role in this scenario.When r = 1 , the bifurcation of the origin occurs at D = 1∕2 .At this value of r, the eigenval- ues of the Jacobian matrix evaluated at the origin are and the corresponding eigenvectors are given by Consider the linear change of variables given by the eigenvectors.
Notice that the new variables y 1 and y 2 represent the differ- ence of population between the patches and the total population, respectively.Negative values of y 1 represent initial conditions with more individuals in the second patch.Conversely, positive values mean that the first patch is more populated.These variables where used in Holt (1985) for different purposes.The change of variables casts the original vector field (with r = 1 ) to: This vector field will be used later to study how the heteroclinic connection depends on the normalized carrying capacity k.For the time being, we focus on the symmetric case.If we chose k = 1 , then, some terms vanish and the resulting ODE is Let us analyse this system more carefully.In first place, we notice that it has only ecological meaning for trajectories within the region y 2 > 0 and −y 1 ≤ y 2 ≤ y 1 .That is, when the total population is positive, and the difference of the subpopulations is less than the total one.It is evident that this region is positively invariant by the flow.In these coordinates, the four equilibria are given by The Jacobian matrix in these coordinates reads P 0 = (0, 0), It is easy to see that the determinant of the matrix for both P 1 and P 2 is 4D 2 − 1 which is negative for D < 1∕2 .This means that these equilibria are always saddle points whenever they exist.The points P 1 and P 2 lie outside the region, and therefore are not ecologically meaningful.When D = 1∕2 , P 1 and P 2 merge with the origin in a pitchfork bifurcation.
The vertical axis, {y 1 = 0} , is invariant and it is, for each value of D, an heteroclinic connection between P 0 and P 3 .Moreover, it is the only heteroclinic connection for D > 1∕2 .Notice that {y 1 = 0} corresponds to the line {x 1 = x 2 } in the original coordinates.The first equation of system (7) captures the horizontal flow of the system.Let us rewrite it as: When 0 < D < 1∕2 and y 2 < (1 − 2D) the horizontal flow is positive for y 1 > 0 and negative for y 1 < 0 , meaning that the vertical axis is repelling for y 2 < (1 − 2D) .This is sketched in Fig. 4. If D > 1∕2 , the vertical axis is always attracting.
This global property of the phase space is relevant in the capacity of the system to exploit the population of one patch to recover population of the other one, at least, at a short time scale.In the region of the phase space in which the vertical axis is repelling, the trajectories of the system that start close to the lines {y 2 = ±y 1 } remain close to them (see Fig. 4).This suggests that the recovering of a patch is much slower when the origin is a source.A natural conjecture, to be tested in "Role of diffusion in metapopulation robustness Fig. 4 Sketch of the phase space of System (7).The red lines delimit the admissible region (i.e. points outside the admissible region correspond to points of the original system with some negative coordinate).The vertical axis is an heteroclinic connection between P 0 and P 3 , and it is repelling if y 2 < 2D − 1 .The dashed arrows determine this attracting character of the heteroclinic connection.See "On the role of the heteroclinic connection" section for more details to perturbations" section, is that having a saddle is favourable for recovery if one of the patches is populated by few individuals, thus providing more robustness to survival to the metapopulation.The area of the subregion in which the vertical axis is attracting is given by 2D(1 − 2D) and decreases linearly with D until it gets zero at D = 1∕2.
The dependence of the unstable manifold of the origin on k starts at second order.We can study the second term of the Taylor expansion of the manifold by applying the parameterization method (Haro et al. 2016).That is, we consider system (6) and consider a parameterization U of the manifold as Where s is a real parameter and a and b are constants to be determined.These constants can be obtained by imposing invariance of the manifold up to second order.That is, we select a and b so the following equation is fulfilled: Here, Y stands for the vector field of system (6) and = 1 is the unstable eigenvalue related to the origin.Solving the latter equation for a and b results in Looking at the sign of a we can understand how the manifold bends (at second order) when k ≠ 1 .As expected, if k < 1 the manifold bends towards the region y 1 > 0 , where the population of the first patch is larger than the population of the second patch.The situation is opposite if k > 1.
We have seen that when the origin is a source the line x 1 = x 2 is locally repelling and gets attracting as the origin becomes a saddle.To do so, we have used adapted coordinates ( y 1 -y 2 ).This allows us to conjecture that having a saddle benefits the recovering capacity of the system.This fact is studied in next section introducing random population fluctuations.

Role of diffusion in metapopulation robustness to perturbations
In "On the role of the heteroclinic connection" section, we analysed the local dynamics close to the origin and concluded that if D is chosen such that the origin is a saddletype equilibrium point, the heteroclinic connection is locally attracting.In this section, we provide numerical evidence showing how this can be beneficial to the persistence of the metapopualtion under the effect of stochastic perturbations.Specifically, we conjecture that a homoclinic connection that attracts initial conditions nearby may have a positive effect on the persistence of the metapopulation when one of the two subpopulations is close to extinction.
To investigate this phenomenon, we designed a simple numerical experiment simulating random deaths due to some external agent (e.g.sick animals, predators, climatic factors, hunting).The process is as follows: 1. We consider system (7) and fix a maximal integration time of T = 10 system units.This value has been selected experimentally to ensure effective stabilization for all the initial conditions in the unperturbed system.
Here, effective stabilization means that the initial condition is closer to P 3 than 10 −3 in 2 -norm.Notice that, in these symmetric conditions, P 3 does not depend on D. 2. We select a grid of 500 × 250 initial conditions in the region 0 < y 2 < 4∕5 and −y 2 < y 1 < y 2 .Note that 0 < y 2 < 4∕5 are the values of y 2 for which the vertical axis is repelling if we pick D = 0.1.3.For each initial condition, we perform an integration and, at random intervals following a uniform distribution U(0, 0.1) , we subtract a random quantity from the total population y 2 .This quantity also follows a uniform dis- tribution U(0, 1∕50) , i.e. the maximum subtracted value is the total carrying capacity divided by 100. 4. If this random perturbation is large enough to cause y 2 to become negative, the integration is terminated, and we consider the population to be extinguished.If y 1 > y 2 , i.e. x 2 < 0 in the original coordinates, we then rear- range the initial condition to y 1 = y 2 , and the integration continues.A similar procedure is done if y 2 < −y 1 , i.e.
x 1 < 0 .When the integration is completed, any initial conditions that result in one of the subpopulations being smaller than 1/100 are labelled as being at risk of extinction in patch 1 (or 2).
We perform this experiment (steps 1-4) N times for each initial condition, associating an extinction probability with each initial condition.The extinction probability associated to each initial condition is computed by dividing the number of simulations that result to extinction by the total number of simulations N (we have chosen N = 100).
In Fig. 5a, we display the direction of the vector field of system (2) for D = 0.1 and D = 0.9 (under symmetric condi- tions, k = r = 1 ).In each of the pictures there appear three trajectories, the heteroclinic connection {x 1 = x 2 } and the trajectories starting at (10 −2 , 0) and (0, 10 −2 ) , respectively.When the origin is a source ( D = 0.1 ) these trajectories spend some time close to the axes.When D = 0.9 , both tra- jectories are rapidly attracted by the heteroclinic connection (Fig. 5a).These two plots provide an illustration of the main point of "On the role of the heteroclinic connection" section, namely the system is better at recovering a patch which is close to extinction if the origin is a saddle.Figure 5b displays some representative time series with local and global extinctions together with survival dynamics at increasing diffusion simulating the random perturbations described above.The first row displays the dynamics of the population at patch 2 ( D = 0.1 first column, D = 0.9 second column).In the upper panel for D = 0.1 , the local population goes to extinction, but recovery is still possible from the other patch although the population remains at extremely low values.The lower panel for this diffusion value displays another simulation with full extinction of the metapopulation.In the case D = 0.9 , the attracting charac- ter of the heteroclinic connection permits the total population to exit the extinction zone rapidly and, therefore, the population persists.The upper plot clearly shows this effect and the population at patch 2 rapidly increases.The lower panel displays another simulation for the entire metapopulation, which already has large population values at initial times while it progressively grows towards larger values.We notice here that the population continues increasing after time = 10 (results not shown).
Finally, panel (c) displays both local ( P 1,2 ) and global (total, P T ) extinction probabilities for D = 0 (first row), D = 0.1 (second row) and D = 0.9 (third row) in the space of initial conditions (x 1 (0), x 2 (0)) .The case D = 0 is included for completeness.We notice that extinction probabilities shall be interpreted differently for local and total populations.Total extinction involves that both time series x 1 and x 2 abandon the first quadrant simultaneously and the integra- tion is terminated.extinction probability for subpopulations refers to the event the number of individuals at patch i = 1, 2 is below 1/100 after 10 units of time.In this case, the ecological interpretation is: if the subpopulation starts in the danger zone, then the system is not able to take it out from it with a certain probability.We notice that, for D = 0.1 the probability of subpopulation extinction is close to one along the axis when the line {x 1 = x 2 } is not attracting.In the case D = 0.9 the system is always able to recover a patch in danger unless the total population becomes extinct.These analyses also show that at low diffusion values, both local and global extinctions occur for wide ranges of initial conditions.A slight increase in diffusion (from D = 0 to D = 0.1 ) decreases local extinctions although global ones are still found.The simulations done with D = 0.9 indicate that local extinctions only take place for smaller population sizes and global extinctions are also restricted to low population values.As expected, diffusion ensures persistence of the populations at a local level keeping the entire metapopulation in a safe state.

Optimal dispersal rate
A well-studied phenomenon from metapopulation theory is the fact that, under suitable conditions, the total population can exceed the sum of the carrying capacities of the patches (Holt 1985).This is known to happen when the conditions in one of the patches are better than the conditions in the other one.When this hypothesis is fulfilled, the less advantageous patch acts as a source of population.
For fixed values of any pair of values r and k, system (2) has always a coexistence equilibrium point P 3 (D) that is globally asymptotically attracting.Therefore, for any selection of the parameters of the system and any initial condition, the population tends to stabilize at P 3 (D) .That is, for any threshold , any values r, k D and any initial condition (x 1 (0), x 2 (0)) , there exist a time t for which Here, D,r,k t denotes the flow of system (2).Therefore, the behaviour of the coordinates of P 3 (D) with respect to D determines the long-term behaviour of the total population ( x 1 + x 2 ) of the system for any initial condition.In Ruiz- Herrera and Torres ( 2018), the authors define the function They are able to prove the following result.Henceforth, the function Ω is increasing nearby D = 0 if r > 1 and k > 1 or if r < 1 and k < 1 .As Ω(0) > 1 + k , under these conditions, Ω can exceed the sum of the carry- ing capacities of the patches.The behaviour of Ω has been analysed previously for arbitrarily large diffusion values (namely D → ∞ ).For instance, results in Holt (1985); Arditi et al. (2015) show that dispersal is beneficial if k > 1 and r∕k > 1 or k < 1 and r∕k < 1 .These conditions are usually refereed as positive (negative) r − k relation.As it is stated in Ruiz-Herrera and Torres (2018), this last condition is more restrictive than r > 1 and k > 1 or if r < 1 and k < 1 .In Arditi  et al. (2015), the authors derive the following formula: Notice that formula (5) from Remark 1 is a compact version of this one.In particular, using (5) it can be shown that Ω(∞) > 1 + k whenever k is contained between 1 and r (without assuming r > 1 or r < 1).
The qualitative behaviour of the function Ω with respect to D is also analysed in Arditi et al. (2015).It is shown that, for some values of the parameters the function is strictly increasing while for some other, it has a maximum.Asking which dispersal rate maximizes Ω is a natural ques- tion.In this section, we provide an answer and, moreover, recover some insights on the qualitative behaviour of the function Ω.
The first step to tackle the problem is to understand the locus of the equilibrium points of system (2).Particularly, while the system is two dimensional, the equilibrium points are located, for fixed values of k and r, in a closed 1-dimensional manifold given by an ellipse, as the following result shows.Proof We consider once again the system of equations for the equilibria: Adding both equations, we get the following expression.
which is equivalent to: This is, in fact, the equation of an ellipse centred at (1/2, k/2) and with semi-axes √ a and √ b .◻ The previous argument only states that, if (x 1 , x 2 ) is an equilibrium point, then it is contained in the ellipse.The converse is tackled in the following Lemma for some points of the ellipse.
Lemma 4.3 (The points of the ellipse are equilibria) If (x 1 , x 2 ) verify Eq. (8), the following identity holds: Moreover, the point (x 1 , x 2 ) is an equilibrium point of sys- tem (2) for Proof Notice that the second statement is trivial if identity (9) holds.Before proving (9), we observe that, if a and b are defined as in the statement of Lemma 4.2, then: Notice also that, by square completion, If (x 1 , x 2 ) verify Eq. ( 8), then the above quantity is equal to The identity is proven by noticing that and recalling that k = r .Indeed,

◻
We stress the fact that identity (9) does not provide a one-to-one map from the ellipse to the set of equilibria of system (2).For instance, the origin always belongs to both the ellipse and the set of equilibria but identity 9 does not make any sense if x 1 = x 2 = 0 .On the other hand, if x 2 = x 1 , the value of D cannot be recovered from iden- tity (9).For the coexistence equilibrium point, this happens, for instance, in the case r = 1 , k = 1 .In this case, however, the coexistence equilibrium does not depend on D (see "On the role of the heteroclinic connection" section).
To obtain the best dispersal rate, i.e. the one that maximizes the total population, one has to look for the maximum of the function Ω on the ellipse (8).Notice that as the objec- tive function is continuous and the set of feasible solutions a compact, this problem has always a solution.solution may correspond to a negative value of the dispersal rate D. Interestingly enough, the singularity of identity (9) implies that, for some choices of r and k, the population is maximized at the limit D → ∞ .The next result provides the solution of the aforementioned optimization problem, and it is tackled using Lagrangian Multipliers.
Proof In view of Lemma 4.2, the maximal total population can be obtained by maximizing the function Ω on the ellipse (8).This is equivalent to solve the following optimisation problem: Here, b − 1 the restriction func- tion.Let us maximize the Lagrangian function.
for some real parameter .If we look for the roots of ∇L we notice that those are: Since we are looking for a maximum of f, it is clear that we must use the solutions with positive sign (the negatives, which have no ecological meaning, correspond to the global minimum).This optimal point being an equilibrium point comes as a straightforward application of Lemma 4.3.The maximal value of f is given by which corresponds to Ω * if we expand a and b in terms of k and r.Moreover, if we recast (X 1 , X 2 ) to the original coordi- nates, we can use Lemma 4.3 to find the value of D for which the optimum is an equilibrium point: with There are several consequences of Theorem 4.4.Notice first that, if r = 1 , then Ω(D) = 1 + k .This means that, if the two patches have the same growth rate, the total population is maximized at the sum of the carrying capacities.This is a known fact that can be recovered from Eq. ( 10).Moreover, the linear coefficient of the quadratic polynomial inside the square root is which is a function whose minimal value is 2 and it is achieved precisely at r = 1 .That is, Hence, Ω * is always superior (or equal) to the sum of the carrying capacities of the patches.As we stated before, this does not mean that the optimal total population is achieved for positive D. To check when the optimal population is achieved for a positive value of the dispersion rate we have to examine the sign of Eq. ( 11).The parameter space (r, k) can be partitioned according to the behaviour of the total population with respect to the dispersion rate (see Fig. 6).
The following result provides this classification.
Corollary 4.5 (Behaviour of the total population) Ω * is achieved for a positive dispersion rate D * if one of the fol- lowing conditions hold: where Proof From the proof of Theorem 4.4 we know that Ω has only one critical value that is an equilibrium point inside the first quadrant.We also know that the maximal value is achieved at: where, Fig. 6 Behaviour of the total population Ω in the parameter space (r, k).Green areas denote those values of the parameters for which the total population is achieved for a positive value of the dispersal rate.Grey ones are areas where r-k positiveness is fulfilled.Note that green and grey areas overlap.Grey area which does not overlap with green area are those values for which the total population increases monotonically for positive dispersal rate but the optimal value is not achieved.Red regions denote those values of r and k for which the total population decreases for positive dispersal rates.See Corollary 4.5 for more details.
Boxes with letters from a to d correspond to the values of r and k used in Fig. 7 where the total population is shown at increasing values of D Of course, D * is positive whenever the numerator and the denominator have equal sign.The numerator has positive sign if and only if r > 1 .On the other hand, the denominator has positive sign if and only if k > r(3 + r)∕(1 + 3r) .This proves points 1 and 2. If r < 1 and k * (r) < k < 1 , the optimal value is achieved for some negative value of D. However, Ω(D) is increasing near D = 0 (see Lemma 4.1).As there are no other relative optimal values of D, Ω(D) must be an increasing function of D for D > 0 .◻ Remark 5 For k = r(3 + r)∕(1 + 3r) , the optimal value is achieved in perfectly mixing conditions, that is, D * → ∞.
The information provided by Corollary 4.5 is shown in Fig. 6.The horizontal axis represents the value of r and, the vertical, the value of k.Different regions are coloured according to the behaviour of the function Ω with respect to the dispersal rate.The red zone represents values of the parameters r and k for which having positive D is always detrimental in terms of total population.This is because the function Ω is monotonically decreasing with respect to D. Zones which are not coloured in red corresponded to the conditions appearing in Ruiz-Herrera (2018).Therefore, in non-red zones Ω is locally increasing for D near to zero.We coloured in grey regions in which r − k positive- ness condition is fulfilled.In this region, the population at the limit case D → ∞ is larger than the sum of the carry- ing capacities of the patches (see, for instance, Freedman and Waltman (1977)).The line k = r determines the limit case Ω(∞) < 1 + k if k > r and r > 1 or k < r and r < 1 .The region determined by the curves k = r(3 + r)∕(1 + 3r) and k = r correspond to values for which Ω * > 1 + k but Ω(∞) < 1 + k .The green regions represent values of the parameters for which Ω * is achieved for some positive D. Notice that the green region and the grey region overlap (this overlapping is seen as a dark green).For values of the parameters in this overlapped region, D * is larger than the limit case.In the grey region, the function Ω is strictly increasing for positive D as the maximal value is achieved for some negative D.
Figure 7 displays the evolution of Ω as function of D for four different choices of the parameters r and k, labelled from a to d in Fig. 6.The plots have been obtained by implementing a pseudo arc-length continuation to the equation of the equilibrium points starting at the point (1, k) for D = 0 .The program is terminated when the characteristic curve crosses the homotopy level {D = 10} .This program has been also used to check the correctness of the formulas obtained in this section.Figure 7a has been obtained with r = 1.5 and k = 1.5 .At D = 0 the value of Ω is 1 + k = 2.5 .
The characteristic curve increases until some maximal value that can be computed using formula (10) (in this case Ω * ≈ 2.5247 ).After the maximal value is reached, the curve decreases monotonically until a horizontal asymptote.This asymptote can be computed by using the asymptotic formula of Remark 1.In this case Ω decreases asymptotically to the value 2.5 (i.e. the sum of carrying capacities: 1 + k ).
In Fig. 7b we set r = 0.5 and k = 1.5 .These are conditions where the theory predicts that positive dispersal rate is detrimental to the total population.Indeed, the function Ω , in this case, is strictly decreasing until the asymptotic limit Ω = 2.25 .Notice that Ω * is achieved for some D < 0 .In Fig. 7c, we have selected r = 1.5 and k = k * (1.5) , where k * is taken as in the statement of Corollary 4.5.As it can be seen in the plot, the theory predicts that the maximal value of Ω is achieved at the limit case and, therefore, the function Ω is strictly increasing for positive values of D. The asymptotic value Ω is, in this case, 2.25.This value can be predicted with both Theorem 4.4 and Remark 1. Finally, Fig. 7d illustrates a case in which the limit of Ω is below the sum of the carrying capacities.The values of the parameters are r = 1.5 and k = 1.7 .The maximal total population is slightly larger ( ≈ 2.72 ), and it is achieved for some positive D.

Conclusions
The investigation of metapopulation mathematical models is very extensive in the literature.Most of these studies have focused on small metapopulations considering few patches (Hastings 1993;Arditi et al. 2015Arditi et al. , 2016;;Sardanyés and Fontich 2010), or in multi-patch systems (Allen et al. 1993).Discrete-time metapopulation models have focused on inspecting the role of dispersal in the stability of chaotic dynamics and its role in metapopulations' persistence (Hastings 1993;Allen et al. 1993).For instance, Allen and coworkers showed that local chaotic dynamics involved lower extinction probabilities under both intrinsic and global noise (Allen et al. 1993).Additionally, two-patch time-continuous models have been extensively investigated.Arditi and colleagues used two coupled logistic systems to study the total population under arbitrarily large dispersal rates (Arditi et al. 2015).Later on they explored the same system of coupled logistic models using the balanced dispersal model instead of linear diffusion (Arditi et al. 2016).Moreover, the exploration of two-patch models with a generic growth function revealed certain conditions of optimality, indicating that the total population can surpass the sum of carrying capacities of each independent patch (Holt 1985).More recently, a similar system was explored for local populations growing hyperbolically instead of exponentially (Sardanyés and Fontich 2010).
In this manuscript, we have revisited one of the simplest metapopulation models, the one in which a diagonal connectivity matrix couples two logistic differential equations.The investigation of metapopulations with few patches allows developing analytical studies, providing clear information on dynamical phenomena such as equilibrium points and bifurcations that could be conserved in higher dimensions.Unlike the previous works on two-patch, time-continuous metapopulations, we have focused on global aspects of the system.Our work provides two novel main contributions to metapopulation theory.
The first contribution has to do with the robustness of the system to external perturbations.That is, the capacity of the system to recover from a situation in which one of the patches is close to extinction.It has been known for decades that, for all the values of the parameters, an empty patch will be eventually populated by individuals from the other patch.In this process, the populations of both patch tend to stabilize at some coexistence equilibrium point.However, references in literature do not take into account that the stability type of the global extinction equilibrium point (namely, the origin) can play a relevant role in the fragility of the metapopulation in terms of extinction.Concerning this point, we have identified two different scenarios: (i) When the origin is a source, the trajectories with initial conditions close to the axes (i.e.situations in which one of the patches is close to extinction) need some time to separate from them; (ii) when the origin is a saddle, there is a unique, locally attracting, heteroclinic connection between the origin and the coexistence equilibrium that leads low initial condition to achieve a safe state in a short period of time.By including noise in the dynamics, we have shown that recovery is more prone when the origin is a saddle point.Notice that, at first sight, the origin being a source could seem a better situation to prevent from under perturbations.To illustrate this fact we have run a simulation in which the system has stochastic losses of individuals in both patches.
The second contribution is related to a well-known counter-intuitive property of the system: Under suitable conditions, having dispersal rate, leads the total population to stabilize overcoming the sum of the carrying capacities of the patches.This phenomenon was established for the limit case under the k − r positiveness hypothesis (Freedman and Waltman 1977;Holt 1985;Arditi et al. 2015).That is, if we name r 1 and k 1 the growth rate and the carrying capacity of the first patch, and r 2 , k 2 the respec- tive quantities related to the second patch, k − r positiveness means that k 1 > k 2 and r 1 k 2 > r 2 k 1 or k 2 > k 1 and r 2 k 1 > r 1 k 2 .
the other hand, a less restrictive hypothesis was identified in (Ruiz-Herrera and Torres 2018).If k 1 > k 2 and r 1 > r 2 (or k 1 < k 2 and r 1 < r 2 ) and D is small enough, the total popula- tion exceeds also the sum of the carrying capacities.
In this work, we study analytically the problem for all the positive values of the dispersal rate.In particular, we derive a formula for the optimal dispersal rate and a bound (that is fulfilled for the optimal dispersal rate) to the total population.To avoid cumbersome notation, we have used dimensionless units, but we stress that the change of units can be pulled-back to recover our formulas for standard units.The bound to the total population reads as: The main takeaway from the "Optimal dispersal rate" section is a unified framework in which both r-k positiveness and the less restrictive conditions appearing in Ruiz-Herrera (2018) and all positive values of the dispersal rate are included.In particular, there are regions in which the r-k positiveness is fulfilled but the maximal population is achieved for some finite value of D. We derive formulae that can be used to predict both the optimal dispersal rate, the maximal population and the asymptotic behaviour of the function Ω.
The simplicity of the model studied in this work has allowed us to conduct an analytic description of global aspects of the system (namely, global dynamics and behaviour of the total population as a function of dispersal).We remark, however, that the phenomena investigated in this paper, e.g. that changes in local behaviour can impact on global dynamics, may appear in other complex systems.For instance, such results found in low-dimensional metapopulation models may be found in higher dimensions.

Technical details
The programs corresponding to "Role of diffusion in metapopulation robustness to perturbations" and "Optimal dispersal rate" sections have been written in C from the scratch and are available upon request.Library LAPACK (Anderson et al. 1999) has been used for linear algebra and the Taylor package (Jorba and Zou 2005)) has been used to perform the numerical integrations of system (7).

Lemma 2. 1 (
Number of equilibrium points) System(2) has two, three or four equilibrium points.

Lemma 2. 2 (
Bifurcation of the origin) The origin is a repelling node for r < D∕(1 − D) and a saddle for r > D∕(1 − D).

Fig. 1
Fig. 1 Stability of the origin for Eq.(2) in the plane (r, D), where the curve D = r∕(1 + r) separates the source from the saddle behaviour.Several phase portraits are shown at increasing values of diffusion for r = k = 1 and: a D = 0 , b D = 0.1 , and c D = 1 .The arrows indicate

Fig. 2
Fig. 2 Interior equilibrium of Eqs.(2) in the parameter space (r, D) for x 1 (left) and x 2 (right) computed numerically with k = 0.5 and x 1 (0) = 0.45 and x 2 (0) = 0.2 .Overlapped we display the curve separating the source from the saddle of the origin shown in Fig.1

Fig. 5 a
Fig. 5 a Vector field (the colours of the arrows denote the norm of the field) of the system (2) for D = 0.1 (left) and D = 0.9 (right).b Time series with local ( x 2 ) and global populations for D = 0.1 and D = 0.9 .c Extinction probabilities ( P 1,2 and P T denote local and global extinctions, respectively) for initial conditions close to the Lemma 4.1 (Ruiz-Herrera and Torres 2018)

Lemma 4. 2 (
Equilibria lie on an ellipse) The equilibrium points of system (2) are contained on an ellipse with centre (1/2, k/2) and axes √ a and √ b , where: