# Optimal dispersal and its role in avoiding stochastic extinctions in a two-
patch metapopulation model

Marc Jorba-Cuscó(1) Ruth I. Oliva-Zúniga(2) Josep Sardanyés(1) Daniel Pérez-
Palau(3)

(August 2, 2023)

# Optimal dispersal and its role in providing resilience in two-patch
metapopulations

Marc Jorba-Cuscó(1) Ruth I. Oliva-Zúniga(2) Josep Sardanyés(1) Daniel Pérez-
Palau(3)

(August 2, 2023)

# About the role of dispersal in two-patch metapopulations resilience

Marc Jorba-Cuscó(1) Ruth I. Oliva-Zúniga(2) Josep Sardanyés(1) Daniel Pérez-
Palau(3)

(August 2, 2023)

# Optimal dispersal and its role on the resilience
of two-patch metapopulations

Marc Jorba-Cuscó(1) Ruth I. Oliva-Zúniga(2) Josep Sardanyés(1) Daniel Pérez-
Palau(3)

(August 2, 2023)

# Optimal dispersal and its role in two-patch metapopulations
resilience in noisy environments

Marc Jorba-Cuscó(1) Ruth I. Oliva-Zúniga(2) Josep Sardanyés(1) Daniel Pérez-
Palau(3)

(August 2, 2023)

# Dispersal-driven mechanisms of resilience to disturbances
in two-patch metapopulations

Marc Jorba-Cuscó(1) Ruth I. Oliva-Zúniga(2) Josep Sardanyés(1) Daniel Pérez-
Palau(3)

(August 2, 2023)

# Dispersal-enhanced resilience in two-patch
metapopulations: origin’s instability type matters

Marc Jorba-Cuscó(1) Ruth I. Oliva-Zúniga(2) Josep Sardanyés(1) Daniel Pérez-
Palau(3)

(August 2, 2023)

###### Abstract

Many populations of birds, reptiles, mammals, insects, or plants, exhibit a
metapopulation structure with close, spatially-separated subpopulations. The
field of metapopulation theory has made significant advancements since the
influential Levins model. Various modeling approaches, considering diverse
ecological interactions and numbers of patches, have provided valuable
insights to theoretical Ecology. Despite extensive research on metapopulation
models, there are still challenging questions that are difficult to answer
from ecological metapopulational data or multi-patch models. Low-dimension
mathematical models offer a promising avenue to address these questions and
provide analytical information, especially for global dynamics which have been
scarcely investigated. In this study, we investigate a two-patch
metapopulation model with logistic growth and diffusion between patches. By
using analytical and numerical methods, we thoroughly analyze the impact of
diffusion on the dynamics of the metapopulation. We identify the equilibrium
points and assess their local and global stability. Furthermore, we
analytically derive the optimal diffusion rate that leads to the highest
metapopulation values. Our findings demonstrate that increased diffusion plays
a crucial role in the preservation of both subpopulations and the full
metapopulation, especially under the presence of stochastic perturbations.
Specifically, at low diffusion values, the origin is a repeller, causing
orbits starting around it to travel closely parallel to the axes. This
configuration makes the metapopulation less resilient and thus more
susceptible to local and global extinctions. However, as diffusion increases,
the repeller transitions to a saddle point, and orbits starting near the
origin rapidly converge to the unstable manifold of the saddle. This
phenomenon reduces the likelihood of stochastic extinctions and the
metapopulation becomes more resilient due to these changes in the vector field
of the phase space.

  * (1)

Centre de Recerca Matemàtica. Edifici C, Campus de Bellaterra, 08193
Cerdanyola del Vallès, Spain<EMAIL_ADDRESS><EMAIL_ADDRESS>
  * (2)

Universidad Nacional Autónoma de Honduras en el Valle de Sula. UNAH-VS,
Boulevard UNAH-VS 21102 San Pedro Sula, Honduras<EMAIL_ADDRESS>
  * (3)

Escuela Superior de Ingeniería y Tecnología, Universidad Internacional de la
Rioja, Av. La Paz 137, 26006 Logroño, Spain<EMAIL_ADDRESS>
Keywords: Dynamical systems; Bifurcations, Metapopulation resilience,
Theoretical Ecology, Stochastic extinctions

##  1 Introduction

The metapopulation concept, originally formulated by Richard Levins [Levins,
1969, Levins, 1970], has historical foundations dating back to Andrewartha,
MacArthur and Wilson [Andrewartha and Birch, 1954, MacArthur and Wilson,
1967]. The emergence and subsequent development of this concept as a
fundamental aspect of modern ecological theory has been extensively documented
in [Simberloff, 1988] and [Hanski and Simberloff, 1997]. Metapopulation theory
has provided key results into the dynamics of species inhabiting fragmented
populations (patches). Since the seminal work by Levins, metapopulation models
have been widely used to investigate the dynamics and persistence of
fragmented populations under different scenarios. Metapopulations, consisting
of semi-independent subpopulations connected by dispersal, have been a subject
of interest for population ecologists [Clobert et al., 2012]. The dynamics of
population size in these metapopulations have been extensively studied, with a
focus on understanding factors influencing their stability [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 widespsread 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]. Other species such as butterflies, amphibia, or birds can also
distribute in metapopulatons. For instance, 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 examples, among many others, suggest that metapopulation theory can play
a crucial role in delving into the dynamics of spatially distributed
population and in providing a theory of applicability for conservation. In
this sense, simple metapopulation models such as those considering two
subpopulations (e.g., inhabitting two habitats such as wetlands and bogs) may
provide a first initial approach to natural metapopulations. Metapopulations
have been widely studied using discrete- [Dey et al., 2014] and continuous-
time [Pulliam, 1988, Sardanyés, J. and Fontich, E., 2010, Sardanyés, J. 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 resilience under local and global
perturbations in coupled subpopulations modeled with logistic and Ricker maps
[Allen et al., 1993]. Another model considering two-patch 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, suggesting the so-called ghosts as transient-generator
mechanisms [Sardanyés, J. and Fontich, E., 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]. To gain a deeper understanding of
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 stochasticity effects, either intrinsic or extrinsic, in the
overall metapopulation dynamics. As far as we know, few works have explored
the interplay between noise and metapopulations’ persistence, mainly providing
numerical results in discrete-time models [Allen et al., 1993, Sardanyés, J.
et al., 2019].

In this article, we provide novel insights for continuous-time two-patch
metapopulation models by providing a thorough analytical and numerical study.
We focus on how dispersal impacts both local and global populations providing
an analytical derivation of the parameter conditions maximizing the population
of the entire metapopulation. Despite the origin, equilibrium at which both
subpopulations are zero, is unstable in the deterministic approach and thus no
possible extinctions are found with ecologically-meaningful parameters, we
show how increased diffusion can change its linear behaviour affecting
resilience. Specifically, we show that the vector field around the origin
changes its unstable nature from a source to a saddle at increasing diffusion.
This dispersal-induced change of local behaviour enhances the persistence of
local populations thus providing resilience to the metapopulation to face
random fluctuations.

This paper is structured as follows. In Section 2 we first set the model (a
two dimensional differential equation) and discuss its parameters.
Particularly, we provide an adimensionalization of units that is useful to
avoid cumbersome notation. In subsequent subsections, we discuss dynamical
aspects of the model. Namely, 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 takeaway of these sections is that the global dynamics is
affected remarkably by the stability of the origin. In Section 3 we
investigate how global dynamics can be relevant under stochastic
perturbations. There, we simulate random deaths to test the role of diffusion
in the metapopulation resilience. In Section 4 we derive analytically a
formula for the optimal dispersal i.e. the one that maximizes the total
population. Finally, in Section 5 the conclusions of the work are discussed.

##  2 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:

$\dot{x}_{i}=r_{i}\,x_{i}\left(1-\dfrac{x_{i}}{k_{i}}\right)+D(x_{j}-x_{i}),\phantom{m}{\rm
with}\phantom{m}i,j=1,2;\phantom{m}i\neq j.$ (1)

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 and the unit of time such that the rate of birth in the first patch is
equal to 1. This fixes $r_{1}=k_{1}=1$ without loss of generality. We let the
parameters $r_{2}$, $k_{2}$ and $D$ free for the present study. In order to
avoid cumbersome notation, we rename $r_{2}$ as $r$ and $k_{2}$ as $k$. With
these modifications, the system reads as

$\begin{array}[]{rcl}\dot{x}_{1}&=&x_{1}\left(1-x_{1}\right)+D(x_{2}-x_{1}),\\\\[10.0pt]
\dot{x}_{2}&=&rx_{2}\left(1-\dfrac{x_{2}}{k}\right)+D(x_{1}-x_{2}).\end{array}$
(2)

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 Section 3. 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].

###  2.1 Local and global dynamics

We first inspect the number of equilibrium points for Eqs. (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.

######  Lemma 2.1 (Number of equilibrium points).

System (2) has two, three or four equilibrium points.

###### Proof.

Let us assume, first, that $D\neq 0$ (the case $D=0$ can be studied
separately). The nonlinear system of equations verified by the equilibrium
points is:

$\displaystyle x_{1}-x_{1}^{2}+D(x_{2}-x_{1})$ $\displaystyle=0,$
$\displaystyle rx_{2}-\alpha x_{2}^{2}+D(x_{1}-x_{2})$ $\displaystyle=0.$

Here, $\alpha=r/k$. We isolate $x_{2}$ with respect to $x_{1}$ using the first
equation:

$x_{2}=\frac{1}{D}x_{1}^{2}+\frac{D-1}{D}x_{1}.$ (3)

Now, we substitute $x_{2}$ in the second equation, obtaining:

$-\frac{\alpha}{D^{2}}x_{1}^{4}-\frac{2\alpha(D-1)}{D^{2}}x_{1}^{3}+\left(\frac{r}{D}-\frac{\alpha(D-1)^{2}}{D^{2}}-1\right)x_{1}^{2}+\left(1+r-\frac{r}{D}\right)x_{1}=0.$
(4)

Obviously, $x_{1}=0$ is always a solution of this quartic. Substituting
$x_{1}=0$ in (3), we get $x_{2}=0$ and we recover the equilibrium 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 case $D=0$ can be studied separately. The system has four equilibrium
points which can be obtained easily. Those are $(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\to\infty$ 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 Section 4.
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.

######  Remark 1 (Limit case).

If $D\to\infty$, the quartic (4) converges uniformly to the polynomial

$-(\alpha+1)x_{1}^{2}+(1+r)x_{1}.$

Its non-trivial root is

$x_{1}=\frac{1+r}{\alpha+1}=\frac{k+kr}{r+k}.$

On the other hand, Eq. (3) implies that

$\lim_{D\to\infty}x_{2}(D)=\lim_{D\to\infty}x_{1}(D).$

It holds that,

$\lim_{D\to\infty}\left(x_{1}(D)+x_{2}(D)\right)=2\frac{k+kr}{r+k}.$ (5)

Figure 1: Stability of the origin in the plane $(r,D)$. The bifurcating curve
separating both instability types is given by $r/(1+r)$.

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.

######  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)$.

###### Proof.

The Jacobian matrix of the vector field evaluated at the point $(0,0)$ is
given by:

$J(0,0)=\begin{pmatrix}1-D&D\\\ D&r-D\end{pmatrix}.$

An elementary computation shows that the eigenvalues are:

$\lambda_{\pm}=\frac{1}{2}\left(1-2D+r\pm\sqrt{4D^{2}+(1-r)^{2}}\right).$

Notice that $\lambda_{+}$ is always positive. On the other hand, $\lambda_{-}$
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 Figure 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. ∎

Figure 2: Sketch of the region defined by the vertical and horizontal
nullclines with $r<D$ (left) and $r>D$ (right).

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. 2.

######  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, $\frac{d}{dt}x_{1}=0$, is described by the parabola:

$x_{2}=g_{v}(x_{1}):=\frac{1}{D}\left(x_{1}^{2}-(1-D)x_{1}\right).$

The horizontal nullcline, $\frac{d}{dt}x_{2}=0$, is given by:

$x_{1}=g_{h}(x_{2}):=\frac{1}{D}\left(\alpha x_{2}^{2}-(r-D)x_{2}\right).$

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)/\alpha)$. It holds that,

$\frac{d}{dx_{1}}g_{v}(N_{1})=(D-D^{2}),$

and

$\frac{d}{dx_{2}}g_{h}(N_{2})=\frac{D(-D+r)}{\alpha}.$

The two parabolas intersect in a point contained in the first quadrant.
Indeed, if $r>D$, at the vertical axis, the horizontal nullcline $g_{h}$ is
above the vertical one $g_{v}$. Moreover, it behaves as $\sqrt{x_{1}}$ when
$x_{1}$ is large. On the other hand, the vertical nullcline behaves as
$x_{1}^{2}$. A direct application of the Bolzano theorem shows that there
exist a crossing point for $x_{1}>0$ and $x_{2}>0$. If $r<D$, 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 component of the vector field is positive along
the vertical axis for values of $0<x_{2}<(r-D)/\alpha$ while the second
component of the vector field is positive along the horizontal axis for
$x_{1}<(1-D)$. This, together with the fact that the parabolas are nullclines
establishes the positive invariance of the region. If $r<D$ the result follows
from the second component 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. 2). For the next result we shall
focus on the intersection between the region enclosed by the isoclines and the
first quadrant. This intersection, let us name it $\mathcal{I}^{+}$, is also
positively invariant.

######  Remark 4 (Unstable eigenvector of the origin).

The eigenvector associated to the $\lambda_{+}$ eigenvalue always points
towards $\mathcal{I}^{+}$. Indeed, the eigenvector associated to $\lambda_{+}$
is given by:

$v_{+}=\left(\dfrac{1-r+\sqrt{(1-r)^{2}+4D^{2}}}{2D},1\right)$

Observe that the first component of $v_{+}$ is always positive, therefore,
pointing to the first quadrant. If $r>D$, the isoclines crosses the axis for
positive values. Therefore, always point to $\mathcal{I}^{+}$. If $r<D$, then
the slope of the eigenvector must be between the slopes of $g_{v}$ and
$g_{h}$:

$\frac{D-1}{D}<\dfrac{2D}{1-r+\sqrt{(1-r)^{2}+4D^{2}}}<\frac{D}{D-r}.$

Both inequalities are true while

$r<D.$

######  Theorem 2.4 (Global dynamics).

The following sentences hold for system (2):

  1. 1.

There are no periodic orbits.

  2. 2.

It has an equilibrium point ($P_{3}$) inside the first quadrant that is
Globally Asymptotically Stable (GAS) in the first quadrant.

  3. 3.

If $D<r/(1+r)$, there exist infinitely many heteroclinic connections between
the origin ($P_{0}$) and $P_{3}$.

  4. 4.

If $D>r/(1+r)$ there is a unique heteroclinic connection between $P_{0}$ and
$P_{3}$.

###### Proof.

We shall prove the statements by order.

  1. 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. 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 non of the equilibria and, therefore, the
$\omega$-limit for a full measure set of initial conditions is $(1,k)$. For
the case $D\neq 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. 3.

If $P_{0}$ is a source, any initial condition on $\mathcal{I}^{+}$ with
$\alpha$-limit $P_{0}$ must have $P_{3}$, which is GAS, as $\omega$-limit.

  4. 4.

If $P_{0}$ is a saddle, the unstable manifold has $P_{3}$ as $\omega$-limit.
Any other initial condition in $\mathcal{I}^{+}$ has its $\alpha$-limit
outside the first quadrant.

∎

###  2.2 On the role of the heteroclinic 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 resilience 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
eigenvalues of the Jacobian matrix evaluated at the origin are

$\lambda_{1}=1,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \lambda_{2}=1-2D.$

and the corresponding eigenvectors are given by

$v_{1}=(1,1)^{T},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ v_{2}=(-1,1)^{T}.$

Consider the linear change of variables given by the eigenvectors.

$\begin{pmatrix}y_{1}\\\ y_{2}\end{pmatrix}=\begin{pmatrix}1&-1\\\
1&1\end{pmatrix}\begin{pmatrix}x_{1}\\\ x_{2}\end{pmatrix}.$

Notice that the new variables $y_{1}$ and $y_{2}$ represent the difference 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:

$\displaystyle\begin{split}\dot{y}_{1}&=(1-2D)y_{1}+\left(\frac{1}{4k}-\frac{1}{4}\right)\left(y_{1}^{2}+y_{2}^{2}\right)-\left(\frac{1}{2}+\frac{1}{2k}\right)y_{1}y_{2},\\\
\dot{y}_{2}&=y_{2}-\left(\frac{1}{4k}+\frac{1}{4}\right)\left(y_{1}^{2}+y_{2}^{2}\right)-\left(\frac{1}{2}-\frac{1}{2k}\right)y_{1}y_{2}.\end{split}$
(6)

This vector field will be used later to study how the heteroclinic connection
depends on the normalized carrying capacity $k$. For the moment being, we
focus on the symmetric case. If we chose $k=1$, then, some terms vanish and
the resulting ODE is

$\displaystyle\begin{split}\dot{y}_{1}&=(1-2D)y_{1}-y_{1}y_{2},\\\
\dot{y}_{2}&=y_{2}-\frac{1}{2}\left(y_{1}^{2}+y_{2}^{2}\right).\\\
\end{split}$ (7)

Let us analyze 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}\leq y_{2}\leq 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

$\displaystyle P_{0}$ $\displaystyle=(0,0),$ $\displaystyle P_{1}$
$\displaystyle=\left(\sqrt{1-4D^{2}},1-2D\right),$ $\displaystyle P_{2}$
$\displaystyle=\left(-\sqrt{1-4D^{2}},1-2D\right),$ $\displaystyle P_{3}$
$\displaystyle=(0,2).$

The Jacobian matrix in these coordinates reads

$\begin{pmatrix}(1-2D)-y_{2}&-y_{1}\\\ -y_{1}&1-y_{2}\end{pmatrix}.$

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:

$\dot{y}_{1}=y_{1}\left((1-2D)-y_{2}\right).$

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. 3. If $D>1/2$ the
vertical axis is always attracting.

$y_{2}$$y_{1}$$y_{2}=y_{1}$$y_{2}=-y_{1}$$1-2D$$P_{0}$$P_{3}$ Figure 3:
Attracting character of the vertical axis in system 7. The red lines delimit
the admissible region. The vertical axis is an heteroclinic connection between
$P_{0}$ and $P_{3}$ and it is repelling if $y_{2}<2D-1$. See Section 2.2 for
more details.

This global property of the phase space is relevant in the capacity of the
system to use 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}=\pm y_{1}\\}$ remain close to them (see
Fig. 3). This suggests that the recovering of a patch is much slower when the
origin is a source. A natural conjecture, to be tested in Section 3, is that
having a saddle is favorable for recovery if one of the patches is populated
by few individuals, thus providing more resilience to the metapopulation. The
area of the sub-region 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

$\begin{cases}U_{1}(s)=as^{2}+\mathcal{O}(s^{3}),\\\
U_{2}(s)=s+bs^{2}+\mathcal{O}(s^{3}).\end{cases}$

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:

$Y(U(s))=\dot{U}(s)\lambda s+\mathcal{O}(s^{3}),$

Here, $Y$ stands for the vector field of system (6) and $\lambda=1$ is the
unstable eigenvalue related to the origin. Solving the latter equation for $a$
and $b$ results in

$a=\frac{1-k}{4k(1+2D)},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
b=-\frac{1+k}{4k}.$

Looking at the sign of $a$ we can understand how the manifold bends (at second
order) when $k\neq 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$.

##  3 Role of diffusion in metapopulaiton resilience to perturbations

In Section 2.2, we analyzed the local dynamics close to the origin and
concluded that if $D$ is chosen such that the origin is a saddle-type
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
recovering capacity (resilience) of the metapopulation when one of the two
subpopulations is close to extinction.

Figure 4: (a) Vector field of the system for $D=0.1$ (left) and $D=0.9$
(right). The colour of the arrows denote the norm of the field. (b) Different
time series that lead to local (upper) and global (lower) extinctions for
$D=0.1$ Notice that the left-upper panel shows local extinctions which are
recovered by the diffusion from $x_{1}$, while the left-lower shows the
collapse of the metapopulation. The local populations and thus the entire
metapopulations rapidly leave the origin with $D=0.9$, ensuring the
persistence of the metapopulation. (c) Extinction probabilities ($P_{1,2}$ and
$P_{T}$ denote local and global extinctions, respectively) for initial
conditions close to the origin for $D=0$ (first row), $D=0.1$ (second row) and
$D=0.9$ (third row). The dashed red lines indicate the values of the
populations of $x_{1}(0)$ where extinctions are observed. The red diagonal
line indicates the union of these values including $x_{2}(0)$ to ease the
visualization of the changes in extinctions due to increased diffusion. In all
the panels we fix $k=r=1$.

To investigate this phenomenon, we designed a simple numerical experiment
simulating random deaths due to some external agent (e.g., predators, climatic
factors, hunting). The process is as follows:

  1. 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 $\ell_{2}$-norm.

  2. 2.

We select a grid of $500\times 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. 3.

For each initial condition, we perform an integration and, at random intervals
following a uniform distribution $\mathcal{U}(0,0.1)$, we subtract a random
quantity from the total population $y_{2}$. This quantity also follows a
uniform distribution $\mathcal{U}(0,1/50)$. (i.e. the maximum subtracted value
is the total carrying capacity divided by $100$).

  4. 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),
then we rearrange 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 labeled as being at risk
of extinction in patch 1 (or 2).

  5. 5.

We perform this experiment 100 times for each initial condition, associating
an extinction probability with each initial condition.

In Fig. 4(a) we display the direction of the vector field of system (2) for
$D=0.1$ and $D=0.9$ (under symmetric conditions, $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
trajectories are rapidly attracted by the heteroclinic connection (Fig. 4(a)).
These two plots provide an illustration of the main point of Section 2.2,
namely, the system is better at recovering a patch which is close to
extinction if the origin is a saddle.

Panel (b) displays some representative examples of local and global
extinctions together with survival dynamics by means of time series computed
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 metapopulaiton. In the case $D=0.9$, the attracting
character 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
metapopulaion, 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 total populations.
Total extinction involves that both time series $x_{1}$ and $x_{2}$ abandon
the first quadrant simulatenously and the integration is terminated. The
extinction probability for subpopulations refer 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.

##  4 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 $\varepsilon$, any values $r$, $k$ $D$, and any initial
condition $(x_{1}(0),x_{2}(0))$, there exist a time $t$ for which

$\left\|\varphi^{D,r,k}_{t}(x_{1}(0),x_{2}(0))-P_{3}(D)\right\|_{2}<\varepsilon.$

Here, $\varphi^{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

$\Omega(D):=\|P_{3}(D)\|_{1}.$

They are able to prove the following result.

######  Lemma 4.1 ([Ruiz-Herrera and Torres, 2018]).

$\Omega^{\prime}(0)=\left(\frac{1}{r}-1\right)(1-k).$

Henceforth, the function $\Omega$ is increasing nearby $D=0$ if $r>1$ and
$k>1$ or if $r<1$ and $k<1$. As $\Omega(0)>1+k$, under these conditions,
$\Omega$ can exceed the sum of the carrying capacities of the patches. The
behaviour of $\Omega$ has been analyzed previously for arbitrarily large
diffusion values (namely $D\to\infty$). 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:

$\lim_{D\to\infty}\Omega(D)=1+k+(1-k)\frac{k-r}{k+r},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ k\geq
1,\leavevmode\nobreak\ \leavevmode\nobreak\ r/k>1.$

Notice that formula (5) from Remark 1 is a compact version of this one. In
particular, using (5) it can be shown that $\Omega(\infty)>1+k$ whenever $k$
is contained between $1$ and $r$ (without assuming $r>1$ or $r<1$).

The qualitative behaviour of the function $\Omega$ 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 $\Omega$ is a natural question.
In this section we provide an answer and, moreover, recover some insights on
the qualitative behaviour of the function $\Omega$.

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.

######  Lemma 4.2 (The equilibria lie on an ellipse).

The equilibrium points of system (2) are contained on an ellipse with center
$(1/2,k/2)$ and axes $\sqrt{a}$ and $\sqrt{b}$, where:

$\displaystyle a$ $\displaystyle=(1+rk)/4,$ $\displaystyle b$
$\displaystyle=(k^{2}+k/r)/4.$

###### Proof.

We consider once again the system of equations for the equilibria:

$\displaystyle x_{1}-x_{1}^{2}+D(x_{2}-x_{1})$ $\displaystyle=0,$
$\displaystyle rx_{2}-\frac{r}{k}x_{2}^{2}+D(x_{1}-x_{2})$ $\displaystyle=0.$

Adding both equations, we get the following expression.

$x_{1}-x_{1}^{2}+rx_{2}-\frac{r}{k}x_{2}^{2}=0,$

which is equivalent to:

$\frac{(x_{1}-\frac{1}{2})^{2}}{a}+\frac{(x_{2}-\frac{k}{2})^{2}}{b}=1,$ (8)

This is, in fact, the equation of an ellipse centered at $(1/2,k/2)$ and with
semi-axes $\sqrt{a}$ and $\sqrt{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:

$\frac{x_{1}^{2}-x_{1}}{x_{2}-x_{1}}=\frac{\alpha
x_{2}^{2}-rx_{2}}{x_{1}-x_{2}},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \alpha=\frac{r}{k}.$ (9)

Moreover, the point $(x_{1},x_{2})$ is an equilibrium point of system (2) for

$D=\frac{x_{1}^{2}-x_{1}}{x_{2}-x_{1}}.$

###### 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:

$\frac{a}{b}=\frac{1+\alpha k^{2}}{k^{2}+\alpha^{-1}}=\alpha.$

Notice also that, by square completion,

$x_{1}^{2}-x_{1}=(x_{1}-1/2)^{2}-1/4.$

If $(x_{1},x_{2})$ verify Eq. (8), then the above quantity is equal to

$a-(1/4+k^{2}\alpha/4)-(\alpha x_{2}^{2}-\alpha kx_{2}).$

The identity is proven by noticing that

$(1/4+k^{2}\alpha/4)=a,$

and recalling that $\alpha k=r$. Indeed,

$\frac{x_{1}^{2}-x_{1}}{x_{2}-x_{1}}=\frac{-(\alpha x_{2}^{2}-\alpha
kx_{2})}{x_{2}-x_{1}}=\frac{\alpha x_{2}^{2}-rx_{2}}{x_{1}-x_{2}}.$

∎

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 identity (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 Section 2.2).

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 $\Omega$ on the
ellipse (8). Notice that as the objective function is continuous and the set
of feasible solutions a compact, this problem has always a solution. However,
the 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\to\infty$.
The next result provides the solution of the aforementioned optimization
problem and it is tackled using Lagrangian Multipliers.

######  Theorem 4.4 (Optimal dispersal rate).

Let $(\bar{x}_{1}(D),\bar{x}_{2}(D))$ be the coexistence solution of system
(2) and $\Omega(D)=\bar{x}_{1}(D)+\bar{x}_{2}(D)$, then:

$\Omega(D)\leq\frac{1}{2}\left(\sqrt{1+\frac{r^{2}+1}{r}k+k^{2}}+(1+k)\right)=:\Omega^{*}.$
(10)

Moreover, this bound is sharp: There exist $D^{*}$ for which
$\Omega(D^{*})=\Omega^{*}$.

###### Proof.

In view of Lemma 4.2, the maximal total population can be obtained by
maximizing the function $\Omega$ on the ellipse (8). This is equivalent to
solve the following optimisation problem:

maximize $\displaystyle X_{1}+X_{2}+\frac{1}{2}(1+k),$ subject to
$\displaystyle\frac{X_{1}^{2}}{a}+\frac{X_{2}^{2}}{b}=1.$

Here, $X_{1}=x_{1}-1/2$ and $X_{2}=x_{2}-k/2$. To solve the problem, we name
$f(X_{1},X_{2})=X_{1}+X_{2}+\frac{1}{2}(1+k)$ the objective function and
$H(X_{1},X_{2})=\frac{X_{1}^{2}}{a}+\frac{X_{2}^{2}}{b}-1$ the restriction
function. Let us maximize the Lagrangian function.

$\mathcal{L}(X_{1},X_{2})=f(X_{1},X_{2})+\lambda H(X_{1},X_{2}),$

for some real parameter $\lambda$. If we look for the roots of
$\nabla\mathcal{L}$ we notice that those are:

$\left(\pm\frac{a}{\sqrt{a+b}},\pm\frac{b}{\sqrt{a+b}}\right).$

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

$\sqrt{a+b}+\frac{1}{2}(1+k),$

which corresponds to $\Omega^{*}$ if we expand $a$ and $b$ in terms of $k$ and
$r$. Moreover, if we we recast $(X_{1},X_{2})$ to the original coordinates, we
can use Lemma 4.3 to find the value of $D$ for which the optimum is an
equilibrium point:

$D^{*}=\frac{x_{1}^{2}-x_{1}}{x_{2}-x_{1}},$ (11)

with

$x_{1}=\frac{a}{\sqrt{a+b}}+\frac{1}{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
x_{2}=\frac{b}{\sqrt{a+b}}+\frac{k}{2}.$

∎

Figure 5: Behaviour of the total population $\Omega$ 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 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. The red
stars correspond to the values of $r$ and $k$ used in Fig. 6 where the total
population is shown at increasing values of $D$.

There are several consequences of Theorem 4.4. Notice first that, if $r=1$,
then $\Omega(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, of course, a known fact that can be recovered from Eq. (10).
Moreover, the linear coefficient of the quadratic polynomial inside the square
root is

$\frac{r^{2}+1}{r},$

which is a function whose minimal value is $2$ and it is achieved precisely at
$r=1$. In other words,

$\Omega^{*}(D)\geq\frac{1}{2}\left(\sqrt{1+2k+k^{2}}+(1+k)\right)=(1+k).$

Hence, $\Omega^{*}$ 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. 5). The following result provides this
classification.

######  Corollary 4.5 (Behaviour of the total population).

$\Omega^{*}$ is achieved for a positive dispersion rate $D^{*}$ if one of the
following conditions hold:

  1. 1.

$r>1$ and $k>k^{*}(r)$,

  2. 2.

$r<1$ and $k<k^{*}(r)$,

where

$k^{*}(r)=\frac{r(3+r)}{1+3r}.$

Moreover, if $r<1$ and $k^{*}(r)<k<1$ then $\Omega^{\prime}(D)>0$ for $D>0$.
In any other case, $\Omega(D)<1+k$ if $D>0$.

###### Proof.

From the proof of Theorem 4.4 we know that $\Omega$ has only one critical
value that is an equilibrium point inside the first quadrant. We also know
that the maximal value is achieved at:

$D^{*}=\frac{x_{1}^{2}-x_{1}}{x_{2}-x_{1}},$

where,

$x_{1}=\frac{a}{\sqrt{a+b}}+\frac{1}{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
x_{2}=\frac{b}{\sqrt{a+b}}+\frac{k}{2}.$

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, $\Omega(D)$ is increasing
near $D=0$ (see Lemma 4.1). As there are no other relative optimal values of
$D$, $\Omega(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^{*}\to\infty$.

The information provided by Corollary 4.5 is shown in Fig. 5. The horizontal
axis represents the value of $r$ and, the vertical, the value of $k$.
Different regions are colored according the behaviour of the function $\Omega$
with respect to the dispersal rate. The red zone represent 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 $\Omega$ is
monotonically decreasing with respect to $D$. Within the grey regions the
total population is an increasing function of $D$ but the optimal value
$\Omega^{*}$ is never achieved for positive values of $D$. The green regions
represent values of the parameters for which $\Omega^{*}$ is achieved for some
positive $D$. The line $k=r$ determines the limit case $\Omega(\infty)<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 $\Omega^{*}>1+k$
but $\Omega(\infty)<1+k$.

Figure 6: Total population $\Omega$ with respect to diffusion with: (a)
$r=1.5$, $k=1.5$, (b) $r=0.5$, $k=1.5$, (c) $r=1.5$, $k=k^{*}(1.5)\approx
1.227$, and (d) $r=1.5$, $k=1.7$

Figure 6 displays the evolution of $\Omega$ as function of $D$ for four
different choices of the parameters $r$ and $k$, shown in small red stars in
Fig. 5. 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 6(a)
has been obtained with $r=1.5$ and $k=1.5$. At $D=0$ the value of $\Omega$ is
$1+k=2.5$. The characteristic curve increases until some maximal value that
can be computed using formula (10) (in this case $\Omega^{*}\approx 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 $\Omega$ decreases asymptotically to the
value $2.5$ (i.e. the sum of carrying capacities: $1+k$). In Fig. 6(b) we set
$r=0.5$ and $k=1.5$. These are conditions where the theory predicts that
positive positive dispersal rate is detrimental to the total population.
Indeed, the function $\Omega$, in this case, is strictly decreasing until the
asymptotic limit $\Omega=2.25$. Notice that $\Omega^{*}$ is achieved for some
$D<0$. In Fig. 6(c) 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 $\Omega$ is achieved at the
limit case and, therefore, the function $\Omega$ is strictly increasing for
positive values of $D$. The asymptotic value $\Omega$ is, in this case,
$2.25$. This value can be predicted with both Theorem 4.4 and Remark 1.
Finally, Fig. 6(d) illustrates a case in which the limit of $\Omega$ 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 ($\approx 2.72$)
and it is achieved for some positive $D$.

##  5 Conclusions

The investigation of metapopulation mathematical models is very extensive in
the literature. Most of these studies have focused on small metapopulations
considering few pacthes [Hastings, 1993, Arditi et al., 2015, Arditi et al.,
2016, Sardanyés, J. and Fontich, E., 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 co-workers 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,
proposing the so-called ghosts as transient-generator mechanisms [Sardanyés,
J. and Fontich, E., 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 equilibrum 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 recovering capacity of the system.
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 resilience of the metapopulation. 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 extinction under perturbations. To illustrate this fact we have
run a simulation in which the system has stochastic losses of individuals in
both patches. These random perturbations simulate a source of external noise
(sick animals, tough environmental conditions, hunting).

The second contribution is related to a well known counter-intuitive property
of the system: under suitable conditions, having positive dispersal rate,
leads the total population to stabilize overcoming the sum of the carrying
capacities of the patches. These 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
respective 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}$. On 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 the total population 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:

$\overline{\Omega}^{*}=\frac{1}{2}\left(\sqrt{(k_{1}+k_{2})^{2}+\frac{(r_{1}-r_{2})^{2}}{r_{1}r_{2}}k_{1}k_{2}}+k_{1}+k_{2}\right),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \overline{D}^{*}=D^{*}/r_{1}.$

The main takeaway of the above formula is the following: in standard $k-r$
positiveness conditions $\Omega^{*}$ is always larger than the limit case
($D\to\infty$). If $r-k$ positiveness is not fulfilled but the hypothesis
appearing in [Ruiz-Herrera and Torres, 2018] is, $\Omega^{*}$ is superior to
the sum of the carrying capacities but smaller than the total population in
the limit case. If non of the hypotheses are verified, $\Omega^{*}$ exceeds
the total carrying capacity of the system provided that $r_{1}\neq r_{2}$ but
is not achieved for a positive value of the dispersal rate.

## Technical details

The programs corresponding to Sections 3 and 4 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).

## Acknowledgments

This work has been supported by the Spanish State Research Agency (AEI)
through grants PGC2018-100928-B-I00 (DP), PGC2018-100699-B-I00
(MCIU/AEI/FEDER, UE) (MJC), and by Catalan Government grant 2017 SGR 1374
(MJC). JS has been supported by the Ramón y Cajal grant RYC-2017-22243 funded
by MCIN/AEI/10.13039/501100011033 ”FSE invests in your future”, and by the
2020-2021 Biodiversa+ and Water JPI joint call under the BiodivRestore ERA-NET
Cofund (GA N°101003777) project MPA4Sustainability with funding organizations:
Innovation Fund Denmark (IFD), Agence Nationale de la Recherche (ANR),
Fundação para a Ciência e a Tecnologia (FCT), Swedish Environmental Protection
Agency (SEPA), and grant PCI2022-132926 funded by
MCIN/AEI/10.13039/501100011033 and by the European Union
NextGenerationEU/PRTR. This work has been also funded through the Severo Ochoa
and María de Maeztu Program for Centers and Units of Excellence in R&D
(CEX2020-001084-M). We thank CERCA Programme/Generalitat de Catalunya for
institutional support. DP has been supported by the Research Program 2021-2023
of Universidad Internacional de la Rioja. RO thanks to the Faculty for the
Future program of the Schlumberger Foundation for the scholarship.

## References

  * [Abbott, 2011] Abbott, K. (2011).  A dispersal-induced paradox: synchrony and stability in stochastic metapopulations.  Ecol. Lett., 14.
  * [Allen et al., 1993] Allen, J., Schaffer, W., and Rosko, D. (1993).  Chaos reduces species extinction by amplifying local population noise.  Nature, 364:229–232.
  * [Anderson et al., 1999] Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., and Sorensen, D. (1999).  LAPACK Users’ Guide.  Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition.
  * [Andrewartha and Birch, 1954] Andrewartha, H. and Birch, L. (1954).  The distribution and abundance of animals.  Oxford University Press, Chicago: University of Chicago Press.
  * [Angelis et al., 1979] Angelis, D. D., Travis, C., and Post, W. (1979).  Persistence and stability of seed-dispersed species in a patchy environment.  Theor. Pop. Biol., 16(2):107–125.
  * [Arditi et al., 2016] Arditi, R., Bersier, L.-F., and Rohr, R. (2016).  The perfect mixing paradox and the logistic equation: Verhulst vs. lotka: Comment.  Ecosphere, 7(11):479–500.
  * [Arditi et al., 2015] Arditi, R., Lobry, C., and Sari, T. (2015).  Is dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation.  Theor. Pop. Biol., 106:45–59.
  * [Bonsall et al., 2002] Bonsall, M., French, D., and Hassell, M. (2002).  Metapopulation structures affect persistence of predator–prey interactions.  J. Anim. Ecol., 71:1075–1084.
  * [Clobert et al., 2012] Clobert, J., Baguette, M., Bemton, T., and Bullock, J. e. (2012).  Dispersal Ecology and Evolution.  Oxford University Press, Great Clarendon Street, Oxford, OX2 6DP.
  * [Dey et al., 2014] Dey, S., Goswami, B., and Joshi, A. (2014).  Effects of symmetric and asymmetric dispersal on the dynamics of heterogeneous metapopulations: Two-patch systems revisited.  J. Theor. Biol., 345:52–60.  Article.
  * [Dey and Joshi, 2011] Dey, S. and Joshi, A. (2011).  Stability via asynchrony in drosophila metapopulations with low migration rates.  Science, 312:434–436.
  * [Doebeli, 1995] Doebeli, M. (1995).  Dispersal and dynamics.  Theoretical Population Biology, 47(1):82–106.
  * [Ellner et al., 2011] Ellner, S., McCauley, E., Kendall, B., Briggs, C., Hosseini, P., Wood, S., Janssen, A., Sabelis, M., Turchin, P.and Nisbet, R., and Murdoch, W. (2011).  Habitat structure and population persistence in an experimental community.  Nature, 412:538–543.
  * [Fahrig and Merriam, 1985] Fahrig, L. and Merriam, G. (1985).  Habitat patch connectivity and population survival.  Ecology, 66(6):1762–1768.
  * [Fang et al., 2020] Fang, M., Wang, Y., Chen, M., and DeAngelis, D. (2020).  Asymptotic population abundance of a two-patch system with asymmetric diffusion.  Discr. Cont. Dyn. Syst., 40(6):3411–3425.
  * [Franklin et al., 2000] Franklin, A., Anderson, D., Gutierrez, R., Burnham, K., and Graham, D. (2000).  Recognizing and managing uncertainty in the use of short-term time-series data in natural resource management. In Decision Making in Natural Resource Management.  Springer, Boston, MA.
  * [Freedman and Waltman, 1977] Freedman, H. I. and Waltman, P. (1977).  Mathematical models of population interactions with dispersal. i: Stability of two habitats with an without a predator.  SIAM Journal on Applied Mathematics, 32(3):631–648.
  * [Gonzalez-Andujar and Perry, 1993] Gonzalez-Andujar, J. and Perry, J. N. (1993).  Chaos, metapopulations and dispersal.  Ecological Modelling, 65(3):255–263.
  * [Gyllenberg et al., 1993] Gyllenberg, M., Soderbacka, G., and Ericsson, S. (1993).  Does migration stabilize local-population dynamics - analysis of a discrete metapopulation model.  Mathematical Biosciences, 118(1):25–49.
  * [Hanski and Gilpin, 1994] Hanski, I. and Gilpin, M. (1994).  Metapopulation dynamics: brief history and conceptual domain.  Biol. J. Linn. Soc., 51:187–208.
  * [Hanski and Ovaskainen, 2003] Hanski, I. and Ovaskainen, O. (2003).  Metapopulation theory for fragmented landscapes.  Theor. Pop. Biol., 64:119–127.
  * [Hanski and Simberloff, 1997] Hanski, I. and Simberloff, D. (1997).  The metapopulation approach, its history, conceptual domain and application to conservation. In Metapopulation biology: Ecology, genetics, and evolution (I. Hanski and M.E. Gilpin, eds., pp. 5–26).  Academic Press, San Diego CA.
  * [Haro et al., 2016] Haro, A., Canadell, M., Luque, A., Mondelo, J.-M., and Figueras, J.-L. (2016).  The Parameterization Method for Invariant Manifolds. From Rigorous Results to Effective Computations, volume 195 of Appl. Math. Sci. Springer-Verlag.
  * [Hastings, 1993] Hastings, A. (1993).  Complex interactions between dispersal and dynamics - lessons from coupled logistic equations.  Ecology, 74(5):1362–1372.
  * [Holt, 1985] Holt, R. (1985).  Population-dynamics in 2-patch environments - some anomalous consequences of an optimal habitat distribution.  Theor. Pop. Biol., 28(2):181–208.
  * [Honnay et al., 2002] Honnay, O., Hermy, M., and Coppin, P. (2002).  Effects of area, age and diversity of forest patches in belgium on plant species richness, and implications for conservation and reforestation.  Biol. Cons., 99:3–12.
  * [Jorba and Zou, 2005] Jorba, À. and Zou, M. (2005).  A software package for the numerical integration of ODEs by means of high-order Taylor methods.  Exper. Math., 14(1):99–117.
  * [Kendall and Fox, 1998] Kendall, B. E. and Fox, G. A. (1998).  Spatial structure, environmental heterogeneity, and population dynamics: Analysis of the coupled logistic map.  Theoretical Population Biology, 54(1):11–37.
  * [Kritzer and Sale, 2006] Kritzer, J. and Sale, P. (2006).  Marine metapopulations.  Elsevier Academic Press, US.
  * [Levins, 1969] Levins, R. (1969).  Some demographic and genetic consequences of environmental heterogeneity for biological control.  Bull. Entomol. Soc. Am., 15:237–240.
  * [Levins, 1970] Levins, R. (1970).  Extinction. In Some mathematical problems in biology (M. Desternhaber, ed.).  American Mathematical Society.Oxford University Press, Providence, RI: American Mathematical Society.
  * [Lloyd, 1995] Lloyd, A. L. (1995).  The coupled logistic map: a simple model for the effects of spatial heterogeneity on population dynamics.  Journal of Theoretical Biology, 173(3):217–230.
  * [MacArthur and Wilson, 1967] MacArthur, R. and Wilson, E. (1967).  The theory of island biogeography. Princeton University Press., Princeton, NJ.
  * [Oro et al., 2023] Oro, D., Alsedà, L., Hastings, A., Genovart, M., and Sardanyés (2023).  Social copying drives a tipping point for nonlinear population collapse.  Proc. Natl. Acad. Sci., 120:e2214055120.
  * [Pulliam, 1988] Pulliam, H. (1988).  Sources, sinks, and population regulation.  Am. Nat., 132(5):652–661.  Article.
  * [Ruiz-Herrera, 2018] Ruiz-Herrera, A. (2018).  Metapopulation dynamics and total biomass: Understanding the effects of diffusion in complex networks.  Theor. Pop. Biol., 121:1–11.  Article.
  * [Ruiz-Herrera and Torres, 2018] Ruiz-Herrera, A. and Torres, P. (2018).  Effects of diffusion on total biomass in simple metacommunities.  Journal of Theoretical Biology, 447:12–24.
  * [Sardanyés, J. and Fontich, E., 2010] Sardanyés, J. and Fontich, E. (2010).  On the metapopulation dynamics of autocatalysis: extinction transients related to ghosts.  Int. J .Bifurc. Chaos, 20:1261–1268.
  * [Sardanyés, J. et al., 2019] Sardanyés, J., Piñero, J., and Solé, R. (2019).  Habitat loss-induced tipping points in metapopulations with facilitation.  Pop. Ecol., 61:436–449.
  * [Simberloff, 1988] Simberloff, D. (1988).  The contribution of population and community biology to conservation science.  Annu. Rev. Ecol. Syst., 19:473–511.
  * [Smith, 2022] Smith, J. (2022).  Coherence, Dynamics and Stability in Spatially Structured and Unstructured Populations: Effects of Dispersal and Crowding.  Phd thesis, Jawaharlal Nehru Centre for Advanced Scientific Research, Bangalore, India.
  * [Ylikarjula et al., 2000] Ylikarjula, J., Alaja, S., Laasko, J., and Tesar, D. (2000).  Effects of patch number and dispersal patterns on population dynamics and synchrony.  Journal of Theoretical Biology, 207(3):377–387.


which the verification of either $(i)$ or $(ii)$ becomes redundant. To prove
$(i)$, we shall adapt the argument in [8]. Note that since any measurable
function $f$ is the pointwise $dm$-a.e limit of its truncation by height
$f_{N}:=f1_{\left\lvert f\right\rvert\leq N}$, it suffices to approximate any
essentially bounded functions on $S$. To this end, we shall regard the
sequence of Taylor polynomials $\\{T_{n}\\}_{n}$ as continuous linear
operators from $X$ into balls of radius $M>0$ in $L^{\infty}(S)$, which
equipped with the weak-star topology are metrizable (for instance, see [38]).
We denote the ball by $B_{M}(L^{\infty}(A))$ and recall that the universality
criterion in THM:Univ1 assures that it suffices to check that for any pair
$(f,g)\in X\oplus B_{M}(L^{\infty}(A))$ there exists analytic polynomials
$\\{Q_{k}\\}$ and a subsequence $\\{T_{n_{k}}\\}_{k}$ such that

$(Q_{k},T_{n_{k}}(Q_{k}))\to(f,g)$

in $X\oplus B_{M}(L^{\infty}(A))$. However, since $T_{n}(Q_{k})=Q_{k}$ for any
$n$ sufficiently large, the universality criterion actually follows from the
assumption that $S$ satisfies the SA-property wrt $X$. By means of repeat the
same argument for any ball $B_{M}(L^{\infty}(A))$, the claim follows.

_Step 2: Statement $(ii)$:_

Let $(f,f^{*})\in\mathcal{P}_{A}(X)$ and assume that $\\{T_{n_{k}}f\\}_{k}$ is
a subsequence of Taylor polynomials of $f$, which converges pointwise $dm$-a.e
to a measurable function $g$. Pick a sequence of analytic polynomials
$(Q_{j})_{j}$ such that $(Q_{j},Q_{j})$ converges to $(f,f^{*})$ in $X\oplus
L^{\infty}(A)$. Fix an $\varepsilon>0$ and write

$m(\\{A:\left\lvert g-Q_{j}\right\rvert>\varepsilon\\})\leq m(\\{A:\left\lvert
g-T_{n_{k}}(f)\right\rvert>\varepsilon/2\\})+\frac{2}{\varepsilon}\sum_{m=0}^{n_{k}}\left\lvert
f_{m}-(Q_{j})_{m}\right\rvert.$

Using the assumption that $Q_{j}\to f$ in
$X\hookrightarrow\text{Hol}(\mathbb{D})$, it follows by standard arguments
involving uniform convergence on compacts, that the Taylor coefficients
$(Q_{j})_{n}$ of $Q_{j}$ converge to the Taylor coefficients $f_{n}$ of $f$.
Sending $j$ to infinity, we obtain

$\limsup_{j\to\infty}m(\\{A:\left\lvert
g-Q_{j}\right\rvert>\varepsilon\\})\leq m(\\{A:\left\lvert
g-T_{n_{k}}(f)\right\rvert>\varepsilon/2\\}),\qquad\forall k,\varepsilon>0.$

Now since $T_{n_{k}}(f)$ converges pointwise $dm$-a.e to $g$ on $A$, they also
converge in Lebesgue measure on $A$, hence we can send $k\to\infty$ and
conclude that $Q_{j}\to g$ in Lebesgue measure on $A$. By means of passing to
a subsequence and using the assumption that $Q_{j}\to f^{*}$ weak-star in
$L^{\infty}(A)$, we conclude that $g=f^{*}$ on $A$, which proves $(ii)$. The
proof is now complete. ∎

Recall that examples of Holomorphic Banach spaces $X$ with SA-sets of positive
Lebesgue measure are provided by THM:SARegX, hence the above theorem has
particularly interesting applications for such spaces. In particular, it
applies to weighted Bloch spaces $\mathcal{B}_{0}(w)$ with majorants $w$
satisfying condition (3). It was asked in [8] (see Remark 2.3 therein) whether
there exists Lebesgue measurable subsets $E\subset\partial\mathbb{D}$ which
induce Menshov universality in $\mathcal{B}_{0}$. An application of THM:Univ1
in conjunction with THM:SAsets, answers the question affirmatively in a strong
sense.

######  Corollary 8.3.

THM:Univ2 Let $w$ be a majorant satisfying

$\int_{0}^{1}\frac{w^{2}(t)}{t}dt=+\infty.$

Then for any number $0<\delta<1$, there exists a set $E_{\delta}$ of Lebesgue
measure $\delta$ and a corresponding subset $\mathcal{U}_{w,\delta}$ of
$\mathcal{B}_{0}(w)$, generic in the sense of Baire categories, such that
every element $f\in\mathcal{U}_{w,\delta}$ gives rise to a universal Menshov
pair $(f,E_{\delta})$ from $\mathcal{B}_{0}(w)$.

It was further asked in [8], which conditions on such sets $E$ in give rise to
the Mensov universality in $\mathcal{B}_{0}$. Our THM:CHAR1 gives a function
theoretical description of such sets, while THM:CHARREM and COR:Hwinv give
some necessary metric and/or geometric conditions.

## References

  * [1] L. Ahlfors and A. Beurling, Conformal invariants and function-theoretic null-sets, (1950).
  * [2] A. B. Aleksandrov, Invariant subspaces of shift operators. An axiomatic approach, Zap. Nauchn. Sem. Leningrad. Otdel. Mat. Inst. Steklov. (LOMI), 113 (1981), pp. 7–26, 264.
  * [3] A. B. Aleksandrov, J. Anderson, and A. Nicolau, Inner functions, Bloch spaces and symmetric measures, Proceedings of the London Mathematical Society, 79 (1999), pp. 318–352.
  * [4] A. Aleman and B. Malman, Density of disk algebra functions in de Branges–Rovnyak spaces, C. R. Math. Acad. Sci. Paris, 355 (2017), pp. 871–875.
  * [5] , Hilbert spaces of analytic functions with a contractive backward shift, Journal of Functional Analysis, 277 (2019), pp. 157–199.
  * [6] A. Aleman, S. Richter, and C. Sundberg, Nontangential limits in ${\mathcal{P}}^{t}(\mu)$-spaces and the index of invariant subspaces, Annals of mathematics, (2009), pp. 449–490.
  * [7] J. Anderson, J. Fernandez, and A. L. Shields, Inner functions and cyclic vectors in the Bloch space, Transactions of the American Mathematical Society, 323 (1991), pp. 429–448.
  * [8] H.-P. Beise and J. Müller, Generic boundary behaviour of Taylor series in Hardy and Bergman spaces, Mathematische Zeitschrift, 284 (2016), pp. 1185–1197.
  * [9] C. Bénéteau and B. Korenblum, Jensen type inequalities and radial null sets, Analysis, 21 (2001), pp. 99–106.
  * [10] S. Bloom and G. S. De Souza, Atomic decomposition of generalized Lipschitz spaces, Illinois Journal of Mathematics, 33 (1989), pp. 181–209.
  * [11] L. Carleson, Sets of uniqueness for functions regular in the unit circle, Acta mathematica, 87 (1952), pp. 325–345.
  * [12] J. Cima, A. Matheson, and W. Ross, The Cauchy transform, vol. 125 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, 2006.
  * [13] J. A. Cima and W. T. Ross, The backward shift on the Hardy space, no. 79, American Mathematical Soc., 2000.
  * [14] J. B. Conway, The theory of subnormal operators, no. 36, American Mathematical Soc., 1991.
  * [15] C. Costara and T. Ransford, Which de Branges-Rovnyak spaces are Dirichlet spaces (and vice versa)?, J. Funct. Anal., 265 (2013), pp. 3204–3218.
  * [16] G. S. De Souza et al., Fourier series and the maximal operator on the weighted special atom spaces, International Journal of Mathematics and Mathematical Sciences, 12 (1989), pp. 579–582.
  * [17] G. S. de Souza, R. O’Neil, and G. Sampson, Several characterizations for the special atom spaces with applications., Revista Matematica Iberoamericana, 2 (1986), pp. 333–355.
  * [18] P. L. Duren, Theory of H p Spaces, Academic press, 1970.
  * [19] J. Garnett, Bounded analytic functions, vol. 236, Springer Science & Business Media, 2007.
  * [20] J. B. Garnett and D. E. Marshall, Harmonic measure, vol. 2, Cambridge University Press, 2005.
  * [21] K.-G. Grosse-Erdmann, Universal families and hypercyclic operators, Bulletin of the American Mathematical Society, 36 (1999), pp. 345–381.
  * [22] V. P. Havin and B. Jöricke, The uncertainty principle in harmonic analysis, vol. 72 of Encyclopaedia Math. Sci., Springer, Berlin, 1995.
  * [23] L. I. Hedberg, Removable singularities and condenser capacities, Arkiv för Matematik, 12 (1974), pp. 181–201.
  * [24] H. Hedenmalm, B. Korenblum, and K. Zhu, Theory of Bergman spaces, vol. 199 of Graduate Texts in Mathematics, Springer-Verlag, New York, 2000.
  * [25] F. Holland and D. Walsh, Growth estimates for functions in the Besov spaces, in Proceedings of the Royal Irish Academy. Section A: Mathematical and Physical Sciences, JSTOR, 1988, pp. 1–18.
  * [26] S. Khrushchev, A continuous function with universal Fourier series on a given closed set of Lebesgue measure zero, Journal of Approximation Theory, 252 (2020), p. 105361.
  * [27] S. V. Khrushchev, The problem of simultaneous approximation and of removal of the singularities of Cauchy type integrals, Trudy Matematicheskogo Instituta imeni VA Steklova, 130 (1978), pp. 124–195.
  * [28] A. Limani, Shift invariant subspaces in growth spaces and sets of finite entropy, arXiv preprint arXiv:2304.09081, (2023).
  * [29] A. Limani and B. Malman, On the problem of smooth approximations in $\mathcal{H}(b)$ and connections with subnormal operators, preprint, (2021).
  * [30] A. Limani and B. Malman, An abstract approach to approximation in spaces of pseudocontinuable functions, Proceedings of the American Mathematical Society, 150 (2022), pp. 2509–2519.
  * [31] , Constructions of some families of smooth Cauchy transforms, Canadian Journal of Mathematics, 76 (2024), pp. 319–344.
  * [32] A. Limani and A. Nicolau, Shift invariant subspaces in the Bloch space, (2023).
  * [33] N. G. Makarov, On a class of exceptional sets in the theory of conformal mappings, Matematicheskii Sbornik, 180 (1989), pp. 1171–1182.
  * [34] B. Malman, Thomson decompositions of measures in the disk, Transactions of the American Mathematical Society, 376 (2023), pp. 8529–8552.
  * [35] P.-O. Parisé, Densité de l’algebre du disque dans les espaces de de Branges–Rovnyak, (2019).
  * [36] K.-M. Perfekt, On m-ideals and oo type spaces, Mathematica scandinavica, (2017), pp. 151–160.
  * [37] C. Pommerenke, Boundary behaviour of conformal maps, vol. 299, Springer Science & Business Media, 2013.
  * [38] W. Rudin, Functional Analysis . MacGraw-Hill Book Company, International Edition, (1973).
  * [39] N. A. Shirokov, Zero sets for functions from $\Lambda_{\omega}$, Zapiski Nauchnykh Seminarov POMI, 107 (1982), pp. 178–188.
  * [40] W. Smith, Inner functions in the hyperbolic little Bloch class., Michigan Mathematical Journal, 45 (1998), pp. 103–114.
  * [41] B. A. Taylor and D. L. Williams, Ideals in rings of analytic functions with smooth boundary values, Canadian Journal of Mathematics, 22 (1970), pp. 1266–1283.
  * [42] J. E. Thomson, Approximation in the mean by polynomials, Annals of Mathematics, 133 (1991), pp. 477–507.