# Slow passage through a Hopf-like bifurcation in piecewise linear systems:
application to elliptic bursting

J. Penalva<EMAIL_ADDRESS>Departament de Matemàtiques i Informàtica & IAC3,
Universitat de les Illes Balears, Spain  M. Desroches
<EMAIL_ADDRESS>MathNeuro Team, Inria Sophia Antipolis Méditerranée
Research Centre, France  A. E. Teruel<EMAIL_ADDRESS>Departament de
Matemàtiques i Informàtica & IAC3, Universitat de les Illes Balears, Spain  C.
Vich<EMAIL_ADDRESS>Departament de Matemàtiques i Informàtica & IAC3,
Universitat de les Illes Balears, Spain

###### Abstract

The phenomenon of slow passage through a Hopf bifurcation is ubiquitous in
multiple-timescale dynamical systems, where a slowly-varying quantity
replacing a static parameter induces the solutions of the resulting slow-fast
system to feel the effect of a Hopf bifurcation with a delay. This phenomenon
is well understood in the context of smooth slow-fast dynamical systems. In
the present work, we study for the first time this phenomenon in piecewise
linear (PWL) slow-fast systems. This special class of systems is indeed known
to reproduce all features of their smooth counterpart while being more
amenable to quantitative analysis and offering some level of simplification,
in particular through the existence of canonical (linear) slow manifolds. We
provide conditions for a PWL slow-fast system to exhibit a slow passage
through a Hopf-like bifurcation, in link with the number of linearity zones
considered in the system and possible connections between canonical attracting
and repelling slow manifolds. In doing so, we fully describe the so-called
way-in/way-out function. Finally, we investigate this slow passage effect in
the Doi-Kumagai model, a neuronal PWL model exhibiting elliptic bursting
oscillations.

Keywords— Slow-fast dynamical systems, Piecewise linear dynamics, Delayed loss
of stability, Hopf-like bifurcation, Bursting oscillations.

##  1 Introduction

Differential systems involving multiple timescales are used in a wide range of
applications to describe the dynamical behavior of a given set of variables.
Let us consider a system of the form

$\begin{split}\dot{\mathbf{x}}&=F(\mathbf{x},z),\\\
\dot{z}&=\varepsilon,\end{split}$ (1)

where $\mathbf{x}\in\mathbb{R}^{q}$, ${z}\in\mathbb{R},$ $F$ is a sufficient
smooth function, and $0\leq\varepsilon\ll 1$ is a small parameter. The
presence of $\varepsilon$ in (1) induces a timescale separation such that the
variables $\mathbf{x}$ are said to be fast while $z$ is slow. If we consider
$\varepsilon=0$ in (1), which yields the so-called fast subsystem or layer
problem, then $z$ becomes a parameter for the differential equation
$\dot{\mathbf{x}}=F(\mathbf{x},z)$. An important object is the so-called
critical manifold
$\mathcal{S}_{0}=\\{(\mathbf{x},z):F(\mathbf{x},z)=\mathbf{0}\\}$, which
corresponds to the $z$-dependent family of fast subsystem equilibria
associated with system (1). In the present work, we consider cases where
$\mathcal{S}_{0}$ possesses a Hopf bifurcation at a certain value $z=z_{H}$
where the fast subsystem equilibrium loses its stability. Without loss of
generality, we further assume that the critical manifold $\mathcal{S}_{0}$ is
a graph over $z$, for $z$ within an interval $I$ that contains $z_{H}$, hence
we can write $\mathcal{S}_{0}=\\{(S_{0}(z),z):z\in I\\}$.

In the full system (1), $z$ drifts slowly in time, hence there is no Hopf
bifurcation anymore. Instead, the system undergoes a dynamic (or delayed) Hopf
bifurcation at $z=z_{H}$, also referred to as a slow passage through the Hopf
bifurcation at $z_{H}$. One may expect that orbits of system (1) with initial
conditions $(\mathbf{x}_{0},z_{0})$, where $z_{0}<z_{H}$ and $\mathbf{x}_{0}$
close to $S_{0}(z_{0})$, stay close to the critical manifold $\mathcal{S}_{0}$
for $0<t<\frac{z_{H}-z_{0}}{\varepsilon}:=t_{H}$ and move away from it as soon
as $t>t_{H}$. Instead, a surprising delay phenomenon occurs, where orbits
remain close to $\mathcal{S}_{0}$ for a long period of time $\tau^{*}>t_{H}$,
that is, past the fast subsystem bifurcation point, and then they move away
from $\mathcal{S}_{0}$ only after $z^{*}=z(\tau^{*})$ [27, 28]. This value
$z^{*}$ (and the corresponding time $\tau^{*}$) marks therefore a delay in the
response of the full system (1) to the fast subsystem instability occurring at
$z_{H}$. This delay depends on $\varepsilon$ and is expected to tend to the
bifurcation value $z_{H}$ as $\varepsilon$ tends to zero, given that, when
$\varepsilon=0$, the change of stability occurs exactly at $z_{H}$.
Nevertheless, different limiting behaviors can be obtained for
$z^{*}(\varepsilon)$, mainly depending on the regularity of the vector field;
see [26] for details. In particular there can be: no delay; a fixed maximal
delay in the interval $(z_{H},+\infty)$, which is then called a buffer point;
or the delay may be arbitrarily large.

The phenomenon of slow passage through a Hopf bifurcation (delayed loss of
stability) occurs in models coming from a vast range of application areas; for
instance, in Neuroscience [1], Physics [29], Population Biology [24] or
Chemistry [17], to name but a few. Slow-fast dynamical systems in which a
delayed loss of stability happens often display complex oscillatory behaviors
such as bursting [18, 19, 32, 4] or mixed-mode oscillations (MMOs) [1, 11, 16,
20]. In this context, delayed Hopf bifurcation can be related to the canard
phenomenon [2, 9, 11, 25].

It is well known that continuous piecewise linear (PWL) systems provide a
simple and minimal framework to study nonlinear dynamics. More specifically,
one can reproduce a wide range of multiple-timescale dynamics, in particular
in the canard regime, with PWL slow-fast systems [3, 7, 8, 10, 14, 30, 31].
Over the past few years, we have indeed shown that such systems are able to
retain all salient features of their smooth counterpart while allowing for
substantial simplification: canonical linear slow manifolds, explicit
expressions for solutions and access to quantitative information related to
various flight times along solutions.

In the present work, we analyze the slow passage through a Hopf bifurcation in
the context of PWL systems, where this bifurcation is known as Hopf-like
bifurcation. We perform this analysis using a series of minimal PWL models
presenting delayed loss of stability near a fast subsystem Hopf-like
bifurcation, and we study the limiting behavior of their respective delay as
the singular parameter $\varepsilon$ tends to zero. We then apply our results
to a PWL model of elliptic bursting oscillations as this type of bursting
relies upon a slow passage through a Hopf bifurcation [21, 34, 22].

The rest of the manuscript is organized as follows. In section 2 we analyze a
minimal PWL system with two linear regions exhibiting delayed loss of
stability due to a slow passage through a Hopf-like bifurcation, such that the
delay tends to the bifurcation value $z_{H}$ as $\varepsilon$ tends to zero.
By adding a third region allowing for the connection between the attracting
and the repelling canonical slow manifolds, we build up in section 3 a minimal
PWL model also exhibiting delayed loss of stability, but with an arbitrarily
large delay. Then, by adding an equilibrium point on the repelling canonical
slow manifold, we analyze the existence of a buffer point and the associated
maximal delay. Finally, in section 4 we apply locally all the previous results
to the Doi-Kumagai PWL elliptic bursting model, which we revisit by
considering first two and then three linearity zones near the fast subsystem
Hopf-like bifurcation point, respectively.

##  2 Two-regions system

Let us consider system

$\displaystyle\dot{x}$ $\displaystyle=f(x)-y,$ $\displaystyle\dot{y}$
$\displaystyle=x-z,$ (2) $\displaystyle\dot{z}$ $\displaystyle=\varepsilon,$

where $\varepsilon$ is a small parameter, i.e. $0\leq\varepsilon\ll 1$, and
$f$ is the continuous piecewise linear function

$f(x)=\left\\{\begin{array}[]{rl}-mx&x\leq 0,\\\ kx&x\geq
0,\end{array}\right.$ (3)

with $0<k,m<2$. It is therefore a 2-fast 1-slow PWL system of the form given
by (1), whose critical manifold $\mathcal{S}_{0}$ is the polygonal line

$\mathcal{S}_{0}=\left\\{\begin{array}[]{ll}(x,-mx,x)&\text{if }x\leq 0,\\\
(x,kx,x)&\text{if }x>0.\end{array}\right.$ (4)

When $\varepsilon=0$, the fast subsystem of (2) corresponds to the planar PWL
system

$\begin{array}[]{rcl}\dot{x}&=&f(x)-y,\\\ \dot{y}&=&x-z,\end{array}$ (5)

with $z$ acting as a parameter. This system exhibits a unique equilibrium
point at $(z,f(z))$ and, when $z\neq 0$, the Jacobian matrix evaluated at the
equilibrium point is

$\begin{pmatrix}f^{\prime}(z)&-1\\\ 1&0\end{pmatrix},$

where $f^{\prime}(z)$ is the derivative of the PWL function $f$ at $z$. This
equilibrium point changes stability at $z=0$, passing from a stable focus when
$z<0$ to an unstable focus when $z>0$. Following Theorem 5 in [15] the system
can exhibit either a supercritical or a subcritical Hopf-like bifurcation at
$z=0$, where a family of stable (resp. unstable) limit cycles emerges and
exists for all $z>0$ (resp., for all $z<0$). The amplitude of the limit cycles
along both these families grows linearly with $\lvert z\rvert$; see also [33].

In order to study the slow passage phenomenon through a Hopf-like bifurcation,
we consider the full system (2), in which $z$ is a slow variable, driving the
dynamics of the fast subsystem (5) through the Hopf-like bifurcation at $z=0$.
In fact, from (2), we only consider $z$ as a slow drift.

In the next result we describe the invariant manifolds of system (2), both
before and after perturbation in $\varepsilon$. In particular, we study the
normally hyperbolic branches of the critical manifold. We recall that normally
hyperbolic manifolds are those where the normal component of the flow
dominates the tangential component of the flow.

######  Lemma 1.

Let us consider system (2)-(3) with $m\neq k$. For $\varepsilon=0$, the system
exhibits two invariant sets. One invariant set is the critical manifold
$\mathcal{S}_{0}$, which has a normally hyperbolic attracting branch,
$S_{0}^{a}$, defined for $x<0$ and a normally hyperbolic repelling one,
$S_{0}^{r}$, defined for $x>0$. Another invariant set is a cone,
$\mathcal{C}$, with vertex at the origin and foliated by periodic orbits,
which is stable when $k<m$ and unstable when $k>m$ (see fig. 1(a)).

For small enough $\varepsilon>0$, the attracting branch of the critical
manifold $\mathcal{S}^{a}_{0}$ perturbs to a canonical attracting slow
manifold $\mathcal{S}_{\varepsilon}^{a}$, and the repelling branch of the
critical manifold $\mathcal{S}_{0}^{r}$ perturbs to a canonical repelling slow
manifold $\mathcal{S}_{\varepsilon}^{r}$ where

$\begin{array}[]{ll}\mathcal{S}_{\varepsilon}^{a}=(x,-mx-\varepsilon,x+m\varepsilon)&\text{if
}x\leq 0,\\\
\mathcal{S}_{\varepsilon}^{r}=(x,kx-\varepsilon,x-k\varepsilon)&\text{if
}x\geq 0.\end{array}$ (6)

These manifolds intersect the plane $\\{x=0\\}$ at points
$\mathbf{p}^{a}=(0,-\varepsilon,m\varepsilon)$ and
$\mathbf{p}^{r}=(0,-\varepsilon,-k\varepsilon)$, respectively (see fig. 1(b)).

###### Proof.

System (2)-(3) can be written as the 3-dimensional PWL system

$\dot{\mathbf{u}}=\left\\{\begin{array}[]{ll}A_{-}\mathbf{u}+\varepsilon\,\mathbf{e}_{3}&\text{if
}x\leq 0,\\\ A_{+}\mathbf{u}+\varepsilon\,\mathbf{e}_{3}&\text{if }x\geq
0,\end{array}\right.$

where $\mathbf{u}=(x,y,z)^{T}$, $\mathbf{e}_{3}=(0,0,1)^{T}$ and

$A_{-}=\begin{pmatrix}-m&-1&0\\\ 1&0&-1\\\ 0&0&0\end{pmatrix},\quad
A_{+}=\begin{pmatrix}k&-1&0\\\ 1&0&-1\\\ 0&0&0\end{pmatrix}.$

The critical manifold $\mathcal{S}_{0}$ decomposes as
$\mathcal{S}_{0}^{a}\cup\\{0\\}\cup\mathcal{S}_{0}^{r}$, $\mathcal{S}_{0}^{a}$
being the subset contained in the half-space $\\{x<0\\}$ and
$\mathcal{S}_{0}^{r}$ the subset contained in $\\{x>0\\}$.

The stability of the branches $\mathcal{S}_{0}^{a}$ and $\mathcal{S}_{0}^{r}$
depends on the sign of the real part of the non-null eigenvalues of $A_{-}$
and $A_{+}$. These eigenvalues are
$-\frac{m}{2}\pm\frac{\sqrt{4-m^{2}}}{2}\rm{i}$ and
$\frac{k}{2}\pm\frac{\sqrt{4-k^{2}}}{2}\rm{i}$, respectively. The real part of
the first ones is negative hence $\mathcal{S}_{0}^{a}$ is a normally
hyperbolic attracting branch, while the real part of the second ones is
positive and hence $\mathcal{S}_{0}^{r}$ is a normally hyperbolic repelling
branch. Notice that the origin is an equilibrium point located at the
switching plane, where the Jacobian matrix is not defined. Therefore, the
origin is not a normally hyperbolic point.

On the other hand, the existence of the stable (resp. unstable) cone
$\mathcal{C}$ is a consequence of the supercritical (resp. subcritical) Hopf-
like bifurcation exhibited by the fast subsystem (5) when $0<k<m<2$ (resp.
$0<m<k<2$) at $z=0$; see Theorem 5(d) in [15], where we consider
$\gamma_{L}=-\frac{m}{\sqrt{4-m^{2}}}$, $\gamma_{R}=\frac{k}{\sqrt{4-k^{2}}}$.
We point out that in the Hopf-like bifurcation, the amplitude of the limit
cycle grows linearly with $\lvert z\rvert$, which guarantees the conic shape.

When $\varepsilon>0$, both rays in expression (6) are invariant under the flow
of the system (2)-(3). Moreover, each of these rays is at distance of order
$\varepsilon$ to the respective branch of the critical manifold
$\mathcal{S}_{0}$. Following [31] we conclude that
$\mathcal{S}^{a}_{\varepsilon}$ is a canonical attracting slow manifold and
$\mathcal{S}^{r}_{\varepsilon}$ is a canonical repelling slow manifold.

The expression of the intersection point $\mathbf{p}^{a}$ (resp.
$\mathbf{p}^{r}$) of the canonical slow manifold
$\mathcal{S}^{a}_{\varepsilon}$ (resp. $\mathcal{S}^{r}_{\varepsilon}$) with
the switching plane $\\{x=0\\}$, follows straightforwardly. ∎

(a)$x$$y$$z$$\mathcal{S}_{0}^{a}$$\mathcal{S}_{0}^{r}$$\mathcal{C}$(b)$x$$y$$z$$\mathcal{S}_{\varepsilon}^{a}$$\mathcal{S}_{\varepsilon}^{r}$$\mathbf{p}^{r}$$\mathbf{p}^{a}$

Figure 1: Invariant objects before and after perturbation in
$\boldsymbol{\varepsilon}$: (a) Attracting, $\mathcal{S}_{0}^{a}$, and
repelling, $\mathcal{S}_{0}^{r}$, branches of the critical manifold, and
stable invariant cone $\mathcal{C}$ obtained after a supercritical Hopf-like
bifurcation in the fast subsystem, for $k<m$ and $\varepsilon=0$. (b)
Attracting canonical slow manifold $\mathcal{S}_{\varepsilon}^{a}$ and
repelling canonical slow manifold $\mathcal{S}_{\varepsilon}^{r}$ appearing
after perturbation of the critical manifold $\mathcal{S}_{0}$ and their
intersection points, $\mathbf{p}^{a}=(0,-\varepsilon,m\varepsilon)$ and
$\mathbf{p}^{r}=(0,-\varepsilon,-k\varepsilon)$, with the switching plane
$\\{x=0\\}$. The perturbation of the cone is not represented in the figure.

Note that both canonical slow manifolds, $\mathcal{S}_{\varepsilon}^{a}$ and
$\mathcal{S}_{\varepsilon}^{r}$ (referred as slow manifolds if no confusion
arises), do not connect on the plane $\\{x=0\\}$, since from (6) the
respective intersection points, $\mathbf{p}^{a}$ and $\mathbf{p}^{r}$, remain
at a distance $(m+k)\varepsilon$ (see lemma 1).

In theorem 1 we will see that this phenomenon forces the delayed loss of
stability to behave in a trivial way. Indeed, consider the local expression of
the flow $\mathbf{u}(t;0,(x_{0},y_{0},z_{0}))$ with respect to the initial
conditions $(x_{0},y_{0},z_{0})$ given by equations (33) and (34) in appendix
A, depending on $x_{0}<0$ and $x_{0}>0$, respectively. For any $t>0$, let
$\pi(z_{0},t)$ be the plane $\\{z=z_{0}+\varepsilon t\\}$. Given an initial
condition $(0,y_{0},z_{0})$ close to $\mathbf{p}^{r}$, we define the distance
along the plane $\pi(z_{0},t)$ from the solution
$\mathbf{u}(t;0,(0,y_{0},z_{0}))$ to the slow manifold
$\mathcal{S}_{\varepsilon}^{r}$ as

$d_{R}(t,(0,y_{0},z_{0}))=\|\mathbf{u}(t;0,(0,y_{0},z_{0}))-(\hat{x}(t),\hat{y}(t),\hat{z}(t))\|,$
(7)

where $t>0$ and $(\hat{x}(t),\hat{y}(t),\hat{z}(t))$ is the intersection
between $\mathcal{S}_{\varepsilon}^{r}$ and $\pi(z_{0},t)$. Similarly, we can
define the distance along the plane $\pi(z_{0},t)$ from the solution
$\mathbf{u}(t;0,(0,y_{0},z_{0}))$ to the slow manifold
$\mathcal{S}_{\varepsilon}^{a}$ as

$d_{L}(t,(0,y_{0},z_{0}))=\|\mathbf{u}(t;0,(0,y_{0},z_{0}))-(\hat{x}(t),\hat{y}(t),\hat{z}(t))\|,$
(8)

where $t<0$ and $(\hat{x}(t),\hat{y}(t),\hat{z}(t))$ is the intersection
between $\mathcal{S}_{\varepsilon}^{a}$ and $\pi(z_{0},t)$.

######  Lemma 2.

Let $\mathbf{u}(t;0,(0,y_{0},z_{0}))$ be the solution of system (2)-(3) with
initial condition $(0,y_{0},z_{0})$. For $t>0$, as
$\mathbf{u}(t;0,(0,y_{0},z_{0}))$ remains in the half-space $\\{x>0\\}$, the
distance along the plane $\pi(z_{0},t)$ from $\mathbf{u}(t;0,(0,y_{0},z_{0}))$
to $\mathcal{S}_{\varepsilon}^{r}$ satisfies

$d_{R}(t,(0,y_{0},z_{0}))=d_{R}(0,(0,y_{0},z_{0}))\mathrm{e}^{\frac{k}{2}t}\frac{\sqrt{C(t)}}{\sqrt{C(0)}},$
(9)

where $\frac{1}{16}<|C(t)|<2$.

For $t<0$, as $\mathbf{u}(t;0,(0,y_{0},z_{0}))$ remains in the half-space
$\\{x<0\\}$, the distance along the plane $\pi(z_{0},t)$ from
$\mathbf{u}(t;0,(0,y_{0},z_{0}))$ to $\mathcal{S}_{\varepsilon}^{a}$ satisfies

$d_{L}(t,(0,y_{0},z_{0}))=d_{L}(0,(0,y_{0},z_{0}))\mathrm{e}^{-\frac{m}{2}t}\frac{\sqrt{D(t)}}{\sqrt{D(0)}},$
(10)

where $\frac{1}{16}<|D(t)|<2$.

###### Proof.

From expression (34), we write the expression of
$\mathbf{u}(t;0,(0,y_{0},z_{0}))$ on the half-space $\\{x>0\\}$ as

$\begin{pmatrix}x\\\ y\\\ z\end{pmatrix}=\begin{pmatrix}k\varepsilon\\\
k^{2}\varepsilon-\varepsilon\\\ 0\end{pmatrix}+z\begin{pmatrix}1\\\ k\\\
1\end{pmatrix}+\vartheta\mathrm{e}^{\frac{k}{2}t}\begin{pmatrix}\cos(\theta_{1}+\xi_{k}t)\\\
\cos(\theta_{2}+\xi_{k}t)\\\ 0\end{pmatrix},$ (11)

where

$\displaystyle\vartheta^{2}$
$\displaystyle=\frac{4}{4-k^{2}}((y_{0}+\varepsilon)^{2}+(z_{0}+k\varepsilon)(z_{0}-ky_{0})),$
$\displaystyle\tan\theta_{1}$
$\displaystyle=\frac{k(z_{0}+k\varepsilon)-2(y_{0}+\varepsilon)}{\sqrt{4-k^{2}}(z_{0}+k\varepsilon)},$
$\displaystyle\tan\theta_{2}$
$\displaystyle=\frac{(z_{0}+k\varepsilon)(k^{2}-2)-k(y_{0}+\varepsilon)}{\sqrt{4-k^{2}}(k(z_{0}+k\varepsilon)-(y_{0}+\varepsilon))},$

and $\xi_{k}$ is given in appendix A.

From (11), the evolution of an orbit in the region $\\{x>0\\}$ can be obtained
as the evolution of a point $(\hat{x},\hat{y},z)$, with
$\hat{x}=k\varepsilon+z$ and $\hat{y}=k^{2}\varepsilon-\varepsilon+kz$, over
the repelling branch of the slow manifold $S_{\varepsilon}^{r}$, plus an
oscillatory term with increasing amplitude. Hence, the evolution of the
distance from the orbit to the slow manifold is computed as

$\sqrt{(x-\hat{x})^{2}+(y-\hat{y})^{2}}=\vartheta\mathrm{e}^{\frac{k}{2}t}\sqrt{\cos^{2}(\theta_{1}+\xi_{k}t)+\cos^{2}(\theta_{2}+\xi_{k}t)}.$

Taking $C(t)=\cos^{2}(\theta_{1}+\xi_{k}t)+\cos^{2}(\theta_{2}+\xi_{k}t)$ and
considering different trigonometric identities, we can rewrite $C(t)$ as

$C(t)=1+\frac{1}{2}\sqrt{\frac{k^{2}+3}{2}}\sin\left(\sqrt{4-k^{2}}t+\frac{\pi}{2}+\varphi\right),$

where

$\tan\varphi=\tan\frac{\theta_{1}+\theta_{2}}{2}.$

Hence, given that the sine function lies between $-1$ and $1$ and that
$0<k<2$, we obtain $\frac{1}{16}<|C(t)|<2$. Finally, we take the initial
distance $d_{R}(0,(0,y_{0},z_{0}))=\vartheta\sqrt{C(0)}$, which proves (9).
The expression of the distance from the solution
$\mathbf{u}(t;0,(0,y_{0},z_{0}))$ along the plane $\pi(z_{0},t)$ to the
attracting slow manifold $\mathcal{S}_{\varepsilon}^{a}$ given in (10) follows
in a similar way. This concludes the proof of the lemma. ∎

From (9)-(10), both distances $d_{L}(t,(0,y_{0},z_{0}))$ and
$d_{R}(t,(0,y_{0},z_{0}))$ depend both on time $t$ and on the initial
condition $\mathbf{p}=(0,y_{0},z_{0})$. As the initial condition gets closer
to $\mathbf{p}^{r}$ along a radial segment (see fig. 2), the initial distance
becomes smaller, and the orbit through $\mathbf{p}$ takes more time to exit
from a tubular neighbourhood of $\mathcal{S}_{\varepsilon}^{r}$. The magnitude
$\delta$ of the radius of this tubular neighbourhood is not relevant as long
as it remains invariant throughout the analysis, so without loss of generality
we consider $\delta=1$. Let $t^{*}$ be the exit time of the orbit starting at
$\mathbf{p}$. Then, the orbit leaves the neighbourhood of
$\mathcal{S}_{\varepsilon}^{r}$ at point $\mathbf{p}_{o}$ with $z$-coordinate
being $z_{o}=z_{0}+\varepsilon t^{*}$. In the particular case where
$\mathbf{p}=\mathbf{p}^{a}$, we write $z_{d}$ for $z_{o}$.

$\\{x=0\\}$$\mathbf{p}^{a}$$\mathbf{p}^{r}$$\mathbf{p}$$\mathbf{p}_{o}$$\mathbf{p}_{i}$$\mathcal{R}$$\mathcal{S}_{\varepsilon}^{a}$$\mathcal{S}_{\varepsilon}^{r}$

Figure 2: Tubular neighbourhoods of the slow manifolds. Magnification of fig.
1(b) and representation of the tubular neighbourhood of
$\mathcal{S}_{\varepsilon}^{a}$ (in blue) and that of
$\mathcal{S}_{\varepsilon}^{r}$ (in red). The orbit through the initial
condition $\mathbf{p}=(0,y_{0},z_{0})$ it is also represented in backward and
forward time, together with the entry point $\mathbf{p}_{i}$ and the exit one
$\mathbf{p}_{o}$ (dashed black curve through $\mathbf{p}$). Solid black line
named $\mathcal{R}$ represents a segment parallel to $S_{\varepsilon}^{a}$
lying on the tubular neighbourhood. Solid spiral lying on the plane
$\\{x=0\\}$ represents the image by the flow of the segment $\mathcal{R}$.

On the other hand, an inverse process takes place in the half-space
$\\{x<0\\}$, where orbits tend to the attracting branch of the slow manifold,
reaching the tubular neighbourhood of radius $\delta=1$ around
$\mathcal{S}_{\varepsilon}^{a}$ at some point denoted by $\mathbf{p}_{i}$ with
$z$-coordinate $z_{i}$. From this point, the orbit continues tending to
$\mathcal{S}_{\varepsilon}^{a}$ and reaching the plane $\\{x=0\\}$ at
$\mathbf{p}$ (see fig. 2). If $z_{i}$ is negative enough, then the orbit
approaches the switching plane $\\{x=0\\}$ very close to $\mathbf{p}^{a}$.

By construction, if we consider a segment $\mathcal{R}$ along the tubular
neighbourhood of $\mathcal{S}_{\varepsilon}^{a}$, and we continue the flow
generated by each point in $\mathcal{R}$, then we obtain a smooth manifold
that intersects the switching plane $\\{x=0\\}$ at a curve spiralling towards
$\mathbf{p}^{a}$ (see fig. 2). In fig. 3(a), we compute this curve by taking
different initial conditions on $\mathcal{R}$, which is at distance $\delta=1$
to $\mathcal{S}_{\varepsilon}^{a}$. We note that, the closer a point is on the
spiral to $\mathbf{p}^{a}$, the smaller $z_{i}$ is and, hence, the closer
$z_{o}$ is to $z_{d}$. In fig. 3(b), we represent the relation between
$|z_{i}|$ and $z_{o}$ obtained by considering different initial conditions on
$\mathcal{R}$. We call this relation the way-in/way-out function since its
behavior and properties are compatible with classical way-in/way-out functions
considered in smooth slow-fast systems, see e.g. [23]. Note that the way-
in/way-out function asymptotically tends to the constant value $z_{d}$, which
can be considered as the maximal delay. Since $z_{d}>0$, this implies that the
PWL system that we are considering exhibits a delayed loss of stability. Next,
we describe how the maximal delay $z_{d}$ behaves as $\varepsilon$ tends to
zero.

|
---|---
(a) | (b)
Figure 3: Graph of the way-in/way-out function. Panel (a) displays the spiral
formed by the intersection points with the plane $\\{x=0\\}$ of the orbits
starting in the segment
$\mathcal{R}=\\{(z-m\varepsilon,\varepsilon-m(z-m\varepsilon)+1,z):\ z<0\\}$
and represented in fig. 2. The endpoint of the spiral is given by the point
$\mathbf{p}^{a}$, plotted in blue. The red dot corresponds to the point
$\mathbf{p}^{r}$. In panel (b), we solve the way-in/way-out function for
points on the spiral, with initial conditions on the segment $\mathcal{R}$.
The blue dotted line corresponds to the maximal delay $z_{d}$. The values of
the remaining parameters are $k=0.1$, $m=1$, and $\varepsilon=0.25$.

######  Theorem 1.

Consider system (2)-(3). Then, the maximal delay $z_{d}$ satisfies

$-\frac{2\varepsilon}{k}\ln\left(\frac{2\sqrt{2}(m+k)}{\sqrt{4-k^{2}}}\varepsilon\right)+m\varepsilon<z_{d}<-\frac{2\varepsilon}{k}\ln\left(\frac{(m+k)}{2\sqrt{4-k^{2}}}\varepsilon\right)+m\varepsilon.$
(12)

In particular, $z_{d}$ tends to zero as $\varepsilon$ tends to zero.

###### Proof.

Consider $\gamma$ to be the orbit through
$\mathbf{p}^{a}=(0,-\varepsilon,m\varepsilon)$. To compute the escape point of
$\gamma$ from the $\mathcal{S}_{\varepsilon}^{r}$ tubular neighbourhood, we
equate equation (9) to $1$, which is the radius of the tubular neighbourhood.
Since $\frac{1}{16}<C(t)<2$, we obtain that
$\frac{1}{16}\vartheta^{2}\mathrm{e}^{kt_{d}}<1<2\vartheta^{2}\mathrm{e}^{kt_{d}}$.
Therefore, the time of flight from $\mathbf{p}^{a}$ to the escape point
satisfies

$-\frac{2}{k}\ln\left(\sqrt{2}\vartheta\right)<t_{d}<-\frac{2}{k}\ln\left(\frac{1}{4}\vartheta\right).$

Since the maximal delay is given by $z_{d}=\varepsilon t_{d}+m\varepsilon$,
expression (12) holds. Moreover, by the squeeze theorem, $z_{d}$ tends to zero
as $\varepsilon$ tends to zero. ∎

Therefore, lemma 1 implies that slow-fast PWL systems with two zones do not
allow to establish a connection between canonical attracting and repelling
slow manifolds, that is, between $\mathcal{S}_{\varepsilon}^{a}$ and
$\mathcal{S}_{\varepsilon}^{r}$. Hence, these systems present a delayed loss
of the stability with the maximal delay tending to zero with $\varepsilon$, as
we have proved in theorem 1. We next modify function $f$ in system (2) in
order to allow for the connection, which presents a behavior closer to what
happens in smooth slow-fast systems.

##  3 Three-regions system

In this section, we introduce a minimal PWL system also presenting a slow-
passage through a Hopf-like bifurcation, but for which the maximal delay does
not tend to zero with $\varepsilon$. For this purpose, we add a small
linearity zone to replace the point of non-normal hyperbolicity at the origin.
Hence, we add an extra segment in the critical manifold (3). This strategy of
considering a 3-piece critical manifold instead of a 2-piece one has proven
useful when studying fine slow-fast phenomena related to non-normal
hyperbolicity, in particular to approximate a smooth slow-fast system near a
quadratic fold of the critical manifold; see [10, 14] for more details on this
approach in the context of canard-explosive systems. Given that delayed Hopf
bifurcation is akin to a canard phenomenon (with different characteristics
though), this further motivates us to add this extra segment to the critical
manifold of system (2). This segment defines a new central linear region with
boundaries $\rho<0<\mu$, chosen such that the orbit following the attracting
branch $\mathcal{S}_{\varepsilon}^{a}$ at the switching plane $\\{x=\rho\\}$
approaches the switching plane $\\{x=\mu\\}$ at the same point that
$\mathcal{S}_{\varepsilon}^{r}$ does. Thus, the connection between both slow
manifolds happens, hence the slow-passage behavior.

Let us consider system (2) where the function $f(x)$ is now given by

$f(x)=\left\\{\begin{array}[]{cc}-mx&x<\rho,\\\ lx+n&\rho\leq x<\mu,\\\
kx&x\geq\mu,\end{array}\right.$ (13)

with parameters $l$ and $n$ fixed as

$l=\frac{m\rho+k\mu}{\mu-\rho},\quad n=-\rho\mu\frac{k+m}{\mu-\rho},$ (14)

to ensure the continuity of the vector field. The slope $l$ is the weighted
mean of the slopes $-m$ and $k$ with weights $\frac{-\rho}{\mu-\rho}$ and
$\frac{\mu}{\mu-\rho}$, respectively, and then $-2<-m<l<k<2$. Notice that,
since parameters $\rho<0<\mu$ organize the boundaries between the different
linear regions, in the particular case where these regions are centered, that
is when $\mu=-\rho$, then $l$ is the arithmetic mean of the slopes $-m$ and
$k$.

As solutions cross the plane $\\{x=\rho\\}$, their local expression changes
from the left-side expression (33) to the central one (35). Similarly, as
solutions cross the plane $\\{x=\mu\\}$, their local expression changes from
the central to the right-side expression (34).

As in the two regions case, slow manifolds perturb from the critical manifold
when $\varepsilon>0$ is small enough, see lemma 1. Then, the segments
$\\{(x,-mx-\varepsilon,x+m\varepsilon),\,x<\rho\\}$ and
$\\{(x,kx-\varepsilon,x-k\varepsilon),\,x>\mu\\}$ are part of the attracting
and of the repelling canonical slow manifolds, $\mathcal{S}_{\varepsilon}^{a}$
and $\mathcal{S}_{\varepsilon}^{r}$, respectively. These segments intersect
with the planes $\\{x=\rho\\}$ and $\\{x=\mu\\}$ at points
$\tilde{\mathbf{p}}^{a}=(\rho,-m\rho-\varepsilon,\rho+m\varepsilon)$ and
$\tilde{\mathbf{p}}^{r}=(\mu,k\mu-\varepsilon,\mu-k\varepsilon)$,
respectively. On the other hand, the segment
$\\{(x,xl+n-\varepsilon,x-\varepsilon l),\,\rho<x<\mu\\}$ is invariant under
the flow of the system defined in the central region, and it is either part of
the attracting ($l<0$) or the repelling ($l>0$) slow manifold, respectively.

To be able to connect both lateral slow manifolds through the central region,
we impose that the solution in the central region (35) passing through
$\tilde{\mathbf{p}}^{a}$ reaches the boundary $\\{x=\mu\\}$ at
$\tilde{\mathbf{p}}^{r}$. This yields

$\mathbf{u}_{C}(\hat{t};0,\tilde{\mathbf{p}}^{a})=\tilde{\mathbf{p}}^{r}.$
(15)

This boundary condition provides three different equations, corresponding to
the three coordinates of the solution in the central region (35), with six
unknowns, $(\hat{t},m,k,\rho,\mu,\varepsilon)$. Next, we solve this system of
equations by obtaining the first three unknowns in terms of the others.

######  Lemma 3.

Consider system (2) with $f$ given by (13) and $\rho<0<\mu$. Assuming

$m=\frac{\mu-\rho}{\varepsilon}-k-\sqrt{\ln^{2}\left\lvert\frac{\rho}{\mu}\right\rvert+\pi^{2}}$
(16)

and

$k=\frac{-1}{\mu-\rho}\left(\frac{2\ln\left\lvert\frac{\rho}{\mu}\right\rvert(\rho-\mu)}{\sqrt{\ln^{2}\left\lvert\frac{\rho}{\mu}\right\rvert+\pi^{2}}}+\frac{\rho}{\varepsilon}\big{(}\mu-\rho-\varepsilon\sqrt{\ln^{2}\left\lvert\frac{\rho}{\mu}\right\rvert+\pi^{2}}\big{)}\right),$
(17)

then equation (15) is satisfied with the time of flight of the solution in the
central region being equal to

$\hat{t}=\frac{\mu-\rho}{\varepsilon}-(k+m).$ (18)

###### Proof.

The time of flight (18) from $x=\rho$ to $x=\mu$ can be obtained from the
third equation in (15) (corresponding to the $z$ component), where
$z(\hat{t})=m\varepsilon+\rho+\varepsilon\hat{t}=-k\varepsilon+\mu$. Moreover,
the first and second equations in (15) yield a linear system in terms of
$\mathrm{e}^{\frac{l}{2}\hat{t}}\cos(\xi_{l}\hat{t})$ and
$\mathrm{e}^{\frac{l}{2}\hat{t}}\sin(\xi_{l}\hat{t})$, whose solution is given
by

$\begin{pmatrix}\mathrm{e}^{\frac{l}{2}\hat{t}}\cos(\xi_{l}\hat{t})\\\
\mathrm{e}^{\frac{l}{2}\hat{t}}\sin(\xi_{l}\hat{t})\end{pmatrix}=\frac{\rho}{\mu}\begin{pmatrix}1\\\
0\end{pmatrix}.$

These expressions are equivalent to

$l=\frac{2}{\hat{t}}\ln\left\lvert\frac{\rho}{\mu}\right\rvert\quad\text{and}\quad\frac{\sqrt{4-l^{2}}}{2}\hat{t}=\pi+2\pi
q,\text{ with}\quad q\in\mathbb{Z}.$ (19)

Considering $q=0$ and merging the first equation in (19) with the second one,
we obtain that the latter can be written as

$\hat{t}^{2}-\ln^{2}\left\lvert\frac{\rho}{\mu}\right\rvert=\pi^{2}.$

Then, replacing $\hat{t}$ in the last expression by its value given by (18),
and isolating $m$ in the resulting expression, we obtain equation (16).

Similarly, equation (17) is obtained by equating the first expressions in (14)
and in (19), replacing therein the value of $m$ obtained in equation (16), and
then isolating $k$.

To ensure that the local solution $\mathbf{u}_{C}(t;0,\tilde{\mathbf{p}}^{a})$
is contained in the central region, that is $\rho\leq
x_{C}(t;0,\tilde{\mathbf{p}}^{a})\leq\mu$, for all $t\in[0,\hat{t}]$, we study
the sign of the derivative

$x_{C}^{\prime}(t;0,\tilde{\mathbf{p}}^{a})=\varepsilon\Big{(}1+\frac{m+l}{\xi_{l}}\sin(\xi_{l}t)e^{\frac{l}{2}t}\Big{)}.$

Then, since $m+l>0$, this function is positive for all $\sin(\xi_{l}t)\geq 0$,
that is for all $t\in[0,\frac{\pi}{\xi_{l}}]=[0,\hat{t}]$. Thus,
$x_{C}(t;0,\tilde{\mathbf{p}}^{a})$ is a monotone function in this interval
with $x_{C}(0;0,{\tilde{\mathbf{p}}}^{a})=\rho$ and
$x_{C}(\hat{t};0,{\tilde{\mathbf{p}}}^{a})=\mu$. ∎

######  Remark 1.

In some contexts, it is usual to set the slope in each region rather than
their boundaries, as it is done in lemma 3. In order to set the slopes, we
need to study the Jacobian matrix of the function
$\left(m(\rho,\mu),k(\rho,\mu)\right)$. Then, by the Inverse Function Theorem,
we can invert the previous functions if the expression

$\frac{2\pi^{2}}{\frac{-\rho}{\mu}A}+\frac{1}{\varepsilon(\frac{-\rho}{\mu}+1)^{2}}\Big{(}(\mu-\rho)\sqrt{A}-\varepsilon
A\Big{)},$

with $A=ln^{2}\left\lvert\frac{\rho}{\mu}\right\rvert+\pi^{2}$, is different
from $0$. However, all its terms are positive except the last one. Hence, it
may happen that for some values we cannot guarantee the existence of the
inverse function. In the particular case where parameters satisfy

$(\mu-\rho)\sqrt{A}>\varepsilon A,$

we can guarantee the existence of the inverse in a neighbourhood of
$(\rho,\mu)$. In fact, if the system is centered, $\mu=-\rho$, the previous
inequality becomes $\frac{2\mu}{\pi}>\varepsilon$, which amounts to ensure
that $m$ and $k$ are positives, see (16)-(17).

lemma 3 provides conditions on the parameters $m$ and $k$ to ensure the
connection between the attracting and repelling slow manifolds. These
conditions are sufficient, but not necessary. Indeed, for some parameter
choices, we can obtain some loops in the central region by considering other
values for $q$ in (19) (see fig. 4(a)), or cross the boundary $x=\mu$ at least
one time before the connection is made (see fig. 4(b)). These situations
provide different ways to connect the slow manifolds. Nevertheless, the full
study of how to make the connection goes beyond the scope of the present work,
given that our aim is to find a minimal system exhibiting any simple
connection between the two slow manifolds.

|
---|---
(a) | (b)
Figure 4: Connecting the slow manifolds. Other possible connections that can
appear between the attracting (blue line) and the repelling (red line) slow
manifolds. The connection can be given; (a) with a loop in the central region;
(b) crossing $x=\mu$ at least one time before connecting with the repelling
slow manifold.

######  Theorem 2.

Consider system (2) with $f$ given by (13) and let $\rho<0<\mu$. Let $m$ and
$k$ satisfy (16) and (17), respectively. Then, the maximal delay tends to
infinity as $\varepsilon$ tends to zero. Moreover, the way-in/way-out function
is asymptotic to the line
$z_{\mathrm{out}}=\frac{m}{k}z_{\mathrm{in}}-\frac{2\varepsilon}{k}\ln\left\lvert\frac{\rho}{\mu}\right\rvert$.

###### Proof.

The assertion about the behavior of the maximal delay as $\varepsilon$ tends
to zero follows immediately from lemma 3, since under these conditions both
slow manifolds $\mathcal{S}_{\varepsilon}^{a,r}$ connect, which implies that
the orbit along the attracting slow manifold continues along the repelling
slow manifold and the delay can be considered infinite.

On the other hand, the slope of the way-in/way-out function can be
approximated as follows. Consider a tubular neighbourhood of radius $\delta$
around the attracting slow manifold, and a similar neighbourhood around the
repelling one; without loss of generality, we can take $\delta=1$ on both
sides. Let $\mathbf{p}_{i}=(x_{i},y_{i},z_{i})$ be an initial condition (also
called entry point) located on the left tubular neighbourhood around the
attracting slow manifold, that is, at a distance $\delta$ to
$\mathcal{S}_{\varepsilon}^{a}$, see (10). Let
$\mathbf{p}_{\rho}=(x_{\rho},y_{\rho},z_{\rho})$ be the point where the orbit
through $\mathbf{p}_{i}$ crosses the boundary $\\{x=\rho\\}$ and
$\mathbf{p}_{\mu}=(x_{\mu},y_{\mu},z_{\mu})$ the one crossing the
$\\{x=\mu\\}$ boundary. Finally, consider $\mathbf{p}_{o}=(x_{o},y_{o},z_{o})$
be the point of the orbit located on the right tubular neighbourhood, that is,
around the repelling slow manifold (also called exit point). The slope of the
way-in/way-out function is given by the ratio of the $z$-component of the
entry and exit points, respectively, that is,
$(z_{o}-z_{\mu})/(z_{\rho}-z_{i})$. By definition, $z_{i}=\varepsilon
t_{i}+z_{\rho}$ and $z_{o}=\varepsilon t_{o}+z_{\mu}$, where $-t_{i}>0$ is the
time necessary to reach $\mathbf{p}_{\rho}$ from the initial condition
$\mathbf{p}_{i}$ and $t_{o}$ is the time required to reach $\mathbf{p}_{o}$
starting at $\mathbf{p}_{\mu}$. Hence,

$z_{o}-z_{\mu}=\frac{t_{o}}{-t_{i}}(z_{\rho}-z_{i}).$ (20)

On the other hand, the contraction along the attracting slow manifold at time
$t_{i}$ yields $\delta=d_{L}(t_{i})=d_{L}(0)\mathrm{e}^{-\frac{m}{2}t_{i}}$,
see equation (10), while the expansion along the repelling slow manifold at
time $t_{o}$ is given by
$\delta=d_{R}(t_{o})=d_{R}(0)\mathrm{e}^{\frac{k}{2}t_{o}}$, see equation (9).
Hence, combining both expressions, it follows
$d_{L}(0)\mathrm{e}^{-\frac{m}{2}t_{i}}=d_{R}(0)\mathrm{e}^{\frac{k}{2}t_{o}}$.
On the other hand, since $d_{R}(0)\approx d_{L}(0)e^{\frac{l}{2}\hat{t}}$,
where $\frac{l\hat{t}}{2}$ can be obtained from (19), it follows
$d_{R}(0)\approx d_{L}(0)\left\lvert\frac{\rho}{\mu}\right\rvert$. Therefore,
we obtain that
$\frac{t_{o}}{-t_{i}}\approx\frac{m}{k}+\frac{2}{kt_{i}}\ln\left\lvert\frac{\rho}{\mu}\right\rvert$
and so, from $\eqref{eq.z1_z2}$ and $-\varepsilon t_{i}=z_{\rho}-z_{i}$ we
have that

$z_{o}-z_{\mu}\approx\frac{m}{k}(z_{\rho}-z_{i})-\frac{2\varepsilon}{k}\ln\left\lvert\frac{\rho}{\mu}\right\rvert,$

which proves the second statement of the theorem. ∎

Given $\mu$ and $\rho$, and under assumptions of Lemmalemma 3, if we perturb
the boundaries, then the slopes $m$, $k$ and $l$ will change. The relative
position of the boundaries imposes both the sign of the slope $l$ of the
critical manifold in the central region, as well as the relative size of the
slopes in the other two regions, as shown in the following lemma.

######  Lemma 4.

Consider the system (2) with $f$ given by (13). Suppose that equations (16)
and (17) are satisfied. Then, the following statements hold:

  * (a)

if $\mu<-\rho$, then $l>0$ and $k>m$,

  * (b)

if $\mu>-\rho$, then $l<0$ and $k<m$,

  * (c)

if $\mu=-\rho$, then $l=0$ and $k=m$.

###### Proof.

Since the first equation in (19) holds, if $\mu<-\rho$, then $l$ is positive.
Consequently, from the first equation in (14), we have that $\mu k+\rho m>0$
and thus

$k>\frac{-\rho}{\mu}m>m.$

Following the same procedure, if we consider $\mu>-\rho$, then $l<0$, $\mu
k+\rho m<0$ and

$k<\frac{-\rho}{\mu}m<m.$

Finally, in the specific case that $\mu=-\rho$, from equation (19) we obtain
that $l=0$ and thus $\mu k+\rho m=0$. Hence, $k=m$. ∎

Take the minimal PWL slow-fast system (2) with $f$ given by (13), and consider
the fast subsystem at $\varepsilon=0$. Then, following [15], the existence of
a Hopf-like bifurcation (supercritical or subcritical) in the fast subsystem
(5) can be obtained according to a magnitude ratio between the slopes of the
$x$-nullcline. Our next result provides conditions to guarantee the existence,
location and criticality of this bifurcation.

######  Lemma 5.

Consider system (5) with $f$ given by (13). Then, the following assertions
hold:

  * (a)

If $l>0$, a Hopf-like bifurcation occurs at $z=\rho$; it is a subcritical if
$l-m>0$ and supercritical otherwise.

  * (b)

If $l<0$, a Hopf-like bifurcation occurs at $z=\mu$; it is a subcritical if
$l+k>0$ and supercritical otherwise.

###### Proof.

Consider the 2-regions system formed by the central and left regions of
(5)-(13). Then, the change of variables given by $x=\hat{x}+\rho$ and
$y=\hat{y}-m\rho$ allows us to rewrite the system in the following form which
is analyzed in [15]

$\begin{array}[]{l}\dot{\hat{x}}=\hat{f}(\hat{x})-\hat{y},\\\
\dot{\hat{y}}=\hat{x}+a,\end{array}$ (21)

with, in the present case, $a=\rho-z$, $\hat{f}(\hat{x})=-m\hat{x}$ if
$\hat{x}<0$ and $\hat{f}(\hat{x})=l\hat{x}$ otherwise. Let $\text{det}_{L}=1$
and $t_{L}=-m$ be the determinant and the trace of the matrix defining the
linear subsystem (21) in $x<\rho$, and $\text{det}_{C}=1$ and $t_{C}=l$ be the
determinant and the trace of the matrix in $x>\rho$. Consider also the
weighted traces

$\gamma_{L}=\frac{t_{L}}{\sqrt{4\text{det}_{L}-t_{L}^{2}}}=\frac{-m}{\sqrt{4-m^{2}}}\quad\text{and}\quad\gamma_{R}=\frac{t_{C}}{\sqrt{4\text{det}_{C}-t_{C}^{2}}}=\frac{l}{\sqrt{4-l^{2}}}.$

Since $-2<-m<l<2$, it follows that

$D_{L}=\text{det}_{L}-\frac{t_{L}^{2}}{4}=\frac{4-m^{2}}{4}>0,\quad\text{and}\quad
D_{R}=\text{det}_{C}-\frac{t_{C}^{2}}{4}=\frac{4-l^{2}}{4}>0.$

In the conditions of statement (a), it follows that $t_{L}\,t_{C}=-ml<0$ and
$\text{sign}(\gamma_{L})+\text{sign}(\gamma_{R})=\text{sign}(l-m)$. According
to the sign of $l-m$, we obtain different situations following Theorem 5 in
[15], which are explained below.

If $l-m>0$, for $z<\rho$ we have $a>0$ and so we have an asymptotically stable
equilibrium which is surrounded by a unique unstable limit cycle (case (a)).
Otherwise, for $z>\rho$ we have $a<0$ and so the equilibrium becomes unstable
and no limit cycles exist (case (c)). Consequently, a supercritical Hopf-like
bifurcation appears at $z=\rho$.

If $l-m<0$, for $z<\rho$ we have $a>0$ and so we have an asymptotically stable
equilibrium and the system do not present limit cycles (case (b)). Otherwise,
for $z>\rho$ we have $a<0$ and so the equilibrium becomes unstable and it is
surrounded by a unique stable limit cycle (case (d)). Hence, a subcritical
Hopf-like bifurcation appears at $z=\rho$.

Otherwise, if $l=m$, notice that in this case no criterion is established
since Theorem 5 only applies when $\text{sign}(a)=\text{sign}(\gamma_{C})$. In
such a case, the continuity of the vector field implies that the linear
systems in the regions $L$ and $C$ are identical and the minimal PWL model is
equivalent to the one analyzed in section 2, but here with the boundary
translated to $x=\mu$. This finishes the proof of statement (a).

If we consider the 2-regions system given by the central and right regions and
proceeding in a similar way that in the proof of the statement (a), the proof
of statement (b) follows. ∎

Let us consider the specific case $\mu=-\rho$, where the central region is
centered at the origin. Under this new assumption, lemma 4 implies that $l=0$
and $m=k$. Hence, from equations (16)-(17), we obtain that
$k=m=\frac{\mu}{\varepsilon}-\frac{\pi}{2}$. In particular, to ensure that
$0<k=m<2$, we need to consider that
$\frac{\pi}{2}\varepsilon<\mu<\frac{\pi+4}{2}\varepsilon$, which forces the
width of the central region to be of order $\varepsilon$. Notice that, for
$l=0$, lemma 5 does not ensure the existence of the Hopf-like bifurcation.
Nonetheless, given $\varepsilon$ and $\mu=-\rho$ lying in the specified
interval, there is a slight perturbation of $\mu$ and $\rho$ such that $m$ and
$k$ remain in the interval $(0,2)$. Hence, by lemma 4, if $\mu<-\rho$, then
$k>m$, and if $\mu>-\rho$, then $k<m$. Therefore, by lemma 5, a Hopf-like
bifurcation appears.

This fact is illustrated in fig. 5. Panel (a) shows the case $\mu=-\rho=0.15$
and $\varepsilon=0.05$, where no Hopf-like bifurcation exists, as can be
concluded from the exponential growth of the amplitude of the oscillations.
However, if we slightly perturb $\rho$, taking for instance $\rho=-0.085$ (see
Panel (b)), we obtain $k<m$ and a supercritical Hopf-like bifurcation appears.

|
---|---
(a) | (b)
Figure 5: Slow passage phenomenon through a Hopf-like bifurcation, with
initial conditions near the attracting slow manifold (blue line). The
repelling slow manifold is depicted in red. The parameters are $\mu=0.15$,
$\varepsilon=0.05$ and (a) $\rho=-0.15$, in this case, no Hopf bifurcation
occurs; (b) $\rho=-0.085$, in this case a supercritical Hopf-like bifurcation
takes place defining an attracting manifold, in blue in the figure. The values
of $m$ and $k$ are given by expressions (16) and (17), respectively.

In fig. 6 we plot the way-in/way-out function for $\mu=0.15$ and
$\rho=-0.085$; here, $k$ and $m$ are fixed so that the relations in lemma 3
are satisfied and guarantee the existence of a connection between the lateral
slow manifolds. In the figure, we can see that the $z$ coordinate of the
escape point from the tubular neighbourhood of the repelling slow manifold is
proportional to the $z$ coordinate of the entry point to the tubular
neighbourhood of the attracting slow manifold. Due to the connection between
both slow manifolds, this relation is expected to persist along the whole
domain of the function, as shown in theorem 2. Nevertheless, in the
simulations, this relation is not satisfied at all and the way-in/way-out
function remains constant to a specific value. This unexpected behavior is due
to the numerical precision in the simulations, which does not provide an exact
connection. Therefore the plateau do not correspond to a buffer point (see
fig. 6).

Figure 6: Way-in/way-out function for a delayed supercritical Hopf-like
bifurcation. Plot of the way-in/way-out function in the case $\rho=-0.085$,
$\mu=0.15$, $\varepsilon=0.05$, for the parameters $k\approx 0.189$, $m\approx
1.318$ and $l\approx-0.356$ computed according to equations (16), (17) and
(14). The way-in/way-out function has been calculated for different numerical
precision. The dashed curve corresponds to the asymptotic line
$z_{\mathrm{out}}=\frac{m}{k}z_{\mathrm{in}}-\frac{2\varepsilon}{k}\ln{\left\lvert\frac{\rho}{\mu}\right\rvert}$
according to theorem 2.

To support this claim, we study the relationship between computational
accuracy and the asymptotic value of the graph at which it saturates. table 1
relates the error in the connection, $\vartheta$, needed to obtain the
observed saturation point in the way-in/way-out function, with the
computational error of the numerical integrator. Using lemma 2 and equation
(9), we can bound $\vartheta$ by the extreme values of $C(t)$, that is,

$\sqrt{\frac{1}{2\mathrm{e}^{kt}}}<\vartheta<\sqrt{\frac{16}{\mathrm{e}^{kt}}}.$

Hence, using the flying time from $x=\mu$ to the escape time, we obtain the
minimum and maximum values of $\vartheta$ given in table 1. As we can see, the
computational precision has approximately the same order as $\vartheta$.
Therefore, the saturation value that can be seen in fig. 6 agrees with the
computational precision.

Computational Precision | min $\vartheta$ | max $\vartheta$
---|---|---
$10^{-12}$ | $2.78\cdot 10^{-12}$ | $1.57\cdot 10^{-11}$
$10^{-9}$ | $1.43\cdot 10^{-10}$ | $8.09\cdot 10^{-10}$
$10^{-6}$ | $1.24\cdot 10^{-6}$ | $7.03\cdot 10^{-6}$
Table 1: Comparison between round-off error and error at the connection.
Comparison between the precision of the numerical integrator and the value of
$\vartheta$, which provides an approximate value of the initial distance
between the attracting and the repelling slow manifolds on the plane
$\\{x=\mu\\}$. To compute the second and the third column, we use the way-
in/way-out function evaluated on the attracting slow manifold at $x=\rho$,
with a tubular neighbourhood of radius $1$ and the same parameters as in fig.
6(a).

###  3.1 Adding a buffer point

As described in theorem 2, the connection between the slow manifolds yields a
maximal delay that is not finite, which implies that the way-in/way-out
function is unbounded. However, the presence of an equilibrium point on the
repelling slow manifold can modify this behavior. In fact, such an equilibrium
blocks the flow not only on the repelling slow manifold that connects to it,
but also for orbits sufficiently close to this slow manifold. Hence, it forces
the $z$-coordinate of the exit point of each such orbits to remain constant
(see [11, 26]). This constant value is called a buffer point.

In order to analyze the slow passage through a Hopf-like bifurcation in the
presence of a buffer point, we extend the previous PWL slow-fast system and
now consider

$\begin{array}[]{ccl}\dot{x}&=&f(x)-y,\\\ \dot{y}&=&x-z,\\\
\dot{z}&=&\varepsilon(a-x),\end{array}$ (22)

where $f(x)$ is still taken as in (13).

We start by establishing conditions on the parameters in order to guarantee a
connection between the attracting and the repelling slow manifolds. As the
expressions of the eigenvalues become more complicated, we now assume that the
matrix of the system in each region has its spectrum formed by one real
eigenvalue $\lambda_{j}$ and a pair of complex conjugated eigenvalues
$\alpha_{j}\pm\beta_{j}\,\rm{i}$, with $j\in\\{L,C,R\\}$; the subscripts $L,C$
and $R$ stand for the left half-space $\\{x<\rho\\}$, the central strip
$\\{\rho<x<\mu\\}$ and the right half-space $\\{x>\mu\\}$, respectively. To
simplify the notations, when $\lambda_{j}\alpha_{j}<0$ we refer to this
configuration as of saddle-focus type, when $\lambda_{j}\alpha_{j}>0$ we refer
to it as of node-focus type and when $\alpha_{j}=0$ we refer to it as of
saddle-center type. In particular, we consider a configuration of saddle-
center type in the central region, that is, $\alpha_{C}=0$. Therefore, the
complex eigenvalues of the matrix of the linear system defined in the central
strip

$\mathbf{\dot{u}}=\begin{pmatrix}l&-1&0\\\ 1&0&-1\\\
-\varepsilon&0&0\end{pmatrix}\mathbf{u}+\begin{pmatrix}n\\\ 0\\\ \varepsilon
a\end{pmatrix},$ (23)

are purely imaginary. We also assume the size of this central region to be
order $O(\varepsilon)$. Following [31], the attracting and repelling slow
manifolds are given by the line segments each formed by the equilibrium point
and eigenvector associated with the slow eigenvalue, contained respectively in
L or R, i.e.

$\begin{array}[]{l}\mathcal{S}_{\varepsilon}^{a}=\Big{\\{}(a,-ma,a)+r_{L}(1,-m-\lambda_{L},-\varepsilon/\lambda_{L});\;r_{L}\leq\rho-a\Big{\\}},\\\
\mathcal{S}_{\varepsilon}^{r}=\Big{\\{}(a,ka,a)+r_{R}(1,k-\lambda_{R},-\varepsilon/\lambda_{R});\;r_{R}\geq\mu-a\Big{\\}},\end{array}$
(24)

respectively. In the next lemma, we establish conditions on the parameters to
guarantee the connection between both slow manifolds. These conditions are
given by power series of $\varepsilon$ obtained using the method of the
undetermined coefficients.

######  Lemma 6.

Consider system (22), with $f$ given by (13), and $\varepsilon>0$ small
enough. Assume that the behavior in the left half space is of saddle-focus
type with $\lambda_{L}\alpha_{L}<0$, in the right half space is of node-focus
type with $\lambda_{R}\alpha_{R}>0$ and the behavior in the central strip is
of saddle-center type, i.e. $\alpha_{C}=0$. Setting
$\rho=-k\varepsilon-\varepsilon^{2}$ and $\mu=m\varepsilon-\varepsilon^{2}$,
there exist positive functions $m(a,\varepsilon)$ and $k(a,\varepsilon)$ such
that, when $m=m(a,\varepsilon)$ and $k=k(a,\varepsilon)$, the slow manifolds
$\mathcal{S}_{\varepsilon}^{a}$ and $\mathcal{S}_{\varepsilon}^{r}$ connect.
Moreover, $m(a,\varepsilon)$ and $k(a,\varepsilon)$ can be expressed in series
expansion of $\varepsilon$ as

$\displaystyle m(a,\varepsilon)$
$\displaystyle=\frac{a\pi}{2(1-a)}+O(\varepsilon),$ $\displaystyle
k(a,\varepsilon)$ $\displaystyle=\frac{a\pi}{2(1-a)}+O(\varepsilon).$

###### Proof.

We first ensure that all the conditions stated in the hypotheses are
compatible, in particular, with the continuity of the system, i.e., that
relation (14) is satisfied.

The characteristic polynomial of the system (23) is
$p_{C}(\lambda)=-\lambda^{3}+l\lambda^{2}-\lambda-\varepsilon$. Since we
assume a configuration of saddle-center type, i.e.
$p_{C}(\lambda)=-(\lambda-\lambda_{C})(\lambda^{2}+\beta_{C}^{2}),$ by
equalling coefficients, we get $l=\lambda_{C}=-\varepsilon$ and $\beta_{C}=1$.
Hence equation (14) writes

$-\varepsilon=\frac{m\rho+k\mu}{\mu-\rho},$ (25)

and from the expressions of $\mu$ and $\rho$ given in the statement of the
lemma, we conclude that equation (14) is satisfied.

Proceeding similarly with the lateral systems, we can relate the eigenvalues
in terms of the parameters by equalling the coefficients of the characteristic
polynomial. From this, we get the following sets of relations

$\left\\{\begin{array}[]{ccl}-m&=&2\alpha_{L}+\lambda_{L},\\\
1&=&2\alpha_{L}\lambda_{L}+(\alpha_{L}^{2}+\beta_{L}^{2}),\\\
-\varepsilon&=&\lambda_{L}(\alpha_{L}^{2}+\beta_{L}^{2}),\end{array}\right.\quad\text{and}\quad\left\\{\begin{array}[]{ccl}k&=&2\alpha_{R}+\lambda_{R},\\\
1&=&2\alpha_{R}\lambda_{R}+(\alpha_{R}^{2}+\beta_{R}^{2}),\\\
-\varepsilon&=&\lambda_{R}(\alpha_{R}^{2}+\beta_{R}^{2}).\end{array}\right.$

We note that these conditions are compatible with $m>0$, $\alpha_{L}<0$,
$\lambda_{L}<0$, and $k>0$, $\alpha_{R}>0$ and $\lambda_{R}<0$.

From (24), the intersection points $\mathbf{p}_{\varepsilon}^{a}$ and
$\mathbf{p}_{\varepsilon}^{r}$, of the attracting and repelling slow manifold
with the planes $\\{x=\rho\\}$ and $\\{x=\mu\\}$ are

$\begin{array}[]{l}\mathbf{p}_{\varepsilon}^{a}=\left(\rho,-ma-(\rho-a)(m+\lambda_{L}),a-(\rho-a)\frac{\varepsilon}{\lambda_{L}}\right)\\\
\\\
\mathbf{p}_{\varepsilon}^{r}=\left(\mu,ka+(\mu-a)(k-\lambda_{R}),a-(\mu-a)\frac{\varepsilon}{\lambda_{R}}\right),\\\
\end{array}$ (26)

respectively. Therefore, the condition for the connection between both slow
manifolds can be written as
$\mathbf{u}_{C}(\tau;0,\mathbf{p}_{\varepsilon}^{a})=\mathbf{p}_{\varepsilon}^{r}$,
where $\mathbf{u}_{C}$ is given in (36). Expressing the solution
$\mathbf{u}_{C}$ in terms of the coordinates, the connection condition becomes
in a system of three equations and five unknowns $(a,m,k,\tau,\varepsilon)$,

$\begin{array}[]{rcl}x_{C}(\tau;0,\mathbf{p}_{\varepsilon}^{a})&=&\mu,\\\
y_{C}(\tau;0,\mathbf{p}_{\varepsilon}^{a})&=&ka+(\mu-a)(k-\lambda_{R}),\\\
z_{C}(\tau;0,\mathbf{p}_{\varepsilon}^{a})&=&a-(\mu-a)\frac{\varepsilon}{\lambda_{R}},\end{array}$
(27)

which is equivalent to the system obtained by dividing by $\varepsilon$ every
equation,

$\begin{array}[]{rcl}\frac{1}{\varepsilon}(x_{C}(\tau;0,\mathbf{p}_{\varepsilon}^{a})-\mu)&=&0,\\\
\frac{1}{\varepsilon}(y_{C}(\tau;0,\mathbf{p}_{\varepsilon}^{a})-ka-(\mu-a)(k-\lambda_{R}))&=&0,\\\
\frac{1}{\varepsilon}(z_{C}(\tau;0,\mathbf{p}_{\varepsilon}^{a})-a+(\mu-a)\frac{\varepsilon}{\lambda_{R}})&=&0.\end{array}$
(28)

Following the expression of $x_{C},y_{C}$ and $z_{C}$ in (36) it can be
concluded that for every $a\not\in\\{0,1\\}$, the vector
$\mathbf{s}=(a,\frac{a\pi}{2(1-a)},\frac{a\pi}{2(1-a)},\pi,0)$ is a solution
of system (28). Moreover, computing the derivative of the system with respect
to the variables $m,k,$ and $\tau$ evaluated at the previous solution we
obtain that

$D_{m,k,\tau}(\mathbf{s})=\begin{pmatrix}2a-1&1&a\\\
0&0&a\frac{a\pi}{2(a-1)}\\\ a-1&a-1&a\end{pmatrix},$

having determinant $\det(D_{m,k,\tau}(\mathbf{s}))=-a^{2}\pi(a-1)\neq 0$. From
the Implicit Function Theorem we conclude the existence of functions
$m(a,\varepsilon)$, $k(a,\varepsilon)$ and $\tau(a,\varepsilon)$, analytic as
functions of $\varepsilon$, and such that the parameters
$(a,m(a,\varepsilon),k(a,\varepsilon),\tau(a,\varepsilon),\varepsilon)$ are
solutions of the system (28). ∎

######  Theorem 3.

Under the conditions to allow the connection of the slow manifolds given in
lemma 6, the maximal delay for system (22) is attained at the buffer point,
i.e. $z_{d}=a-\mu$.

###### Proof.

Under the conditions of lemma 6 both slow manifolds connect and the
equilibrium point at $(a,f(a),a)$ is the limit of the solution on the
repelling side. Therefore, the $z$ coordinate of the exit point of any orbit
coincides with that of the buffer point and then, the distance from the right
switching line to the exit point is $a-\mu$. ∎

theorem 3 describes the behavior of the delay in the presence of a buffer
point. In fig. 7, we show this behavior through the plot of the way-in/way-out
function. This plot has a theoretical maximum value at $a-\mu$. However, for
the simulations, we consider a truncated expression of the series appearing in
lemma 6, and consequently the slow manifolds do not exactly connect and the
coordinate of the exit point from the tubular neighbourhood tends to a value
lower than the theoretical one. The curve in fig. 7 was obtained by computing
4 terms in the series expansion of $m(a,\varepsilon)$ and $k(a,\varepsilon)$.
If we add new terms in the series expansion then we will get the asymptotic
value to be closer to this theoretical value. This is in contrast to fig. 6,
where the way-in/way-out function grows indefinitely as we improve the
numerical precision.

Figure 7: Way-in/way-out function in the presence of a buffer point, defined
by considering a tubular neighbourhood of radius $1$ for $a=0.2$,
$\varepsilon=0.05$, $k\approx 0.3915$ and $m\approx 0.3990$, which are given
by the first four terms of the series expansion in lemma 6. The green line
corresponds to the distance from the buffer point to the boundary at $x=\mu$,
i.e. $a-\mu$; the dashed red line corresponds to the distance from the exit
coordinate of the repelling slow manifold to the boundary; the dashed blue
line corresponds to distance from the exit coordinate of the attracting slow
manifold to the boundary.

Another illustration of the lack of precision in the connection between the
slow manifolds is provided in fig. 8. Under the connection conditions given by
lemma 6 and satisfied only up to the fourth term, the orbit through the
attracting slow manifold leaves the repelling slow manifold instead of
approaching the equilibrium.

As mentioned above, this lack of precision is the cause of the discrepancy
between theorem 3 and the results presented in fig. 7. What is more, this
effect is amplified because of the stiffness of the problem.

Figure 8: Solution of the system with a buffer point. Plot of the solution (in
black), the attracting (blue) and repelling (red) slow manifolds and the
buffer point (purple dot), with $a=0.2$, $\varepsilon=0.05$ and the other
parameters $k\approx 0.3915$ and $m\approx 0.3990$ are fixed for the
connection. We take an initial condition on the attracting slow manifold.

##  4 Application to the Doi-Kumagai neural burster

The Doi-Kumagai (DK) model [12] is a very simplified neural burster obtained
by the piecewise linearization of a 3-dimensional extension of the Bonhoefer-
van der Pol oscillator; see [19] for details. Its equations read

$\begin{split}\dot{x}&=f(x)-y-z+I,\\\ \dot{y}&=\eta(x-ay),\\\
\dot{z}&=\varepsilon(x-bz),\end{split}$ (29)

where

$f(x)=-x+\lvert x+1\rvert-\lvert x-1\rvert,$

and $\varepsilon$ is a small parameter, i.e., $0<\varepsilon\ll 1$.

As shown in [12, 13], the DK model exhibits complex oscillations of elliptic
bursting type. This periodic behavior can be understood through the slow
dynamics of variable $z$, which drives the fast variables $(x,y)$ between
quiescence (resting states of the fast subsystem) and burst regimes (spiking
states of the fast subsystem). In particular, in the DK model the slow
variable $z$ drives the system from quiescence to burst through a delayed
subcritical Hopf bifurcation, and it organises the reverse transition through
a delayed fold bifurcation of limit cycles.

fig. 9(a) displays a bursting periodic orbit of the DK model, superimposed
onto the bifurcation diagram of its fast subsystem. Therefore we also show the
critical manifold as family of equilibria of the fast subsystem, together with
the subcritical Hopf bifurcation point at the beginning of the oscillatory
regime. The fast subsystem possesses two families of limit cycles, which we
represent as transparent surfaces: the invariant cone (red) formed by the
unstable limit cycles born at the subcritical Hopf bifurcation, as well as the
invariant cylinder (blue) of stable limit cycles born through a restabilizing
fold of cycles bifurcation, which ends the burst phase in the full system.

The bursting phenomenon in the DK model is obtained by setting parameters in
such a way that the system presents a unique equilibrium point. This
equilibrium point is an attracting node-focus when it is in the lateral
regions, and it becomes a saddle-focus when it is in the central region, which
is the case for the bursting cycle presented in fig. 9.

From [12] we consider the bursting regime with parameter values $a=0.8$,
$\eta=0.5$ and $b=0.5$. Parameter $I$ is chosen such that $-2.25<I<2.25$, in
order for the DK model to have two virtual equilibrium points in the lateral
regions and a real equilibrium point
$\mathbf{e}=-\frac{I}{1-1/a-1/b}(1,\frac{1}{a},\frac{1}{b})$ in the central
one.

|
---|---
(a) | (b)
Figure 9: Eliptic bursting in the DK model. (a) A 3 dimensional view of a
bursting cycle in the DK model. We can see the canonical slow manifolds (the
repelling slow manifold in red and the attracting one in blue), the unstable
cone (in red) and the stable cylinder (in blue) of fast subsystem limit
cycles, together with a bursting orbit. (b) Time-series of the solution
plotted in (a); $x(t)$ in blue, $y(t)$ in green and $z(t)$ in red. We take
$I=-1.5$ and $\varepsilon=10^{-3}$.

The critical manifold of the DK model is the polygonal curve

$\mathcal{S}_{0}=\Big{\\{}(x,y,z):f(x)-y-z+I=0,\;x-ay=0\Big{\\}}.$

It is normally hyperbolic everywhere except at the points
$(-1,-\frac{1}{a},-1+\frac{1}{a}+I)$ and $(1,\frac{1}{a},1-\frac{1}{a}+I)$,
located at the intersection between $\mathcal{S}_{0}$ and the switching planes
$\\{x=\pm 1\\}$, at which the equilibrium points change stability.

Owing to the extension of Fenichel Theory to PWL slow-fast systems developed
in [31], the attracting and repelling branches of the critical manifold will
perturb for $\varepsilon>0$ small enough, to attracting and repelling slow
manifolds, respectively. In fig. 9(a) we show the canonical attracting slow
manifold (in blue), contained in the half-space $L=\\{x<-1\\}$ (left region),
as well as the canonical repelling slow manifold (in red), contained in the
strip $M=\\{-1<x<1\\}$ (middle region). In the following result we give
conditions for the slow manifolds to connect.

######  Lemma 7.

Consider the DK system (29). For $\varepsilon$ small enough, the attracting
slow manifold (in the half-space $L$) and the repelling slow manifold (in the
middle strip $M$) connect if, and only if, the equilibrium point $\mathbf{e}$
is located on the boundary $\\{x=-1\\}$.

###### Proof.

For $\varepsilon$ small enough, we can consider $\eta a\neq\varepsilon b$. In
this case, the attracting and the repelling slow manifolds are given by the
line segments parameterised by $\tau$

$\mathbf{u}_{j}(\tau)=\Big{\\{}\mathbf{u}_{0,j}+\tau\mathbf{v}_{j}\Big{\\}},\quad\mathbf{v}_{j}=\left((\lambda_{j}+\eta
a)(\lambda_{j}+\varepsilon b),\eta(\lambda_{j}+\varepsilon
b),\varepsilon(\lambda_{j}+\eta a)\right)^{T}$ (30)

where $j\in\\{L,M\\}$, $\mathbf{u}_{0,j}=(x_{0,j},y_{0,j},z_{0,j})$ is the
equilibrium point in the $j$ region, $\lambda_{j}$ is the slow eigenvalue in
the $j$ region and $\mathbf{v}_{j}$ the eigenvector associated with
$\lambda_{j}$. Parameter $\tau$ then satisfies:
$x_{L}(\tau)=x_{0,L}+\tau(\lambda_{L}+\eta a)(\lambda_{L}+\varepsilon b)<-1$
and $|x_{M}(\tau)|<1$; see [31] for details.

Suppose that the slow manifolds connect on the boundary $x=-1$. Therefore,
there exist $\tau_{L}$ and $\tau_{M}$ such that (i)
$\mathbf{u}_{L}(\tau_{L})=\mathbf{u}_{M}(\tau_{M})$ and (ii) their first
coordinates $x_{L}(\tau_{L})=x_{M}(\tau_{M})=-1$. Hence, from (ii) we obtain
that

$\tau_{j}=-\frac{1+x_{0,j}}{(\lambda_{j}+\eta a)(\lambda_{j}+\varepsilon
b)},\quad j\in{L,M}.$

|
---|---
(a) | (b)
Figure 10: Way-in/way-out function and maximal delay for the DK model,
computed for the parameter values $a=0.8$, $\eta=0.5$, $b=0.5$, $I=2$, and;
(a) $\varepsilon=10^{-3}$; (b) $\varepsilon=10^{-5}$. After an initial
decrease, the function increases and stabilises around the value of the
maximal delay. Observe that, as $\varepsilon$ is divided by $10^{2}$, so the
maximal delay is divided by a similar amount.

From (ii), equalling $y_{L}(\tau_{L})=y_{M}(\tau_{M})$ and
$z_{L}(\tau_{L})=z_{M}(\tau_{M})$ in expression (30) and considering also the
$\tau_{j}$ expressions obtained, we obtain the following linear system
depending on the $x$-coordinates of the equilibrium points,

$\begin{array}[]{l}\displaystyle\frac{x_{0,L}}{a}-\frac{(1+x_{0,L})\eta}{(\lambda_{L}+\eta
a)}=\frac{x_{0,M}}{a}-\frac{(1+x_{0,M})\eta}{(\lambda_{M}+\eta a)},\\\
\displaystyle\frac{x_{0,L}}{b}-\frac{(1+x_{0,L})\varepsilon}{(\lambda_{L}+\varepsilon
b)}=\frac{x_{0,M}}{b}-\frac{(1+x_{0,M})\varepsilon}{(\lambda_{M}+\varepsilon
b)}.\end{array}$ (31)

Notice that the determinant of the matrix associated with each system is given
by

$\det(\Lambda)=\frac{-\lambda_{L}\lambda_{M}(\lambda_{M}-\lambda_{L})(\varepsilon
b-\eta a)}{ab(\lambda_{L}+\eta a)(\lambda_{M}+\eta a)(\lambda_{L}+\varepsilon
b)(\lambda_{M}+\varepsilon b)}.$

In order to see if system (31) has a unique solution, we have to determine
whether the previous determinant is zero or not. To do so, we consider the
characteristic polynomials of system (29) in each region, which are given by

$p_{j}(\lambda)=-\lambda^{3}+(f^{\prime}-\varepsilon b-\eta
a)\lambda^{2}-(\eta+\varepsilon-\varepsilon bf^{\prime}+\varepsilon b\eta
a-\eta af^{\prime})\lambda-\eta\varepsilon(b+a-baf^{\prime})$

where $f^{\prime}$ stands for the derivative of $f(x)$ in the corresponding
region $j\in\\{L,M\\}$. Since the determinant is zero when
$\lambda_{L}=\lambda_{M}$, this leads to
$p_{L}(\lambda_{L})-p_{M}(\lambda_{L})=0$. Thus,

$\lambda_{L}^{2}+(\varepsilon b+\eta a)\lambda_{L}+\varepsilon b\eta a=0,$

which implies that either $\lambda_{L}=-\varepsilon b$ or $\lambda_{L}=-\eta
a$. In any case, $p_{j}(-\eta a)=\eta(\eta a-\varepsilon b)$,
$p_{j}(-\varepsilon b)=\varepsilon(\varepsilon b-\eta a)$ and both imply that
$\varepsilon b=\eta a$. Therefore, $\det(\Lambda)\neq 0$ and so the system has
a unique solution at $x_{0,L}=x_{0,M}=-1$. In particular, this solution is the
equilibrium point of the differential equation located in between the half-
space $L$ and the middle strip $M$. ∎

Following lemma 7, when the equilibrium point is not on the boundary, the
attracting and the repelling slow manifolds do not connect. Locally, this
configuration is similar to those analyzed in section 2, theorem 1. We
concluded that in this no-connection scenario, the maximal delay decreases to
zero as $\varepsilon$ tends to zero. This behavior of the maximal delay can
hence be observed in the DK model; see fig. 10 where the way-in/way-out
function is plotted for two different values of $\varepsilon$. In panel (a) we
take $\varepsilon=10^{-3}$, whereas in panel (b) $\varepsilon=10^{-5}$. In the
latter case we observe that the maximal delay is also divided by approximately
$10^{2}$, which suggests a dependence on $\varepsilon$. theorem 1 also
provides, in the minimal model, the dependence of the maximal delay on
$\varepsilon$, for $\varepsilon$ small enough. In particular, the maximal
delay decreases like $u_{1}\varepsilon+u_{2}\varepsilon\ln(\varepsilon)$. In
fig. 11, we represent the maximal delay in the DK model as a function of
$\varepsilon$ for two values of parameter $I$. For suitable values of the
constants $u_{1}$ and $u_{2}$ we also trace the curve
$u_{1}\varepsilon+u_{2}\varepsilon\ln(\varepsilon)$ which shows that it fits
well the predicted behavior of the maximal delay. On the other hand, as
$\varepsilon$ increases, the maximal delay tends to the value $1$ or $1/2$,
which corresponds with the coordinate of the equilibrium point for the chosen
$I$ value, that is, when $I=2$ or $I=0$, respectively.

Figure 11: Maximal delay as a function of $\varepsilon$ in the DK model,
computed over an interval of $\varepsilon$-values, for two different values of
parameter $I$. These values are obtained by setting initial conditions on the
attracting slow manifold at $\\{x=-1\\}$. The remaining parameter values are
$a=0.8$, $\eta=0.5$ and $b=0.5$, and the radius of the tubular neighbourhood
is equal to $1$. For suitable values of the constants $u_{1}$ and $u_{2}$, we
also plot the curve $u_{1}\varepsilon+u_{2}\varepsilon\ln(\varepsilon)$
(dashed), which fits well with the computed maximal delay curves.

Next, we slightly modify the DK model in order to allow a connection between
the attracting slow manifold (contained in the half-space L) and the repelling
one (in the middle strip M). That is, we add an extra linearity zone in
between these two regions. Since we maintain the equilibrium in the middle
strip $\\{|x|\leq 1\\}$, the expected behavior of the slow passage is the one
described in section 3.1 about the presence of a buffer point, in particular
in theorem 3. Therefore, let us consider system (29), where we rewrite $f$ as

$f(x)=-x-2+(1+s)\frac{\lvert x-\rho\rvert+(x-\rho)}{2}+(1-s)\frac{\lvert
x-\mu\rvert+(x-\mu)}{2}-\lvert x-1\rvert-(x-1).$

This extended version of the DK model has 3 new parameters, namely $\rho$ and
$\mu$, which are the left and the right boundary of the new region,
respectively, and parameter $s$, which is the slope of the $x$-nullcline in
this new region. Considering $\mu=\frac{\rho(s+1)+2}{s-1}$, we set one of
these parameters in such a way that, for $x\not\in(\rho,\mu)$, the vector
field is exactly the same as the classical DK model.

The slow manifolds in the modified DK are given by equation (30). Hence, the
intersection points with these manifolds are given by

$\displaystyle\tilde{\mathbf{p}}^{a}=\mathbf{u}_{0,L}+\rho((\lambda_{L}+\eta
a)(\lambda_{L}+\varepsilon b),\eta(\lambda_{L}+\varepsilon
b),\varepsilon(\lambda_{L}+\eta a)),$
$\displaystyle\tilde{\mathbf{p}}^{r}=\mathbf{u}_{0,R}+\mu((\lambda_{R}+\eta
a)(\lambda_{R}+\varepsilon b),\eta(\lambda_{R}+\varepsilon
b),\varepsilon(\lambda_{R}+\eta a)).$

In order to guarantee the connection between the two slow manifolds, we solve
the following system

$\varphi(t;0,\tilde{\mathbf{p}}^{a})=\tilde{\mathbf{p}}^{r}$

by applying the shooting method considering $\varepsilon,\rho$ and $s$ as
unknowns. This procedure allows us to obtain the aforementioned connection up
to a prescribed error.

|
---|---
(a) | (b)
Figure 12: Eliptic bursting in the modified DK model. (a) 3 dimensional view
of the bursting cycle in the extended DK model, with a fourth region between
the L and M regions. Also shown are the invariant branches of the critical
manifold (the repelling branch in red and the attracting one in blue), the
unstable cone (in red) and the stable cylinder (in blue) of fast subsystem
limit cycles. (b) Time-series of the solution plotted in (a); $x(t)$ in blue,
$y(t)$ in green and $z(t)$ in red. We take $I=-1.4$ and $\varepsilon=10^{-3}$.

Under suitable conditions for the connection to occur, fig. 12 shows a
bursting cycle in the modified DK model, whose way-in/way-out function takes
the form depicted in fig. 13. Following theorem 3, the way-in/way-out function
has to reach the same plateau value for every value of $\varepsilon$, this
plateau value being related to the position of the equilibrium point. However,
after comparing panels (a) and (b) in fig. 13, the maximal delay seems to
decrease when $\varepsilon$ gets smaller, hinting at the fact that the
connection has not occurred. This discrepancy is not only related to the
tolerance of the shooting method we use to approach the connection conditions,
but also affected by the strong repulsion of the repelling slow manifold.

To end this section, we study the latter effect through the analysis of the
ratio between the real eigenvalue $\lambda_{M}$ and the real part $\alpha_{M}$
of the complex one, which organises the linear dynamics of the system in the
strip $\\{\mu<x<1\\}$. We have

|
---|---
(a) | (b)
Figure 13: Way-in/way-out function in the modified DK model. Fixed parameter
values are $a=0.8$, $\eta=0.5$, $b=0.5$, and $I=-1.4$ (a)
$\varepsilon=10^{-3}$, $s=0.4552$ and $\rho=-1.0067$, with the radius of the
tubular neighbourhood as $\delta=10^{-2}$, (b) $\varepsilon=10^{-5}$,
$s=0.461129$ and $\rho=-1.000068$, with $\delta=10^{-2}$. The paramaters $s$
and $\rho$ are chosen to guarantee a connection between the attracting and the
repelling slow manifolds. The red line in panel (a) corresponds to the buffer
point value, which is not shown in (b) for a clearer representation.

$\displaystyle\lambda_{M}=\frac{1}{a+1}(ab-
a-b)\varepsilon+O(\varepsilon^{2}),$
$\displaystyle\alpha_{M}=\frac{1-a\eta}{2}+\frac{a}{2(a+1)}(1-2b)\varepsilon+O(\varepsilon^{2}).$

As can be seen, $\lambda_{M}$ is of order 1 in $\varepsilon$, whereas
$\alpha_{M}$ is of order 0 in $\varepsilon$. As a result, orbits escape very
fast from the repelling slow manifold, increasing the stiffness of the problem
as $\varepsilon$ tends to zero. In order to avoid this effect, we set
parameter $\eta$ such that

$\eta=1/a+\eta_{1}\varepsilon\quad\text{with}\quad\eta_{1}<0,$ (32)

which allows for $\lambda_{M}$ and $\alpha_{M}$ to be of the same order in
$\varepsilon$. In fig. 14, we show the way-in/way-out function of the modified
Doi-Kumagai model, after imposing the relation (32). As can be observed, the
maximal delay behaves in the way described in theorem 3, which predicts that
the maximal delay will be reached at the buffer point, located at $a-\mu$,
independently of the chosen $\varepsilon$ value.

|
---|---
(a) | (b)
Figure 14: Way-in/way-out function in the modified DK model, with eigenvalues
of the same order. Fixed parameter values are $\eta=1/a-10\varepsilon$,
$a=0.8$, $b=0.5$ and $I=-1.4$. The red lines correspond to the buffer point.
(a) $\varepsilon=10^{-3}$, $s=0.6165$ and $\rho=-1.0018$, with
$\delta=10^{-2}$, (b) $\varepsilon=10^{-5}$, $s=0.625649$ and $\rho=1.000018$,
with $\delta=10^{-2}$. Parameters $s$ and $\rho$ are chosen to guarantee a
connection between the attracting and the repelling slow manifolds.

##  5 Discussion

In this paper, we have analyzed the phenomenon of slow passage through a Hopf
bifurcation in the context of piecewise linear slow-fast dynamical systems,
both qualitatively and quantitatively. Delayed loss of stability is a common
behavior naturally present in smooth systems and we have shown that the PWL
framework allows to reproduce this delay. However, different scenarios can
arise depending on the minimal model we use, as described in all cases through
the way-in/way-out function and the maximal delay’s asymptotic value.

We first considered a minimal PWL model exhibiting a delayed loss of
stability, in such a way that the maximal delay $z_{d}$ depends on
$\varepsilon$, and tends to zero with $\varepsilon$. We obtained an explicit
expression for the dependence of $z_{d}$ on $\varepsilon$ by taking advantage
of the PWL framework. The main ingredient for such a behavior of the maximal
delay is the $O(\varepsilon)$ distance between the attracting and the
repelling slow manifolds.

We then extended this minimal PWL model by introducing a third linearity
region allowing the connection between the canonical slow manifolds. This
strategy has proven useful to fully analyse PWL versions of smooth slow-fast
dynamics near a quadratic fold of the critical manifold, that is, recovering
also the canard regime; see [10, 14]. Hence, the connection of the canonical
slow manifolds allowed the maximal delay to tend to infinity. This behavior is
more similar to what happens in smooth systems. However, because in the PWL
context the divergence along the repelling slow manifold remains constant from
the moment it crosses the separation plane, the slow passage phenomenon has a
very unstable character. Hence, given that the conditions for a connection to
occur are not, in practice, entirely satisfied, and due to round-off errors
during the simulations, the computed way-in/way-out function always behaves as
if the maximal delay has a finite value that cannot be exceeded, similar to a
buffer point.

In order to study the presence of a proper buffer point, we also considered
the case of a three-region system having an equilibrium point in the central
region. Results show the appearance of a boundary for the solutions, which
does not depend on the previously observed drawbacks, revealing the existence
of the buffer point.

Since the slow passage phenomenon is a key element to generate bursting
dynamics, we checked the validity of our results on a PWL elliptic bursting
model, namely, the DK model. This adds to the existing literature on bursting
dynamics in the PWL context [5, 6]. However, the slow passage phenomenon in
this model appears to be too weak since, as in section 2, the distance between
the attracting and the repelling slow manifold is linear with $\varepsilon$.
To make it stronger, we added to the DK model a new linearity region allowing
for the connection between both the attracting and the repelling slow
manifolds. However, even though we could connect them, the problem appeared to
be so unstable that we could not observe any significant difference with the
case where the slow manifolds do not connect. Thus, even when theoretically
the maximal delay must approach the buffer point independently of
$\varepsilon$, in practice, the numerical precision cannot afford the
stiffness of the problem and the maximal delay decreases with $\varepsilon$.

To overcome the stiffness problem, we finally studied a suitable parameter set
of the modified DK system, not only allowing for the connection, but also
providing the real eigenvalue and the real part of the complex eigenvalues of
the linear system in $M$ to be of the same order in $\varepsilon$. This
reduced the stiffness of the system and maintained the maximal delay high,
irrespective the value of $\varepsilon$. However, further efforts would be
needed to find suitable parameter sets allowing to exhibit bursting
oscillations for this new scenario as well. This is an interesting question
for future work.

##  6 Acknowledgments

JP, AET and CV are partially supported by the Ministerio de Ciencia,
Innovación y Universidades (MCIU) project PID2020-118726GB-I00. AET and CV are
also partially supported by the Ministerio de Economia y Competitividad
through the project MTM2017-83568-P (AEI/ERDF,EU).

## References

  * [1] S. Baer, T. Erneux, and J. Rinzel.  The slow passage through a hopf bifurcation: Delay, memory effects, and resonance.  SIAM Journal on Applied Mathematics, 49(1):55–71, 1989.
  * [2] E. Benoît, J.-L. Callot, F. Diener, and M. Diener.  Chasse au canard.  Collectanea Mathematica, 32(1-2):37–119, 1981.
  * [3] V. Carmona, S. Fernández-García, and A. E. Teruel.  Saddle-node of limit cycles in planar piecewise linear systems and applications.  Discrete & Continuous Dynamical Systems, 39(9):5275–5299, 2019\.
  * [4] S. Coombes and P. Bressloff.  Bursting: The Genesis of Rhythm in the Nervous System.  World Scientific, 2005.
  * [5] B. Deng.  Conceptual circuit models of neurons.  Journal of Integrative Neuroscience, 8(03):255–297, 2009.
  * [6] M. Desroches, S. Fernández-García, and M. Krupa.  Canards in a minimal piecewise-linear square-wave burster.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 26(7):073111, 2016.
  * [7] M. Desroches, S. Fernández-García, M. Krupa, R. Prohens, and A. E. Teruel.  Piecewise-linear (pwl) canard dynamics : Simplifying singular perturbation theory in the canard regime using piecewise-linear systems.  In Nonlinear Systems, volume 1 of Mathematical Theory and Computational Methods. Springer, 2018.
  * [8] M. Desroches, E. Freire, S. J. Hogan, E. Ponce, and P. Thota.  Canards in piecewise-linear systems: explosions and super-explosions.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 469(2154):20120603, 2013.
  * [9] M. Desroches, J. Guckenheimer, B. Krauskopf, C. Kuehn, H. M. Osinga, and M. Wechselberger.  Mixed-mode oscillations with multiple time scales.  SIAM Review, 54(2):211–288, 2012.
  * [10] M. Desroches, A. Guillamon, E. Ponce, R. Prohens, S. Rodrigues, and A. E. Teruel.  Canards, folded nodes and mixed-mode oscillations in piecewise-linear systems.  SIAM Review, 58(4):653–691, 2016.
  * [11] F. Diener and M. Diener.  Maximal delay.  In E. Benoît, editor, Dynamic Bifurcations, volume 1493 of Lecture Notes in Mathematics, pages 71–86. Springer, Berlin, Heidelberg, 1991.
  * [12] S. Doi and S. Kumagai.  Complicated slow oscillations with simple switching dynamics in piecewise linear neuronal model.  In The 2004 47th Midwest Symposium on Circuits and Systems, 2004\. MWSCAS’04., volume 2, pages II–II. IEEE, 2004.
  * [13] S. Doi and S. Kumagai.  Generation of very slow neuronal rhythms and chaos near the hopf bifurcation in single neuron models.  Journal of Computational Neuroscience, 19(3):325–356, 2005.
  * [14] S. Fernández-García, M. Desroches, M. Krupa, and A. E. Teruel.  Canard solutions in planar piecewise linear systems with three zones.  Dynamical Systems, 31(2):173–197, 2016.
  * [15] E. Freire, E. Ponce, and F. Torres.  Hopf-like bifurcations in planar piecewise linear systems.  Publicacions Matemàtiques, 41(1):135–148, 1997.
  * [16] X. Han and Q. Bi.  Slow passage through canard explosion and mixed-mode oscillations in the forced van der pol’s equation.  Nonlinear Dynamics, 68(1):275–283, 2012.
  * [17] L. Holden and T. Erneux.  Slow passage through a hopf bifurcation: From oscillatory to steady state solutions.  SIAM Journal on Applied Mathematics, 53(4):1045–1058, 1993.
  * [18] L. Holden and T. Erneux.  Understanding bursting oscillations as periodic slow passages through bifurcation and limit points.  Journal of Mathematical Biology, 31(4):351–365, 1993.
  * [19] J. Honerkamp, G. Mutschler, and R. Seitz.  Coupling of a slow and a fast oscillator can generate bursting.  Bulletin of Mathematical Biology, 47(1):1–21, 1985.
  * [20] W. Huagan, Y. Ye, M. Chen, X. Quan, and B. Bao.  Extremely slow passages in low-pass filter-based memristive oscillator.  Nonlinear Dynamics, 97:2339–2353, 2019.
  * [21] E. M. Izhikevich.  Neural excitability, spiking and bursting.  International Journal of Bifurcation and Chaos, 10(06):1171–1266, 2000.
  * [22] E. M. Izhikevich.  Subcritical elliptic bursting of bautin type.  SIAM Journal on Applied Mathematics, 60(2):503–535, 2000.
  * [23] M. Krupa and M. Wechselberger.  Local analysis near a folded saddle-node singularity.  Journal of Differential Equations, 248(12):2841–2888, 2010.
  * [24] Yu. A. Kuznetsov and S. Rinaldi.  Remarks on food chain dynamics.  Mathematical biosciences, 134(1):1–33, 1996.
  * [25] C. Lobry.  Dynamic bifurcations.  In E. Benoît, editor, Dynamic Bifurcations, volume 1493 of Lecture Notes in Mathematics, pages 1–13. Springer, Berlin, Heidelberg, 1991.
  * [26] A. Neishtadt.  On stability loss delay for dynamical bifurcations.  Discrete & Continuous Dynamical Systems - Series S, 2(4):897, 2009\.
  * [27] A. I. Neishtadt.  Persistence of stability loss for dynamical bifurcations I.  Differential Equations, 23:1385–1391, 1987.
  * [28] A. I. Neishtadt.  Persistence of stability loss for dynamical bifurcations II.  Differential Equations, 24:171–176, 1988.
  * [29] D. Premraj, K. Suresh, T. Banerjee, and K. Thamilmaran.  An experimental study of slow passage through hopf and pitchfork bifurcations in a parametrically driven nonlinear oscillator.  Communications in Nonlinear Science and Numerical Simulation, 37:212–221, 2016.
  * [30] R. Prohens and A. E. Teruel.  Canard trajectories in 3d piecewise linear systems.  Discrete and Continuous Dynamical Systems - Series A, 33(10):4595–4611, 2013.
  * [31] R. Prohens, A. E. Teruel, and C. Vich.  Slow–fast n-dimensional piecewise linear differential systems.  Journal of Differential Equations, 260(2):1865–1892, 2016.
  * [32] J. Rinzel.  A formal classification of bursting mechanisms in excitable systems.  In International Congress of Mathematicians, Berkeley, California, USA, August 3-11, 1986, volume II, pages 1578–1593. American Mathematical Society, 1987.
  * [33] D. J. W. Simpson.  A compendium of hopf-like bifurcations in piecewise-smooth dynamical systems.  Physics Letters A, 382(35):2439–2444, 2018.
  * [34] J. Su, J. Rubin, and D. Terman.  Effects of noise on elliptic bursters.  Nonlinearity, 17(1):133–157, 2003.

##  Appendix A Solutions of the systems

Considering the two-region system (2)-(3), its local solution through the
point $\mathbf{p}=(x_{0},y_{0},z_{0})$, which we denote as
$\mathbf{u}(t;0,\mathbf{p})$, has the first and second coordinates given by

$\begin{split}x_{L}(t;0,\mathbf{p})&=z(t;0,\mathbf{p})-m\varepsilon+\mathrm{e}^{-\frac{m}{2}t}\Big{(}(x_{0}-z_{0}+m\varepsilon)\cos\left(\xi_{m}t\right)\\\
&+\frac{-m(z_{0}-m\varepsilon+x_{0})-2(y_{0}+\varepsilon)}{\sqrt{4-m^{2}}}\sin\left(\xi_{m}t\right)\Big{)},\\\
y_{L}(t;0,\mathbf{p})&=-m(z(t;0,\mathbf{p})-m\varepsilon)-\varepsilon\\\
&+\mathrm{e}^{-\frac{m}{2}t}\Big{(}(m(z_{0}-m\varepsilon)+(y_{0}+\varepsilon))\cos\left(\xi_{m}t\right)\\\
&+\frac{m^{2}(z_{0}-m\varepsilon)+m(y_{0}+\varepsilon)+2(x_{0}-z_{0}+m\varepsilon)}{\sqrt{4-m^{2}}}\sin\left(\xi_{m}t\right)\Big{)},\end{split}$
(33)

if it is contained in the left region, and

$\begin{split}x_{R}(t;0,\mathbf{p})&=z(t;0,\mathbf{p})+k\varepsilon+\mathrm{e}^{\frac{k}{2}t}\Big{(}(x_{0}-z_{0}-k\varepsilon)\cos\left(\xi_{k}t\right)\\\
&+\frac{k(z_{0}+k\varepsilon+x_{0})-2(y_{0}+\varepsilon)}{\sqrt{4-k^{2}}}\sin\left(\xi_{k}t\right)\Big{)},\\\
y_{R}(t;0,\mathbf{p})&=k(z(t;0,\mathbf{p})+k\varepsilon)-\varepsilon\\\
&+\mathrm{e}^{\frac{k}{2}t}\Big{(}(-k(z_{0}+k\varepsilon)+(y_{0}+\varepsilon))\cos\left(\xi_{k}t\right)\\\
&+\frac{k^{2}(z_{0}+k\varepsilon)-k(y_{0}+\varepsilon)+2(x_{0}-z_{0}-k\varepsilon)}{\sqrt{4-k^{2}}}\sin\left(\xi_{k}t\right)\Big{)},\end{split}$
(34)

if it is contained in the right region, where
$\xi_{m}=\frac{\sqrt{4-m^{2}}}{2}$, $\xi_{k}=\frac{\sqrt{4-k^{2}}}{2}$. The
third coordinate remains similar in all regions and it is the slow drift given
by $z(t;0,\mathbf{p})=z_{0}+\varepsilon t$.

In the case of the three-region system (2)-(13), the local solution through
$\mathbf{p}$ is given by (33) and (34), if it is contained in the left or
right regions and by

$\begin{split}x_{C}(t;0,\mathbf{p})&=z(t;0,\mathbf{p})+l\varepsilon+\mathrm{e}^{\frac{l}{2}t}\Big{(}(x_{0}-(z_{0}+l\varepsilon))\cos\left(\xi_{l}t\right)\\\
&+\frac{1}{\xi_{l}}(\frac{l}{2}(x_{0}+z_{0}+l\varepsilon)+n-(y_{0}+\varepsilon))\sin\left(\xi_{l}t\right)\Big{)},\\\
y_{C}(t;0,\mathbf{p})&=l(z(t;0,\mathbf{p})+l\varepsilon)+n-\varepsilon+\mathrm{e}^{\frac{l}{2}t}\Big{(}(y_{0}-(l(z_{0}+l\varepsilon)+n-\varepsilon))\cos\left(\xi_{l}t\right)\\\
&+\frac{1}{2\xi_{l}}(l^{2}(z_{0}+l\varepsilon)+2(x_{0}-(z_{0}+l\varepsilon))+l(n-(y_{0}+\varepsilon)))\sin\left(\xi_{l}t\right)\Big{)},\end{split}$
(35)

if it is contained in the central region, where
$\xi_{l}=\frac{\sqrt{4-l^{2}}}{2}$. The third coordinate is also
$z(t;0,\mathbf{p})=z_{0}+\varepsilon t$.

The local expression of the solution in the central region of the buffer-point
system (22)-(13) is given by:

$\mathbf{u}_{C}(t;0,\mathbf{p})=P\mathrm{e}^{Jt}P^{-1}\left(\mathbf{p}+A^{-1}\mathbf{a}_{0}\right)-A^{-1}\mathbf{a}_{0},$
(36)

where

$A=\begin{pmatrix}-\varepsilon&-1&0\\\ 1&0&-1\\\
-\varepsilon&0&0\end{pmatrix},\quad P=\begin{pmatrix}1&1&0\\\
0&-\varepsilon&-1\\\
1&0&\varepsilon\end{pmatrix},\quad\mathbf{a}_{0}=\begin{pmatrix}n\\\ 0\\\
\varepsilon\end{pmatrix},$
$\mathrm{e}^{Jt}=\begin{pmatrix}\mathrm{e}^{-\varepsilon t}&0&0\\\
0&\cos(t)&\sin(t)\\\ 0&-\sin(t)&\cos(t)\end{pmatrix}.$

##  Appendix B Way-in/way-out function

In this section we provide the pseudo-code algorithm used to compute the way-
in/way-out function given a tubular neighborhood of radius $\delta$. When a
two-region system is considered, we use $\mu=0$ while in the three-region
systems, we consider $\mu$ being the $x$ boundary between the central region
and the right one.

1: Define vector of $z$-coordinates $z_{1},\ldots,z_{N}$

2: Let $\delta$ be the radius of the neighbourhood

3: for each $i=1,\ldots,N$ do

4: Let
$\mathbf{p}_{i}=(z_{i}-m\varepsilon,\varepsilon-m(z_{i}-m\varepsilon)+\delta,z_{i})$

5: Solve the ODE system
$\mathbf{u}(t_{i,j};0,\mathbf{p}_{i})=(x(t_{i,j}),y(t_{i,j}),z(t_{i,j}))$

6: Find the time $t_{i,s}$ such that $x(t_{i,s})=\mu$

7: Take $j=s$ and $d_{i,j}=\delta$

8: while $d_{i,j}\leq\delta$ do

9: Find the plane $\pi(z_{0},t_{i,j})$ generated by the complex eigenvectors

10: Find the intersection point
$\hat{\mathbf{u}}(t_{i,j})=(\hat{x}(t_{i,j}),\hat{y}(t_{i,j}),\hat{z}(t_{i,j}))=\mathcal{S}_{\varepsilon}^{r}\cap\pi(z_{0},t_{i,j})$

11: Find the distance
$d_{i,j}=d(\hat{\mathbf{u}}(t_{i,j}),\mathbf{u}(t_{i,j};0,\mathbf{p}_{i}))$

12: if $d_{i,j}>\delta$ then

13: out: $\hat{z}_{i}=\hat{z}(t_{i,j})$

14: else

15: $j=j+1$

16: end if

17: end while

18: end for

19: return pairs $(z_{i},\hat{z}_{i})$

Algorithm 1 Way-in/way-out function