# Resilience of dynamical systems

Hana Krakovská  ,  Christian Kuehn  and  Iacopo P. Longo
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Slovak Academy of
Sciences, Institute of Measurement Science, Department of Theoretical Methods,
Dúbravská cesta 9, 841 04 Bratislava 4, Slovakia.  Technische Universität
München, Forschungseinheit Dynamics, Zentrum Mathematik, M8, Boltzmannstraße
3, 85748 Garching bei München, Germany.

###### Abstract.

Stability is among the most important concepts in dynamical systems. Local
stability is well-studied, whereas determining how “globally stable” a
nonlinear system is very challenging. Over the last few decades, many
different ideas have been developed to address this question, primarily driven
by concrete applications. In particular, several disciplines suggested a web
of concepts under the headline “resilience”. Unfortunately, there are many
different variants and explanations of resilience, and often the definitions
are left relatively vague, sometimes even deliberately. Yet, to allow for a
structural development of a mathematical theory of resilience that can be used
across different areas, one has to ensure precise starting definitions and
provide a mathematical comparison of different resilience measures. In this
work, we provide a systematic review of the most relevant indicators of
resilience in the context of continuous dynamical systems, grouped according
to their mathematical features. The indicators are also generalized to be
applicable to any attractor. These steps are important to ensure a more
reliable, quantitatively comparable and reproducible study of resilience in
dynamical systems. Furthermore, we also develop a new concept of resilience
against certain non-autonomous perturbations to demonstrate, how one can
naturally extend our framework. All the indicators are finally compared via
the analysis of a classic scalar model from population dynamics to show that
direct quantitative application-based comparisons are an immediate consequence
of a detailed mathematical analysis.

###### Key words and phrases:

Resilience, dynamical systems, attractors, stability.

###### 2020 Mathematics Subject Classification:

37C20 34D10 34D05 37G35 37N99

HK was supported by the Slovak Grant Agency for Science under the grant No.
2/0081/19. CK was partly supported by a Lichtenberg Professorship of the
VolkswagenStiftung. CK also acknowledges partial support of the EU within the
TiPES project funded the European Unions Horizon 2020 research and innovation
programme under grant agreement No. 820970. IPL was partly supported by
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No 754462, by MICIIN/FEDER project
RTI2018-096523-B-100, by TUM International Graduate School of Science and
Engineering (IGSSE) and by the University of Valladolid under project PIP-
TCESC-2020.

##  1\. Introduction

In 1973, Holling [28] prompted the mathematical community in ecology to
investigate a new concept of “nonlocal” stability to capture the “persistence
of relationships within a system, […] the ability to absorb changes of state
variables, driving variables, and parameters, and still persist”, and by doing
so, complement the classic Lyapunov stability. He named this property
_resilience_. The theory of dynamical systems still proved to be the most
promising theoretical framework to treat the question rigorously (see Gruemm
[22]). Nevertheless, the newly introduced concept elicited a tremendous, and
yet mostly disorganized, research effort in several areas of applied science.
The result is a plethora of different indicators of resilience, in setups
which are often not directly comparable (Meyer [49]).

It is possible to categorize most of the available indicators within two main
branches: _engineering resilience_ and _ecological resilience_ (Peterson
[58]). The first aims at answering the question “how long will it take for a
system that suffers a perturbation to readjust to its original state?” and can
be traced back to Maynard Smith [70] and to May [44]. The second aims at
answering the question “which are the features of the maximal perturbations
that a system can endeavour and still return to its original state?” and it is
directly related to Holling [28]. The indicators treated in this paper are
listed in Table 1 according to this classic subdivision.

Engineering resilience | Ecological resilience
---|---
Characteristic return time§ |  $T_{R}$ |  3.1 |  Latitude in width§ |  $L_{w}$ |  4.1
Reactivity§ |  $R_{0}$ |  3.3 |  Distance to threshold§ |  $DT$ |  4.3
Max amplification§ |  $\rho_{max}$ |  3.6 |  Precariousness†‡ |  $P_{\mathcal{A}}$ |  4.7
Max amplification time§ |  $t_{max}$ |  3.6 |  Latitude in volume§ |  $L_{v}$ |  4.9
Stochastic invariability† |  $\mathcal{I}_{S}$ |  3.8 |  Basin stability§ |  $S_{\mathcal{B}}$ |  4.11
Deterministic invariability† |  $\mathcal{I}_{D}$ |  3.10 |  Resistance (potential)§ |  $W$ |  5.3
Return Time§ |  $T_{R}$ |  5.1 |  Flow-kick resilience§† |  – |  5.5
Resistance (Harrison)†‡ |  $R$ |  6.2 |  Intensity of attraction† |  $\mathcal{I}$ |  5.6
Elasticity†‡ |  $E$ |  6.3 |  Expected escape times† |  $\tau$ |  5.7
|  |  |  Persistence†‡ |  $P$ |  6.5
|  |  |  Rate-induced tipping†‡ |  – |  6.7
Table 1. Summary of the indicators reviewed in this work according to the
ample categories of engineering and ecological resilience (see Peterson [58]).
The presentation in the paper privileges a subdivision based on the
mathematical techniques which they entail. The notation and number of
definition are displayed next to the name of each indicator. The superscripts
indicate the type of perturbation for which each indicator is designed: ($\S$)
perturbations of initial conditions, ($\dagger$) time-dependent perturbations,
($\ddagger$) perturbations of parameters.

In this paper we have three main objectives. Firstly, we carry out an
extensive review of indicators of resilience within the framework of
continuous dynamical systems theory. While doing so, we provide a systematic
organization based on the mathematical techniques which they entail. To our
knowledge, such a comprehensive effort does not exist in the literature, nor
for dimension, nor in the organizational spirit. As an immediate effect, this
paper provides an easier access to a very vast and fragmented literature.
Secondly, we generalize these definitions (wherever possible) to general
attractors, whereas they are usually designed only for stable equilibria and
sometimes periodic orbits. However, equilibria and periodic orbits, while
playing a fundamental role in the understanding of many phenomena, are far
from exhausting the dynamical possibilities in higher dimensional and complex
systems. Therefore, we believe that our effort largely increases the
applicability of this theory. In particular, this applies to the local theory
based on the first order approximation of the dynamics near a hyperbolic
trajectory. In this context, it is natural to use of the notion of exponential
dichotomy from the theory of nonautonomous linear systems. Incidentally, we
note that this fact encourages the investigation of resilience also in
nonautonomous nonlinear dynamical systems, although, for the sake of
simplicity, we limit our presentation to the autonomous case. Lastly, we
enrich the discussion by bringing together recent accomplishments in the
theory of non-autonomous/rate-induced tipping phenomena, where the time-
dependent variation of parameters of some systems has to be generalized beyond
classical autonomous bifurcations to the non-autonomous setting.

The paper is organized as follows. In Section 2, we set the notation and
recall some notions on continuous dynamical systems. Each of the following
four sections addresses a class of indicators of resilience grouped according
to the most relevant mathematical technique upon which they are based: Section
3 contains the local indicators obtained via the linearization in the
neighbourhood of a hyperbolic trajectory. Section 4 deals with the indicators
focusing on geometric features of the basins of attraction. In Section 5, we
present all the indicators related to the transient behaviour of a system
during and after a perturbation. The variation of parameters and the
consequent indicators are treated in Section 6. Finally, in Section 7, we
carry out an analysis of resilience of a model from population dynamics, using
and comparing the indicators introduced in the previous sections.

##  2\. Notation, assumptions and preliminary definitions

Although the content of this paper applies to any metric space $(X,d)$, for
the sake of clarity and accessibility for applications, we will restrict to
the $N$-dimensional Euclidean space $\mathbb{R}^{N}$ with norm $|\cdot|$. When
$N=1$, we will simply write $\mathbb{R}$ and the symbol $\mathbb{R}^{+}$ will
denote the set of nonnegative real numbers. The symbol $\mathbb{R}^{N\times
M}$ represents the set of matrices with $N$ rows and $M$ columns, and given
$A\in\mathbb{R}^{N\times M}$, $A^{\top}$ will denote its transpose. As a rule,
a singleton $\\{x\\}$ in $\mathbb{R}^{N}$ will be identified with the element
$x$ itself and, given $E\subset\mathbb{R}^{N}$, the symbol $\overline{E}$ will
denote the closure of $E$ with respect to the topology induced by the
Euclidean norm. Moreover, by $B_{r}(p)$, we denote the open ball of
$\mathbb{R}^{N}$ centered at $p$ and with radius $r$. If $p$ is the origin of
$\mathbb{R}^{N}$, we simplify the notation to $B_{r}$. On the other hand,
given any set $E\subset\mathbb{R}^{N}$, $B_{r}(E)$ will denote the set of
points $x\in\mathbb{R}^{N}$ such that $x\in B_{r}(p)$ for some $p\in E$.
Furthermore, for any $U\subseteq\mathbb{R}^{M}$ and any
$W\subset\mathbb{R}^{N}$, $\mathcal{C}(U,W)$ will denote the space of
continuous functions from $U$ to $W$ endowed with the usual supremum norm
$\|\cdot\|_{\mathcal{C}(U,W)}$.

Unless otherwise stated, we will deal with a continuous dynamical system on
$\mathbb{R}^{N}$ which we will identify with its continuous (local) _flow_

$\phi:\mathcal{U}\subset\mathbb{R}\times\mathbb{R}^{N}\to\mathbb{R}^{N},\quad(t,x_{0})\mapsto\phi(t,x_{0}).$

Moreover, given $E\subset\mathbb{R}^{N}$ and $t>0$, we use $\phi(t,E)$ to
denote the set

$\phi(t,E)=\big{\\{}y\in\mathbb{R}^{N}\mid y=\phi(t,x)\text{ for some }x\in
E\big{\\}}.$

In particular, we recall that a set $E\subset\mathbb{R}^{N}$ is called
_forward invariant_ under the flow $\phi$ if $\phi(t,E)\subset E$ for all
$t>0$ and _invariant_ if $\phi(t,E)=E$ for all $t\in\mathbb{R}$.

Where necessary, we will assume that the continuous flow $\phi$ is induced by
an autonomous ordinary differential equation of the form

$\dot{x}=f(x),\quad x\in\mathbb{R}^{N},$ (2.1)

where $f\colon\mathbb{R}^{N}\to\mathbb{R}^{N}$ is regular enough so that the
initial value problem $x(0)=x_{0}\in\mathbb{R}^{N}$ admits a unique solution
$x(\cdot,x_{0})\in C(I_{x_{0}},\mathbb{R}^{N})$, with $I_{x_{0}}$ being its
maximal interval of definition. It is well-known that (2.1) induces a
continuous (local) flow on $\mathbb{R}^{N}$ via the relation

$\phi(t,x_{0})=x(t,x_{0}).$

We will also assume that the considered dynamical system has a local attractor
$\mathcal{A}\subset\mathbb{R}^{N}$ in the sense of the definition below (which
is taken from Kloeden and Rasmussen [35]).

######  Definition 2.1 (Local attractor and basin of attraction).

A compact set $\mathcal{A}\subset\mathbb{R}^{N}$ invariant under a flow $\phi$
on $\mathbb{R}^{N}$ is called a _local attractor_ if there exists $\eta>0$
such that for all $x_{0}\in B_{\eta}(\mathcal{A})$, $x(\cdot,x_{0})$ is
defined for all $t>0$ and

$\lim_{t\to\infty}\mathrm{dist}\big{(}\phi(t,B_{\eta}(\mathcal{A})),\mathcal{A}\big{)}=0,$

where for all $A,B\subset\mathbb{R}^{N}$, $\mathrm{dist}(A,B)$ denotes the
Hausdorff semi-distance between $A$ and $B$. We call the basin of attraction
of $\mathcal{A}$ the set

$\mathcal{B(A)}=\\{x_{0}\in\mathbb{R}^{N}:\emptyset\neq\omega(x_{0})\subset\mathcal{A}\\}.$

where, given $E\subset\mathbb{R}^{N}$, $\omega(E)$ represents the _omega limit
set of $E$_ under the flow $\phi$, i.e. the set

$\omega(E):=\bigcap_{t\geq 0}\overline{\bigcup_{s\geq t}\phi(s,E)}.$

######  Remark 2.2.

There is not a complete agreement in the literature about the definition of
local attractor. For example, the above definition is sometimes completed by
the condition of minimality: “there is no proper subset of $\mathcal{A}$ which
satisfies the properties in Definition 2.1”. In order to understand the
implication of such condition, consider the two-dimensional differential
system in polar coordinates $(r,\varphi)\in[0,\infty)\times[0,2\pi)$ where $0$
is identified with $2\pi$,

$\dot{r}=r(r-1)(r-3),\qquad\dot{\varphi}=1.$ (2.2)

Figure 1. Sketch of the phase space for (2.2). According to Definition 2.1,
one can choose either the whole closed ball of radius one, the origin and the
periodic orbit of radius one, or just the latter as a local attractor for the
induced dynamical system. Under the minimality condition in Remark 2.2, only
the periodic orbit of radius one qualifies for being a local attractor and the
origin belongs to the boundary of the basin of attraction.

It is easy to show that (2.2) has an unstable equilibrium at the origin and
two hyperbolic periodic orbits of radius $r=1$ (stable) and radius $r=3$
(unstable), respectively (see Figure 1 for a sketch of the phase space). If
one employs the minimality condition, the local attractor of the system
consists of the points in the periodic orbit of radius $r=1$ and its basin of
attraction is the open ball of radius $r=3$ with the exception of the origin.
In spite of belonging to the boundary of the basin of attraction, the origin
is merely an isolated point. Every solution starting in a sufficiently small
ball centered at the origin, except for the origin itself, will eventually
converge to the periodic orbit of radius $r=1$. If the minimality condition is
not taken into account, we are allowed to choose as a local attractor the
subset of $\mathbb{R}^{2}$ made of the points in the periodic orbit of radius
$r=1$ and also the origin. Its basin of attraction is now the whole open ball
of radius $r=3$. It is worth noting that an alternative approach is given by
the Milnor’s definition of local attractor [52], where the asymptotic
behaviour of solutions starting in negligible sets (with respect to the
Lebesgue measure in this context) are disregarded. In this case, the local
attractor would be made of the points in the sole periodic orbit of radius
$r=1$. In order to avoid measure theoretic arguments as much as possible, and
in view of the nonlocal nature of the concept of ecological resilience, in
this work we privilege the flexibility provided by Definition 2.1. Note also
that, although in this example the boundary of the basin is an unstable
periodic orbit, this is not always the case (see Figure 2).

\begin{overpic}[width=281.85034pt]{Duffing.pdf} \put(97.0,40.0){$x$}
\put(55.0,63.0){$y$} \end{overpic} Figure 2. The unforced Duffing oscillator
is a classic example given by the equations $\dot{x}=y$,
$\dot{y}=x-x^{3}-\delta y$, where $\delta\geq 0$. The system has three
equilibria, namely $(0,0)$ and $(\pm 1,0)$. Easy calculations show that for
$\delta>0$, the origin is a saddle and the fixed points at $(\pm 1,0)$ are
asymptotically stable. In particular, the boundary of the basin of attraction
of each one of the latter is given by the stable manifold of the saddle at the
origin.

######  Definition 2.3 (Global Attractor).

A compact set $\mathcal{A}\subset\mathbb{R}^{N}$ invariant under the flow
$\phi$ induced by (2.1) is called a _global attractor_ if $\mathcal{A}$
attracts each bounded subset of $\mathbb{R}^{N}$.

It is easy to show that the global attractor for $\phi$ is the minimal compact
set that attracts each bounded subset of $\mathbb{R}^{N}$, and also that it is
the maximal closed and bounded invariant set (see Lemma 1.6 in [8]). As a
matter of fact, the global attractor is composed of all the points of
$\mathbb{R}^{N}$ which belong to a bounded global solution for (2.1) (see
Theorem 1.7 in [8]).

We also recall that, under the considered assumptions, the basin of attraction
of a local attractor is an open set. The proof is a consequence of simple
topological arguments and it is therefore omitted.

######  Proposition 2.4.

Assume the evolution operator $\phi$ is continuous. The basin of attraction
$\mathcal{B(A)}$ of a local attractor $\mathcal{A}$ is an open set.

We conclude this section by recalling the definitions of exponential dichotomy
and of hyperbolic solutions.

######  Definition 2.5 (Exponential dichotomy).

Given a locally integrable function $A:\mathbb{R}\to\mathbb{R}^{N\times N}$,
the linear system

$\dot{y}=A(t)y$ (2.3)

is said to have an _exponential dichotomy_ on an interval $I\subset\mathbb{R}$
if there are a projection $P$ (i.e. a matrix $P\in\mathbb{R}^{N\times N}$
satisfying $P^{2}=P$), and constants $\alpha>0$, $K\geq 1$ such that

$\begin{split}|Y(t)PY^{-1}(s)|\leq Ke^{-\alpha(t-s)},\quad\text{for all
}t,s\in I,\,t\geq s,&\text{ and}\\\ |Y(t)(Id-P)Y^{-1}(s)|\leq
Ke^{\alpha(t-s)},\quad\text{for all }t,s\in I,\,t\leq s,&\end{split}$

where $Y:\mathbb{R}\to\mathbb{R}^{N\times N}$ is a fundamental matrix solution
of (2.3), and Id is the identity matrix on $\mathbb{R}^{N\times N}$.

The exponential dichotomy on $\mathbb{R}$ fulfils the role that eigenvalues
with nonzero real part play in the study of stability of hyperbolic equilibria
or periodic orbits. However, it also allows to treat more complex trajectories
(if they exist) which have the equivalent role of determining the asymptotic
behaviour of a system but involve a nontrivial dependence on time. For the
sake of consistency and to avoid introducing further notation, we will present
the notion of hyperbolic solutions only for autonomous ordinary differential
equations although it is defined for general nonautonomous systems [35].

######  Definition 2.6 (Hyperbolic solution).

Consider a continuous dynamical system induced by an ordinary differential
equation $\dot{x}=f(x)$, with $x\in\mathbb{R}^{N}$ and a globally defined
solution $\widetilde{x}:\mathbb{R}\to\mathbb{R}^{N}$ such that $f$ is
continuously differentiable on $\widetilde{x}(t)$ for all $t\in\mathbb{R}$.
Then, $\widetilde{x}$ is said to be hyperbolic if the corresponding
variational equation
$\dot{y}={\textnormal{D}}f\big{(}\widetilde{x}(t)\big{)}y$ has an _exponential
dichotomy_ on $\mathbb{R}$. In particular, we will call an hyperbolic solution
$\widetilde{x}$ _locally attractive_ if the exponential dichotomy has
projector the identity $P={\textnormal{Id}}$.

######  Remark 2.7.

The reader who is interested in a deeper understanding of the extent of such
generalization, can refer for example to [31, Section 1.4] and the references
cited therein (see for example the combined role of [31, Proposition 1.56] and
[24, Theorem III.2.4]). On the other hand, the reader who deals with
attractors made of equilibria and/or periodic orbits can intend the term
hyperbolic in the classic sense.

##  3\. Local Indicators

Linear stability analysis is one of the classic tools in the study of
dynamical systems. It allows to infer the asymptotic dynamics of a system in
the surrounding of a reference trajectory by looking at its linear
approximation and the dominant Lyapunov exponent of a fundamental matrix
solution. If the dominant Lyapunov exponent’s real part is nonzero, the sign
qualifies the reference trajectory as stable or unstable. Its absolute value
measures the asymptotic speed of convergence or divergence after a small
perturbation. When the reference trajectory is a stationary state or a
periodic orbit, this reduces to the calculation of the eigenvalues of the
Jacobian of the initial vector field evaluated on such orbits.

Such ideas have been used to determine the resilience of an attractor via the
rate of convergence of the nearby solutions. Due to the inherent local nature
of linear stability analysis, these indicators of resilience overlook the
topological structure of the phase space away from the considered attractor
and are not designed to determine the highest possible perturbation which a
system can absorb before tipping to a different state.

The section contains five subsections and six indicators. In subsection 3.1 we
present the classic _characteristic return time_ which relates the resilience
of a system in the nearby of an attractor to the asymptotic rate of
convergence of the solutions. Subsections 3.2 and 3.3 address respectively the
question of local resilience in the short-term horizon and in the transient
after a perturbation of the initial conditions. The first contains the
indicator _reactivity_ whereas the latter the indicators _maximal
amplification_ and _maximal amplification time_. Subsection 3.4 addresses the
question of local resilience against time-dependent (random or deterministic)
perturbations and contains the indicators _intrinsic stochastic variability_
and _intrinsic deterministic variability_. Finally, subsection 3.5 contains a
short discussion of the relations among the previously introduced indicators.

Besides the assumptions in Section 2, we shall also consider the following
assumption.

  1. (H1)

The function $f:\mathbb{R}^{N}\to\mathbb{R}^{N}$ in (2.1) is assumed to be
continuously differentiable and, for every $x\in\mathbb{R}^{N}$,
${\textnormal{D}}f(x)$ will denote the Jacobian of $f$ calculated at $x$.
Moreover, assume that $\mathcal{A}=\\{\widetilde{x}(t)\mid t\in\mathbb{R}\\}$,
and $\widetilde{x}$ is a locally attractive hyperbolic solution for
$\dot{x}=f(x)$.

###  3.1. Characteristic Return Time

The notion of characteristic return time in the context of resilience for
ecological systems already dates back to May [44]. A very commonly used
version is due to Beddington et al. [4] for discrete dynamical systems and
then to Pimm and Lawton [59] for the continuous case. It is indistinctly used
under different names, e.g., return time [59], characteristic return time
[61], engineering resilience [23, 29, 71], resilience [55]. A qualitative
description of the underlying idea is provided by Pimm [61] as “how fast the
variables return towards their equilibrium following a perturbation”, or, more
specifically, as the “time taken for a perturbation to return to $1/e$ of its
initial value”. The definition is motivated by the fact that a trajectory
starting in the nearby of a locally stable equilibrium $x^{*}$ will approach
it in a time which is proportional to the reciprocal of the eigenvalue with
largest real part for the system linearization at $x^{*}$. The presentation
below is however given for the more general case of a locally attractive
hyperbolic solution (see Definition 2.6).

######  Definition 3.1 (Characteristic return time).

Consider $f:\mathbb{R}^{N}\to\mathbb{R}^{N}$ satisfying assumption (H1). The
characteristic return time $T_{R}$ of the system $\dot{x}=f(x)$ for the
attractor $\mathcal{A}$ is defined as

$T_{R}(\mathcal{A})=\frac{1}{\widehat{\alpha}},$

where

$\widehat{\alpha}:=\inf\big{\\{}\alpha>0\ \mid\ |Y(t)Y^{-1}(s)|\leq
Ke^{-\alpha(t-s)},\quad\text{for all }t\geq s\big{\\}},$

and $Y:\mathbb{R}\to\mathbb{R}^{N\times N}$ is a fundamental matrix solution
of $\dot{y}={\textnormal{D}}f\big{(}\widetilde{x}(t)\big{)}y$.

######  Remark 3.2.

If the considered local attractor $\mathcal{A}=x^{*}$ is an attractive
hyperbolic fixed point, then $\widehat{\alpha}$ coincides with the opposite of
the real part of the dominant eigenvalue of ${\textnormal{D}}f(x^{*})$, i.e.
$\widehat{\alpha}=-Re\big{(}\lambda_{dom}({\textnormal{D}}f(x^{*}))\big{)}$ .

The definition of characteristic return time motivated the introduction of the
following indicator of resilience for a stable hyperbolic equilibrium as the
rate of decay

$EV(\mathcal{A})=\widehat{\alpha}.$

This indicator has been widely used and studied for both continuous and
discrete systems in the case that $\mathcal{A}=x^{*}$ is hyperbolic and
attracting (see for example Arnoldi et al. [2], DeAngelis et al. [14, 15],
Harwell and Ragsdale [27], Pimm and Lawton [60], Rooney et al. [64], Van Nes
and Scheffer [73], Vincent and Anderson [74]).

###  3.2. Reactivity

Neubert and Caswell [55] proposed and studied different indicators with the
aim of capturing the transient behavior of a trajectory starting in the
neighbourhood of a stable equilibrium as it may substantially differ from the
asymptotic one. Specifically, the reactivity corresponds to the maximum
instantaneous rate at which an asymptotically stable linear homogeneous system
responds if the initial condition is taken outside the origin. The reactivity
of a nonlinear system $\dot{x}=f(x)$ in the neighbourhood of a stable
hyperbolic equilibrium $x^{*}$ is obtained through its linearization at
$x^{*}$.

######  Definition 3.3 (Reactivity).

Let $\dot{y}=A(t)y$, with $A:\mathbb{R}\to\mathbb{R}^{N\times N}$ locally
integrable, be an asymptotically stable linear homogeneous system, and denote
by $y(\cdot,t_{0},y_{0})$ its unique solution satisfying
$y(t_{0},t_{0},y_{0})=y_{0}$. We shall call the reactivity of the system at
time $t_{0}\in\mathbb{R}$, the quantity

$R_{t_{0}}=\max_{|y_{0}|\neq
0}\left(\frac{1}{|y(t,t_{0},y_{0})|}\frac{d|y(t,t_{0},y_{0})|}{dt}\right)\bigg{|}_{t=t_{0}}.$
(3.1)

The system $\dot{y}=A(t)y$ is called _reactive_ if there is
$t_{0}\in\mathbb{R}$ such that $R_{t_{0}}>0$ and nonreactive otherwise. A
nonlinear system $\dot{x}=f(x)$ satisfying Assumption (H1) is called reactive
if there exists a neighbourhood of a locally attractive hyperbolic solution
$\widetilde{x}$ such that
$\dot{y}={\textnormal{D}}f\big{(}\widetilde{x}(t)\big{)}y$ is reactive.

If a system is reactive in a neighbourhood of a locally attractive hyperbolic
solution $\widetilde{x}$, some trajectories starting in a neighborhood of
$\widetilde{x}$ may initially move away from $\widetilde{x}$, before
converging to it. In other words, the finite-time Lyapunov exponents for
$\widetilde{x}$ can be positive. When $\mathcal{A}=x^{*}$, a sufficient
condition preventing this scenario is that ${\textnormal{D}}f(x^{*})$ has
orthogonal eigenvectors. In such a case the eigenvalues of
${\textnormal{D}}f(x^{*})$ determine both the asymptotic and the transient
behaviour of the system. Neubert and Caswell [55] show that the reactivity of
a linear homogeneous system $\dot{y}=Ay$, $A\in\mathbb{R}^{N\times N}$, is
equal to the dominant eigenvalue of the symmetric part of the Toeplitz
decomposition of $A$ (notably, recall that every real symmetric matrix is
Hermitian and therefore its eigenvalues are real),

$R_{0}=\lambda_{dom}\big{(}H(A)\big{)}\in\mathbb{R},\qquad\text{where
}H(A)=\frac{A+A^{\top}}{2}.$

######  Example 3.4.

In this example we compare three planar asymptotically stable linear systems
$\dot{x}=A_{i}x,$ where $i=\\{1,2,3\\}:$

$A_{1}=\begin{pmatrix}-2&0\\\ 5&-1\end{pmatrix},\ \
A_{2}=\begin{pmatrix}-2&1\\\
\sqrt{26}+\sqrt{50}&-\sqrt{26}-\sqrt{50}-1\end{pmatrix},\ \
A_{3}=\begin{pmatrix}-2&0\\\ 1&-1\end{pmatrix}.$

All three systems have the same value of the dominant eigenvalue and
consequently the same characteristic return time given as $T_{R}(0)=1$. The
systems $A_{1}$ and $A_{2}$ are reactive with a positive value of reactivity
given by $R_{0}=(\sqrt{26}-3)/2\approx 1.05.$ The system $A_{3}$ is
nonreactive with $R_{0}=(\sqrt{2}-3)/2\approx-0.79.$ In the first row of
Figure 3, the plots show the time evolution of the magnitudes of trajectories
starting from different perturbed initial conditions around the origin.
Whereas the trajectories of the reactive systems $A_{1}$ and $A_{2}$ can
exhibit a transient growth in their magnitude, this is not possible for the
nonreactive system $A_{3}$. Furthermore, this example shows that despite
having the same characteristic return time and reactivity indicator values,
the transient behaviour of the trajectories in systems $A_{1}$ and $A_{2}$ may
vary substantially. This motivates the introduction of the so-called
_amplification envelope_ (see Subsection 3.3). The amplification envelopes of
three considered systems are plotted in the second row of Figure 3.

(a) $A_{1}$

(b) $A_{2}$

(c) $A_{3}$

(d) $A_{1}$

(e) $A_{2}$

(f) $A_{3}$

Figure 3. The transient behaviour of three systems $A_{1},A_{2},A_{3}$ from
Example 3.4 is investigated. In the first row, the time evolution of the
magnitude of four trajectories starting from different initial conditions
around the origin is depicted. In the second row, we see the respective
amplification envelopes of each system. Even though the systems have the same
characteristic return time, and additionally, the systems $A_{1}$ and $A_{2}$
have also the same value of reactivity, the indicator amplification envelope
and its characteristics $\rho_{max}$ and $t_{max}$ are able to capture
different transient behaviour.

###  3.3. Maximal amplification and maximal amplification time

Example 3.4 shows that reactive systems with same reactivity and
characteristic return time do not necessarily share the same transient
behaviour. This fact led Neubert and Caswell [55] to the definition of
amplification envelope, which records the maximal deviation from the attractor
for trajectories starting at perturbed initial conditions of fixed norm.

######  Definition 3.5 (Amplification envelope).

The amplification envelope for an asymptotically stable linear system
$\dot{y}=A(t)y$, with $A:\mathbb{R}\to\mathbb{R}^{N\times N}$ locally
integrable, is defined as the continuous function

$\displaystyle\rho:\mathbb{R}^{+}\times\mathbb{R}\to\mathbb{R}^{+},\quad(t,t_{0})\mapsto\rho(t,t_{0})=\sup_{|y_{0}|\neq
0}\frac{|y(t+t_{0},t_{0},y_{0})|}{|y_{0}|},$ (3.2)

where $y(\cdot,t_{0},y_{0})$ is the unique solution of $\dot{y}=A(t)y$
satisfying $y(t_{0},t_{0},y_{0})=y_{0}$. Therefore, $\rho(t,t_{0})$ is in fact
equal to the operator norm (induced by the Euclidean norm) of the principal
matrix solution $Y(t+t_{0},t_{0})$ of $\dot{y}=A(t)y$ for the initial time
$t_{0}$, that is

$\rho(t,t_{0})=\big{\|}Y(t+t_{0},t_{0})\big{\|}_{op}.$

For a nonlinear system $\dot{x}=f(x)$ such that $f$ satisfies Assumption (H1),
the amplification envelope in a neighbourhood of a hyperbolic solution
$\widetilde{x}$ is defined as the amplification envelope of the linear system
$\dot{y}={\textnormal{D}}f(\widetilde{x}(t))y$.

Note in particular that if $A(t)=A\in\mathbb{R}^{N\times N}$ for all
$t\in\mathbb{R}$, the amplification envelope depends only on the variable
$t\in\mathbb{R}^{+}$, i.e.

$\rho(t,t_{0})=\rho(t)=\big{\|}e^{At}\big{\|}_{op}.$

From the amplification envelope we can derive the two following indicators of
resilience. The maximal amplification corresponds to the maximal magnification
over time for trajectories starting in the nearby of an attractor, whereas the
maximal amplification time records the occurrence of the maximal
amplification.

######  Definition 3.6 (Maximal amplification and maximal amplification time).

Consider the setting of Definition 3.5. The maximal amplification and maximal
amplification time are respectively defined as:

$\rho_{max}(t_{0})=\max_{t\geq 0}\rho(t,t_{0}),\qquad
t_{max}(t_{0})=\operatorname*{argmax}_{t\geq 0}\rho(t,t_{0}).$

######  Remark 3.7.

The characteristic return time (see Definition 3.1) and the Reactivity (see
Definition 3.3) of a linear system are also calculable from the amplification
envelope. Indeed, for every fixed $t_{0}\in\mathbb{R}$, the first corresponds
to the slope of $ln(\rho(\widetilde{x},t,t_{0}))$ as $t\to+\infty$, while the
second is equal to the slope of $ln(\rho(\widetilde{x},t,t_{0}))$ as $t\to
t_{0}$.

###  3.4. Intrinsic stochastic and deterministic invariability

The amplification envelope provides an effective tool to describe the local
transient behaviour of a system close to an attractor if it is affected by an
“isolated and impulsive” perturbation. Arnoldi et al. [2] investigate the same
transient behaviour when the system is subjected to a time-dependent (random
or deterministic) forcing. In the spirit of this paper, we generalize the
presentation of [2] to hyperbolic solutions (see Definition 2.6) up to the
point where this makes sense.

Consider a continuous and bounded function $A:\mathbb{R}\to\mathbb{R}^{N\times
N}$ such that the linear homogeneous system

$\dot{y}=A(t)y$ (3.3)

has an exponential dichotomy on $\mathbb{R}$ with projector the identity (see
Definition 2.5), i.e. denoted by $U(t,s)$ the principal matrix solution of
(3.3) at time $s\in\mathbb{R}$, one has that there are constants $\alpha>0$
and $K\geq 1$ such that

$|U(t,s)|\leq Ke^{-\alpha(t-s)},\quad\text{for all }t\geq s.$ (3.4)

Firstly, consider the stochastic differential problem with additive white
noise

$dy=A(t)y\,dt+S\,dW(t)$ (3.5)

where $S\in\mathbb{R}^{N\times M}$, and $W=(W_{1},\dots W_{M})^{\top}$ is a
vector of $M$ independent Brownian motions. As $t\to\infty$ the distribution
of each strong solution of (3.5) converges to a stationary Gaussian
distribution centered at the origin. Moreover, the covariance matrix of the
system at time $t\in\mathbb{R}$ and initial data at $t_{0}\in\mathbb{R}$ is
given by

$C(t,t_{0})=\int_{t_{0}}^{t}U(t,s)\Sigma U(t,s)^{\top}\,ds,$ (3.6)

where $\Sigma:=SS^{\top}$. Notice also that the differentiation of
$C(t,t_{0})$ with respect to $t$ shows that it solves the matrix differential
equation

$\frac{d}{dt}C=A(t)C+CA(t)^{\top}+\Sigma.$ (3.7)

In particular, $C_{*}(\Sigma):=C(t,-\infty)$ is the only bounded solution of
(3.7) over the whole real line, and (3.4) and (3.6) together imply that (3.7)
has a local attractor (see for example [35, Theorem 1.23]) that must thus
coincide with $C_{*}(\Sigma)$. $C_{*}(\Sigma)$ shall be called the _stationary
covariance_ of the system.

In order to construct an indicator of intrinsic resilience against stochastic
perturbations, one looks at the largest stationary response among the possible
perturbations of given “magnitude”.

######  Definition 3.8 (Intrinsic stochastic invariability).

Consider a continuous dynamical system induced by an ordinary differential
equation $\dot{x}=f(x)$, $x\in\mathbb{R}^{N}$, satisfying (H1) and the
assumptions in Section 2. Moreover, assume that $\widetilde{x}$ is a locally
attractive hyperbolic solution and consider the local attractor
$\mathcal{A}=\\{\widetilde{x}(t)\mid t\in\mathbb{R}\\}$. Fixed a norm
$\|\cdot\|$ on the matrix space $\mathbb{R}^{N\times N}$ (e.g. the operatorial
norm, the Froebenius norm), the _stochastic variability_
$\mathcal{V}_{S}(\mathcal{A})$ and the _intrinsic stochastic invariability_ of
the attractor $\mathcal{A}$ with respect to $\|\cdot\|$ are respectively
defined as

$\mathcal{V}_{S}(\mathcal{A})=\sup_{\Sigma\geq
0,\,\|\Sigma\|=1}\big{\\{}\|c\|\mid c\in
C_{*}(\Sigma)\big{\\}}\qquad\text{and}\qquad\mathcal{I}_{S}(\mathcal{A})=\frac{1}{2\mathcal{V}_{S}(\mathcal{A})},$

where $C_{*}(\Sigma)$ is the stationary covariance of (3.5) with
$A(t)={\textnormal{D}}f\big{(}\widetilde{x}(t)\big{)}$.

######  Remark 3.9.

If $\mathcal{A}=x^{*}$, then $A={\textnormal{D}}f(x^{*})$, and (3.7) does not
depend on time. In this case, one can show that the linear operator
$\widehat{A}(C):=AC+CA^{\top}$ has $N^{2}$ eigenvalues of the form
$\lambda_{j}+\lambda_{i}$, where $\lambda_{i},\lambda_{j}$ are eigenvalues of
$A$, and the stationary covariance matrix $C_{*}(\Sigma)$ is given by the
unique solution of the Lyapunov equation $AC+CA^{\top}+\Sigma=0$ [6, Lemma
5.1.2]. Therefore, one has that

$\mathcal{V}_{S}(\mathcal{A})=\sup_{\Sigma\geq
0,\,\|\Sigma\|=1}\|-\widehat{A}^{-1}(\Sigma)\|.$

On the other hand, one can consider a deterministic forcing of (3.3), that is,

$\dot{y}=A(t)y+g(t),$ (3.8)

where $g:\mathbb{R}\to\mathbb{R}^{N}$ is a bounded function. Denoted by $Y(t)$
a fundamental matrix solution of (3.3), note that

$x(t,g)=Y(t)\int_{-\infty}^{t}Y^{-1}(s)g(s)\,ds,\quad t\in\mathbb{R},$

is a solution of (3.8), which can be regarded as the _stationary system
response_ , and its _mean square deviation_ from the origin (which is the
global attractor of the homogeneous problem) is

$m(g):=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}|x(s,g)|^{2}\,ds.$

Arnoldi et al. [2] suggest to use the largest mean square deviation of
$x(t,g)$ from the origin over the possible perturbations $g$ of given norm, as
an indicator of resilience of the attractor.

######  Definition 3.10 (Intrinsic deterministic invariability).

Consider a continuous dynamical system induced by an ordinary differential
equation $\dot{x}=f(x)$, $x\in\mathbb{R}^{N}$, satisfying (H1) and the
assumptions in Section 2. Moreover, assume that $\widetilde{x}$ is a locally
attractive hyperbolic solution and consider the local attractor
$\mathcal{A}=\\{\widetilde{x}(t)\mid t\in\mathbb{R}\\}$. The _deterministic
variability_ $\mathcal{V}_{D}(\mathcal{A})$ and the _intrinsic deterministic
invariability_ of the attractor $\mathcal{A}$ are respectively defined as

$\mathcal{V}_{D}(\mathcal{A})=\sup_{\|g\|_{\infty}=1}2\sqrt{m(g)},\qquad\text{and}\qquad\mathcal{I}_{D}(\mathcal{A})=\frac{1}{\mathcal{V}_{D}(\mathcal{A})},$
(3.9)

where $m(g)$ is the mean square deviation of the stationary response of (3.8)
with $A(t)={\textnormal{D}}f\big{(}\widetilde{x}(t)\big{)}$.

######  Remark 3.11.

If $\mathcal{A}=x^{*}$, then $A={\textnormal{D}}f(x^{*})$, and one can use
standard frequency analysis to improve the information given by (3.9). This is
particularly true when one limits the possible deterministic perturbation to
the so-called “wide-sense stationary signals”. Indeed, since any deterministic
signal can be developed into a sum of harmonic terms, or Fourier modes, and
due to the fact that in the linear approximation, the system response to a
general perturbation is equal to the sum of the system response to the single-
frequency modes, a convexity argument yields that the perturbation generating
the largest system response is a single-frequency mode [2]. Therefore, when
$g$ is a single-frequency periodic forcing one obtains that

$\mathcal{V}_{D}(\mathcal{A})=\sup_{\omega\in\mathbb{R}}\|\textnormal{i}\omega-A^{-1}\|,$

where $\omega$ is the forcing frequency of $g$, i is the imaginary unit, and
$\|\cdot\|$ is the induced matrix operator norm from the norm on
$\mathbb{R}^{N}$.

When $\mathcal{A}=x^{*}$, an important result in [2] establishes a chain of
order for some of the indicators of local resilience presented in this
section.

######  Proposition 3.12.

If $\mathcal{A}=x^{*}$ is a hyperbolic stable equilibrium, the following chain
of inequalities holds true:

$-R_{0}\leq\mathcal{I}_{S}(\mathcal{A})\leq\mathcal{I}_{D}(\mathcal{A})\leq
EV(\mathcal{A}).$

###  3.5. Discussion

  * •

Local indicators are intrinsically robust against “small” perturbations. The
reason lies in the property of persistence of the exponential dichotomy which
goes under the name of _roughness_ (see Coppel [11]). Note that this is a
further reason justifying the notion of exponential dichotomy. Also in the
case in which our original system admits only equilibria a sufficiently small
but general time-dependent forcing perturbs the equilibria into trajectories
whose properties of stability persist but, in general, cannot be analyzed via
the sign of the eigenvalues of the variational problem anymore. In other
words, the Lyapunov spectrum of an equilibrium is perturbed into the dichotomy
spectrum [65].

  * •

Lundström [42] provides calculable definitions to approximate the recovery
rate (reciprocal of the characteristic return time) and the slowest return
time while at the same time estimating the distance to threshold (see
Definition 4.3) and the volume of the basin of attraction (see Definition
4.9).

  * •

Ives and Carpenter [30] propose the idea that the dynamics in the nearby of
the boundary of a basin of attraction can be important to infer a measure of
nonlocal stability for the attractor. Recalling the unforced Duffing
oscillator (see Figure 2) one can intuitively see that if a system lingers
close to the boundary after a first perturbation, it is possibly more
susceptible to tip once a new instantaneous perturbation takes place (see also
Definition 4.7). The linearization of the model along the boundary is
therefore suggested as a valuable tool. In practical terms, it is important to
point out that identifying the boundary may be as difficult, if not more
difficult, compared to identifying the attractor itself.

##  4\. Basin Shape Indicators

This section contains the core indicators of the classic research in
ecological resilience which dates back to the groundbreaking work by Holling
[28]. The underlying assumption is that the considered phenomenon, and its
mathematical model, admits multiple stable states and thus the phase space of
the induced dynamical system can be partitioned into basins of attraction.
Consequently, one aims to estimate the minimal perturbation of initial
conditions which can drive a system lying on an attractor, outside of its
basin of attraction. The study of the geometrical features of the basins
becomes the central focus, aiming to relate such characteristics to inherent
properties of stability of the respective attractors.

This section contains three subsections and five indicators. In subsection
4.1, we present _latitude in width_ , _distance to threshold_ and
_precariousness_. In subsection 4.2, we present _latitude in volume_ and
_basin stability_. In subsection 4.3 we briefly discuss the relations between
the indicators presented above.

All the indicators in this section are defined for a continuous dynamical
system on $\mathbb{R}^{N}$ with local flow

$\phi:\mathcal{U}\subset\mathbb{R}\times\mathbb{R}^{N}\to\mathbb{R}^{N},\quad(t,x_{0})\mapsto\phi(t,x_{0}).$

###  4.1. Latitude in width, distance to threshold and precariousness

The notion of “width” of a basin of attraction as an indicator of resilience
dates back directly to the seminal work by Holling [28]. Loosely speaking, the
width of a basin of attraction corresponds to the length of the “minimal”
segment crossing through the attractor and intersecting (at its extrema) the
boundary of the basin of attraction, if it applies. Although a precise
mathematical formulation seems hard to find in the literature, the geometrical
representations used in many works (for example Peterson et al. [58] and
Walker et al. [75]) permit to precisely formalize the underlying idea.

######  Definition 4.1 (Latitude in width).

Consider the (possibly empty) set

$S=\big{\\{}(y,z)\in\mathbb{R}^{2N}\,\big{|}\,y,z\in\partial\mathcal{B(A)},\text{
and }\alpha y+(1-\alpha)z\in\mathcal{A}\text{ for some
}\alpha\in(0,1)\big{\\}}.$ (4.1)

The latitude in width of an attractor $\mathcal{A}\subset\mathbb{R}^{N}$ for a
continuous flow $\phi:\mathbb{R}\times\mathbb{R}^{N}\to\mathbb{R}^{N}$ is
defined as

$L_{w}(\mathcal{A})=\begin{cases}\infty,&\text{if }S=\emptyset\\\
\inf_{(y,z)\in S}|y-z|,&\text{otherwise}.\end{cases}$

If $L_{w}(\mathcal{A})<\infty$, the inferior in the previous formula is in
fact a minimum as shown in the next result.

######  Proposition 4.2.

Under the notation and assumptions of Definition 4.1 and if $L_{w}<+\infty$,
we have:

$L_{w}(\mathcal{A})=\inf_{(y,z)\in S}|y-z|=\min_{(y,z)\in S}|y-z|.$

###### Proof.

By definition $\mathcal{A}$ is a compact subset of $\mathbb{R}^{N}$.
Therefore, there is $\rho>L_{w}(\mathcal{A})>0$ such that for all
$a_{1},a_{2}\in\mathcal{A}$, $|a_{1}-a_{2}|<\rho$. Moreover, notice that
$\partial\mathcal{B(A)}$ is nonempty since $L_{w}(\mathcal{A})<\rho<\infty$.
By the compactness of $\mathcal{A}$, there are
$a_{1},\dots,a_{n}\in\mathcal{A}$ such that

$\mathcal{A}\subset\bigcup_{i=1}^{n}\overline{B}_{\rho}(a_{i})=:P.$

In particular, notice that if $y,z\in S$ and $|y-z|<\rho$, then $(y,z)\in
P\times P$, which is a compact subset of $\mathbb{R}^{2N}$. Therefore, by the
continuity of the norm in $\mathbb{R}^{N}$ and Weierstrass Theorem, we
immediately obtain the result. ∎

Already in Peterson et al. [58], while still employing the idea of width of
the basin of attraction (under the name of ecological resilience), a focus was
aimed at understanding the smallest perturbation of initial conditions able to
drive the system away from the original attractor (or equivalently the largest
perturbation of initial conditions for which the systems returns to the
attractor). This approach, qualitatively, yet unequivocally presented by
Beisner et al. [5], emphasizes the importance of the points of minimal
distance between an attractor and the boundary of its basin of attraction (in
some cases called threshold or separatrix). The definition has the advantage
of revealing that an attractor with infinite latitude in width (see Definition
4.1) is not necessarily “very resilient” (see Figure 4). This idea has been
used consistently ever since for example in Lundström and Adainpää [43],
Lundström [42], Kerswell et al. [33], Klinshov et al. [34], Mitra et al. [54],
and Fassoni and Yang [19] where it is called precariousness (on the subject
see also Definition 4.7 and Proposition 4.8 below).

######  Definition 4.3 (Distance to threshold).

Consider the (possibly empty) set

$S=\\{(a,y)\in\mathbb{R}^{2N}\mid
a\in\mathcal{A},y\in\partial\mathcal{B(A)}\\}.$ (4.2)

The distance to threshold for an attractor $\mathcal{A}\subset\mathbb{R}^{N}$
for a continuous flow $\phi:\mathbb{R}\times\mathbb{R}^{N}\to\mathbb{R}^{N}$
is defined as

$DT(\mathcal{A})=\begin{cases}\infty,&\text{if }S=\emptyset\\\ \inf_{(a,y)\in
S}|a-y|,&\text{otherwise}.\end{cases}$

######  Remark 4.4.

Reasoning as for the proof of Proposition 4.2 one can easily show that if
$DT(\mathcal{A})<\infty$, then it is attained at a point of minimum. Moreover,
notice that $DT(\mathcal{A})=+\infty$ only if
$\partial\mathcal{B(A)}=\varnothing$ which means that $\mathcal{A}$ is the
global attractor (see Definition 2.3) for the system.

\begin{overpic}[trim=42.67912pt 284.52756pt 42.67912pt
284.52756pt,clip,width=433.62pt]{Fig4p.pdf} \put(50.0,0.0){$x_{1}$}
\put(4.0,27.0){$x_{2}$} \put(57.0,29.0){$X_{s}$} \put(11.0,17.0){$X_{0}$}
\put(74.0,31.5){$X_{1}$} \end{overpic} Figure 4. Representation of the
dynamics induced by the planar system $\dot{x}_{1}=-x_{1}+10x_{2}$,
$\dot{x}_{2}=x_{2}(10\exp(\frac{-x_{1}^{2}}{100})-x_{2})(x_{2}-1)$. The system
has three equilibria. The stable manifold of the saddle node in $X_{S}$ is the
separatrix between the basins of attraction of $X_{0}$ and $X_{1}$. Both
stable equilibria have an infinite latitude in width, while their distance to
threshold can be made as a small as wished through a suitable scaling. The
resilience of this system has been thoroughly analyzed by Kerswell et al.
[33].

When dealing with a concrete model of a real phenomenon, it is possible to
tune the previous definitions according to the available information on the
phenomenon. In particular, if one can anticipate the set of possible perturbed
initial conditions, the calculation of the latitude in width and the distance
to threshold can be restricted to a suitable subset $C$ of the phase space
called _region of interest_. For example, such an adjustment on the distance
to threshold has been explicitly presented in [49] where the set $S$ in (4.2)
is substituted by

$S_{C}=\\{(a,y)\in\mathbb{R}^{2N}\mid
a\in\mathcal{A},y\in\partial\mathcal{B(A)}\cap C\\},$ (4.3)

for some $C\subset\mathbb{R}^{N}$. Likewise, the latitude in width can be
modified with respect to the region of interest. Then, for some fixed
$C\subset\mathbb{R}^{N}$, the set $S$ (4.1) is changed for

$S_{C}=\\{(y,z):(y,z)\in S\text{ and }y\in C\\}.$

######  Example 4.5.

Figure 5. Sketch of the phase space for the dynamical system induced by (4.4)
and of the distance to threshold and latitude in width of the attractor
$\mathcal{A}=\\{0\\}$. Additionally, we can see that when we specify the
anticipated perturbations as $\mathbb{R}^{+},$ the distance to threshold
changes, whereas, the latitude in width remains the same (for details, see
Example 4.5).

$\displaystyle\dot{x}$
$\displaystyle=x(x-1/\varepsilon)(x+\varepsilon)\quad\quad\text{where }\quad
x\in\mathbb{R},\varepsilon\in(0,1]$ (4.4)

The dynamical system induced by (4.4) has a stable equilibrium at $x=0$, which
we consider as the relevant attractor, and two unstable fixed points at
$-\varepsilon$ and $1/\varepsilon.$ See Figure 5 for a sketch of the phase
space. Notice that when $\varepsilon$ becomes smaller, the distance to
threshold $DT(0)=\varepsilon$ decreases, but the latitude in width
$L_{w}(0)=1/\varepsilon+\varepsilon$ increases. If the set of anticipated
perturbed initial conditions $C=\mathbb{R}^{+}$ is considered, the distance to
threshold becomes $DT_{\mathbb{R}^{+}}(0)=1/\varepsilon,$ since we have the
set $S_{C}=\\{(0,1/\varepsilon)\\}$ (see (4.3)), while the latitude in width
$L_{w,\mathbb{R}^{+}}(0)$ remains the same. On the other hand, if we fix
$C=(0,\varepsilon),$ the set $S_{C}$ becomes empty for both distance to
threshold and latitude in width and therefore
$DT_{(0,\varepsilon)}(0)=L_{w,(0,\varepsilon)}(0)=+\infty.$

######  Example 4.6.

Let us consider the example presented in Remark 2.2 (see also Figure 6).
Depending on the different choice of the local attractor we may obtain
different values of the distance to threshold $DT.$ If we consider only the
limit cycle of radius $1$ as the attractor $\mathcal{A}_{1},$ the unstable
equilibrium at the origin is part of the basin’s boundary
$\partial\mathcal{B}(\mathcal{A}_{1}).$ Therefore, the distance to threshold
$DT(\mathcal{A}_{1})=1$. However, if we consider as a local attractor the set
of points in the limit cycle and the origin,
$\mathcal{A}_{2}=\mathcal{A}_{1}\cup\\{(0,0)\\},$ or the whole closed ball of
radius one $\mathcal{A}_{3}=B_{1}$, we obtain
$DT(\mathcal{A}_{2})=DT(\mathcal{A}_{3})=2.$ It is obvious that in terms of
actual resilience of the system’s state against a perturbation the origin
plays a negligible role.

Figure 6.  Sketch of the phase space for the dynamical system induced by
(2.2), and representation of the distance to threshold depending on the choice
of the local attractor (see Remark 2.2). For the local attractor
$\mathcal{A}_{1}$ consisting of the points in the limit cycle of radius $1$,
$DT(\mathcal{A}_{1})=1$. For the local attractor
$\mathcal{A}_{2}=\mathcal{A}_{1}\cup\\{(0,0)\\}$, $DT(\mathcal{A}_{2})=2.$ The
choice of the attractor $\mathcal{A}_{2}$ admits to consider only distances to
the “significant” parts of the boundary.

Instead of focusing on an attractor and its basin, in practice it is
convenient to look at the current state of the system within the basin. The
idea of studying the minimal distance from the boundary of the basin has
qualitatively presented by Walker et al. [75], where it is described as only
one of four components of ecological resilience together with the latitude in
width (see Definition 4.1), the resistance of a gradient system (see
Definition 5.3) and the cross-scale interaction among the previous components
and eventual subsystems. We push the idea in [75] a bit forward and propose a
mathematical definition which probes also the points outside the basin of
attraction. This is not only a reasonable generalization, as it also allows to
quantify the difficulty for a point outside the basin to enter it after a
given perturbation, but also permits to establish a connection with a recently
presented indicator (see Definition 5.5).

######  Definition 4.7 (Precariousness).

The function $P_{\mathcal{A}}:\mathbb{R}^{N}\to\mathbb{R}\cup\\{\infty\\}$
defined by

$x_{0}\mapsto P_{\mathcal{A}}(x_{0})=\begin{cases}\infty&\text{if
}\partial\mathcal{B(A)}=\varnothing,\\\
\displaystyle\inf_{y\in\partial\mathcal{B(A)}}|x_{0}-y|&\text{if
}\partial\mathcal{B(A)}\neq\varnothing\text{ and }x_{0}\in\mathcal{B(A)}\\\
\displaystyle-\inf_{y\in\partial\mathcal{B(A)}}|x_{0}-y|&\text{if
}\partial\mathcal{B(A)}\neq\varnothing\text{ and }x_{0}\notin\mathcal{B(A)}\\\
\end{cases}$

is called the _precariousness_ of the point $x_{0}\in\mathbb{R}^{N}$ with
respect to the basin of attraction $\mathcal{B(A)}$.

The following simple proposition establishes a relation between the distance
to threshold and the asymptotic value of precariousness.

######  Proposition 4.8.

Consider a continuous dynamical system $\phi$ on $\mathbb{R}^{N}$. If
$\partial\mathcal{B(A)}$ is nonempty, then

$\lim_{t\to\infty}P_{\mathcal{A}}\big{(}\phi(t,x_{0})\big{)}\geq
DT(\mathcal{A}),\quad\text{for all }x_{0}\in\mathcal{B(A)},$

and the equality holds if every trajectory in $\mathcal{A}$ is dense.

###### Proof.

The first inequality is a direct consequence of the definition of basin of
attraction (see Definition 2.1) and of distance to threshold (see Definition
4.3) i.e., for any $x_{0}\in\mathcal{B(A)}$,

$\begin{split}\lim_{t\to\infty}P_{\mathcal{A}}\big{(}\phi(t,x_{0})\big{)}&=\lim_{t\to\infty}\inf_{y\in\partial\mathcal{B(A)}}|\phi(t,x_{0})-y|=\lim_{t\to\infty}\inf_{\begin{subarray}{c}x\in\omega(x_{0}),\\\
y\in\partial\mathcal{B(A)}\end{subarray}}|\phi(t,x_{0})\pm x-y|\\\
&\geq\inf_{\begin{subarray}{c}x\in\omega(x_{0}),\\\
y\in\partial\mathcal{B(A)}\end{subarray}}|x-y|\geq\inf_{\begin{subarray}{c}a\in\mathcal{A},\\\
y\in\partial\mathcal{B(A)}\end{subarray}}|a-y|=DT(\mathcal{A}).\end{split}$

If all trajectories in $\mathcal{A}$ are dense, the equality is a direct
consequence of the fact that $\omega(x_{0})\subset\mathcal{A}$ is invariant
for the flow and thus it is dense in $\mathcal{A}$. ∎

###  4.2. Latitude in volume and Basin stability

For one-dimensional dynamical systems, the latitude in width of a basin of
attraction (see Definition 4.1) coincides, in fact, with its Lebesgue measure
$\mu(\mathcal{B(A)})$ (recall that a basin of attraction is always an open set
and therefore Lebesgue-measurable; see Proposition 2.4). The notion acquires
special interest when the study can be restricted to a measurable set
$C\subset\mathbb{R}^{N}$ of finite measure since one can interpret the
“relative volume” of the basin of attraction $\mathcal{B(A)}$ (or the portion
of $\mathcal{B(A)}$ which lies in $C$) in terms of the probability that a
perturbed initial condition still belongs to $\mathcal{B(A)}$.

According to Grümm [22, Page 4], Holling and collaborators are the first to
use this approach. Grümm [22, Page 7], proposes to measure an attractor’s
resilience as the volume of its basin of attraction. Wiley et al. [78], and
Fassoni and Yang [19] suggest to measure the volume in a relative sense to
some bounded subset, as follows.

######  Definition 4.9 (Latitude in volume).

Let $C\subset\mathbb{R}^{N}$ be Lebesgue-measurable and such that
$0<\mu(C)<\infty$. We define the latitude in volume of the attractor
$\mathcal{A}$ with respect to the _region of interest_ $C$ as:

$L_{v}(\mathcal{A},C)=\frac{\mu(\mathcal{B(A)}\cap C)}{\mu(C)}.$

$L_{v}$ returns a nondimensional value in $[0,1]$. Particularly,
$L_{v}(\mathcal{A},C)=0$ if and only if all the trajectories starting in $C$
(except at most those starting in a negligible subset of $C$) do not converge
to $\mathcal{A}$ as $t\to\infty$. On the other hand, $L_{v}(\mathcal{A},C)=1$
if and only if all the trajectories starting in $C$ (except at most those
starting in a negligible subset of $C$) converge to $\mathcal{A}$ as
$t\to\infty$.

######  Remark 4.10 (Region of interest).

The “region of interest” $C$ implicitly reveals the characteristics of the
models to whom the indicator of latitude in volume should be applied. The
system at hand should be at most subjected to isolated impulsive perturbations
of bounded magnitude $\delta>0$, where the term isolated means that the
interval of time between two occurrences should be longer than the transient
required by the system to reach an “indistinguishable proximity” to the
attractor $\mathcal{A}$ from any point in the basin of attraction. Then, the
subset $C$ is defined as the set of possible perturbed initial conditions from
the attractor [47], i.e. $B_{\delta}(\mathcal{A})$. If additional information
on the probability distribution of such perturbations is available, the
indicator can be refined (see Definition 4.11). On the other hand, this
indicator assumes additional meaning in the case of multi-stable systems [78,
19]. In this case the basin’s volume can be interpreted as a measure of the
“attractor’s relevance” in the phase space, where it competes against other
attractors and their basins over the phase space. Thus, the region of interest
should include all the (relevant) attractors $\mathcal{A}_{i}$,
$i=1,\dots,\nu$ (assuming that they are included in a set of bounded measure)
as well as the possible perturbed initial conditions of magnitude $\delta$
from them, i.e. $\bigcup_{i=1}^{\nu}B_{\delta}(\mathcal{A}_{i})\subset C$.
Even so, there is no consensus in the literature as how to choose the region
$C$ itself. For example, in [19] the smallest $N$-dimensional interval $I$
containing $\bigcup_{i=1}^{\nu}B_{\delta}(\mathcal{A}_{i})$ is considered (see
the conceptual sketch in Figure 7(b)). This approach privileges the immediate
calculation of the measure of $I$.

(a)

(b)

Figure 7. Depiction of two different approaches in choosing the region of
interest $C$ depending on the available information. The relevant attractor is
sketched as a black dot and identified by the symbol $\mathcal{A}$, whilst its
basin of attraction corresponds to the gray area denoted by the symbols
$\mathcal{B(A)}$. The region of interest $C$ is shown as the area encircled by
a black dashed line and denoted by $C$. (A) The latitude in volume of the
attractor $\mathcal{A}$ is calculated with respect to the region of interest
$C$ corresponding to the set of expected perturbations of the initial
conditions (red region). (B) The latitude in volume of the attractor
$\mathcal{A}$ is calculated with respect to the $n-$dimensional interval $C$
containing also the further attractors $\mathcal{A}_{1},\mathcal{A}_{2}$ as
well as the set of each local attractor’s expected perturbed initial
conditions (red regions).

Menck et al. [47, 46] propose a variation of the latitude in volume (see
Definition 4.9) where a weight over the phase space is considered, according
to an a priori chosen probability density $\rho$ of the expected perturbations
of initial conditions from the attractor. This indicator and its numerical
calculation has been investigated in several works including for example
Lundström [42], Lundström and Adainpää [43], Evans and Swartz [18], and Mitra
et al. [53].

######  Definition 4.11 (Basin stability).

Let $\rho:\mathbb{R}^{N}\to\mathbb{R}^{+}$ be a Lebesgue integrable function
such that $\int_{\mathbb{R}^{N}}\rho(x)\,dx=1$. The basin stability
$S_{\mathcal{B(A)}}$ of an attractor $\mathcal{A}$ with respect to $\rho$ is
defined as:

$S_{\mathcal{B(A)}}(\rho)=\int_{\mathbb{R}^{N}}\chi_{\mathcal{B(A)}}(x)\rho(x)\,dx,$
(4.5)

where $\chi_{Z}(x)$ is the characteristic function of
$Z\subset\mathbb{R}^{N}$. $S_{\mathcal{B(A)}}$ returns a nondimensional value
in $[0,1]$.

######  Remark 4.12.

Given $C\subset\mathbb{R}^{N}$ Lebesgue-measurable and such that
$0<\mu(C)<\infty$ and considered

$\rho_{C}(x)=\begin{cases}1/\mu(C)&\text{for }x\in C,\\\ 0&\text{else,}\\\
\end{cases}$

one immediately has that $S_{\mathcal{B(A)}}(\rho_{C})=L_{v}(\mathcal{A},C)$.

###  4.3. Discussion

  * •

Basin shape indicators, while being a useful intuitive tool, should be treated
very carefully. Indeed, often a real-world system is inherently “open”,
meaning that it is constantly subjected to a time-dependent forcing
(deterministic and/or random). Hence, the boundary of the basin of attraction
of the unperturbed problem may not have a dynamical meaning and the indicator
can become misleading. A concrete example of this fact can be recognized in
the phenomenon of tipping induced by non-autonomous terms as treated in
Section 6. Furthermore, even if this does not apply and the system is subject
to only instantaneous isolated perturbations, the estimation of a basin’s
boundary is usually a difficult task. For example Kerswell et al. [33] use a
variational method employing the maximal amplification (see Definition 3.6) to
estimate the distance to threshold. On the other hand, basins of attraction
with fractal-like boundaries (or intermingled and riddled basins, which do not
arise for trapped attractors but are possible for the Milnor’s definition) are
discussed in the work of Schultz et al. [68].

  * •

The following result establishes a direct relation between latitude in width,
latitude in volume and distance to threshold. It is intuitively clear that the
latter delivers the safest information (see Figure 8). Nevertheless, we
include a brief proof as it seems to not be anywhere in the literature.

Figure 8. Depiction of the phase space of the system induced by
$\dot{r}=r(r-cos(7\varphi)-(1+\varepsilon))$, $\dot{\varphi}=0$, where
$(r,\varphi)\in[0,\infty)\times[0,2\pi)$ and $\varepsilon>0$. The problem has
an asymptotically stable equilibrium at the origin. The basin of attraction is
in a shape resembling a “flower”, and has a relatively big volume. However,
the distance to threshold is only $DT=\varepsilon.$

######  Proposition 4.13.

Consider a continuous dynamical system $\phi$ on $\mathbb{R}^{N}$. Then, for a
given attractor $\mathcal{A}$ and for any $C\subset\mathbb{R}^{N}$ measurable
with $0<\mu(C)<\infty$ we have that

$2DT(\mathcal{A})\leq L_{w}(\mathcal{A}),\quad\text{and}\quad
L_{v}(\mathcal{A},C)\leq
L_{v}\big{(}\mathcal{A},B_{DT}(\mathcal{A})\big{)}=1.$

###### Proof.

We let $S_{w}$ and $S_{DT}$ denote the sets introduced in Definition 4.1 and
4.3, respectively. Firstly, notice that if $DT(\mathcal{A})=+\infty$, then
$S_{DT}=\emptyset$. Therefore, the basin’s boundary $\partial B(A)=\emptyset$
and thus also $S_{w}=\emptyset$ and $L_{w}=+\infty$. On the other hand, if
$DT(\mathcal{A})<\infty$ and $L_{w}(\mathcal{A})=+\infty$ the aimed inequality
is trivially true. Assume now that $DT(A)<+\infty$ and $L_{w}(A)<+\infty$. Due
to Proposition 4.2 and Remark 4.4 there is $(y,z)\in S_{w}$ such that
$L_{w}(\mathcal{A})=|y-z|$ and $(u,v)\in S_{DT}$ such that
$DT(\mathcal{A})=|u-v|$. Let $\alpha^{*}\in(0,1)$ be such that
$x=\alpha^{*}y+(1-\alpha^{*})z\in\mathcal{A}$. Then, we have that
$|u-v|\leq\min\\{|y-x|,|x-z|\\}$. Therefore,

$2DT(\mathcal{A})=2|u-v|\leq|y-x|+|x-z|=|y-z|=L_{w}(\mathcal{A}),$

which is the first aimed inequality. In regards to the second inequality, an
argument by contradiction allows to easily show that
$B_{DT}(\mathcal{A})\subset\mathcal{B(A)}$. Thus, we have

$L_{v}(\mathcal{A},B_{DT}(\mathcal{A}))=\frac{\mu(\mathcal{B(A)}\cap
B_{DT}(\mathcal{A}))}{\mu(B_{DT}(\mathcal{A}))}=\frac{\mu(B_{DT}(\mathcal{A}))}{\mu(B_{DT}(\mathcal{A}))}=1.$

On the other hand,

$1=\max\\{L_{v}(\mathcal{A},C)\mid C\subset\mathbb{R}^{N},\text{ measurable
with }0<\mu(C)<\infty\\}$

which immediately gives the second inequality. ∎

##  5\. nonlinear transient dynamics

The indicators of resilience introduced in the previous sections focus either
on the properties of the linearized flow, or on the geometrical features of
the basin of attraction. However, they do not take into account the nonlinear
transient dynamics within the basin, which may play a crucial role. In order
to better exemplifies this idea let us briefly present an example from Meyer
and McGehee [51]. Consider the differential equations

$\dot{x}=f(x)=\begin{cases}x&\text{if }x<0,\\\ \pi^{-1}\sin(\pi x)&\text{if
}0\leq x\leq 1,\\\ 1-x&\text{if }x\geq
1,\end{cases}\quad\text{and}\quad\dot{x}=g(x)=-x(x-1).$

The dynamical systems induced by both equations have an unstable equilibrium
at $x=0$ and a stable equilibrium at $x=1$ . Therefore the distance to
threshold for $\mathcal{A}=1$ is the same for both systems, namely
$DT(\mathcal{A})=1$. Moreover, the characteristic return times are also the
same $T_{R}(\mathcal{A})=1$. Nonetheless, the two functions $f$ and $g$ attain
different maxima, respectively $\mu=\pi^{-1}$ and $\widehat{\mu}=1/4$. For a
given uncertainty $\widehat{\mu}<\varepsilon<\mu$ on both vector fields, it
can happen that the induced dynamical systems cease to be topologically
equivalent (see Figure 9). In particular the dynamical system induced by
$\dot{x}=g(x)-\varepsilon$ does not have any bounded solution anymore.

\begin{overpic}[width=433.62pt,trim=0.0pt 398.33858pt 0.0pt
0.0pt,clip]{images/Fig9.pdf} \put(0.0,22.0){(a)} \put(35.0,22.0){(b)}
\put(70.0,22.0){(c)} \put(3.0,18.0){{\color[rgb]{0,0,1}$\mu$}}
\put(3.0,15.0){{\color[rgb]{1,0,0}$\widehat{\mu}$}}
\put(12.0,19.0){{\color[rgb]{0,0,1}$f(x)$}}
\put(47.0,19.0){{\color[rgb]{0,0,1}$f(x)$}}
\put(12.0,13.5){{\color[rgb]{1,0,0}$g(x)$}}
\put(82.0,18.0){{\color[rgb]{1,0,0}$g(x)$}} \put(27.0,13.5){$x$}
\put(23.5,8.5){$1$} \put(59.7,8.5){$\varepsilon$}
\put(59.7,13.0){$\varepsilon$} \put(94.5,8.5){$\varepsilon$}
\put(94.5,13.0){$\varepsilon$} \end{overpic} Figure 9. (a) Two continuous
vector fields $f$ and $g$ induce ordinary differential equations sharing the
same distance to threshold and characteristic return time for the attractor
$\mathcal{A}=1$. However, a common uncertainty $\varepsilon>0$ on both problem
may lead to dynamical systems which are not topologically equivalent, since
(b) the first maintains always a nonempty set of bounded solutions, (a)
whereas it may occur that the latter has no bounded solutions. The example is
thoroughly described by Meyer and McGehee [51].

This section contains five indicators. _Return time_ in subsection 5.1,
_resistance of a gradient system_ in subsection 5.2, _resilience boundary_ in
subsection 5.3, _intensity of attraction_ in subsection 5.4 and _expected
escape times_ in subsection 5.5.

Except for the return time, which is defined for a general continuous flow on
$\mathbb{R}^{N}$, all the other indicators are defined for dynamical systems
induced by ordinary differential equations of the type $\dot{x}=f(x)$,
$x\in\mathbb{R}^{N}$, for which we assume that the assumptions in Section 2
are satisfied.

###  5.1. Return Time

Between the ’70s and the ’90s of the last century, the notion of return time
as an indicator of resilience (firstly introduced by O’Neill [56]) gathered
attention and stimulated a fair amount of research (e.g. Cottingham and
Carpenter [12], DeAngelis et al. [14, 15], Harte [26]). The reason lies in the
fact that this definition was able to capture the transient behaviour of the
system after an initial perturbation in contrast with the more classical
characteristic return time (see Definition 3.1). Although there are several
variations in the literature (mostly differing only for the constant of
normalization), they all are similar to the definition below [55, 56].

######  Definition 5.1 (Return time).

Consider a continuous dynamical system $\phi$ on $\mathbb{R}^{N}$, and assume
that $\mathcal{A}$ is an attractor for $\phi$ and $\mathcal{B(A)}$ its basin
of attraction. The return time for the system from
$x_{p}\in\mathcal{B(A)}\setminus\mathcal{A}$ to $\mathcal{A}$ is defined as:

$T_{R}(\mathcal{A},x_{p})=\frac{1}{\mathrm{dist}(x_{p},\mathcal{A})}\int_{0}^{\infty}\mathrm{dist}(\phi(t,x_{p}),\mathcal{A})\,dt,$
(5.1)

where $\mathrm{dist}(A,B)$ denotes the Hausdorff semi-distance between the
sets $A,B\subset\mathbb{R}^{N}$. Moreover, for a fixed measurable set
$C\subseteq\mathcal{B(A)}\setminus\mathcal{A}$, the maximal return time and
average return time from $C$ to $\mathcal{A}$ are respectively defined as

$T_{R}^{sup}(\mathcal{A},C)=\sup_{x_{p}\in
C}T_{R}(\mathcal{A},x_{p}),\quad\text{and}\quad
T_{R}^{mean}(\mathcal{A},C)=\frac{1}{\mu(C)}\int_{C}T_{R}(\mathcal{A},x)\
d\mu(x).$

The definition of return time is a valuable theoretical tool, but the
indefinite integral in (5.1) is a clear obstacle in practical applications.
Therefore, more calculable indicators, still accounting for the transient
dynamics (see Neubert and Caswell [55]), have been developed and were
frequently preferred in applications.

The following result shows (as expected) that the integral in (5.1) is finite
when the attractor consists of a hyperbolic solution (see Definition 2.6).

######  Proposition 5.2.

Consider a continuous dynamical system on $\mathbb{R}^{N}$ induced by an
ordinary differential equation $\dot{x}=f(x)$, with $f$ continuously
differentiable. Furthermore, assume that the attractor $\mathcal{A}$ is made
of the points of a locally attractive hyperbolic solution $\widetilde{x}$ (see
Definition 2.6). Then,

$T_{R}(\mathcal{A},x_{p})<+\infty,\qquad\text{for all
}x_{p}\in\mathcal{B(A)}.$

###### Proof.

Due to the First Approximation Theorem (see [24, Theorem III.2.4]), there is
$\delta>0$ such that if $\overline{x}\in\mathcal{B(A)}$ satisfies
$|\overline{x}-\widetilde{x}(s)|\leq\delta$ for some $s\in\mathbb{R}$, then
for all $t\geq s$, $|\phi(t,\overline{x})-\widetilde{x}(t)|<Ke^{-\beta(t-s)}$
for some $K>1$ and $\beta>0$. On the other hand, by definition, for every
$x_{p}\in\mathcal{B(A)}\setminus\mathcal{A}$ there is $s>0$ such that
$\mathrm{dist}(\phi(s,x_{p}),\mathcal{A})<\delta$. Then the result is a
consequence of the linearity of the integral in (5.1). ∎

###  5.2. Resistance of a gradient system

Stability landscapes are a commonly used tool in mathematical ecology as they
provide immediate intuitive information on the dynamical features of
equilibria (when they exist) of low-dimensional systems. Limitations and
warnings of cautious employment of this instrument have been pointed out by
many (see [49] for example). In mathematical terms a stability landscape can
be seen as the potential $V$ of a gradient system (and related to the notion
of Lyapunov function [10]). This means that the phenomenon can be represented
via an ordinary differential equation of the type $\dot{x}=-\nabla V(x)$,
where $V$ is a twice continuously differentiable real valued function. Then,
one can consider the minimal work required to bring the system from an
attractor $\mathcal{A}$ to a perturbed initial condition outside the basin of
attraction $\mathcal{B(A)}$ as an indicator of resilience (see for example
Lundström [42]).

######  Definition 5.3 (Resistance of a gradient system).

Consider a continuous dynamical system induced by an ordinary differential
equation $\dot{x}=-\nabla V(x)$, where $V:\mathbb{R}^{N}\to\mathbb{R}$ is a
twice continuously differentiable. The _resistance (potential) of a gradient
system_ with respect to the attractor $\mathcal{A}$ is

$W(\mathcal{A})=\inf_{x_{0}\notin\mathcal{B(A)},\,a\in\mathcal{A}}\big{(}V(x_{0})-V(a)\big{)}.$

######  Remark 5.4.

Walker et al. [75] call _resistance of an attractor_ $\mathcal{A}$ the
quantity defined by

$R_{W}(\mathcal{A})=\frac{W(\mathcal{A})}{L_{w}(\mathcal{A})},$

where $L_{w}(\mathcal{A})$ is the latitude in width of $\mathcal{A}$. As
pointed out also by Meyer [49], care should be taken in the interpretation of
this indicator. In particular, if the timescales of disturbance and internal
dynamics are considerably different, the slope of the potential function at
the state of the system does not necessarily relate to the capability of the
system to recover.

###  5.3. Resilience Boundary

Meyer et al. [50] suggest to measure the resilience of an attractor by
introducing and applying the so-called “flow-kick framework”, a sequence of
repeated “impulsive” disturbances. Given a continuous dynamical system induced
by an ordinary differential equation $\dot{x}=f(x)$, $x\in\mathbb{R}^{N}$,
assume that solutions are defined for all $t\in\mathbb{R}^{+}$. For fixed
$\tau>0$, $\kappa\in\mathbb{R}^{N}$, and $a\in\mathcal{A}$ a _flow-kick
trajectory_ $\widetilde{x}_{\tau,\kappa}(\cdot,a)$ starting in
$a\in\mathbb{R}^{N}$ is constructed by considering a partition of the positive
real line in intervals $[j\tau,(j+1)\tau)$, $j\in\mathbb{N}$, and a sequence
of functions $\widetilde{x}_{j}:[j\tau,(j+1)\tau]\to\mathbb{R}^{N}$ such that
$\widetilde{x}_{j}$ solves the Cauchy problem

$\dot{x}=f(x),\qquad\widetilde{x}_{j}(j\tau)=\widetilde{x}_{j-1}(j\tau)+\kappa,\quad\text{
and }\quad\widetilde{x}_{0}(0)=a\in\mathcal{A}.$

One may interpret $\kappa$ as a _kick_ , while $\tau$ represents the _flow
time_ , and a pair $(\tau,\kappa)\in(0,\infty)\times\mathbb{R}^{N}$ is called
a _disturbance pattern_. Notice that except for the trivial case where
$\kappa$ is the origin in $\mathbb{R}^{N}$, a flow-kick trajectory is not a
solution of the original problem.

######  Definition 5.5 (Flow-kick resilience set and resilience boundary).

Under the previously stated assumptions and notation, a disturbance pattern
$(\tau,\kappa)\in(0,\infty)\times\mathbb{R}^{N}$ is said to be in the _flow-
kick resilience set_ if

$\inf_{t>0}P_{\mathcal{A}}\big{(}\widetilde{x}_{\tau,\kappa}(t,a)\big{)}>0\quad\text{for
all }a\in\mathcal{A},$

where $P_{\mathcal{A}}(y)$ is the precariousness of $y\in\mathbb{R}^{N}$ (see
Definition 4.7). Moreover, the boundary of the flow-kick resilience set in
$(0,\infty)\times\mathbb{R}^{N}$ is called the _(flow-kick) resilience
boundary_ for the basin of attraction $\mathcal{B(A)}$.

If $N>1$ the resilience boundary should be considered as a precautionary
threshold. Indeed, it allows to identify the flow-kick trajectories which
leave $\mathcal{B(A)}$ at some point, but does not allow to distinguish them
from those who eventually return to $\mathcal{B(A)}$ [50]. On the other hand,
if $N=1$ and $\widetilde{x}_{\tau,\kappa}(\overline{t},a)\notin\mathcal{B(A)}$
for some $a\in\mathcal{A}$ and disturbance pattern $(\tau,\kappa)$, the same
is true for all $t\geq\overline{t}$. In fact, under appropriate assumptions,
the scalar case allows to associate a quantitative indicator based on the area
of the region in between the flow-kick resilience set and the line
$|\kappa|=DT(\mathcal{A})$; see also Figure 10, where $DT(\mathcal{A})$
denotes the distance to threshold for $\mathcal{A}$. This indicator has the
merit of establishing a relation with distance to threshold (Definition 4.3)
and distance to bifurcation (Definition 6.1).

Figure 10. Example of flow-kick resilience and resilience boundary for scalar
models in population dynamics. The upper left-hand side panel portrays two
populations. Population 1 modelled by
$\dot{x}=x\left(1-\frac{x}{100kt}\right)\left(\frac{x}{20kt}-1\right),$ (solid
blue line) where $kt$ stands for kilo-tonnes, and population 2 by
$\dot{x}=x\left(1-\frac{x}{100kt}\right)\left(\frac{x}{20kt}-1\right)(0.0002x^{2}-0.024x+1.4)$
(dashed red line). Both systems have an attracting fixed point at
$\mathcal{A}=\\{100kt\\}$ whose basin of attraction is marked in grey. The two
systems have same distance to threshold and characteristic return time. Their
resilience boundaries are shown in the upper right-hand side panel where also
three kick flow patterns are highlighted. In the lower panels, two pairs of
flow-kick trajectories are shown for population 1 on the left-hand side and
population 2 on the right-hand side. In black the flow-kick trajectory
relative to the disturbance pattern $(6\text{ months},-24kt)$ , in red the one
relative to $(12\text{ months},-48kt)$ and in green the one relative to
$(18\text{ months},-30kt)$. The disturbance pattern $(6\text{ months},-24kt)$
lies in the resilience set for population 1, but outside the resilience set of
population 2. We refer the reader to [50] for further details.

###  5.4. Intensity of Attraction

The intensity of attraction aims to measure the strength of the attraction of
an attractor. It was first developed by McGehee [45] in the context of
continuous maps acting on a compact metric space, and it has been recently
extended by Meyer [48] and Meyer and McGehee [51] to the continuous-time
dynamics induced from an ordinary differential equation on $\mathbb{R}^{N}$.
For the sake of consistency with the set-up of this work we shall present the
latter.

Consider a continuous dynamical system induced by an ordinary differential
equation $\dot{x}=f(x)$, with $f:U\subset\mathbb{R}^{N}\to\mathbb{R}^{N}$, $U$
open, and $f$ bounded and globally Lipschitz continuous on $U$. Moreover,
consider an interval $I\subset\mathbb{R}$ containing zero and an essentially
bounded function $g\in L^{\infty}(I,\mathbb{R}^{N})$, that is,
$g:I\to\mathbb{R}^{n}$ such that $\inf\\{c\geq
0:\mu(\\{t\in\mathbb{R}:|g(t)|>c\\})=0\\}<\infty$ and the perturbed system

$\dot{x}=f(x)+g(t).$ (5.2)

Equation (5.2) is well-defined (it is in fact a special type of so-called
Carathéodory differential equations) and it satisfies the conditions of
existence and uniqueness for the associated Cauchy problems [9]. In fact, it
is also assumed that for any $x_{0}\in\mathbb{R}^{N}$ the solution
$x_{g}(\cdot,x_{0})$ of (5.2) with $x(0)=x_{0}$ can be continuously extended
to $I$. Notice also that the function

$\phi:I\times\mathbb{R}^{N}\times
L^{\infty}\to\mathbb{R}^{N},\quad(t,x_{0},g)\mapsto\phi(t,x_{0},g)=x_{g}(t,x_{0})$

vary continuously in $I\times\mathbb{R}^{N}\times L^{\infty}$ [51]. As a
matter of fact, under the given assumptions, the problem also induces a
continuous skew-product flow [39, 40]. The fundamental idea behind the concept
of intensity is to conceive each bounded set of perturbing functions $g\in
L^{\infty}$ as a class of controls, and to investigate which parts of the
phase space can be reached through them. In this way one is able to relate the
transient dynamics of a system to the minimal class of controls, which are
able to drive the points of an attractor away from it for all future time. In
order to state the mathematical definition, let us recall some notation. Given
$S\subset\mathbb{R}^{n}$, by $\phi(t,S,g)$ we shall denote the set
$\bigcup_{x_{0}\in S}\phi(t,x_{0},g)$ and we shall call the reachable set from
$S$ via controls of norm $r>0$, the set

$P_{r}(S)=\bigcup_{t\geq 0}\bigcup_{g\in F_{r}}\phi(t,S,g),$ (5.3)

where $F_{r}(0)=\\{g\in L^{\infty}(I,\mathbb{R}^{n})\colon\left\lVert
g\right\rVert_{\infty}<r\\}$.

######  Definition 5.6 (Intensity of attraction).

The intensity of attraction $\mathcal{I}(\mathcal{A})$ of an attractor
$\mathcal{A}$ is defined as:

$\mathcal{I}(\mathcal{A})=\sup\\{r\geq 0\mid P_{r}(\mathcal{A})\subset
K\subset\mathcal{B(A)}\text{ for some compact }K\subset\mathbb{R}^{N}\\}.$
(5.4)

Interestingly, the intensity of attraction provides information on the
persistence of an attractor not only against a time-dependent perturbation.
Meyer et al. prove that for any vector field, which differs in norm from $f$
less than its intensity of attraction, a compact positively invariant subset
of $\mathbb{R}^{N}$ exists which contains the attractors of both systems (see
[51, Theorem 6.12]).

###  5.5. Expected escape times

It is a well-known fact that the presence of a stochastic perturbation can
have a non-negligible impact on the dynamics of a system and possibly trigger
a critical transition [36]. Dennis et al. [16] who, in turn, refer to Gardiner
[20] for the main idea, suggest to judge the resilience of an attractor for a
system subjected to white noise via the notion of mean first passage time,
also known as expected escape time.

######  Definition 5.7 (Expected escape times).

Consider a scalar stochastic differential equation

$dX=f(X)\,dt+\sqrt{\nu(X)}\,dW(t),$ (5.5)

where $\nu:\mathbb{R}\to\mathbb{R}^{+}$ is continuous, and $W$ is a standard
Brownian motion. If (5.5) has a stationary distribution with probability
density $p(x)$, the _expected escape times_ are defined as

$\begin{split}\tau_{1}(x)&=2\int_{x_{0}}^{x}\frac{\int_{0}^{y}p(z)\,dz}{\nu(x)p(y)}\,dy\quad\text{for
}x_{0}<x,\\\
\tau_{2}(x)&=2\int^{x_{0}}_{x}\frac{\int_{y}^{\infty}p(z)\,dz}{\nu(x)p(y)}\,dy\quad\text{for
}x_{0}>x.\end{split}$

The expected escape times represent the mean amount of time necessary for the
process $X$ to attain the value $x$ starting from $x_{0}<x$ and $x_{0}>x$,
respectively. As explained in [16], one or both these integrals can be
infinite. For example if $x=0$ acts like an absorbing value and the system
lies close to it, there is a positive probability that the first escape time
takes an infinite value and so $p(x)$ is not integrable near $0$. An analogous
reasoning holds for the second integral if $\infty$ acts as an absorbing
boundary. In this case the right tail of $p$ is not integrable.

###  5.6. Discussion

In 1979 Harrison [25] introduced the concept of _stress period_ for a
differential equation, that is, a compact interval of time in which the
problem is perturbed and one can measure the ability of the system to resist
displacement, and the rate of recovery to the original condition once the
perturbation is over. In the original formulation (presented in subsection
6.2) the perturbation consists only of a time dependent piece-wise-constant
change of a parameter. However, it can be easily generalized as exposed in
Remark 6.6 and included among the indicators of resilience focusing on
transient dynamics.

##  6\. Variation of parameters

In addition to the previous notions of resilience, one can also view
parameters as static phase space variables, which are changed infinitely
slowly. Due to their special role, one can then ask, how one may define
resilience in parameter space. Therefore, the term resilience acquires a wider
significance as it does not limit to the capability of a system to endure a
perturbation and revert to a desirable state but to the overall persistence of
the dynamical relations which determine the qualitative behaviour of all the
trajectories. Let us set some common notation for this section. We shall
consider a family of continuous dynamical system induced by a parametric
family of ordinary differential equations

$\dot{x}=f(x,\lambda),\quad
x\in\mathbb{R}^{N},\,\lambda\in\Lambda\subset\mathbb{R}^{M},$ (6.1)

with $f\in C(\mathbb{R}^{N}\times\mathbb{R}^{M},\mathbb{R}^{N})$ satisfying
the assumptions in Section 2 for all $\lambda\in\Lambda$. When we fix a
parameter $\lambda_{0}\in\mathbb{R}^{M}$, we will refer to the respective
equation by the symbol (6.1)${}_{\lambda_{0}}$ and denote by
$\phi_{\lambda_{0}}$ the local flow on $\mathbb{R}^{N}$ induced by
(6.1)${}_{\lambda_{0}}$.

This section contains three subsections. In subsection 6.1, we present the
indicator _distance to bifurcation_. In subsection 6.2, we present three
indicators: _resistance_ , _elasticity_ and _persistence_. Lastly, the
relation between resilience and tipping points in non-autonomous systems is
explored in subsection 6.3.

###  6.1. Distance to Bifurcation

Bifurcation theory is a classical branch of nonlinear dynamical systems which
studies the lack of topological equivalence in phase space upon (smooth)
change of one or more of the system’s parameters. The classic theory for
autonomous dynamical systems encompasses the bifurcations of local steady
states, periodic orbits, as well as homoclinic and heteroclinic orbits [38].
For its characteristics, a bifurcation point entails a profound change in the
dynamics of the system, which may include the change of stability or the
disappearance of a considered attractor. In practical terms, a bifurcation
point can sometimes be undesirable, especially when the system exhibits
hysteresis or irreversibility [7]. It is therefore natural that the distance
to a bifurcation point has been suggested as a measure of resilience [17, 30].

######  Definition 6.1 (Distance to bifurcation).

Let $\mathcal{A}$ be an attractor for the dynamical system induced by
(6.1)${}_{\lambda_{0}}$, and denote by
$\Lambda_{bif}(\mathcal{A})\subset\mathbb{R}^{M}$ the set of bifurcation
points for $\mathcal{A}$. The distance to bifurcation for $\mathcal{A}$ is
defined as

$D_{\rm{bif}}(\mathcal{A})=\inf_{\lambda^{*}\in\Lambda_{bif}}|\lambda_{0}-\lambda^{*}|.$

###  6.2. Resistance, elasticity and persistence

Harrison [25] examines the transient behaviour of a system which undergoes an
instantaneous change of parameters and returns to the original values, just as
instantly, after an interval of time. The system’s response during and after
the perturbation are suggested as indicators of the system’s capacity to
withstand a change and to recover. Harrison names these features respectively
resistance and resilience (although, to avoid confusion, we will use Orians’s
nomenclature [57] and call Harrison’s resilience elasticity). Using Webster’s
words [76], the first describes “the ability to resist displacement”, whereas
the latter “the rate of recovery to the original condition”. With the notation
set at the beginning of the section we have the following definitions:

######  Definition 6.2 (Resistance).

Given an attractor $\mathcal{A}$ for (6.1)${}_{\lambda_{0}}$, and fixed $T>0$
such that for all $\lambda\in\Lambda$, and $a\in\mathcal{A}$,
$\phi_{\lambda}(\cdot,a)$ is defined at least in $[0,T]$, the _resistance_ of
system (6.1)${}_{\lambda_{0}}$ at $\mathcal{A}$ with respect to $\Lambda$ and
the _stress period_ $[0,T]$ is defined as:

$R(\mathcal{A},\Lambda,T)=\sup_{\begin{subarray}{c}t\in[0,T],\,a\in\mathcal{A}_{\lambda_{0}}\\\
\lambda\in\Lambda\end{subarray}}|\phi_{\lambda}(t,a)-\phi_{\lambda_{0}}(t,a)|.$

######  Definition 6.3 (Elasticity).

Given an attractor $\mathcal{A}$ for (6.1)${}_{\lambda_{0}}$, and fixed $T>0$
such that $\phi_{\lambda}(T,\mathcal{A})\subset\mathcal{B}(\mathcal{A})$ for
all $\lambda\in\Lambda$, consider the continuous function
$\Phi_{\lambda_{0}}:[T,\infty)\times\mathcal{A}\times\Lambda\to\mathbb{R}^{N}$
defined by

$\Phi_{\lambda_{0}}(t,a,\lambda):=|\phi_{\lambda_{0}}\big{(}t-T,\phi_{\lambda}(T,a)\big{)}-\phi_{\lambda_{0}}(t,a)|.$

The elasticity of the system (6.1)${}_{\lambda_{0}}$ at $\mathcal{A}$ with
respect to $\Lambda$ and the _stress period_ $[0,T]$ is defined as:

$E(\mathcal{A},\Lambda,T)=\sup_{\begin{subarray}{c}t\in[T,\infty),\,a\in\mathcal{A}_{\lambda_{0}}\\\
\lambda\in\Lambda\end{subarray}}\left(\frac{1}{\Phi_{\lambda_{0}}(t,a,\lambda)}\frac{d\,\Phi_{\lambda_{0}}(t,a,\lambda)}{dt}\right).$

Figure 11. Effect of a stress period of Harrison’s type on two populations
following the logistic model $\dot{x}=rx\left(1-x/K\right)$ with carrying
capacity $K=1$, and growth rates respectively equal to $r_{1}=1$ and to
$r_{2}=2$. A stress period $[0,1]$ is applied to both models during which the
carrying capacity is diminished of the $20\%$ (i.e. $\Lambda=\\{0.8\\})$ . The
resistance of the attractor $x^{*}=1$ is lower for the first species than the
second. Conversely, the elasticity of the attractor $x^{*}=1$ is higher for
the first species than the second.

######  Remark 6.4.

The definitions of resistance and elasticity presented above admit a natural
generalization. Instead of using the original trajectory inside the attractor
as a reference, in practice it may be sufficient to record the displacement
from the attractor or the rate of convergence to the attractor. This would
lead to the following weaker definitions

$R^{w}(\mathcal{A},\Lambda,T)=\sup_{\begin{subarray}{c}t\in[0,T],\,a,a^{\prime}\in\mathcal{A}_{\lambda_{0}}\\\
\lambda\in\Lambda\end{subarray}}|\phi_{\lambda}(t,a)-a^{\prime}|,$

and

$\Phi^{w}_{\lambda_{0}}(t,a,\lambda):=\inf_{a^{\prime}\in\mathcal{A}}|\phi_{\lambda_{0}}\big{(}t-T,\phi_{\lambda}(T,a)\big{)}-a^{\prime}|,$
$E^{w}(\mathcal{A},\Lambda,T)=\sup_{\begin{subarray}{c}t\in[T,\infty),\,a\in\mathcal{A}_{\lambda_{0}}\\\
\lambda\in\Lambda\end{subarray}}\left(\frac{1}{\Phi^{w}_{\lambda_{0}}(t,a,\lambda)}\frac{d\,\Phi^{w}_{\lambda_{0}}(t,a,\lambda)}{dt}\right).$

Ratajczak et al. [63] use a setup similar to Harrison’s to address the problem
of the interaction between the duration and the “intensity” of a stress period
in determining a tipping of the system from the original basin of attraction
to a different one (if it applies). In accordance with Harrison’s
nomenclature, we shall call this property persistence.

######  Definition 6.5 (Persistence).

The _persistence_ of an attractor $\mathcal{A}$ for (6.1)${}_{\lambda_{0}}$
with respect to a _stress intensity_ $\lambda\in\Lambda$ is given by

$P_{\lambda}(\mathcal{A})=\sup\\{T>0\mid\phi_{\lambda}(t,\mathcal{A})\subset\mathcal{B}(\mathcal{A})\text{
for all }0\leq t\leq T\\}.$

The _persistence_ of an attractor $\mathcal{A}$ for (6.1)${}_{\lambda_{0}}$
with respect to a _stress duration_ $T>0$ is given by

$P_{T}(\mathcal{A})=\sup\\{\rho>0\mid\phi_{\lambda}(t,\mathcal{A})\subset\mathcal{B}(\mathcal{A})\text{
for all }0\leq t\leq T\text{ and }|\lambda-\lambda_{0}|\leq\rho\\}.$

######  Remark 6.6.

It is clear that Harrison’s formulations of resistance and elasticity as well
as of persistence generalize to any perturbation on a finite time interval
rather than just a piece-wise constant change of parameters. Consider
(6.1)${}_{\lambda_{0}}$ for some $\lambda_{0}\in\mathbb{R}^{M}$, $T>0$, and a
set $\Lambda$ of functions from $\mathbb{R}\times\mathbb{R}^{N}$ into
$\mathbb{R}^{M}$ such that if $\lambda\in\Lambda$, then for every compact set
$K\subset\mathbb{R}^{N}$, $\sup_{x\in K}\lambda(t,x)$ is locally integrable
and $\lambda(t,x)=\lambda_{0}$ for all $x\in\mathbb{R}^{N}$ and
$t\notin[0,T]$. Hence, one can consider the problem

$\dot{x}=f\big{(}x,\lambda(t,x)\big{)},$

and still define resistance, elasticity and persistence as above.

###  6.3. Resilience against rate-induced tipping

Sizable changes in the output of a system upon a negligible variation of the
input are referred to as _critical transitions_ or _tipping points_. Motivated
by current challenges in nature and society [21, 66], the study of the several
mechanisms leading to a critical transition has experienced significant
mathematical interest as well [36, 37]. In recent years for example, the
phenomenon of rate-induced tipping [77] has been proposed as an alternative
mechanism for a critical transition, with respect to the more classical
autonomous bifurcations and noise-induced tipping [3]. Rate-induced tipping
can be seen as a special type of nonautonomous bifurcation for the
differential equation

$\dot{x}=f\big{(}x,\gamma(rt)\big{)},\qquad
x\in\mathbb{R}^{N},t\in\mathbb{R},r>0$ (6.2)

obtained from (6.1) through a time-dependent variation of parameters
$\gamma:\mathbb{R}\to\mathbb{R}^{M}$ which is assumed to be bounded,
continuous and asymptotically constant. In particular, the parameter $r>0$
represents the rate at which the time-dependent variation takes place.
Furthermore, assume that for all $\lambda$ in the image of $\gamma$, the
dynamical systems induced by (6.1)λ, are all topologically equivalent. One
might expect that for any fixed attractor $\mathcal{A}_{0}$ of the _past
limit-problem_ , and denoted by $\mathcal{A}_{f}$ the “corresponding”
attractor of the _future limit-problem_ , for all $r>0$ there is an attractor
of (6.2)r which approaches $\mathcal{A}_{0}$ as $t\to-\infty$ and
$\mathcal{A}_{f}$ as $t\to\infty$. However, this is not always the case and,
depending on $r$, a local attractor of (6.2)r limiting at $\mathcal{A}_{0}$ as
$t\to-\infty$ may fail to track $\mathcal{A}_{f}$ as $t\to\infty$.

In order to be more precise and provide a rigorous presentation, we will need
to require further additional assumptions with respect to those stated at the
beginning of the section. In particular, we assume that $f\in
C^{2}(\mathbb{R}^{N}\times\mathbb{R}^{M},\mathbb{R}^{N})$ and that
$\Lambda\subset\mathbb{R}^{M}$ is compact and connected. Moreover, we will
assume that for all $\lambda\in\Lambda$ the dynamical systems induced by (6.1)
are topologically equivalent and their dynamics is completely determined by
hyperbolic fixed points (for attractors which are not fixed points see for
example [1, 32, 41]). Moreover, fix $\lambda_{0},\lambda_{\infty}\in\Lambda$
and a twice continuously differentiable function $\gamma:\mathbb{R}\to\Lambda$
such that

$\lim_{t\to-\infty}\gamma(t)=\lambda_{0},\quad\lim_{t\to\infty}\gamma(t)=\lambda_{\infty}\quad\text{and}\quad\lim_{t\to\pm\infty}\gamma^{\prime}(t)=0.$

Let $x_{0}^{s}\in\mathbb{R}^{N}$ be any fixed stable hyperbolic equilibrium
for (6.1)${}_{\lambda_{0}}$ and $x_{\lambda}^{s}\in\mathbb{R}^{N}$ be the
“corresponding” stable hyperbolic equilibrium for (6.1)λ. This means that if
$\\{h_{\lambda}:\mathbb{R}^{N}\to\mathbb{R}^{N}\\}$ is a $\lambda$-continuous
family of time-preserving homeomorphisms guaranteeing the topological
equivalence, then $x_{\lambda}^{s}=h_{\lambda}(x_{0}^{s})$. Thanks to Theorem
2.2 in [3], for every $r>0$, there is a unique maximal solution of (6.2)r
$\widetilde{x}_{\gamma,r}:(-\infty,\beta_{\gamma,r})\to\mathbb{R}^{N}$, with
$\beta_{\gamma,r}\in\mathbb{R}\cup\\{\infty\\}$, such that
$\lim_{t\to-\infty}\widetilde{x}_{\gamma,r}(t)=x_{0}^{s}$. Additionally, this
solution is _locally pullback attractive_ (a weak form of attractivity which
holds only for the past and is specific of nonautonomous systems; see [35] for
more information). Notice also that, thanks to Lemma 2.3 in [3], there is
$r_{0}>0$ such that for all $0<r<r_{0}$, $\beta_{\gamma,r}=\infty$ and
$\lim_{t\to\infty}\widetilde{x}_{\gamma,r}(t)=x_{\lambda_{\infty}}^{s}$. We
can now give the definition of rate-induced tipping.

######  Definition 6.7.

Under the introduced notation and assumptions, we say that the locally
pullback attracting solution $\widetilde{x}_{\gamma,r}$ _undergoes a rate-
induced tipping at $r_{\gamma}^{*}(x_{0}^{s})>0$_ if $\beta_{\gamma,r}=\infty$
and $\lim_{t\to\infty}\widetilde{x}_{\gamma,r}(t)=x_{\lambda_{\infty}}^{s}$
for all $r\in(0,r^{*}_{\gamma}(x_{0}^{s}))$, and either
$\beta_{\gamma,r_{\gamma}^{\scaleto{*}{3pt}}(x_{0}^{s})}<\infty$ or
$\lim_{t\to\infty}\widetilde{x}_{\gamma,r_{\gamma}^{\scaleto{*}{3pt}}(x_{0}^{s})}(t)\neq
x_{\lambda_{\infty}}^{s}$. Moreover, we will call a rate-induced tipping
_transversal_ if there is $\delta>0$ such that for all
$r\in[r_{\gamma}^{*}(x_{0}^{s}),r_{\gamma}^{*}(x_{0}^{s})+\delta]$,
$\beta_{\gamma,r}<\infty$ or $\lim_{t\to\infty}\widetilde{x}_{\gamma,r}(t)\neq
x_{\lambda_{\infty}}^{s}$.

It seems therefore natural that the quantity

$r_{\gamma}^{*}(x_{0}^{s})=\sup\left\\{\rho>0\ \bigg{|}\
\beta_{\gamma,\rho}=\infty\text{ and
}\lim_{t\to\infty}\widetilde{x}_{\gamma,\rho}(t)=x_{\lambda_{\infty}}^{s}\text{
for all }\rho\leq r\right\\}\in\mathbb{R}\cup\\{\infty\\}$

determines a threshold of resilience for the equilibrium $x_{0}^{s}$ of
(6.1)${}_{\lambda_{0}}$ against the time-dependent change of parameters
$\gamma$, in the sense that for all $r<r_{\gamma}^{*}(x_{0}^{s})$, the
solution $\widetilde{x}_{\gamma,r}(t)$ of (6.2) is able to _end-point track
the curve of quasi-static equilibria_ $x^{s}_{\gamma(t)}$. On the other hand
the inferior of $r_{\gamma}^{*}(\cdot)$ over the stable equilibria of
(6.1)${}_{\lambda_{0}}$ can be regarded as a threshold of resilience for the
whole system (6.1)${}_{\lambda_{0}}$ against the time-dependent change of
parameters $\gamma$.

This statement assumes further significance if we look at rate-induced tipping
as of a _nonautonomous bifurcation_. Then, the bifurcation point
$r_{\gamma}^{*}$ can be used to construct a distance to bifurcation analogous
to the one in Definition 6.1. There is, however, a substantial difference. The
tipping point $r_{\gamma}^{*}$ depends on the function $\gamma$ that realizes
the variation of parameters. Therefore, it is beneficial to investigate the
continuous dependence of $r^{*}_{\gamma}$ with respect to the variation of
$\gamma$.

We hereby propose a first result guaranteeing the lower semi-continuity (and
if it applies the continuity) of $r^{*}_{\gamma}$ with respect to $\gamma$ in
a specific set of functions: fixed $t_{0}\in\mathbb{R}$, we consider the set
$\Omega(\gamma,t_{0})$ of twice continuously differentiable functions
$\omega:\mathbb{R}\to\Lambda$ such that $\omega(t)=\gamma(t)$ for all $t\leq
t_{0}$ and $\lim_{t\to\infty}\omega(t)=\lambda\in\Lambda$.

######  Proposition 6.8.

Assume that $f$ is bounded and uniformly continuous, and that there are
$\varepsilon>0$ and $\overline{\delta}>0$ such that for every $\overline{r}>0$
a $T(\gamma,\overline{r})>0$ exists for which if
$\omega\in\Omega(\gamma,t_{0})$ and
$\|\gamma-\omega\|_{\mathcal{C}(\mathbb{R},\Lambda)}<\overline{\delta}$, and
if $\widetilde{x}_{\omega,\overline{r}}$ is a solution of
$\dot{x}=f(x,\omega(\overline{r}t))$ defined on
$[T(\gamma,\overline{r}),\infty)$ and
$y^{\prime}={\textnormal{D}}f\big{(}\widetilde{x}_{\omega,\overline{r}}(t),\omega(\overline{r}t)\big{)}$
has an exponential dichotomy with projector the identity on
$[T(\gamma,\overline{r}),\infty)$, then the distance between
$\widetilde{x}_{\omega,\overline{r}}$ and any solution of
$\dot{x}=f(x,\omega(\overline{r}t))$ starting at time $s\geq
T(\gamma,\overline{r})$ at $x_{0}\in\mathbb{R}^{N}$ with
$|\widetilde{x}_{\omega,\overline{r}}(s)-x_{0}|<\varepsilon$, converges to
zero as $t\to\infty$. Moreover, assume that if $\widetilde{x}_{\gamma,r}$ has
a rate-induced tipping point $r^{*}_{\gamma}$, such tipping point is unique
and transversal. The following statements hold true.

  * 1.

For all $r>r^{*}_{\gamma}$ there is $\delta_{1}>0$ such that if
$\omega\in\Omega(\gamma,t_{0})$ and
$\|\gamma-\omega\|_{\mathcal{C}(\mathbb{R},\Lambda)}<\delta_{1}$ then $r\geq
r^{*}_{\omega}$.

  * 2.

Assume additionally that for all $\omega\in\Omega(\gamma,t_{0})$, if
$\widetilde{x}_{\omega,r}$ undergoes a rate-induced tipping, the tipping point
$r^{*}_{\omega}$ is unique and transversal. Then, for all $r<r^{*}_{\gamma}$
there is $\delta_{2}>0$ such that if
$\|\gamma-\omega\|_{\mathcal{C}(\mathbb{R},\Lambda)}<\delta_{2}$ then $r\leq
r^{*}_{\omega}$.

###### Proof.

In order to prove 1 let us reason by contradiction. Assume that there is
$\overline{r}>r^{*}_{\gamma}$ and a sequence $(\omega_{n})_{n\in\mathbb{N}}$
in $\Omega(\gamma,t_{0})$ such that
$\|\gamma-\omega_{n}\|_{\mathcal{C}(\mathbb{R},\Lambda)}<1/n$ and
$\overline{r}<r^{*}_{\omega_{n}}$ for all $n\in\mathbb{N}$. Now, note that the
sequence of functions $f\big{(}x,\omega_{n}(t)\big{)}$ converges to
$f\big{(}x,\gamma(t)\big{)}$ in the compact-open topology. Therefore, thanks
to Kamke’s lemma (see Sell [69]) the sequence of solutions
$(\widetilde{x}_{\omega_{n},\overline{r}})_{n\in\mathbb{N}}$ converges, up to
a sub-sequence, to $\widetilde{x}_{\gamma,\overline{r}}$ uniformly on compact
intervals contained in the maximal interval of definition of the latter. In
fact, due to construction,
$\widetilde{x}_{\omega_{n},\overline{r}}(t)=\widetilde{x}_{\gamma,\overline{r}}(t)$
for all all ${n\in\mathbb{N}}$ and $t\leq t_{0}$, and thus the uniform
convergence holds on all intervals of the form $(-\infty,T]$ with $T$ in the
maximal interval of definition of $\widetilde{x}_{\gamma,\overline{r}}$.

Now, for all ${n\in\mathbb{N}}$ let $\lambda_{n}\in\Lambda$ be such that
$\lim_{t\to\infty}\omega_{n}(t)=\lambda_{n}$. Consequently,
$\lim_{n\to\infty}\lambda_{n}=\lambda_{\infty}$. Moreover, denoted by
$x^{s}_{\lambda_{n}}$ the stable hyperbolic equilibrium of
(6.1)${}_{\lambda_{n}}$ such that
$\lim_{t\to\infty}\widetilde{x}_{\omega_{n},\overline{r}}(t)=x^{s}_{\lambda_{n}}$,
since $\overline{r}<r^{*}_{\omega_{n}}$, we have that

$\lim_{n\to\infty}x^{s}_{\lambda_{n}}=\lim_{n\to\infty}h_{\lambda_{n}}(x^{s}_{0})=h_{\lambda_{\infty}}(x^{s}_{0})=x_{\lambda_{\infty}}^{s}.$

Therefore, the sequence
$(\widetilde{x}_{\omega_{n},\overline{r}})_{n\in\mathbb{N}}$ is in fact
uniformly bounded on $\mathbb{R}$ and thus it must be
$\beta_{\gamma,\overline{r}}=\infty$, and furthermore,
$\widetilde{x}_{\gamma,\overline{r}}$ is bounded. Now, since the
$\omega$-limit set of $\widetilde{x}_{\gamma,\overline{r}}$ is nonempty,
thanks to the cocycle property, the continuity of the skew-product flow
induced by (6.2)${}_{\overline{r}}$ (see [69]), and the fact that the dynamics
of (6.1)${}_{\lambda_{\infty}}$ is completely determined by hyperbolic fixed
points, there is a fixed point $\overline{x}_{\lambda_{\infty}}$ of
(6.1)${}_{\lambda_{\infty}}$ such that
$\lim_{t\to\infty}\widetilde{x}_{\gamma,\overline{r}}(t)=\overline{x}_{\lambda_{\infty}}$.
Moreover, by assumption, it must be
$|\overline{x}_{\lambda_{\infty}}-x_{\lambda_{\infty}}^{s}|>\varepsilon>0$. On
the other hand, notice that, for the $\varepsilon>0$ given by the assumptions,
since $|\gamma(\overline{r}t)-\lambda_{\infty}|$ tends to zero as
$t\to\infty$, there are $t_{1}=t_{1}(\gamma,\overline{r},\varepsilon)>0$, and
due to the persistence of hyperbolic solutions (see Theorem 3.8 in [62]), a
solution $\overline{\varphi}$ of (6.2)${}_{\overline{r}}$ such that
$\overline{\varphi}$ is defined at least in $[t_{1},\infty)$, and
$\|\overline{x}_{\lambda_{\infty}}-\overline{\varphi}\|_{\mathcal{C}([t_{1},\infty),\mathbb{R}^{N})}<\varepsilon/5$.
Furthermore, since
$\lim_{t\to\infty}\widetilde{x}_{\gamma,\overline{r}}(t)=\overline{x}_{\lambda_{\infty}}$,
there is $t_{2}=t_{2}(\gamma,\overline{r},\varepsilon)\geq t_{1}>0$ such that
$\|\widetilde{x}_{\gamma,\overline{r}}-\overline{\varphi}\|_{\mathcal{C}([t_{2},\infty),\mathbb{R}^{N})}<\varepsilon/5$.
Now, consider $s>\max\\{t_{2},T(\overline{r})\\}$. Thanks to the continuous
variation of the solutions,
$|\widetilde{x}_{\omega_{n},\overline{r}}(s)-\widetilde{x}_{\gamma,\overline{r}}(s)|<\varepsilon/5$
for all $n\in\mathbb{N}$ sufficiently big. However, again for the persistence
of the hyperbolic solutions, there is $n_{0}\in\mathbb{N}$ such that for all
$n>n_{0}$ a solution $\overline{\psi}_{n}$ of
$\dot{x}=f(x,\omega_{n}(\overline{r}t))$ exists which is defined at least in
$[t_{2},\infty)$ and
$\|\overline{\psi}_{n}-\overline{\varphi}\|_{\mathcal{C}([t_{2},\infty),\mathbb{R}^{N})}<\varepsilon/5$
and
$\lim_{t\to\infty}|\widetilde{x}_{\omega_{n},\overline{r}}(t)-\overline{\psi}_{n}(t)|=0$.
However, at least for one $n>n_{0}$ this is in contradiction with the fact
that
$\lim_{t\to\infty}\widetilde{x}_{\omega_{n},\overline{r}}(t)=x^{s}_{\lambda_{n}}$
because for all $t>\max\\{t_{2},T(\overline{r})\\}$

$\begin{split}0<\varepsilon<|\overline{x}_{\lambda_{\infty}}-x^{s}_{\lambda_{\infty}}|<&|\overline{x}_{\lambda_{\infty}}-\overline{\phi}(t)|+|\overline{\phi}(t)-\overline{\psi}_{n}(t)|+\\\
+&|\overline{\psi}_{n}(t)-\widetilde{x}_{\omega_{n},\overline{r}}(t)|+|\widetilde{x}_{\omega_{n},\overline{r}}(t)-x^{s}_{\lambda_{n}}|+|x^{s}_{\lambda_{n}}-x^{s}_{\lambda_{\infty}}|<\varepsilon.\end{split}$

Therefore, 1 must be true.

In order to prove 2, we shall reason similarly; once more by contradiction.
Assume that there is $\overline{r}<r^{*}_{\gamma}$ and a sequence
$(\omega_{n})_{n\in\mathbb{N}}$ in $\Omega(\gamma,t_{0})$ such that
$\|\gamma-\omega_{n}\|_{\mathcal{C}(\mathbb{R},\Lambda)}<1/n$ and
$\overline{r}>r^{*}_{\omega_{n}}$ for all $n\in\mathbb{N}$. Therefore, thanks
to Kamke’s lemma (see Sell [69]) the sequence of solutions
$(\widetilde{x}_{\omega_{n},\overline{r}})_{n\in\mathbb{N}}$ converges, up to
a sub-sequence, to $\widetilde{x}_{\gamma,\overline{r}}$ uniformly on compact
intervals contained in the maximal interval of definition of the latter.
Additionally, as for 1, for the $\varepsilon>0$ given by the assumptions,
since $|\gamma(\overline{r}t)-\lambda_{\infty}|$ tends to zero as
$t\to\infty$, there are $t_{1}=t_{1}(\gamma,\overline{r},\varepsilon)>0$ and a
solution $\overline{\varphi}$ of (6.2)${}_{\overline{r}}$ such that
$\overline{\varphi}$ is defined at least in $[t_{1},\infty)$, and
$\|x_{\lambda_{\infty}}^{s}-\overline{\varphi}\|_{\mathcal{C}([t_{1},\infty),\mathbb{R}^{N})}<\varepsilon/3$.
Moreover, since $\overline{r}<r^{*}_{\gamma}$, there is also $t_{2}\geq t_{1}$
such that
$\|\widetilde{x}_{\gamma,\overline{r}}-\overline{\varphi}\|_{\mathcal{C}([t_{2},\infty),\mathbb{R}^{N})}<\varepsilon/3$.
On the other hand, the persistence of hyperbolic solutions guarantees also
that there is $n_{0}\in\mathbb{N}$ such that for all $n>n_{0}$ a solution
$\overline{\psi}_{n}$ of $\dot{x}=f(x,\omega_{n}(\overline{r}t))$ exists which
is defined at least in $[t_{2},\infty)$ and
$\|\overline{\psi}_{n}-\overline{\varphi}\|_{\mathcal{C}([t_{2},\infty),\mathbb{R}^{N})}<\varepsilon/3$.
In particular, thanks to the assumptions, for all $n>n_{0}$, it must be
$\lim_{t\to\infty}|\overline{\psi}_{n}(t)-x_{\lambda_{n}}^{s}|=0$, where
$x_{\lambda_{n}}^{s}=h_{\lambda_{n}}(x_{\lambda_{\infty}}^{s})$. Now, consider
$t>t_{2}$. From the continuous variation of the solutions, we obtain the
contradiction. Indeed, there must be at least one $n>n_{0}$ such that

$0\leq|\widetilde{x}_{\omega_{n},\overline{r}}(t)-\overline{\psi}_{n}(t)|\leq|\widetilde{x}_{\omega_{n},\overline{r}}(t)-\widetilde{x}_{\gamma,\overline{r}}(t)|+|\widetilde{x}_{\gamma,\overline{r}}(t)-\overline{\varphi}(t)|+|\overline{\varphi}(t)-\overline{\psi}_{n}(t)|,$

and the left-hand side can be made smaller than $\varepsilon$ by choosing $n$
sufficiently large. In this case, however, the assumptions guarantee that
$\lim_{t\to\infty}|\widetilde{x}_{\omega_{n},\overline{r}}(t)-x_{\lambda_{n}}^{s}|=0$
which is against the fact that $\overline{r}>r^{*}_{\omega_{n}}$ for all
$n\in\mathbb{N}$. Therefore,we obtain the thesis. ∎

##  7\. Example

In this section, we consider a classic one-dimensional population model with
an Allee effect, i.e. with declining individual fitness for low population
densities (see Courchamp et al. [13]). Firstly, we derive the values of the
indicators presented in Sections 3 to 6. Then, we perform a parametric study
of their behaviour, and finally compare the resilience of chosen populations
with concrete parameter settings.

### The model

We shall study the parametric differential problem

$\dot{x}=f(x,r,L)=rx\left(1-\frac{x}{K}\right)\left(\frac{x}{L}-1\right),$
(7.1)

where $x\in\mathbb{R}$ represents the population size, $r>0$ the intrinsic
growth rate, $K>0$ the carrying capacity of the environment, and $0<L\leq K$
an Allee threshold. For the sake of convenience, we will set the carrying
capacity to $K=1$. For $0<L<1$ the problem has three equilibria, at $x_{0}=0$,
$x^{*}=1,$ and $x_{1}=L$. The first two are stable and the latter unstable.
Therefore, $L$ determines a threshold below which the population can not
persist. On the other hand, when $L=1$, the system undergoes a transcritical
bifurcation and the equilibria $x^{*}$ and $x_{1}$ collide. As the attractor
of interest, we shall consider $x^{*}=1$.

###  7.1. Resilience of the system

In this subsection, we apply the presented indicators to the considered model
$\dot{x}=f(x,r,L)$. When the calculations can not be performed explicitly, the
software Matlab is used for numerical computation and simulation. In
particular, we use the function ode45 with the options on the relative and
absolute tolerance respectively set to RelTol=1e-12 and AbsTol=1e-12.

#### Local indicators

The characteristic return time of $x^{*}=1$ (see Definition 3.1 and Remark
3.2) is easily obtained as the opposite of the reciprocal of the derivative of
$f$ at $x^{*}$: $T_{R}(1)=-1/{\textnormal{D}}f(1)=L/(r-rL).$ On the other
hand, notice that for scalar systems, the amplification envelope is trivially
determined—in our case $\rho(1,t)=e^{{\textnormal{D}}f(1)t}$—and does not
provide further information. In fact, thanks to Proposition 3.12, we
immediately have that
$-R_{0}=I_{\mathcal{S}}=I_{\mathcal{D}}=-{\textnormal{D}}f(1)$.

#### Basin shape indicators

The distance to threshold (Definition 4.3) for the attractor is
$DT(x^{*})=1-L$. Note also that considering the set $C=[0,1],$ the latitude in
volume (Definition 4.9) and the basin stability (Definition 4.11) (with a
uniform density $\rho(x)=\chi_{[0,1]}$) with respect to $C$ both yield the
same value as the distance to threshold. Furthermore, the latitude in width
(Definition 4.1) has an infinite value since the basin is unbounded on the
right-hand side.

#### Nonlinear transient dynamics

The return time (Definition 5.1) does not admit an explicit closed form.
Therefore, we employ a routine to estimate it numerically. In Figure 12 we
show an approximation of the mean return time from the set $C=(L,1)$ estimated
via a Monte Carlo simulation (uniform sampling on the set $(L+10^{-7},1)$ with
$10000$ points) depending on the parameters $r\in[0.01,0.5]$ and
$L\in[0.5,1]$. For all points, the return time is approximated via a finite-
time integration ending when the considered trajectory achieves $1-10^{-10}$.
The considered problem can also be interpreted as a gradient system. Its
resistance $W$ (Definition 5.3) is obtained as

$W(x^{*})=\int_{L}^{1}f(x,r,L)\,dx=-\frac{(L-1)^{3}(L+1)r}{12L}.$

The intensity of attraction $\mathcal{I}$ (Definition 5.6) is obtained by
calculating the maximum of $f$ on the interval $[L,1].$ We show the behaviour
of $\mathcal{I}(x^{*})$ depending on $r\in[0.01,0.5]$ and $L\in[0.5,1]$ in
Figure 13.

#### Variation of parameters

The distance to bifurcation (Definition 6.1) for this problem is $D_{\rm
bif}(x^{*})=1-L$ and therefore it coincides with the distance to threshold. In
order to study the resistance and elasticity of $x^{*}=\\{1\\}$, we will
consider a perturbation of the system’s environmental capacity. Particularly,
we aim to record the system’s response to a reduction of the carrying capacity
to $K_{p}=0.9,$ for an interval of time $[0,T]$, with $T=10$. If $K_{p}<L,$
the perturbed model is not ecologically meaningful. Hence, we will restrict
the analysis only to $L<0.9.$ Since the problem at hand is scalar, the
resistance (Definition 6.2) can be calculated simply as
$R(x^{*},K_{p},T)=1-\phi_{K_{p}}(T,x^{*})$. On the other hand, the elasticity
(Definition 6.3) corresponds to the supremum of the weighted vector field
$-f(x)/(1-x)$ over the interval $[\phi_{K_{p}}(T,x^{*}),1),$ which is attained
at $\phi_{K_{p}}(T,x^{*}).$ In regards to persistence (Definition 6.5), fixed
$K_{p}>L,$ the persistence of $x^{*}$ with respect to the stress intensity
$K_{p}$ is $P_{K_{p}}=+\infty,$ since the perturbed trajectories can not be
driven out of the basin of attraction $(K_{p},+\infty)\subset(L,+\infty).$ On
the other hand, the persistence of $x^{*}$ with respect to any stress duration
$T>0$ is $P_{T}(x^{*})=1-L$, which again corresponds to the distance to
threshold but also to the distance to bifurcation.

Figure 12. Comparison of the indicators for the attractor $x^{*}=1$ of the
model (7.1) upon the variation of the parameters $r\in[0.01,0.5]$ and
$L\in[0.5,0.95]$. From top to bottom: dominant eigenvalue
$EV=\frac{1}{T_{R}}$, distance to threshold $DT$, and mean return time
$T_{R}^{mean}(1,(L+10^{-7},1))$.

####  7.1.1. Parametric study

We perform a parametric analysis of (7.1), with the aim of portraying the
behaviour of some of the considered indicators as the transcritical
bifurcation point $L=1$ is approached (results are summarized in Figures 12
and 13).

We focus on the parameter subspace given by $r\in[0.01,0.5]$,
$L\in[0.5,1-0.05]$ and calculate seven indicators (see previous Subsection 7.1
for details): reciprocal of characteristic return time $EV=\frac{1}{T_{R}}$,
distance to threshold $DT$, mean return time $T_{R}^{mean}(1,(L+10^{-7},1)),$
resistance $W,$ intensity of attraction $\mathcal{I}$, resistance $R$ and
elasticity $E$. The last two indicators are calculated for the perturbation of
the environmental capacity $K_{p}=0.9$ for a stress period $T=10.$

Figure 13. Continuation of the comparison of the indicators for different
parameter choices of the model (7.1): resistance (potential) $W,$ intensity of
attraction $\mathcal{I},$ resistance $R$ and elasticity $E$ (please note that
in order to have higher resilience corresponding to a higher numerical value,
we worked with the reciprocal of $R,E$). We consider perturbation of the
environmental capacity to $K=0.9,$ for time $t=[0,10].$ On the left, we see
the scaled indicator values for a given parametric subspace
$[0.01,0.5]=r\times L=[0.5,0.89].$ In the middle there is a heat map for the
parametric subspace, where yellow indicates the highest and blue the lowest
values of the indicators (see the color scale on the right). Black lines
correspond to the contour lines, which mark the parameter combinations with
the same values of the indicator. We see that resistance and elasticity behave
in an opposite manner, elasticity being higher for lower values of $\\{r,L\\}$
while resistance being higher.

For a better comparison of the indicators, the estimates shown in Figures 12
and 13 are scaled to lie in $[0,1]$ by $(x-x_{\min})/(x_{\max}-x_{\min}),$
where $x$ represents the indicator, and minimum $x_{\min}$ and maximum
$x_{\max}$ corresponds to the local minimum and maximum of the indicator on
the considered subset of parameters. For the same reason, instead of
$T_{R},R,E$, we represent their reciprocal so that the value $1$ always
corresponds to the highest estimate of resilience. Additionally, for return
time, resistance and elasticity we further restrict the parameter space so to
guarantee that these indicators are well-defined.

It is possible to appreciate that upon fixing a value $r>0$ and increasing
$L$, all the considered indicators capture a loss of resilience for the system
as the bifurcation point approaches. The only exception is given by the
resistance which increases since the rate of exponential decay towards the
attractor decreases in the nearby of the bifurcation point. This phenomenon is
also directly related to the so-called slowing-down effect [72, 67]. On the
other hand, there are also parameter combinations, indicated by the black
contour lines, where the indicator stays constant.

####  7.1.2. Comparison of five species

Finally, we aim to showcase the behaviour of the indicators with respect to
five different species growth strategies for the Allee effect model (7.1). The
different values of the parameters as well as the respective graphs of
$f(x,r,L)$ for $x\in[0,1]$ are shown in Figure 14. Given the previous
considerations, we shall focus on the following indicators of resilience: the
opposite of the dominant eigenvalue $EV=1/T_{R},$ distance to threshold $DT,$
mean return time $T_{R}(1,(L+10^{-7},1)),$ resistance (potential) of a
gradient system $W,$ intensity of attraction $\mathcal{I},$ resistance $R$ and
elasticity $E.$

_Species_ | _r_ | _L_
---|---|---
Species $1$ | $0.5$ | $0.2$
Species $2$ | $1.3$ | $0.3$
Species $3$ | $2.5$ | $0.4$
Species $4$ | $5$ | $0.6$
Species $5$ | $10$ | $0.7$

Figure 14. Five population strategies with different choices of parameters in
equation (7.1). The table on the left-hand side, contains the parameters for
each species. The plot on the right-hand side, showcases the graphs of
$f(x,r,L)$ for $x\in[0,1]$ depending on the chosen values of parameters $r$
and $L$.

The numerical values for these indicators are shown in Figure 15. In all
cases, $x^{*}=1$ is the reference attractor. Each row corresponds to a single
indicator, each column to a species. The chromatic scale applies row-wise as
the values of the indicators have not been normalized. Darker tones of blue
correspond to higher values of resilience. It is apparent how answering the
question ”Which species is the most resilient?” is not trivial. The ranking of
the resilience of each species varies depending on the considered indicator.

In particular, it is interesting to notice that, not only there is no overall
agreement on the resilience ranking, but some indicators ($EV$ and $DT$ for
example) provide antipodal responses. It is possible to single out three main
clusters: a) $EV$ and $1/T_{R}^{mean}$, b) $1/E$, $I$ and $W$, c) $1/R$ and
$DT$. With the exception of species ”4” all the other species are considered
as the most resilient for at least one of the represented indicators.

Figure 15. Each row represents one indicator and each column one species. The
chromatic scale applies row-wise as the values of the indicators have not been
normalized. Darker tones of blue correspond to higher values of resilience.
With the exception of species 4 all the other species are considered as the
most resilient for at least one of the represented indicators.

Figure 16. On the left-hand side, the resilience boundaries of five species
(see Figure 14) modelled through (7.1). On the right-hand side, first row, the
numerical values of the areas between the resilience boundary and the
respective distance to threshold calculated for each species. The second row
contains the same values now divided by the respective distances to threshold
$DT(x^{*})=1-L$. The color gradient should be intended row-wise. Darker shades
of blue correspond to a higher resilience.

Finally, in Figure 16 we can see the resilience boundaries for the considered
five species. As studied in [50], it is possible to extrapolate a measure of
resilience for the attractor by calculating the area between the resilience
boundary and the line $\kappa=DT(x^{*})=1-L$ (see Subsection 5.3). However,
since the distances to threshold for the five species are different it seems
natural to ask if the obtained values are representative and comparable. The
right panel of Figure 16 shows the numerical values obtained for the indicator
before and after dividing them by the respective distances to threshold.

## References

  * [1] H. M. Alkhayuon and P. Ashwin, Rate-induced tipping from periodic attractors: Partial tipping and connecting orbits, Chaos: An Interdisciplinary Journal of Nonlinear Science, 28 (2018), p. 033608.
  * [2] J.-F. Arnoldi, M. Loreau, and B. Haegeman, Resilience, reactivity and variability: A mathematical comparison of ecological stability measures, Journal of theoretical biology, 389 (2016), pp. 47–59.
  * [3] P. Ashwin, C. Perryman, and S. Wieczorek, Parameter shifts for nonautonomous systems in low dimension: bifurcation-and rate-induced tipping, Nonlinearity, 30 (2017), p. 2185.
  * [4] J. Beddington, C. Free, and J. Lawton, Concepts of stability and resilience in predator-prey models, The Journal of Animal Ecology, (1976), pp. 791–816.
  * [5] B. E. Beisner, C. L. Dent, and S. R. Carpenter, Variability of lakes on the landscape: roles of phosphorus, food webs, and dissolved organic carbon, Ecology, 84 (2003), pp. 1563–1575.
  * [6] N. Berglund and B. Gentz, Noise-induced phenomena in slow-fast dynamical systems: a sample-paths approach, Springer Science & Business Media, 2006.
  * [7] S. R. Carpenter, D. Ludwig, and W. A. Brock, Management of eutrophication for lakes subject to potentially irreversible change, Ecological Applications, 9 (1999), pp. 751–771.
  * [8] A. Carvalho, J. A. Langa, and J. Robinson, Attractors for infinite-dimensional non-autonomous dynamical systems, vol. 182, Springer Science & Business Media, 2012.
  * [9] E. A. Coddington and N. Levinson, Theory of Ordinary Differential Equations, Tata McGraw-Hill Education, 1955.
  * [10] C. Conley, The gradient structure of a flow: I, Ergodic Theory and Dynamical Systems, 8 (1988), pp. 11–26.
  * [11] W. A. Coppel, Dichotomies in stability theory, vol. 629, Springer, 2006\.
  * [12] K. L. Cottingham and S. R. Carpenter, Predictive indices of ecosystem resilience in models of north temperate lakes, Ecology, 75 (1994), pp. 2127–2138.
  * [13] F. Courchamp, L. Berec, and J. Gascoigne, Allee Effects in Ecology and Conservation, Oxford University Press, 2008.
  * [14] D. Deangelis, S. Bartell, and A. Brenkert, Effects of nutrient recycling and food-chain length on resilience, American Naturalist, 134 (1989), pp. 778–805.
  * [15] D. L. DeAngelis, P. J. Mulholland, A. V. Palumbo, A. D. Steinman, M. A. Huston, and J. W. Elwood, Nutrient dynamics and food-web stability, Annual Review of Ecology and Systematics, 20 (1989), pp. 71–95.
  * [16] B. Dennis, L. Assas, S. Elaydi, E. Kwessi, and G. Livadiotis, Allee effects and resilience in stochastic populations, Theoretical Ecology, 9 (2016), pp. 323–335.
  * [17] I. Dobson, Computing a closest bifurcation instability in multidimensional parameter space, Journal of nonlinear science, 3 (1993), pp. 307–327.
  * [18] M. Evans and T. Swartz, Approximating Integrals via Monte Carlo and Deterministic Methods, vol. 20, OUP Oxford, 2000.
  * [19] A. C. Fassoni and H. M. Yang, An ecological resilience perspective on cancer: insights from a toy model, Ecological Complexity, 30 (2017), pp. 34–46.
  * [20] C. W. Gardiner et al., Handbook of stochastic methods, vol. 3, springer Berlin, 1985.
  * [21] M. Gladwell, The tipping point: How little things can make a big difference, Little, Brown, 2006.
  * [22] H.-R. Gruemm, Definitions of resilience, International Institute for Applied Systems Analysis research report RR-76-005, (1976).
  * [23] L. H. Gunderson, Ecological resilience—in theory and application, Annual review of ecology and systematics, 31 (2000), pp. 425–439.
  * [24] J. K. Hale, Ordinary differential equations, re krieger pub, Co., Huntington, NY, 46 (1980).
  * [25] G. W. Harrison, Stability under environmental stress: resistance, resilience, persistence, and variability, The American Naturalist, 113 (1979), pp. 659–669.
  * [26] J. Harte, Ecosystem stability and the distribution of community matrix eigenvalues, in Theoretical Systems Ecology. Advances and Case Studies, Academic Press New York, 1979, pp. 453–466.
  * [27] M. A. Harwell and H. L. Ragsdale, Eigengroup analyses of linear ecosystem models, Ecological Modelling, 7 (1979), pp. 239 – 255.
  * [28] C. S. Holling, Resilience and stability of ecological systems, Annual review of ecology and systematics, 4 (1973), pp. 1–23.
  * [29] C. S. Holling, Engineering resilience versus ecological resilience, Engineering within ecological constraints, 31 (1996), p. 32.
  * [30] A. R. Ives and S. R. Carpenter, Stability and diversity of ecosystems, science, 317 (2007), pp. 58–62.
  * [31] R. Johnson, R. Obaya, S. Novo, C. Núñez, and R. Fabbri, Nonautonomous Linear Hamiltonian Systems: Oscillation, Spectral Theory and Control, Springer, 2016.
  * [32] B. Kaszás, U. Feudel, and T. Tél, Tipping phenomena in typical dynamical systems subjected to parameter drift, Scientific reports, 9 (2019), pp. 1–12.
  * [33] R. Kerswell, C. C. Pringle, and A. Willis, An optimization approach for analysing nonlinear stability with transition to turbulence in fluids as an exemplar, Reports on Progress in Physics, 77 (2014), p. 085901.
  * [34] V. V. Klinshov, V. I. Nekorkin, and J. Kurths, Stability threshold approach for complex dynamical systems, New Journal of Physics, 18 (2015), p. 013004.
  * [35] P. E. Kloeden and M. Rasmussen, Nonautonomous Dynamical Systems, no. 176, American Mathematical Soc., 2011.
  * [36] C. Kuehn, A mathematical framework for critical transitions: Bifurcations, fast–slow systems and stochastic dynamics, Physica D: Nonlinear Phenomena, 240 (2011), pp. 1020–1035.
  * [37] , A mathematical framework for critical transitions: normal forms, variance and applications, Journal of Nonlinear Science, 23 (2013), pp. 457–510.
  * [38] Y. A. Kuznetsov, Elements of applied bifurcation theory, vol. 112, Springer Science & Business Media, 2013.
  * [39] I. P. Longo, S. Novo, and R. Obaya, Topologies of $L_{loc}^{p}$ type for Carathéodory functions with applications in non-autonomous differential equations, Journal of Differential Equations, 263 (2017), pp. 7187–7220.
  * [40] , Weak topologies for Carathéodory differential equations: continuous dependence, exponential dichotomy and attractors, Journal of Dynamics and Differential Equations, 31 (2019), pp. 1617–1651.
  * [41] I. P. Longo, C. Núnez, R. Obaya, and M. Rasmussen, Rate-induced tipping and saddle-node bifurcation for quadratic differential equations with nonautonomous asymptotic dynamics, SIAM Journal on Applied Dynamical Systems, 20 (2021), pp. 500–540.
  * [42] N. L. Lundström, How to find simple nonlocal stability and resilience measures, Nonlinear dynamics, 93 (2018), pp. 887–908.
  * [43] N. L. Lundström and J.-O. Aidanpää, Dynamic consequences of electromagnetic pull due to deviations in generator shape, Journal of Sound and Vibration, 301 (2007), pp. 207–225.
  * [44] R. May, Stability and Complexity in Model Ecosystems, Princeton Landmarks in Biology, Princeton University Press, 2019.
  * [45] R. P. McGehee, Some metric properties of attractors with applications to computer simulations of dynamical systems, preprint, 290 (1988).
  * [46] P. J. Menck, J. Heitzig, J. Kurths, and H. J. Schellnhuber, How dead ends undermine power grid stability, Nature communications, 5 (2014), p. 3969.
  * [47] P. J. Menck, J. Heitzig, N. Marwan, and J. Kurths, How basin stability complements the linear-stability paradigm, Nature physics, 9 (2013), p. 89.
  * [48] K. Meyer, Metric properties of attractors for vector fields via bounded, nonautonomous control., Ph.D. thesis, (2019).
  * [49] K. Meyer, A mathematical review of resilience in ecology, Natural Resource Modeling, 29 (2016), pp. 339–352.
  * [50] K. Meyer, A. Hoyer-Leitzel, S. Iams, I. Klasky, V. Lee, S. Ligtenberg, E. Bussmann, and M. L. Zeeman, Quantifying resilience to recurrent ecosystem disturbances using flow–kick dynamics, Nature Sustainability, 1 (2018), p. 671.
  * [51] K. J. Meyer and R. P. McGehee, Intensity—a metric approach to quantifying attractor robustness in odes, arXiv preprint arXiv:2012.10786, (2020).
  * [52] J. Milnor, On the concept of attractor, in The Theory of Chaotic Attractors, Springer, 1985, pp. 243–264.
  * [53] C. Mitra, A. Choudhary, S. Sinha, J. Kurths, and R. V. Donner, Multiple-node basin stability in complex dynamical networks, Physical Review E, 95 (2017), p. 032317.
  * [54] C. Mitra, J. Kurths, and R. Donner, An integrative quantifier of multistability in complex systems based on ecological resilience, Scientific Reports, 5 (2015).
  * [55] M. G. Neubert and H. Caswell, Alternatives to resilience for measuring the responses of ecological systems to perturbations, Ecology, 78 (1997), pp. 653–665.
  * [56] R. O’Neill, Ecosystem persistence and heterotrophic regulation, Ecology, 57 (1976), pp. 1244–1253.
  * [57] G. H. Orians, Diversity, stability and maturity in natural ecosystems, in Unifying concepts in ecology, Springer, 1975, pp. 139–150.
  * [58] G. Peterson, C. R. Allen, and C. S. Holling, Ecological resilience, biodiversity, and scale, Ecosystems, 1 (1998), pp. 6–18.
  * [59] S. Pimm and J. Lawton, Number of trophic levels in ecological communities, Nature, 268 (1977), p. 329.
  * [60] S. Pimm and J. H. Lawton, On feeding on more than one trophic level, Nature, 275 (1978), pp. 542–544.
  * [61] S. L. Pimm, The complexity and stability of ecosystems, Nature, 307 (1984), p. 321.
  * [62] C. Pötzsche, Nonautonomous continuation of bounded solutions, Commun. Pure Appl. Anal, 10 (2011), pp. 937–961.
  * [63] Z. Ratajczak, P. D’Odorico, S. L. Collins, B. T. Bestelmeyer, F. I. Isbell, and J. B. Nippert, The interactive effects of press/pulse intensity and duration on regime shifts at multiple scales, Ecological Monographs, 87 (2017), pp. 198–218.
  * [64] N. Rooney, K. McCann, G. Gellner, and J. C. Moore, Structural asymmetry and the stability of diverse food webs, Nature, 442 (2006), pp. 265–269.
  * [65] R. J. Sacker and G. R. Sell, A spectral theory for linear differential systems, Journal of Differential Equations, 27 (1978), pp. 320–358.
  * [66] M. Scheffer, Critical transitions in nature and society, vol. 16, Princeton University Press, 2020.
  * [67] M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. Van Nes, M. Rietkerk, and G. Sugihara, Early-warning signals for critical transitions, Nature, 461 (2009), pp. 53–59.
  * [68] P. Schultz, P. J. Menck, J. Heitzig, and J. Kurths, Potentials and limits to basin stability estimation, New Journal of Physics, 19 (2017), p. 023005.