On the transition between autonomous and nonautonomous systems: The case of FitzHugh–Nagumo’s model

Longo, I. P.; Queirolo, E.; Kuehn, C.

doi:10.1063/5.0234833

This work deals with a parametric linear interpolation between an autonomous FitzHugh–Nagumo model and a nonautonomous skewed problem with the same fundamental structure. This paradigmatic example allows us to construct a family of nonautonomous dynamical systems with an attracting integral manifold and a hyperbolic repelling trajectory located within the nonautonomous set enclosed by the integral manifold. Upon the variation of the parameter the integral manifold collapses, the hyperbolic repelling solution disappears and a unique globally attracting hyperbolic solution arises in what could be considered yet another pattern of nonautonomous Hopf bifurcation. Interestingly, the three phenomena do not happen at the same critical value of the parameter, yielding, thus, an example of a nonautonomous bifurcation in two steps. We provide a mathematical description of the dynamical objects at play and analyze the periodically forced case via rigorously validated continuation.

We investigate the transition between autonomous and nonautonomous systems using the FitzHugh–Nagumo model as a paradigmatic example. We analyze the structure of the global attractor of the system while a parameter varies and unveil a nonautonomous bifurcation in which an attracting set collapses, a repelling solution disappears, and a unique globally attracting solution arises. We provide rigorous mathematical justification of the involved dynamical objects and validated numerics of the bifurcation via continuation techniques in the periodic case.

I. INTRODUCTION

Multiple time scale and non-autonomous nonlinear dynamics have been separated in many instances. However, these classes share many physical features, e.g., a slowly drifting variable with constant speed directly leads to a non-autonomous system. Furthermore, a system with time-dependent forcing can be made formally autonomous by augmenting time as a phase space variable. If the forcing is very fast or very slow, a natural multiple time scale structure emerges. This leads to the natural approach that one might want to interpolate between these two system classes in the hope of relating results from their respective analysis. One goal of this work is to follow this approach and uncover not only its advantages but also its substantial challenges. We are going to focus on one of the simplest paradigmatic nonlinear dynamics models and consider a parametric linear interpolation between a FitzHugh–Nagumo model (including the classic Van der Pol system) and a skewed nonautonomous differential problem sharing the same fundamental structure. Specifically, we consider the following nonautonomous planar system of ordinary differential equations:

{\begin{cases} \dot{x} = (1 - δ) y + δ v (t) - \frac{x^{3}}{3} + x, \\ \dot{y} = ε (a - x - b y), \\ v \in L^{\infty} (R, R), \end{cases}

(1.1)

with

0 \leq δ \leq 1, ε, b \geq 0, a \in R

⁠. The first immediate observations concern the family of dynamical systems as

δ

varies in

[0, 1]

⁠. If

δ = 0

⁠, is the classic autonomous FitzHugh–Nagumo model with constant forcing

a

⁠,

{\begin{cases} \dot{x} = y - \frac{x^{3}}{3} + x, \\ \dot{y} = ε (a - x - b y) . \end{cases}

(1.2)

Note that if

b = 0

⁠, then reduces to the Van der Pol model,¹ which is an example of an oscillator with nonlinear damping, energy being dissipated at large amplitudes and generated at low amplitudes. Typically, this feature guarantees the existence of a limit cycle in a suitable region of the parameter space, that is, sustained oscillations around a state at which energy generation and dissipation balance. Van der Pol’s equation is a classic topic of nonlinear dynamics and it has been extensively studied (see, for example, Refs. 2 and 3, and references therein). The FitzHugh–Nagumo model⁴ generalizes the Van der Pol model, enabling the appearance of a larger variety of oscillations and modeling several processes in biophysics in a bit more detail, such as excitable neuron dynamics. Unsurprisingly, the bifurcation analysis of the FitzHugh–Nagumo model is even richer and more involved. We point the interested reader toward the detailed presentation in Ref. 5 for further information. Interestingly, the FitzHugh–Nagumo model is known to showcase dynamical complexity when subjected to nonautonomous forcing.⁶ More in general, the appearance of almost automorphic dynamics and strange nonchaotic attractors in almost periodically forced, damped nonlinear oscillators is a classic field of study.^7–12

Important for our work is the existence of an attractive limit cycle for certain values of the parameters of the autonomous model. This is guaranteed by the occurrence of a supercritical Hopf bifurcation where a stable equilibrium becomes unstable and the attractive periodic orbit appears.

For

δ = 1

⁠, is the skewed nonautonomous problem,

{\begin{cases} \dot{x} = v (t) - \frac{x^{3}}{3} + x, \\ \dot{y} = ε (a - x - b y) . \end{cases}

(1.3)

By the term “skewed,” we mean that is composed of an uncoupled cubic nonautonomous equation and by a further scalar linear problem in

y

that can be easily integrated with the variation of constants formula once the solution of the first equation is known. The dynamical scenarios of the first equation in have been studied by Tineo^13,14 and more recently also by Dueñas et al.¹⁵ under far more general assumptions. In particular, for any

v

bounded and uniformly continuous, there are at least one and at most three hyperbolic solutions of

\dot{x} = v (t) - \frac{x^{3}}{3} + x,

i.e., trajectories that are hyperbolic in the non-autonomous sense as explained in Sec. II B. Additionally, if three hyperbolic solutions exist, then the upper and the lower ones are attractive and the one in between is repelling, whereas if only one hyperbolic solution exists, then it is always attractive. We will work under the assumption that this equation has indeed only one hyperbolic solution and, therefore, it is attracting. Additionally, if

v

is assumed periodic, then we show that an attractive periodic solution exists for .

If $0 < δ < 1$ ⁠, system (1.1) can be understood (bio-)physically as a FitzHugh–Nagumo type of model with time-dependent injected current. More theoretically, (1.1) corresponds to an interpolation between two different FitzHugh–Nagumo models via the parameter $δ$ ⁠. The two limits $δ = 0$ and $δ = 1$ correspond to classical multiple time scale and non-autonomous systems, respectively.

The variation of the parameter $δ$ connects dynamical systems living on different phase spaces: the autonomous problem for $δ = 0$ has phase-space $R^{2}$ ⁠, whereas the nonautonomous problems for $δ > 0$ are defined on the extended phase-space $R \times R^{2}$ ⁠. In fact, for $δ > 0$ ⁠, standard techniques from topological dynamics allow us to construct an autonomous flow on an extended phase space, seen as a bundle of a base and fibers, which is dynamically much more interesting. The base corresponds to a functional metric space on which the time shift (Bebutov flow) produces a well-defined dynamical system. Under certain assumptions on $v$ ⁠, the base of this bundle is, in fact, compact, unlike the real line of the trivial augmented phase space $R \times R^{2}$ ⁠. The fiber, instead, is the original Euclidean space now parametrized over the elements of the base. The resulting flow is called a skew-product flow. A brief summary of these facts is presented in Sec. II A. As an example, if $v$ is periodic, a standard trick allows us to reduce the nonautonomous problems to autonomous ones on the extended phase space $S^{1} \times R^{2}$ —the variable of time can be considered modulo the period of $v (t)$ which is fixed a priori. In order to provide a unified framework, also the phase space of (1.2) is accordingly augmented. This fact, however, has important consequences on the interpretation of the dynamics of (1.2). For example, in the simplest nontrivial case with $v$ nonconstant periodic, an equilibrium point and a periodic orbit of the autonomous FitzHugh–Nagumo (or the Van der Pol) model will, respectively, correspond to a periodic orbit and to the boundary of a two-torus in the augmented phase space $S^{1} \times R^{2}$ ⁠.

Consequently, even in the case where

y

is substituted by the

y

-coordinate of a twin independent Fitzhugh–Nagumo model,

\begin{matrix} {\begin{cases} \dot{x} = v (t) - \frac{x^{3}}{3} + x, \\ \dot{y} = ε (a - x - b y), \end{cases} \\ where (u (t), v (t)) solves {\begin{cases} \dot{u} = v - \frac{u^{3}}{3} + u, \\ \dot{v} = ε (a - u - b v), \end{cases} \end{matrix}

the dynamical scenarios of and will be fundamentally different. In the first case, the global attractor is a two-torus, while in the second case is a single attracting periodic orbit.

We investigate the change of the global attractor and specifically the occurrence of bifurcations upon the variation of the parameter $δ \in [0, 1]$ ⁠. In the periodic case, this also corresponds to the breaking of the torus. The description of the global attractor is carried out by combining analytical techniques and rigorous validated numerics. The former allows us a precise account of the fine structure when $δ$ assumes values close to $0$ and $1$ ⁠. The latter shed light on the bifurcation. The pattern we encounter can be considered yet another version of nonautonomous Hopf bifurcation^16–22 in the sense that the global attractor with positive measure in the fiber existing close to $δ = 0$ discontinuously collapses giving birth to a unique attracting hyperbolic solution. In the process, a hyperbolic completely repelling solution, contained in the original global attractor, also disappears. Unlike an autonomous Hopf bifurcation, the diameter of the global attractor does not shrink close to the bifurcation point but it abruptly vanishes, something not uncommon in the nonautonomous setting.²² Reference 19 contains, to the best of our knowledge, the most similar setup, in that their object of study is a nonautonomous perturbation of the normal form of the autonomous Hopf bifurcation. We, however, do not assume to work in a range of parameters close to an autonomous Hopf bifurcation but rather take the existence of a periodic limit cycle for the autonomous problem as a starting point and vary a parameter, which gradually leads to a nonautonomous problem that is not a perturbation of the autonomous one.

The work is organized as follows. In Sec. II, we present some background notions on nonautonomous dynamical systems, such as the skew-product formalism and the process formalism, the notions of attractors for nonautonomous systems, hyperbolic solutions, and integral manifolds. This section might seem a little technical. However, as we aim to link non-autonomous and autonomous multiple time scale systems, we first have to carefully set up the background for both areas.

Section III deals with the description of the global attractor of (1.1) whose existence for all $δ \in [0, 1]$ is proved in Sec. III A; we distinguish three cases: in Sec. III B, we deal with the case $δ = 0$ employing classical techniques from multiple time scale dynamics²³ to recall the existence of a strongly attracting periodic orbit, which, in turn, induces a trivial attracting integral manifold in the augmented phase space. Its persistence under small perturbations is used to inherit the existence of a nontrivial integral manifold for strictly positive values of $δ$ close to zero. In Sec. III C, we apply recent results on the characterization of hyperbolic solutions for scalar d-concave ordinary differential equations¹⁵ to prove the existence of a unique attracting hyperbolic solution for (1.3). Its persistence under small perturbations is used to guarantee the same dynamical scenario for values of $δ$ close to one. In Sec. III D, we give a rigorous description of the finite-time behavior of the attractor for intermediate values of $δ$ using recent results in singularly perturbed nonautonomous systems.²⁴

Section IV tackles the analysis of the nonautonomous bifurcation leading to the collapse of the integral manifold, the disappearance of the repelling hyperbolic solution, and the appearance of the globally attracting hyperbolic solution. After providing numerical evidence of the two-step bifurcation for the periodically and quasi-periodically forced cases, we focus on a four-dimensional autonomous system where the forcing is generated by a twined FitzHugh–Nagumo system whose $y$ -component is fed into the two-dimensional parametric problem as $v$ ⁠. In this setup, we are able to employ a numerically validated approach to carry out the continuation of the hyperbolic periodic trajectories of the two-dimensional forced problem and study their stability by approximating the Lyapunov spectra for the variational equations along them.

II. BASIC FACTS ON FLOWS, ATTRACTORS, HYPERBOLIC SOLUTIONS, AND INTEGRAL MANIFOLDS

A. Skew-product flows and processes

Let

v : R \to R

be a bounded and uniformly continuous function. A standard technique allows us to construct an autonomous flow for on a compact metric space. Let

H (v)

denote the closure with respect to the compact-open topology of the set of time-translations of

v

⁠, i.e., the set

cls {v_{t} : R \to R ∣ t \in R}

⁠, where

v_{t} (s) = v (t + s)

⁠.

H (v)

is called the hull of

v

⁠. It is well-known that

H (v)

is a compact metric space and the shift map

σ : R \times H (v) \to H (v)

⁠,

(t, v) \mapsto σ_{t} (v) = v_{t}

defines a continuous flow on

H (v)

called the Bebutov flow. This flow is also minimal provided that

v

has certain recurrent behavior in time, such as periodicity, almost periodicity, or other weaker properties of recurrence. If

v

is almost periodic, then the flow on

H (v)

is minimal and almost periodic and, thus, uniquely ergodic. Let us now consider

w \in H (v)

⁠,

X_{0} = (x_{0}, y_{0}) \in R^{2}

⁠, and denote by

X (t, w, X_{0})

the maximal solution of the Cauchy problem,

{\begin{cases} \dot{x} = (1 - δ) y + δ w (t) - \frac{x^{3}}{3} + x, \\ \dot{y} = ε (a - x - b y), \\ x (0) = x_{0}, y (0) = y_{0}, \end{cases}

(2.1)

with

0 \leq δ \leq 1, ε, b \geq 0, a \in R .

The map

\begin{matrix} Φ : R \times H (v) \times R^{2} \to H (v) \times R^{2}, \\ (t, w, X_{0}) \mapsto Φ (t, w, X_{0}) = (w_{t}, X (t, w, X_{0})) \end{matrix}

is a continuous skew-product (autonomous) flow on

H (v) \times R^{2}

on the compact metric base

H (v)

⁠. In the special case that

v

is periodic, the identification

H (v) = S^{1}

holds. If

v

is quasi-periodic,

H (v) = T^{n}

for some

n \in N

⁠, but already for

v

almost periodic,

H (v)

is not a manifold in general. We note that, in fact, the construction of a continuous skew-product flow for holds under substantially weaker assumptions akin to local integrability (see Ref. 25, and references therein). For example, induces a continuous skew-product flow also if

v

is essentially bounded, i.e.,

v \in L^{\infty} (R, R)

⁠, where an integral topology is used in place of the compact-open topology [see Theorem 5.6(i) in Ref. 26]. In this case, a solution is intended in the sense of Carathéodory, that is, a locally absolutely continuous function that solves the associated integral problem. Although the compactness of the base does not hold in general under such mild assumptions of regularity, it does hold in the specific case of (for general characterizations of compactness with respect to strong and weak integral topologies, see Refs. 27 and 28).

Proposition II.1

Let

v \in L^{\infty} (R, R)

⁠. Then,

H (v)

is a compact metric space and induces a continuous skew-product flow on a compact metric base

\begin{matrix} Φ : R \times H (v) \times R^{2} \to H (v) \times R^{2}, \\ (t, w, X_{0}) \mapsto Φ (t, w, X_{0}) = (w_{t}, X (t, w, X_{0})) . \end{matrix}

Proof.

Let

f (t, X)

denote the right-hand side of with

X = (x, y)

⁠, and let

J_{X} f (t, X)

denote the Jacobian of

f

with respect to the variable

X

⁠. Additionally, let

‖ v ‖_{\infty} = V

⁠. For any

r > 0

⁠, there are positive constants

m = m (V, r)

and

l = l (r)

such that for almost every

t \in R

⁠,

| f (t, X) | \leq m

and

\begin{aligned} | f (t, X_{1}) - f (t, X_{2}) | \\ \leq | X_{1} - X_{2} | \int_{0}^{1} | J_{X} f (t, X_{1} + θ (X_{2} - X_{1})) | d θ \leq l | X_{1} - X_{2} | \end{aligned}

for every

X, X_{1}, X_{2} \in B_{r} (0)

⁠. Therefore, in the terminology of Carathéodory functions,

f

has

L_{l o c}^{1}

-uniformly integrable

m

-bounds and

L_{l o c}^{1}

-bounded

l

-bounds—in fact, essentially bounded

m

-bounds and

l

-bounds. Therefore, Theorem 4.1 in Ref. 27 guarantees that

H (v)

is a compact metric space, whereas Theorem 3.5(i) and Theorem 2.33 in Ref. 25 show that the induced skew-product flow is continuous.

Instead of considering a flow on an extended phase space, sometimes it is useful to look at a nonautonomous problem in terms of an evolution operator, which now depends on two parameters, the initial and the final time. Such evolution operator goes under the name of process or cocycle.

Definition II.2

Given a metric space $(E, d)$ ⁠, a process is a family of continuous maps ${S (t, s) ∣ t \geq s} \subset C (E)$ satisfying

$S (t, t) x = x$ for every $t \in R$ and $x \in E$ ⁠,
$S (t, s) = S (t, r) S (r, s)$ ⁠, for every $t \geq r \geq s$ ⁠, and
$(t, s, x) \mapsto S (t, s) x$ is continuous for every $t \geq s$ and $x \in E$ ⁠.

Remark II.3

If

\dot{X} = f (t, X)

is such that for any

X_{0} \in R^{N}

and any

t_{0} \in R

⁠, there is a unique solution

X (\cdot, t_{0}, X_{0})

and it is defined on

[t_{0}, \infty)

⁠, then a process is induced by

S_{f} (t + t_{0}, t_{0}) X_{0} = X (t + t_{0}, t_{0}, X_{0}),

(2.2)

where

t \geq 0

and

t_{0} \in R

⁠.

B. Attractors for nonautonomous differential equations

The double dependence on time of nonautonomous problems gives rise to a new type of attractivity—pullback attractivity—which is independent from the classic notion of forward attractivity. The two notions are indistinguishable for autonomous systems. Next, we recall the definitions of global and pullback attractor for skew-product flows as well as the definition of pullback attractor for a process on $R^{2}$ ⁠. For an extensive treatment on the matter, we point the reader to Refs. 29 and 30, and the references therein.

(Pullback attractor for a process)

Definition II.4

(Pullback attractor for a process)

A family of subsets $A = {A (t) ∣ t \in R}$ of the phase space $R^{2}$ is said to be a pullback attractor for the process $S (\cdot, \cdot)$ if

$A (t)$ is compact for each $t \in R$ ⁠;
$A (\cdot)$ is invariant, that is, $S (t, s) A (s) = A (t)$ for all $t \geq s$ ⁠;
for each $t \in R$ ⁠, $A (t)$ pullback attracts bounded sets at time t, i.e., for any bounded set $B \subset R^{2}$ ⁠, one has

$lim_{s \to - \infty} d (S (t, s) B, A (t)) = 0,$

where $d (A, B)$ is the Hausdorff semi-distance between two nonempty sets $A, B \subset R^{2}$ ⁠, i.e., $d (A, B) := sup_{x \in A} inf_{y \in B} d (x, y)$ ⁠; and
$A$ is the minimal family of closed sets with property (iii).

The pullback attractor is said to be bounded if

⋃_{t \in R} A (t)

is bounded.

(Pullback and global attractors for a skew-product flow)

Definition II.5

(Pullback and global attractors for a skew-product flow)

Assume that for any $w \in H (v)$ and any $X_{0} \in R^{2}$ ⁠, the solution $X (\cdot, w, X_{0})$ of (2.1) is defined on $[0, \infty)$ ⁠, i.e., the induced skew-product flow is defined on $R^{+} \times H (v) \times R^{2}$ ⁠. (This fact will be proved true in Theorem III.1.)

A family $\hat{A} = {A_{w} \subset R^{2} ∣ w \in H (v)}$ of nonempty, compact sets of $R^{2}$ is said to be a pullback attractor for the skew-product flow if it is invariant, i.e.,

$X (t, w, A_{w}) = A_{w_{t}} for each t \geq 0 and w \in H (v),$

(2.3)
and, for every nonempty bounded set $D$ of $R^{2}$ and every $w \in H (v)$ ⁠, one has

$lim_{t \to \infty} d (X (t, w_{- t}, D), A_{w}) = 0,$

(2.4)
where $d (A, B)$ denotes the Hausdorff semi-distance of two nonempty sets $A$ ⁠, $B$ of $R^{2}$ ⁠.

A pullback attractor for the skew-product flow is said to be bounded if

$⋃_{w \in H (v)} A_{w} is bounded.$
A compact set $A$ of $H (v) \times R^{2}$ is said to be a global attractor for the skew-product flow if it is the maximal nonempty compact subset of $H (v) \times R^{2}$ ⁠, which is $Φ$ -invariant, i.e.,

$Φ (t, A) = A for each t \geq 0,$

and attracts all compact subsets $D$ of $H (v) \times R^{2}$ ⁠, i.e.,

$lim_{t \to \infty} d (Φ (t, D), A) = 0,$

where now $d (B, C)$ denotes the Hausdorff semi-distance of two nonempty sets $B$ ⁠, $C$ of $H (v) \times R^{2}$ ⁠.

Remark II.6

We wish to make explicit the relation between the family of sets

A_{w} \subset R^{2}

⁠, with

w \in H (v)

⁠, in the previous definition of the pullback attractor of the skew-product flow

\hat{A} = {A_{w} \subset R^{2} ∣ w \in H (v)}

and the pullback attractors of the processes induced by each

w \in H (v)

⁠. Recalling that for every

w \in H (v)

⁠, one has

S_{w} (s + r, r) x_{0} = x (s, w_{r}, x_{0}),

(2.5)

then, for every nonempty bounded set

D

of

R^{2}

⁠, becomes

lim_{t \to \infty} d (S_{w} (0, - t) D, A_{w}) = 0,

(2.6)

which implies that, for the given process

S_{w} (\cdot, \cdot)

⁠,

A_{w}

pullback attracts bounded sets at time

0

⁠. In turn, for every

τ \in R

⁠, one can write for

w_{τ}

and

A_{w_{τ}}

instead of

w

and

A_{w}

⁠, respectively, i.e.,

lim_{t \to \infty} d (S_{w_{τ}} (0, - t) D, A_{w_{τ}}) = 0.

(2.7)

Nevertheless, using twice, we also have

S_{w_{τ}} (0, - t) D = x (t, w_{τ - t}, D) = S_{w} (τ, τ - t) D,

and, thus, can be written as

lim_{t \to \infty} d (S_{w} (τ, τ - t) D, A_{w_{τ}}) = 0,

(2.8)

which implies that, for the given process

S_{w} (\cdot, \cdot)

⁠,

A_{w_{τ}}

pullback attracts bounded sets at time

τ

⁠. Therefore, as a consequence of the invariance contained in and of the fact that

A_{w}

is taken compact for any

w \in H (v)

⁠, we deduce that the process

S_{w} (\cdot, \cdot)

has a pullback attractor. In particular, if

A_{w} = {A (τ) = A_{w_{τ}} ∣ τ \in R}

is the minimal family of closed sets satisfying , then

A_{w}

is the pullback attractor for the process

S_{w} (\cdot, \cdot)

⁠. On the other hand, if for any

w \in H (v)

⁠, the induced process

S_{w} (\cdot, \cdot)

has a pullback attractor and

A_{w}

denotes the section at time

0

of the pullback attractor of

S_{w} (\cdot, \cdot)

⁠, then one has that

\hat{A} = {A_{w} ∣ w \in H (v)}

is a pullback attractor for the skew-product flow on

H (v) \times R^{2}

⁠.

A special case of nonautonomous (local) attractors is given by hyperbolic attracting solutions, i.e., globally defined solutions which uniformly pullback and forward attract nearby solutions at an exponential rate. Let us be more precise. Let

L : R \to R^{2 \times 2}

be a continuous function, and recall that a linear homogeneous system

\dot{z} = L (t) z

(2.9)

is uniformly asymptotically stable if and only if the trivial solution

x \equiv 0

is exponentially stable (see, for example, Theorem 4.4.2 in Ref. 31), i.e., there are constants

K > 0

and

γ > 0

such that

‖ Ψ (t, s) ‖ \leq K e^{- γ (t - s)} for all t \geq s,

(2.10)

where

Ψ (t, s)

is the principal matrix solution of at

s \in R

⁠. In this case, one also says that is exponentially stable. Since we are in a finite dimensional case and all matrix norms are equivalent, we do not explicitly choose one, unless otherwise stated. On the other hand, is asymptotically uniformly stable in backward time if there are constants

K > 0

and

γ > 0

such that

‖ Ψ (t, s) ‖ \leq K e^{γ (t - s)} for all s \geq t .

(2.11)

Given a nonlinear problem $\dot{x} = f (t, x)$ (possibly of Carathéodory type), where $f$ is continuously differentiable in $x$ for almost every $t \in R$ ⁠, a globally defined solution $\tilde{x}$ is said to be hyperbolic attracting (respectively, hyperbolic completely repelling) if the corresponding variational equation $\dot{z} = J_{f} (t, \tilde{x} (t)) z$ satisfies an inequality of the type (2.10) [respectively, (2.11)]. Note that hyperbolic attracting solutions are both (locally) pullback and forward attracting. Hyperbolic saddle solutions are also possible and this requires a generalization of the previous formulas through the concept of exponential dichotomy. This is not required for this work, and, therefore, we point the interested reader, for example, to Refs. 29 and 32 for the details and further references. An important feature of hyperbolic solutions is that they persist under small perturbations.³²

C. On Floquet’s theory and normally hyperbolic integral manifolds

In this section, we shall briefly relate the Floquet exponents of a periodic solution $γ$ of a planar system to the normal hyperbolicity of the integral manifold associated to $γ$ when the planar system is embedded in the extended phase space $R \times R^{2}$ ⁠.

Let us start by recalling Floquet’s theorem in the planar case. Every fundamental matrix solution $Y : R \to R^{2 \times 2}$ of the $T$ -periodic linear problem $\dot{y} = L (t) y$ ⁠, has the form $Y (t) = P (t) e^{t B}$ ⁠, where $P (t)$ and $B$ are $2 \times 2$ matrices, $P (t + T) = P (t)$ for all $t \in R$ and $B$ is constant.

A monodromy matrix of $\dot{y} = L (t) y$ is a nonsingular matrix $M$ associated with a fundamental matrix solution $Y (t)$ of $\dot{y} = L (t) y$ through the relation $Y (t + T) = Y (t) M$ ⁠. The eigenvalues $μ$ of a monodromy matrix are called the characteristic (or Floquet) multipliers of $\dot{y} = L (t) y$ and any $λ$ such that $μ = e^{λ T}$ is called a characteristic (or Floquet) exponent of $\dot{y} = L (t) y$ ⁠. Note that the characteristic exponents are not uniquely defined, albeit their real parts are, and they do not depend on the chosen monodromy matrix.³³

We shall now turn our attention to non-linear autonomous planar problems subjected to a small time-dependent perturbation. Consider the system of equations

\begin{matrix} \dot{x} = f (x, y) + δ f_{0} (t, x, y), \\ \dot{y} = g (x, y) + δ g_{0} (t, x, y), \end{matrix}

(2.12)

where

f, f_{0}, g, g_{0} : R^{2} \to R

are smooth. At first, consider

δ = 0

and assume that a non-constant T-periodic solution

γ : R \to R^{2}

exists. The Floquet exponents of the variational equation along

γ

are³³

μ_{1} = 1 and μ_{2} = \exp \int_{0}^{T} (\frac{\partial f}{\partial x} (γ (t)) + \frac{\partial g}{\partial y} (γ (t))) d t,

and, thus,

γ

is asymptotically stable if

μ_{2} < 1.

In such a case, a constant

c > 0

exists such that for any

(x_{0}, y_{0})

with

d ((x_{0}, y_{0}), γ (R)) < c

⁠, there is

τ = τ (x_{0}, y_{0})

such that the solution

ϕ_{0} : R \to R^{2}

with initial condition

(x_{0}, y_{0})

satisfies

| ϕ_{0} (t) - γ (t - τ) | \leq c e^{λ_{2} t}

⁠, with

λ_{2} < 0

⁠.

We shall now look at for

δ = 0

in the augmented phase space. In view of the above, the cylinder in

R \times R^{2}

⁠,

S_{0} = {(t, x, y) ∣ (x, y) = γ (θ), θ \in [0, T], t \in R},

and the torus in

S^{1} \times R^{2}

⁠,

T_{0} = {(t, x, y) ∣ (x, y) = γ (θ), θ, t mod T \in [0, T)},

are invariant for the flow induced by with

δ = 0

⁠, and asymptotically stable. In their respective spaces, they are examples of integral manifolds,³³ i.e., invariant sets for the flow induced by with a special structure. Considering a change of coordinates

(x, y) = γ (θ) + Z (θ) ρ

(a moving orthonormal system along

γ

as described in Theorem VI.1.1 in Ref. 33), the solutions of the original planar system can be described in a neighborhood of

S_{0}

(or

T_{0}

⁠) by

\begin{matrix} \dot{θ} = 1 + Θ (θ, ρ), \\ \dot{ρ} = λ_{2} ρ + R (θ, ρ), \end{matrix}

where

θ

is a radial coordinate, whereas

ρ

represents the normal deviation from

γ (θ)

⁠,

Θ (θ, ρ) = O (| ρ |)

⁠, and

R (θ, ρ) = O (| ρ |^{2})

as

ρ \to 0

⁠, and

Θ

and

R

are

T

-periodic in

θ

⁠. Details on this construction can be found in Sec. VII 1 of Ref. 33.

Importantly, such an integral manifold is robust against perturbations: let

δ \neq 0

in , where

f_{0}

and

g_{0}

are bounded in a neighborhood of

S_{0}

uniformly in

t

⁠, and

δ

is sufficiently small. Then, the above change of coordinates applied to the perturbed problem leads to

\begin{matrix} \dot{θ} = 1 + Θ (θ, ρ) + δ Θ_{0} (t, θ, ρ), \\ \dot{ρ} = λ_{2} ρ + R (θ, ρ) + δ R_{0} (t, θ, ρ), \end{matrix}

where

Θ_{0} (t, θ, 0)

and

R_{0} (t, θ, 0)

are bounded. Thus, Corollary VII.2.1 in Ref. 33 ensures the existence of a continuous function

r : R^{2} \times [0, δ_{0}] \to R^{+}

⁠,

(t, θ, δ) \mapsto r (t, θ, δ)

⁠, for some

δ_{0} > 0

⁠, bounded and Lipschitz continuous in

θ

⁠, so that

S_{δ} = {(t, θ, r (t, θ, δ)) ∣ θ \in [0, T], t \in R}

is an integral manifold of in

R \times R^{2}

⁠. Moreover, if the functions

f_{0}, g_{0}

are

τ

-periodic or almost periodic in

t

⁠, the same holds for

r (\cdot, θ, δ)

⁠.

III. DESCRIPTION OF THE GLOBAL ATTRACTOR

A. Existence of a global attractor for all δ ∈ [0, 1]

Our first result proves the dissipative nature of (1.1), which gives rise to a pullback and a global attractor for the system.

Theorem III.1

Let $v : R \to R$ be essentially bounded. Then, the following statements are true for (1.1) for any $δ \in [0, 1]$ ⁠:

there exists a unique bounded pullback attractor ${\hat{A}}^{δ} = {A_{w}^{δ} \subset R^{2} ∣ w \in H (v)}$ for the skew-product flow and it is defined by

$A_{w}^{δ} = ⋂_{τ \geq 0} \bar{⋃_{t \geq τ} x (t, w_{- t}, B_{\bar{r}}; δ)} for each w \in H (v),$

where $\bar{r} = \bar{r} (a, ε, v)$ is a positive constant, and $B_{\bar{r}}$ is the closed ball in $R^{2}$ of radius $\bar{r}$ centered at the origin; and
there is a global attractor of the skew-product flow

$A^{δ} = ⋂_{τ \geq 0} \bar{⋃_{t \geq τ} Φ_{δ} (t, H (v) \times B_{\bar{r}})} = ⋃_{w \in H (v)} {{w} \times A_{w}^{δ}} .$

Proof.

First, note that, thanks to Proposition II.1, (1.1) induces a continuous skew-product flow $Φ_{δ} : U \subset R \times H (v) \times R^{2} \to H (v) \times R^{2}$ on the compact metric space $H (v)$ ⁠. Moreover, since $v$ is essentially bounded, then every element of $H (v)$ shares the same norm in $L^{\infty}$ ⁠.

For any

w \in H (v)

⁠, considering the induced dynamical systems in polar coordinates, we have that the radial equation satisfies

\begin{aligned} \dot{r} & = r^{3} [- \frac{\cos^{4} (θ)}{3} + \frac{δ w (t) + a ε \sin (θ)}{r^{3}} \\ + \frac{\cos^{2} (θ) - ε \sin^{2} (θ) + (1 - δ - ε) \sin (θ) \cos (θ)}{r^{2}}] . \end{aligned}

Hence, there is

\bar{r} = \bar{r} (a, ε, v)

such that

\dot{r} < 0

whenever

r \geq \bar{r}

uniformly on

w \in H (v)

⁠,

δ \in [0, 1]

and

t \in R

⁠. Consequently, the flow

Φ

is uniformly bounded dissipative in the sense that all solutions of are defined on positive half-lines and uniformly ultimately bounded by the sphere of radius

r = \bar{r}

⁠. Therefore, among other references, (i) follows from Theorem 3.20 in Ref. 29. In turn, the existence of a global attractor

A^{δ}

follows from the compactness of

H (v)

and from Theorem 2.2 in Ref. 34 and, as shown in Theorem 16.2 in Ref. 30,

A_{w}^{δ}

is the section of

A^{δ}

over

w

⁠, that is,

A^{δ} = ⋃_{w \in H (v)} {{w} \times A_{w}^{δ}},

which finishes the proof.

In Secs. III B–III D, we aim to give a finer description of the global attractor. We shall first focus on the boundary cases $δ = 0$ and $δ = 1$ and nearby values of $δ$ ⁠.

B. Global attractor for δ = 0 (FitzHugh–Nagumo model) and nearby values

We briefly recall the construction of an attracting limit cycle for the FitzHugh–Nagumo model via geometric singular perturbation theory. Equation is generally referred to as in its fast time scale as opposed to its slow time scale obtained via the change of variable

τ = ε t

⁠,

{\begin{cases} ε \dot{x} = y - \frac{x^{3}}{3} + x, \\ \dot{y} = a - x - b y . \end{cases}

(3.1)

By setting

ε = 0

in , we obtain a reduced (or slow) system with an algebraic constraint,

{\begin{cases} 0 = y - \frac{x^{3}}{3} + x, \\ \dot{y} = a - x - b y . \end{cases}

(3.2)

The algebraic constraint,

S = {(x, y) \in R^{2} ∣ y = x^{3} / 3 - x},

is called the critical manifold of . Note also that the points of

S

are equilibrium points for the so-called layer problem (or fast subsystem), which is obtained by setting

ε = 0

in the fast time scale problem ,

{\begin{cases} \dot{x} = y - \frac{x^{3}}{3} + x, \\ \dot{y} = 0. \end{cases}

(3.3)

Fenichel theory³⁵ guarantees the persistence of any compact connected and normally hyperbolic subset of

S

when

ε > 0

is small enough. Recall that an invariant manifold is called normally hyperbolic if, loosely speaking, the attraction to the manifold in forward and/or backward time is stronger than the dynamics on the manifold itself; see Refs. 2 and 36–39 for the exact definition. In the case of fast-slow system like ours, normal hyperbolicity of a manifold can be verified in a simpler way: if

M

is a subset of

S

⁠, it is sufficient to show that all the points in

M

are hyperbolic equilibria of the layer problem.

Therefore, since

\frac{\partial}{\partial x} f (x, y) = 1 - x^{2} \neq 0

for all

x \in R ∖ {\pm 1}

⁠, the critical manifold is normally hyperbolic except for

p_{\pm} = (\pm 1, \mp 2 / 3)

⁠. Specifically, the sign of

\frac{\partial}{\partial x} f (x, y)

characterizes the branches

S \cap {(x, y) \in R^{2} ∣ | x | > 1}

as attracting and the branch

S \cap {(x, y) \in R^{2} ∣ | x | < 1}

as repelling, regardless of the value of

a \in R

and

b \geq 0

⁠. The so-called slow flow, that is, the dynamics of the reduced system can be obtained via the invariance equation

\begin{aligned} a - x - b y (x) = \frac{d}{d τ} y = \frac{d}{d τ} (\frac{x^{3}}{3} - x) & = (x^{2} - 1) \frac{d}{d τ} x \\ \Rightarrow \dot{x} = \frac{a - x - b (x^{3} / 3 - x)}{x^{2} - 1} . \end{aligned}

The slow flow is well-defined everywhere except for

x = \pm 1

[assuming that

a \neq \pm (1 - 2 b / 3)

]. For fixed parameters satisfying

| a | < 1 - 2 b / 3

⁠, the slow flow is such that

\dot{x} > 0

for

x < 1

and

\dot{x} < 0

for

x > 1

⁠. A desingularization at

x = \pm 1

can be obtained by rescaling time with the factor

(x^{2} - 1)

although at a cost, the reversal of time direction on the repelling branch of the critical manifold. We will momentarily ignore this issue as it has no effect on what comes next. Consider the candidate closed orbit obtained by concatenating orbits on the attracting branches of the critical manifolds with quick horizontal jumps from one of the attracting branches of

S

onto the other attracting branch of

S

at the fold points

p_{\pm} = (\pm 1, \mp 2 / 3)

⁠. In other words, the candidate orbit follows the slow flow on

S

up to a fold point and then it quickly transitions via a trajectory of the layered problem to the opposite attracting branch of

S

and so on. Under some genericity conditions, which for read as

a - 1 + 2 b / 3 \neq 0 and a + 1 - 2 b / 3 \neq 0,

(3.4)

and some assumptions on the slow flow, which for read as

a - x - b (x^{3} / 3 - x) {\begin{cases} > 0 & if x < - 1, \\ < 0 & if x > 1, \end{cases}

(3.5)

Theorem 2.1 in Ref. 23 shows that, for

ε > 0

sufficiently small, the candidate orbit perturbs into a strongly attracting periodic orbit of that lies

O (ε^{2 / 3})

-close to the candidate and approaches it in Hausdorff distance as

ε \to 0

⁠. Since in our analysis

ε > 0

is kept fixed, we shall denote this orbit by

γ_{0} (t) = (x_{0} (t), y_{0} (t))

⁠, where the subindex

0

refers to the fact that we are dealing with the case

δ = 0

⁠. The term strongly attracting refers to the fact that its non-unitary Floquet exponent is negative.

As we have seen in Sec. II C, this means that, whenever $ε > 0$ is small enough, the attracting periodic orbit of (1.2) induces an attracting integral manifold $S_{0}$ for the process associated to it on the augmented phase space $R \times R^{2}$ ⁠. In particular, by construction, it can be shown that $S_{0}$ is homeomorphic to $R \times S^{1}$ since the graph of the periodic limit cycle $γ_{0} (t)$ can be parameterized by the angle $θ \in [0, 2 π]$ modulo $2 π$ ⁠. Figure 1 depicts the strongly attracting periodic orbit of (3.1) for $a = 0$ ⁠, $b = 0.7$ and $ε = 0.1$ and the induced integral manifold in $R \times R^{2}$ ⁠.

FIG. 1.

View large Download slide

On the left-hand side: dynamics of the autonomous FitzHugh–Nagumo model (1.2) for $a = 0$ ⁠, $b = 0.7$ ⁠, and $ε = 0.1$ ⁠. The attracting periodic orbit is depicted in blue. The unstable equilibrium is depicted in red. The black dotted lines represent sample trajectories. On the right-hand side, the system is portrayed in the extended phase space $R \times R^{2}$ giving rise to a trivial integral manifold.

Furthermore, this integral manifold persists under small perturbations.³³ Hence, for

δ > 0

sufficiently small, an attractive integral manifold

S_{δ}

—also homeomorphic to

R \times S^{1}

—exists for . In accordance to Sec. II C, the equation on

S_{δ}

reduces to

\dot{θ} = 1 + δ Θ_{0} (t, θ, δ) .

In the extended notation of the skew-product flow of Sec. II A, we can consider the integral manifold

S_{δ}

parameterized with respect to the elements of

H (v)

⁠, that is,

S_{δ} = {(w, θ, r_{δ} (w, θ)) ∣ θ \in [0, 2 π], w \in H (v)} .

This approach allows to interpret the flow

φ

on

H (v) \times S_{δ}

—that is the restriction of the skew-product flow

Φ

to

H (v) \times S_{δ}

—as a circle extension of the base flow

σ_{t}

on

H (v)

in the sense that

π ° φ (t, w, p) = σ_{t} ° π (w, p),

where

π : S_{δ} \to H (v)

is the natural projection

(w, p) \mapsto π (w, p) = w

⁠. In fact, circle flows arising from almost periodically forced, damped nonlinear oscillators are a well-studied topic.^7–12 The resulting dynamical scenarios can be fairly complicated and certainly exceeding the complexity of the base flow due to the interaction between internal and forcing frequencies. For example, if the model is periodically forced, it follows from the classical Poincaré–Birkhoff–Denjoy theory that a minimal set for

φ

can be periodic, almost periodic or almost automorphic with a Cantor structure. Under almost periodic forcing, the dynamics on

H (v) \times S_{δ}

can be even more complicated. Reference 7 provides extended information in this regard. In particular, even in quasi-periodically forced nonlinear oscillators as simple as Van der Pol, numerical studies have shown the existence of so-called strange non-chaotic attractors, that is, invariant sets, which have a fractal-like geometric structure but admit no positive Lyapunov exponent.^9,40

In order to complete our analysis of the global attractor for (1.1) when $δ$ is close to zero, we note that under the assumptions (3.4) and (3.5), the autonomous FitzHugh–Nagumo model (1.2) has also one hyperbolic completely unstable equilibrium point, which lies inside the periodic orbit in $R^{2}$ ⁠. Note that an unstable hyperbolic equilibrium for an autonomous differential equation is, in particular, a completely unstable hyperbolic solution in the nonautonomous sense presented at the end of Sec. II B. In turn, the (nonautonomous notion) hyperbolicity guarantees robustness against small perturbations.³² Hence, for any $f_{0}, g_{0}$ bounded and $δ \neq 0$ sufficiently small, the hyperbolic equilibrium is perturbed into a hyperbolic completely unstable trajectory for the flow induced by (2.12). Figure 2 shows the persisting integral manifold and hyperbolic repelling solution for $δ = 0.07$ and same values for the other parameters as in Fig. 1.

FIG. 2.

Numerical approximation of the integral manifold of (1.1) for δ = 0.07, a = 0, b = 0.7, ε = 0.1, and quasi-periodic forcing v ( t ) = cos ⁡ ( 2 π t / 30 ) + sin ⁡ ( 2 π t / ( 30 5 ) ). The numerically approximated hyperbolic completely unstable trajectory is depicted in red.

View large Download slide

Numerical approximation of the integral manifold of (1.1) for $δ = 0.07$ ⁠, $a = 0$ ⁠, $b = 0.7$ ⁠, $ε = 0.1$ ⁠, and quasi-periodic forcing $v (t) = \cos (2 π t / 30) + \sin (2 π t / (30 \sqrt{5}))$ ⁠. The numerically approximated hyperbolic completely unstable trajectory is depicted in red.

C. Global attractor for δ = 1 (non-autonomous skewed problem) and nearby values

In this subsection, we show the existence of a hyperbolic attracting solution for the problem . First of all, note that the dynamical scenarios of the first equation in have been exhaustively studied by Dueñas et al.¹⁵ In particular, for any

v \in C (R, R)

bounded and uniformly continuous, there is at least one, and at most two, hyperbolic attracting solutions of

\dot{x} = v (t) - \frac{x^{3}}{3} + x

⁠, which are copies of the base, i.e., if

v

is

τ

-periodic (or almost periodic), the same holds for the hyperbolic attracting solutions. We shall be working under the assumption that there is exactly one hyperbolic attracting solution

ϕ_{1} : R \to R

for

\dot{x} = v (t) - \frac{x^{3}}{3} + x

⁠. In this case, in fact,

ϕ_{1}

is also globally attracting. As explained in Sec. II B, hyperbolicity of

ϕ_{1}

means that there are

K \geq 1

and

α > 0

such that

ψ (t, t_{0}) := \exp (\int_{t_{0}}^{t} (1 - {(ϕ_{1} (s))}^{2}) d s) \leq K e^{- α (t - t_{0})} .

(3.6)

By plugging

ϕ_{1} (t)

into , and integrating the scalar linear differential equation

\dot{y} = ε (a - ϕ_{1} (t) - b y)

⁠, with

b > 0

⁠, for initial datum

y_{0}

at time

t_{0}

⁠, we obtain the following solution to defined for all

t > t_{0}

{\begin{cases} x (t) = ϕ_{1} (t), \\ y (t, t_{0}, y_{0}) = \frac{a}{b} + e^{- b ε (t - t_{0})} (y_{0} - \frac{a}{b}) - ε \int_{t_{0}}^{t} e^{- b ε (t - s)} ϕ_{1} (s) d s . \end{cases}

The pullback limit as

t_{0} \to - \infty

gives the bounded solution

ϕ : R \to R^{2}

defined by

{\begin{cases} x (t) = ϕ_{1} (t), \\ y (t) = \frac{a}{b} - ε \int_{- \infty}^{t} e^{- b ε (t - s)} ϕ_{1} (s) d s =: ϕ_{2} (t) . \end{cases}

We shall show that this is a hyperbolic solution of . The variational equation of along

ϕ

\dot{z} = (\begin{matrix} 1 - {(ϕ_{1} (t))}^{2} & 0 \\ - ε & - b ε \end{matrix}) z

admits the principal matrix solution at

t_{0} \in R

⁠,

Ψ (t, t_{0}) = (\begin{matrix} ψ (t, t_{0}) & 0 \\ - ε \int_{t_{0}}^{t} e^{- b ε (t - s)} ψ (s, t_{0}) d s & e^{- b ε (t - t_{0})} \end{matrix}) .

Now note that, using ,

\begin{aligned} e^{- b ε (t - t_{0})} + ε \int_{t_{0}}^{t} e^{- b ε (t - s)} ψ (s, t_{0}) d s \\ \leq e^{- b ε (t - t_{0})} + \frac{ε K}{b ε - α} (e^{- α (t - t_{0})} - e^{- b ε (t - t_{0})}) \end{aligned}

and for any

b > 0

and

α > 0

⁠, we can choose

ε > 0

sufficiently small so that

b ε - α < 0

and the right-hand side of the previous inequality is smaller than

e^{- b ε (t - t_{0})}

⁠. Therefore,

‖ Ψ (t, t_{0}) ‖_{\infty} \leq K e^{- b ε (t - t_{0})},

which proves the sought-for hyperbolicity of

ϕ (t) = (ϕ_{1} (t), ϕ_{2} (t))

⁠. Moreover, it is immediate to prove that this entire bounded solution is globally attractive. Being hyperbolic, the solution

ϕ

perturbs into a hyperbolic attracting solution of for

δ

sufficiently close to

1

⁠.³²

Figure 3 shows the numerical approximation of the unique hyperbolic attracting solution to (1.3) (upper panel) and its persistence upon decreasing $δ$ (lower panel). The rest of parameters are shown in Figs. 1 and 2.

FIG. 3.

View large Download slide

Upper panel: numerical approximation of the hyperbolic attracting solution $ϕ$ for (1.1) with $δ = 1$ ⁠, $a = 0$ ⁠, $b = 0.7$ ⁠, $ε = 0.1$ ⁠, and almost periodic forcing $v (t) = \cos (2 π t / 30) + \sin (2 π t / (30 \sqrt{5}))$ as in Fig. 2. Lower panel: the same problem with $δ = 0.7$ ⁠. The orbits were obtained by numerically integrating (1.1) for 300 initial conditions starting at uniformly spaced initial times in $[- 200, - 100]$ and then discarding the transient until $t = 0$ ⁠.

Note also that, if

v

is

τ

-periodic, then using the change of variable

s = u + τ

in

\begin{aligned} ϕ_{2} (t + τ) & = \frac{a}{ε} + \int_{- \infty}^{t + τ} e^{- ε (t + τ - s)} ϕ (s) d s \\ = \frac{a}{ε} + \int_{- \infty}^{t} e^{- ε (t - u)} ϕ (u + τ) d u \\ = \frac{a}{ε} + \int_{- \infty}^{t} e^{- ε (t - u)} ϕ (u) d u = ϕ_{2} (t) \end{aligned}

yields the required periodicity for

ϕ

⁠.

D. Global attractor for intermediate values of δ

The case of intermediate values of

δ

away from 0 and 1 is considerably more complicated and highly dependent on the choice of

v

⁠. Nonetheless, we are able to investigate the finite-time behavior of singularly perturbed problems featuring an explicit dependence on the fast time such as thanks to the results developed in Ref. 24. Assume

v

is almost periodic. Then, the flow on the hull of

v

is minimal and uniquely ergodic. Let us also set some further notation. Write in the slow time scale

{\begin{cases} ε \dot{x} = (1 - δ) y + δ v (τ / ε) - \frac{x^{3}}{3} + x, \\ \dot{y} = a - x - b y, \end{cases}

and call

(x_{ε} (τ), y_{ε} (τ))

the solution of the previous problem with initial condition

(x_{0}, y_{0})

at time zero. Moreover, consider the so-called layered problems obtained from for

ε = 0

⁠,

\dot{x} = (1 - δ) y + δ v (t) - \frac{x^{3}}{3} + x

(3.7)

for

y \in [- k, k]

⁠, with

k > 0

and

δ \in [0, 1]

⁠. This is family of scalar nonautonomous d-concave differential equations dependent on the parameters

y \in [- k, k]

and

δ \in [0, 1]

⁠. Let us fix any

δ \in [0, 1]

and remove the dependence in this parameter from the upcoming notation. For each

y \in [- k, k]

⁠, induces a skew-product flow

\begin{matrix} Φ_{y} : R \times H (v) \times R \to H (v) \times R, \\ (t, w, x_{0}) \mapsto Φ_{y} (t, w, x_{0}) = (w_{t}, x (t, w, x_{0}; y)), \end{matrix}

where

x (t, w, x_{0}; y)

denotes the solution at time of

t

of [for a fixed pair

(y, δ)

] with initial condition

x (0) = x_{0}

⁠. Each of these flows admits a global attractor

A^{y} = \cup_{w \in H (v)} {w} \times A_{w}^{y}

whose structure has been described in detail in Ref. 15 and briefly recalled in Sec. III C. The dependence on

δ

for

A^{y}

is left tacit as we will not be varying

δ

within this paragraph. Numerical evidence suggests that for every

y \in [- 3, 3]

⁠, the skew-product flows

Φ_{y}

only have one copy of the base, which corresponds to the global attractor

A^{y}

⁠. In other words, there is a continuous map

η : H (v) \times [- 3, 3] \to R

such that

A_{w}^{y} = η (w, y)

⁠. Consequently, the flow

Φ_{y}

has a unique invariant measure concentrated on

A^{y}

⁠; let us call

μ^{y}

its projection on

R

⁠. Then, Theorem 4.8 in Ref. 24 allows us to draw the following conclusions. Fix

τ_{0} > 0

⁠,

for any sequence $ε_{j} \to 0$ ⁠, there is a subsequence $ε_{k} ↓ 0$ such that $y_{ε_{k}}$ converges uniformly in $[0, τ_{0}]$ to the unique solution $y_{0} (τ)$ of the differential equation

$\dot{y} = \int_{R} (a - x - b y) d μ^{y} (x), y (0) = y_{0} .$

The right-hand side above can be understood as an averaged problem in the sense of Artstein⁴¹ and Sanders et al.⁴²
Given $σ > 0$ ⁠, there exists $T = T (σ, x_{0}, A) > 0$ and $k_{0} = k_{0} (σ, T) > 0$ with $2 T < τ_{0} / ε_{k_{0}}$ such that for every $k \geq k_{0}$ ⁠,

$| x_{ε_{k}} (t) - η (v_{t}, y_{0} (ε_{k} t)) | \leq σ for all t \in [T, τ_{0} / ε_{k}] .$

In other words, as $ε \to 0$ ⁠, the slow variable converges uniformly on $[0, τ_{0}]$ to the unique solution of an averaged scalar problem anchored at the union of the hyperbolic solutions of (3.7) upon the variation of $y \in [- k, k]$ ⁠. The fast variable, in turn, approaches a path determined by the value at time $t \in [0, τ_{0} / ε_{k}]$ of the hyperbolic solution of (3.7) for $y = y_{0} (ε_{k} t)$ ⁠.

IV. A NONAUTONOMOUS BIFURCATION

In Sec. III, we have studied the nature of the global attractor of (1.1) for values of $δ \in [0, 1]$ ⁠. In Sec. III B, we have seen that for values of $δ$ sufficiently close to $δ = 0$ ⁠, there exists a normally hyperbolic integral manifold enclosing the rest of bounded solutions, including a hyperbolic repelling solution for the problem. If $δ$ takes values close to $δ = 1$ ⁠, we showed in Sec. III C that there is a unique bounded solution and it is hyperbolic attracting and a copy of the base. For intermediate values of $δ$ ⁠, we have used a nonautonomous version of Tychonov’s theorem in Sec. III D to provide a local qualitative description of the global attractor that, although rigorous, has practical limitations when a quantitative estimate is desired.

It is clear that somewhere in the transition between the extreme values of $δ$ ⁠, the normally hyperbolic integral manifold breaks up, the hyperbolic completely unstable solution ceases to exist, and a unique hyperbolic attracting solution arises. We can reasonably affirm that a bifurcation takes place.

A natural candidate is a nonautonomous Hopf bifurcation pattern.^16–22 The numerical evidence points at a two-step bifurcation.^20,43,44 At first, the integral manifold collapses giving rise to a hyperbolic stable solution, and the hyperbolic repelling solution disappears only at a higher value of $δ$ ⁠. This is appreciable in Figs. 4 and 5 for periodic and quasi-periodic forcing, respectively. Note, in particular, that at small values of delta, initial conditions in a certain region of the phase space give rise to bounded solutions converging to the hyperbolic repelling trajectory in backward time. Upon varying the parameter some of these solutions start to blow up in finite time while others keep converging to the hyperbolic repelling solution in backward time; this fact points toward the conclusion that the integral manifold does not exist anymore. For even bigger values of the parameter, all solutions cease to converge in backward time and rather blow up in finite time.

FIG. 4.

View large Download slide

A phenomenology of the bifurcation through the appearance of unbounded solutions upon the variation of $δ$ ⁠. In this example, $v (t) = \cos (2 π t / 30)$ ⁠, $a = 0$ ⁠, $b = 0.4$ ⁠, $ε = 0.1$ ⁠. 500 solutions with initial conditions uniformly distributed in $[- 200, 50]$ are used to approximate the attracting invariant set of (1.1) (a transient of at least 50 time steps is ignored and not plotted). The red trajectory represents the hyperbolic repelling solution of (1.1). It is approximated numerically integrating backward in time 20 distinct initial conditions at time $t = 20 100$ ⁠. The blue dotted lines are 50 solutions starting at time $t = 66.7$ on the circle of radius $0.3$ centered at the origin and numerically integrated backward and forward in time. At $δ = 0.1$ ⁠, these solutions are all bounded. At $δ = 0.6$ ⁠, some of these solutions remain bounded whereas others blow up in finite time suggesting that the integral manifold does not exist anymore. The hyperbolic repelling solution seems to persist, as the remaining bounded solutions converge to it in backward time. At $δ = 0.72$ ⁠, all considered solutions in blue blow up in finite time and the hyperbolic repelling solution does not seem to exist anymore.

FIG. 5.

As for Fig. 4, but for quasi-periodic forcing v ( t ) = c o s ( 2 π t / 30 ) + s i n ( 2 π t / ( 30 5 ) ) with a = 0, b = 0.4, ε = 0.1. 500 solutions with initial conditions uniformly distributed in [ − 200 , 50 ] are used to approximate the attracting invariant set of (1.1) (a transient of at least 50 time steps is ignored and not plotted). The red trajectory represents the hyperbolic repelling solution of (1.1). It is approximated numerically integrating backward in time 20 distinct initial conditions at time t = 20 100. The blue dotted lines are 50 solutions starting at time t = 66.7 on the circle of radius 0.3 centered at the origin and numerically integrated backward and forward in time. At δ = 0.1, these solutions are all bounded. At δ = 0.35, some of these solutions remain bounded whereas others blow up in finite time suggesting that the integral manifold does not exist anymore or it has been pierced. The hyperbolic repelling solution seems to persist, as the remaining bounded solutions converge to it in backward time. At δ = 0.4, all considered solutions in blue blow up in finite time and the hyperbolic repelling solution does not seem to exist anymore.

View large Download slide

As for Fig. 4, but for quasi-periodic forcing $v (t) = c o s (2 π t / 30) + s i n (2 π t / (30 \sqrt{5}))$ with $a = 0$ ⁠, $b = 0.4$ ⁠, $ε = 0.1$ ⁠. 500 solutions with initial conditions uniformly distributed in $[- 200, 50]$ are used to approximate the attracting invariant set of (1.1) (a transient of at least 50 time steps is ignored and not plotted). The red trajectory represents the hyperbolic repelling solution of (1.1). It is approximated numerically integrating backward in time 20 distinct initial conditions at time $t = 20 100$ ⁠. The blue dotted lines are 50 solutions starting at time $t = 66.7$ on the circle of radius $0.3$ centered at the origin and numerically integrated backward and forward in time. At $δ = 0.1$ ⁠, these solutions are all bounded. At $δ = 0.35$ ⁠, some of these solutions remain bounded whereas others blow up in finite time suggesting that the integral manifold does not exist anymore or it has been pierced. The hyperbolic repelling solution seems to persist, as the remaining bounded solutions converge to it in backward time. At $δ = 0.4$ ⁠, all considered solutions in blue blow up in finite time and the hyperbolic repelling solution does not seem to exist anymore.

In order to acquire a better understanding of the phenomenon, we further explore the simpler case of periodic forcing. In particular, we consider the special case where the

y

-component of the periodic attracting solution (with phase zero at time zero) of a twin FitzHugh–Nagumo system is fed into as

v

⁠,

\begin{matrix} {\begin{cases} \dot{x} = (1 - δ) y + δ v (t) - \frac{x^{3}}{3} + x, \\ \dot{y} = ε (a - x - b y), \end{cases} \\ where (u (t), v (t)) solves {\begin{cases} \dot{u} = v - \frac{u^{3}}{3} + u, \\ \dot{v} = ε (a - u - b v) . \end{cases} \end{matrix}

(4.1)

In this case, we are able to employ standard continuation techniques on the four-dimensional autonomous problem .

A. Continuation of hyperbolic solutions

Let us first introduce some notation. Considering

ϕ : R \to R

⁠,

t \mapsto ϕ (t)

⁠, a periodic function of period

T

⁠, we can write it in its Fourier expansion as

ϕ (t) = \sum_{k = - \infty}^{\infty} c_{k} e^{i k ω t}, c_{k} \in C,

with

ω = 2 π T^{- 1}

⁠. Thanks to the periodicity in

t

⁠,

c_{k} = \bar{c_{- k}}

⁠. For a given

ν

⁠, we introduce the Banach space

ℓ_{ν}^{1}

of complex sequences with the norm

‖ c ‖_{ℓ_{ν}^{1}} := \sum_{k = - \infty}^{\infty} | c_{k} | ν^{| k |} .

If

ν > 1

⁠, to have a bounded norm, the coefficients

| c_{k} |

need to decrease geometrically away from

k = 0

⁠. It is then numerically reasonable to consider the truncated Fourier sequence

\hat{c} = [c_{- k}, \dots, c_{0}, \dots, c_{k}] \in C^{2 k + 1}

for a given problem dependent

k < \infty

⁠, and consider the embedding into the

ℓ_{ν}^{1}

space.

Let

y = (c_{1}, c_{2}, c_{3}, c_{4}) \in {(ℓ_{ν}^{1})}^{4}

be the Fourier representation of the solution to the ordinary differential equation , written as

\begin{aligned} F (y) = (\begin{matrix} - ω i K c_{1} + c_{2} - \frac{1}{3} c_{1} ⋆ c_{1} ⋆ c_{1} + c_{1} \\ - ω i K c_{2} + ϵ (a - c_{1} - b c_{2}) \\ - ω i K c_{3} + (1 - δ) c_{4} + δ c_{2} - \frac{1}{3} c_{3} ⋆ c_{3} ⋆ c_{3} + c_{3} \\ - ω i K c_{4} + ϵ (a - c_{3} - b c_{4}) \end{matrix}) = 0, \end{aligned}

(4.2)

with

{(i K c)}_{k} = i k c_{k}

and

{(c ⋆ d)}_{k} = \sum_{i + j = k} c_{i} d_{j}

⁠. The latter is a finite sum if either

c

or

d

have finitely many elements.

Considering

δ

variable, the solution space for Eq. is

(ω, δ, y) \in R \times R \times {(ℓ_{ν}^{1})}^{4}

⁠, thus leaving us with two degrees of freedom, one associated with the parameter

δ

⁠, the other with time shifts of periodic orbits, i.e., if

y

is a solution so is

y (t + τ)

for any

τ

⁠. To fix this time shift, a new equation is introduced

g (y) = ⟨ y, y_{N}^{k e r} ⟩ + β = 0,

where

⟨ y, y_{N}^{k e r} ⟩ = \sum_{i = 1}^{4} ⟨ c_{i}, {\dot{c}}_{i} ⟩

⁠,

⟨ c_{i}, {\dot{c}}_{i} ⟩ = \sum_{k} {(c_{i})}_{k} {({\dot{c}}_{i})}_{k}

and

y_{N}^{k e r}

is a numerically computed kernel of the Frechet derivative of

F

⁠. Then,

β

is defined as

⟨ y_{N}, y_{N}^{k e r} ⟩

for

y_{N}

a chosen numerical solution. The choice of this phase condition follows Ref. 45.

Remark IV.1

When considering the four-dimensional problem (4.1), we can view it either as a forced system, where the time-matching of the solution $(x, y)$ and the forcing $(u, v)$ is critical, or as a four-dimensional system of ordinary differential equations, where inherently $x, y, u,$ and $v$ are defined over the same time. When working in Fourier space, this second formalism is the one used, thus removing the necessity to explicitly take care of the time-matching problem.

Once introduced the correct functional spaces and root-finding equations, the solution to

(F (ω, δ, y), g (y)) = 0

is, generically, a one-dimensional solution curve that we want to approximate by a sequence of points

{\hat{x}}_{i}

⁠, starting from a given numerical approximation

{\hat{x}}_{0} = (ω_{0}, δ_{0}, y_{0})

⁠. Let

v_{1}

be a numerically approximated tangent to the solution curve, computed as the kernel of the Frechet derivative of

(F, g)

at

{\hat{x}}_{0}

⁠. Then, an approximation for a new continuation point on the curve is

{\tilde{x}}_{1} = {\hat{x}}_{0} + ϵ v_{1}

⁠, where

ϵ

is an appropriately chosen step size. To refine this approximation, we apply a Newton algorithm of the zero-finding problem

(\begin{matrix} F (x), g (x), h (v, x) \end{matrix}) = 0,

where

h (v_{1}, x) = ⟨ v_{1}, x - {\tilde{x}}_{1} ⟩

ensures not only unicity but also that the solution found by the Newton algorithm can be written as

{\tilde{x}}_{1} + p

⁠, and the perturbation

p

is perpendicular to the tangent, as presented in Fig. 6.

FIG. 6.

View large Download slide

A sketch of the predictor–corrector method.

Remark IV.2

The predictor–corrector here introduced rely on the Newton method to converge to orbits in a given branch of solutions. This makes it particularly well-suited for the continuation of hyperbolic solutions because the Newton method’s convergence is independent of the stability properties of the solution.

Finally, the segment

x_{s}

of analytical solutions connecting

x_{0}

and

x_{1}

and approximated by

x_{0} + s (x_{1} - x_{0})

solves the full zero-finding problem

G_{s} (x) = (\begin{matrix} F (x) \\ g_{s} (x) \\ h_{s} (x) \end{matrix}),

where

s \in [0, 1]

⁠,

g_{s} (x) = ⟨ y, y_{N s}^{k e r} ⟩ - β_{s}

⁠, and

h_{s} (x) = ⟨ v_{s}, x - {\tilde{x}}_{s} ⟩

⁠, where the appendix

s

indicates a linear interpolation on

s

⁠.

Let ${\hat{x}}_{0}$ and ${\hat{x}}_{1}$ be two numerical solutions approximating $x_{0}$ and $x_{1}$ ⁠. We here want to tackle the problem of proving the existence of an analytical branch of solutions in the neighborhood of the numerical segment connecting the two approximations we computed. This is achieved here with the application of the radii polynomial approach, based on the definition of a map $T_{s} (x) = x - A_{s} G_{s} (x)$ ⁠, where $A_{s}$ is an approximation of the inverse of $D G_{s} ({\hat{x}}_{s})$ ⁠, the Frechet derivative of $G_{s}$ ⁠. If $T_{s}$ is a contraction in a neighborhood of ${\hat{x}}_{s}$ for all $s \in [0, 1]$ ⁠, then by the Banach contraction theorem a unique fixed point $T_{s} (x_{s}) = x_{s}$ exists in that neighborhood. This is equivalent to saying that $A_{s} G_{s} (x_{s}) = 0$ ⁠. If $A_{s}$ is non-singular, it then follows that $G_{s}$ has a unique solution $x_{s}$ in the same neighborhood. To make matters more precise, we sketch here the radii polynomial theorem.^46,47

Theorem IV.3

Let

T_{s}

and

{\hat{x}}_{s}

be as previously defined. Let us further define

Y \geq max_{s \in [0, 1]} ‖ T_{s} ({\hat{x}}_{s}) - {\hat{x}}_{s} ‖_{sup}

and

Z (r) \geq max_{s \in [0, 1]} max_{z, w \in B_{1} (0)} ‖ D T_{s} ({\hat{x}}_{s} + z r) w r ‖_{sup},

where the norm is defined as

‖ x ‖_{sup} = ‖ (ω, δ, c_{1}, c_{2}, c_{3}, c_{4}) ‖_{sup} = max {| ω |, | δ |, max_{i = 1, 2, 3, 4} ‖ c_{i} ‖_{ℓ_{ν}^{1}}} .

Then, if for a given

r^{*}

it holds

Y + Z (r^{*}) < r^{*}

⁠,

T_{s}

is a contraction in the neighborhood

B_{r^{*}} ({\hat{x}}_{s})

for all

s \in [0, 1]

and there exists a unique solution

x_{s}

to

G_{s} (x_{s}) = 0

in the same neighborhood.

For the implementation, the BiValVe library⁴⁸ in combination with the Intlab library⁴⁹ has been used to continue and validate the periodic orbits presented. The complete code for the validation can be found on Git.⁵⁰ Figure 7 represents the validated orbits computed. Thanks to the validation, we proved that two locally unique branches of solutions exist, one connecting the blue orbits and one the red ones. The validation of the blue orbit starts at $δ = 0.9$ and succeeds until reaching $δ = 0.2563$ ⁠, while the red orbit starts at $δ = 0.1$ and reaches $δ = 0.6101$ ⁠. Both branches have error bounds lower than $r^{*} = 2.8 \times 10^{- 4}$ ⁠. This result also proves that for all $δ$ in at least the interval $[0.2563 + r^{*}, 0.6101 - r^{*}]$ ⁠, both orbits coexist.

FIG. 7.

View large Download slide

Continuation and numerical validation of hyperbolic periodic orbits of (4.1) with $a = 0$ ⁠, $b = 0.5$ ⁠, and $ε = 0.1$ ⁠. The starting hyperbolic repelling solution is approximated numerically for $δ = 0.1$ and then continued forward. The hyperbolic attracting solution is approximated numerically for $δ = 0.9$ and then continued backward. The periodic orbits are here represented as projected on the $x y$ -plane. It is possible to appreciate that the critical values of $δ$ where they appear/disappear do not coincide.

In order to complete the continuation analysis carried out in Sec. IV A, we numerically investigate the behavior of the Floquet exponents of the variational equation along the continued hyperbolic solutions. Particularly, the Floquet exponents are determined by approximating the Lyapunov spectra. We recall that for the periodic case, and more in general for the almost periodic case, Lyapunov exponents are well-defined thanks to the Birkhoff ergodic theorem and the fact that the flow on the extended phase space $H (v) \times R^{2}$ is uniquely ergodic. We implemented the algorithm described in Sec. 2 of Ref. 51. This method uses the QR factorization of a sequence of fundamental matrix solutions of the variational equation along the periodic orbits. The estimation of the Lyapunov exponents is achieved by truncating the limit for the calculation after ten periods of the periodic solution. The results of the numerical simulation are shown in Fig. 8, where the stability of the blue orbit and the instability of the red one are confirmed.

FIG. 8.

View large Download slide

The Lyapunov spectra of the periodic solutions plotted with respect to $δ$ ⁠. The hyperbolically unstable solution is given in red, and the stable solution is given in blue.

It is also apparent that the two-step nature of the bifurcation cannot be resolved as a crossing of the zero for the spectrum at two different values of the parameters, since the simulation seems to suggest that the periodic hyperbolic orbits conserve their properties of stability until they cease to exist. In this sense, this phenomenon is different from the bifurcation discussed in Ref. 44.

V. CONCLUSIONS

This work explores a natural and yet understudied framework: the parametric transition between autonomous and nonautonomous systems. We use the FitzHugh–Nagumo model as a paradigmatic example. This is a planar singularly perturbed system, which has a periodic limit cycle and an unstable hyperbolic equilibrium in a suitable range of the parameter space. Moreover, the structure of the model allows us to modify the equation for $\dot{x}$ by introducing a parametric linear interpolation between the $y$ -component and a nonautonomous forcing $v (t)$ in a very simple fashion. An interesting case consists of taking $v$ as the $y$ -component of a twin FitzHugh–Nagumo model. In the extended phase space $R \times R^{2}$ ⁠, the autonomous FitzHugh–Nagumo model has a global attractor with positive measure, whose boundary is an integral manifold. By its nature, this attractor persists when $δ > 0$ is sufficiently small. We, therefore, obtain a nonautonomous system with an integral manifold enclosing the set of bounded solutions of the system, including a unique hyperbolic repelling solution. On the other hand, if $δ = 1$ ⁠, we obtain a nonautonomous skewed problem, made of a nonautonomous uncoupled cubic scalar equation and a linear inhomogenous scalar equation that can be explicitly solved once we have a solution to the first equation. Therefore, under suitable assumptions, we can analytically prove the existence of a unique globally attracting hyperbolic solution, which persists for values of $δ$ close to 1. This result is further supported by validated numerics techniques, allowing for a quantitative understanding of the solutions. Hence, we analytically and numerically studied the bifurcation phenomenon that leads to the collapse of the attracting integral manifold, the disappearance of the hyperbolic repelling solution, and the birth of the unique hyperbolic globally attracting solution. The paper offers a rigorous technique to generate similar bifurcation phenomena in nonautonomous problems derived by autonomous equations of Liénard type. A finer description of the bifurcation is, however, auspicable. Techniques from topological dynamics akin to the Diliberto map used in Ref. 19 seem to be promising. In this sense, our work intends to stimulate the further investigation of this, similar and even more complicated scenarios using the transition between the autonomous and the nonautonomous realms as a basis. For example, we have hereby limited the presentation to a case where the FitzHugh–Nagumo model does not feature any particularly complicated dynamics. This is, however, not always the case as this differential problem is classically known to have several (autonomous) bifurcation points and a rather complicated unfolding. The case of higher dimensional systems and more complicated attractors also deserves attention.

ACKNOWLEDGMENTS

The authors wish to deeply thank Courtney Quinn from the University of Tasmania for sharing her invaluable knowledge on the numerical approximation of Lyapunov exponents. This allowed us to rigorously complete the last part of Sec. IV.

I.P.L. acknowledges partial support by the UKRI under Grant Agreement No. EP/X027651/1, the MICIIN/FEDER (Project No. PID2021-125446NB-I00), the TUM International Graduate School of Science and Engineering (IGSSE), and the University of Valladolid under Project No. PIP-TCESC-2020.

E.Q. acknowledges full support by the DFG Walter Benjamin Programme (No. QU 579/1-1).

C.K. acknowledges partial support by a Lichtenberg Professorship of the VolkswagenStiftung and by the DFG Sachbeihilfe (Grant No. 444753754).

C.K. and I.P.L. also acknowledge partial support of the EU within the TiPES project funded by the European Union’s Horizon 2020 research and innovation programme under Grant Agreement No. 820970.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

I. P. Longo: Conceptualization (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). E. Queirolo: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). C. Kuehn: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

REFERENCES

1.

B.

Van Der Pol

, “

VII. Forced oscillations in a circuit with non-linear resistance (reception with reactive triode)

,”

London Edinb. Dublin Philos. Mag. J. Sci.

3

,

65

–

80

(

1927

).

https://doi.org/10.1080/14786440108564176

Google Scholar

Crossref

2.

C.

Kuehn

,

Multiple Time Scale Dynamics

(

Springer

,

2015

), Vol. 191.

Google Scholar

Crossref

3.

J.

Guckenheimer

and

P.

Holmes

,

Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields

(

Springer Science & Business Media

,

2013

), Vol. 42.

Google Scholar

4.

R.

FitzHugh

, “

Mathematical models of threshold phenomena in the nerve membrane

,”

Bull. Math. Biophys.

17

,

257

–

278

(

1955

).

https://doi.org/10.1007/BF02477753

Google Scholar

Crossref

5.

C.

Rocsoreanu

,

A.

Georgescu

, and

N.

Giurgiteanu

,

The FitzHugh-Nagumo Model: Bifurcation and Dynamics

(

Springer Science & Business Media

,

2012

), Vol. 10.

Google Scholar

6.

C.

Kuehn

, “

Quenched noise and nonlinear oscillations in bistable multiscale systems

,”

Europhys. Lett.

120

,

10001

(

2017

).

https://doi.org/10.1209/0295-5075/120/10001

Google Scholar

Crossref

7.

W.

Huang

and

Y.

Yi

, “

Almost periodically forced circle flows

,”

J. Funct. Anal.

257

,

832

–

902

(

2009

).

https://doi.org/10.1016/j.jfa.2008.12.005

Google Scholar

Crossref

8.

Y.

Yi

, “

On almost automorphic oscillations

,”

Fields Inst. Commun.

42

,

75

–

99

(

2004

).

Google Scholar

9.

F. J.

Romeiras

and

E.

Ott

, “

Strange nonchaotic attractors of the damped pendulum with quasiperiodic forcing

,”

Phys. Rev. A

35

,

4404

(

1987

).

https://doi.org/10.1103/PhysRevA.35.4404

Google Scholar

Crossref

10.

A. S.

Pikovsky

,

U.

Feudel

, and

S. P.

Kuznetsov

,

Strange Nonchaotic Attractors: Dynamics Between Order and Chaos in Quasiperiodically Forced Systems

(

World Scientific

,

2006

), Vol. 56.

Google Scholar

11.

T. H.

Jäger

and

G.

Keller

, “

The Denjoy type of argument for quasiperiodically forced circle diffeomorphisms

,”

Ergod. Theory Dyn. Syst.

26

(2),

447

–

465

(

2006

).

https://doi.org/10.1017/S0143385705000477

Google Scholar

Crossref

12.

T. H.

Jäger

and

J.

Stark

, “

Towards a classification for quasiperiodically forced circle homeomorphisms

,”

J. London Math. Soc.

73

(3),

727

–

744

(

2006

).

https://doi.org/10.1112/S0024610706022782

Google Scholar

Crossref

13.

A.

Tineo

, “

A result of Ambrosetti–Prodi type for first-order ODEs with cubic non-linearities. Part I

,”

Ann. Mat. Pura Appl.

182

,

113

–

128

(

2003

).

https://doi.org/10.1007/s10231-002-0038-0

Google Scholar

Crossref

14.

A.

Tineo

, “

A result of Ambrosetti–Prodi type for first-order ODEs with cubic non-linearities. Part II

,”

Ann. Mat. Pura Appl.

182

,

129

–

141

(

2003

).

https://doi.org/10.1007/s10231-002-0055-z

Google Scholar

Crossref

15.

J.

Dueñas

,

C.

Núñez

, and

R.

Obaya

, “Generalized pitchfork bifurcations in d-concave nonautonomous scalar ordinary differential equations,” arXiv:2206.14559 (2022).

16.

B.

Braaksma

and

H. W.

Broer

, “

On a quasi-periodic Hopf bifurcation

,”

Ann. Inst. Henri Poincaré C

4

,

115

–

168

(

1987

).

https://doi.org/10.1016/s0294-1449(16)30370-5

Google Scholar

Crossref

17.

B.

Braaksma

,

H.

Broer

, and

G.

Huitema

, “Toward a quasi-periodic bifurcation theory,” in Unfoldings and Bifurcations of Quasi-Periodic Tori (Memoirs of the American Mathematical Society, 1990), pp. 83–167.

18.

R.

Johnson

and

Y.

Yi

, “

Hopf bifurcation from non-periodic solutions of differential equations, II

,”

J. Differ. Equ.

107

,

310

–

340

(

1994

).

https://doi.org/10.1006/jdeq.1994.1015

Google Scholar

Crossref

19.

M.

Franca

,

R.

Johnson

, and

V.

Muñoz-Villarragut

, “

On the nonautonomous Hopf bifurcation problem

,”

Discrete Contin. Dyn. Syst. Ser. S

9

,

1119

–

1148

(

2016

).

https://doi.org/10.3934/dcdss.2016045

Google Scholar

Crossref

20.

V.

Anagnostopoulou

,

T.

Jäger

, and

G.

Keller

, “

A model for the nonautonomous Hopf bifurcation

,”

Nonlinearity

28

,

2587

(

2015

).

https://doi.org/10.1088/0951-7715/28/7/2587

Google Scholar

Crossref

21.

M.

Franca

and

R.

Johnson

, “

On the non-autonomous Hopf bifurcation problem: Systems with rapidly varying coefficients

,”

Electron. J. Qual. Theory Differ. Equ.

2019

,

1

–

24

(

2019

).

https://doi.org/10.14232/ejqtde.2019.1.56

Google Scholar

Crossref

22.

C.

Núñez

and

R.

Obaya

, “

Li–Yorke chaos in nonautonomous Hopf bifurcation patterns—I

,”

Nonlinearity

32

,

3940

(

2019

).

https://doi.org/10.1088/1361-6544/ab28ab

Google Scholar

Crossref

23.

M.

Krupa

and

P.

Szmolyan

, “

Relaxation oscillation and canard explosion

,”

J. Differ. Equ.

174

,

312

–

368

(

2001

).

https://doi.org/10.1006/jdeq.2000.3929

Google Scholar

Crossref

24.

I. P.

Longo

,

R.

Obaya

, and

A. M.

Sanz

, “Tracking nonautonomous attractors in singularly perturbed systems of ODEs with dependence on the fast time,”

J. Differ. Equ.

414

,

609–644

(

2025

).

https://doi.org/10.1016/j.jde.2024.09.011

Google Scholar

Crossref

25.

I. P.

Longo

, “

Topologies of continuity for Carathéodory differential equations with applications in non-autonomous dynamics,

” Ph.D. thesis (Universidad de Valladolid, 2018).

https://doi.org/10.35376/10324/44162

26.

I. P.

Longo

,

S.

Novo

, and

R.

Obaya

, “

Topologies of

L_{l o c}^{p}

type for Carathéodory functions with applications in non-autonomous differential equations

,”

J. Differ. Equ.

263

,

7187

–

7220

(

2017

).

https://doi.org/10.1016/j.jde.2017.08.006

Google Scholar

Crossref

27.

Z.

Artstein

, “

Topological dynamics of an ordinary differential equation

,”

J. Differ. Equ.

23

,

6

–

2

(

1977

).

https://doi.org/10.1016/0022-0396(77)90127-9

Google Scholar

Crossref

28.

R.

George

, “

Compact sets of nonlinear operators

,”

Funkc. Ekvacioj

11

,

131

–

138

(

1968

).

Google Scholar

29.

P. E.

Kloeden

and

M.

Rasmussen

,

Nonautonomous Dynamical Systems

(

American Mathematical Society

,

2011

), Vol. 176.

Google Scholar

Crossref

30.

A.

Carvalho

,

J. A.

Langa

, and

J.

Robinson

,

Attractors for Infinite-Dimensional Non-Autonomous Dynamical Systems

(

Springer Science & Business Media

,

2012

), Vol. 182.

Google Scholar

31.

L.

Adrianova

, “

Introduction to linear systems of differential equations

,”

Transl. Math. Monogr.

146

,

1

(

1995

).

Google Scholar

32.

C.

Pötzsche

, “

Nonautonomous continuation of bounded solutions

,”

Commun. Pure Appl. Anal.

10

(3),

937–961

(

2011

).

https://doi.org/10.3934/cpaa.2011.10.937

Google Scholar

Crossref

33.

J. K.

Hale

,

Ordinary Differential Equations

(

Wiley-Interscience

,

New York

,

1969

).

Google Scholar

34.

D. N.

Cheban

,

P. E.

Kloeden

, and

B.

Schmalfuß

, “

The relationship between pullback, forwards and global attractors of nonautonomous dynamical systems

,”

Nonlinear Dyn. Syst. Theory

2

,

125

–

144

(

2002

).

Google Scholar

35.

N.

Fenichel

, “

Geometric singular perturbation theory for ordinary differential equations

,”

J. Differ. Equ.

31

,

53

–

98

(

1979

).

https://doi.org/10.1016/0022-0396(79)90152-9

Google Scholar

Crossref

36.

N.

Fenichel

, “

Persistence and smoothness of invariant manifolds for flows

,”

Indiana Univ. Math. J.

21

,

193

–

226

(

1971

).

https://doi.org/10.1512/iumj.1972.21.21017

Google Scholar

Crossref

37.

N.

Fenichel

, “

Asymptotic stability with rate conditions

,”

Indiana Univ. Math. J.

23

,

1109

–

1137

(

1974

).

https://doi.org/10.1512/iumj.1974.23.23090

Google Scholar

Crossref

38.

N.

Fenichel

, “

Asymptotic stability with rate conditions, II

,”

Indiana Univ. Math. J.

26

,

81

–

93

(

1977

).

https://doi.org/10.1512/iumj.1977.26.26006

Google Scholar

Crossref

39.

M. W.

Hirsch

,

C. C.

Pugh

, and

M.

Shub

,

Invariant Manifolds

(

Springer

,

2006

), Vol. 583.

Google Scholar

40.

C.

Grebogi

,

E.

Ott

,

S.

Pelikan

, and

J. A.

Yorke

, “

Strange attractors that are not chaotic

,”

Physica D

13

,

261

–

268

(

1984

).

https://doi.org/10.1016/0167-2789(84)90282-3

Google Scholar

Crossref

41.

Z.

Artstein

, “

Averaging of time-varying differential equations revisited

,”

J. Differ. Equ.

243

,

146

–

167

(

2007

).

https://doi.org/10.1016/j.jde.2007.01.022

Google Scholar

Crossref

42.

J. A.

Sanders

,

F.

Verhulst

, and

J.

Murdock

, Averaging Methods in Nonlinear Dynamical Systems, Applied Mathematical Sciences Vol. 59 (Science+Business Media, LLC, 1985).

43.

L.

Arnold

, “Random dynamical systems,” in Dynamical Systems, Fondazione C.I.M.E. (Springer-Verlag, Berlin, 1994), pp. 1–43.

44.

R. A.

Johnson

,

P.

Kloeden

, and

R.

Pavani

, “

Two-step transition in nonautonomous bifurcations: An explanation

,”

Stochastic Dyn.

2

,

67

–

92

(

2002

).

https://doi.org/10.1142/S0219493702000297

Google Scholar

Crossref

45.

A. R.

Champneys

and

B.

Sandstede

, “

Numerical computation of coherent structures

,” in

Numerical Continuation Methods for Dynamical Systems: Path Following and Boundary Value Problems

, edited by B. Krauskopf, H. M. Osinga, and J. Galán-Vioque (

Springer Netherlands

,

Dordrecht

,

2007

), pp.

331

–

358

.

https://doi.org/10.1007/978-1-4020-6356-5_11

Google Scholar

Crossref

46.

A.

Hungria

,

J.-P.

Lessard

, and

J.

Mireles James

, “

Rigorous numerics for analytic solutions of differential equations: The radii polynomial approach

,”

Math. Comput.

85

,

1427

–

1459

(

2016

).

https://doi.org/10.1090/mcom/3046

Google Scholar

Crossref

47.

J. B.

van den Berg

and

E.

Queirolo

, “

A general framework for validated continuation of periodic orbits in systems of polynomial ODEs

,”

J. Comput. Dyn.

8

,

59

–

97

(

2020

).

https://doi.org/10.3934/jcd.2021004

Google Scholar

Crossref

48.

K.

Church

and

E.

Queirolo

, “

Computer-assisted proofs of Hopf bubbles and degenerate Hopf bifurcations

,”

J. Dyn. Differ. Equ.

36

,

3385

–

3439

(

2024

).

https://doi.org/10.1007/s10884-023-10279-x

Google Scholar

Crossref

PubMed

49.

S. M.

Rump

, “Intlab interval laboratory,” in Developments in Reliable Computing (Springer, 1999), pp. 77–104.

50.

E.

Queirolo

, https://github.com/ (2024).

51.

D.

Breda

and

D.

Liessi

, “

A practical approach to computing Lyapunov exponents of renewal and delay equations

,”

Math. Biosci. Eng.

21

,

1249

–

1269

(

2024

).

https://doi.org/10.3934/mbe.2024053

Google Scholar

Crossref

PubMed

© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/).

Published open access through an agreement with JISC Collections

On the transition between autonomous and nonautonomous systems: The case of FitzHugh–Nagumo’s model

I. INTRODUCTION

II. BASIC FACTS ON FLOWS, ATTRACTORS, HYPERBOLIC SOLUTIONS, AND INTEGRAL MANIFOLDS

A. Skew-product flows and processes

B. Attractors for nonautonomous differential equations

(Pullback attractor for a process)

(Pullback and global attractors for a skew-product flow)

C. On Floquet’s theory and normally hyperbolic integral manifolds

III. DESCRIPTION OF THE GLOBAL ATTRACTOR

A. Existence of a global attractor for all δ ∈ [0, 1]

B. Global attractor for δ = 0 (FitzHugh–Nagumo model) and nearby values

C. Global attractor for δ = 1 (non-autonomous skewed problem) and nearby values

D. Global attractor for intermediate values of δ

IV. A NONAUTONOMOUS BIFURCATION

A. Continuation of hyperbolic solutions

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

On the transition between autonomous and nonautonomous systems: The case of FitzHugh–Nagumo’s model

I. INTRODUCTION

II. BASIC FACTS ON FLOWS, ATTRACTORS, HYPERBOLIC SOLUTIONS, AND INTEGRAL MANIFOLDS

A. Skew-product flows and processes

B. Attractors for nonautonomous differential equations

(Pullback attractor for a process)

(Pullback and global attractors for a skew-product flow)

C. On Floquet’s theory and normally hyperbolic integral manifolds

III. DESCRIPTION OF THE GLOBAL ATTRACTOR

A. Existence of a global attractor for all δ ∈ [0, 1]

B. Global attractor for δ = 0 (FitzHugh–Nagumo model) and nearby values

C. Global attractor for δ = 1 (non-autonomous skewed problem) and nearby values

D. Global attractor for intermediate values of δ

IV. A NONAUTONOMOUS BIFURCATION

A. Continuation of hyperbolic solutions

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only