# Dynamics of a Piecewise-Linear Morris-Lecar Model: Bifurcations and Spike Adding

Jordi Penalva, Mathieu Desroches, Antonio Teruel, Catalina Vich

► To cite this version:

Jordi Penalva, Mathieu Desroches, Antonio Teruel, Catalina Vich. Dynamics of a Piecewise-Linear Morris-Lecar Model: Bifurcations and Spike Adding. Journal of Nonlinear Science, 2024, 34 (3), pp.52. 10.1007/s00332-024-10029-3. hal-04545904

HAL Id: hal-04545904

https://inria.hal.science/hal-04545904

Submitted on 14 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamics of a piecewise-linear Morris-Lecar model: bifurcations and spike adding

J. Penalva¹*, M. Desroches², A. E. Teruel¹ and C. Vich¹

¹*Departament de Matemàtiques i Informàtica & IAC3, Universitat de les Illes Balears, Palma, 07141, Balearic Islands, Spain.

²MathNeuro Team, Inria Centre at Université Côte d'Azur, Sophia Antipolis, 06902, France.

*Corresponding author(s). E-mail(s): j.penalva@uib.cat;

## Abstract

Multiple-timescale systems often display intricate dynamics, yet of great mathematical interest and well suited to model real-world phenomena such as bursting oscillations. In the present work, we construct a piecewise-linear version of the Morris-Lecar neuron model, denoted PWL-ML, and we thoroughly analyse its bifurcation structure with respect to three main parameters. Then, focusing on the homoclinic connection present in our PWL-ML, we study the slow passage through this connection when augmenting the original system with a slow dynamics for one of the parameters, thereby establishing a simplified framework for this slow-passage phenomenon. Our results show that our model exhibits equivalent behaviors to its smooth counterpart. In particular, we identify canard solutions that are part of spike-adding transitions. Focusing on the one-spike and on the two-spike scenarios, we prove their existence in a more straightforward manner than in the smooth context. In doing so, we present several techniques that are specific to the piecewise-linear framework and with the potential to offer new tools for proving the existence of dynamical objects in a wider context.

**Keywords:** Homoclinic bifurcation, piecewise-linear systems, slow-fast dynamics, slow passage, spike adding, bursting oscillations
# Dynamics of a piecewise-linear Morris-Lecar model: bifurcations and spike adding

J. Penalva¹*, M. Desroches², A. E. Teruel¹ and C. Vich¹

¹*Departament de Matemàtiques i Informàtica & IAC3, Universitat de les Illes Balears, Palma, 07141, Balearic Islands, Spain.

²MathNeuro Team, Inria Centre at Université Côte d'Azur, Sophia Antipolis, 06902, France.

*Corresponding author(s). E-mail(s): j.penalva@uib.cat;

## Abstract

Multiple-timescale systems often display intricate dynamics, yet of great mathematical interest and well suited to model real-world phenomena such as bursting oscillations. In the present work, we construct a piecewise-linear version of the Morris-Lecar neuron model, denoted PWL-ML, and we thoroughly analyse its bifurcation structure with respect to three main parameters. Then, focusing on the homoclinic connection present in our PWL-ML, we study the slow passage through this connection when augmenting the original system with a slow dynamics for one of the parameters, thereby establishing a simplified framework for this slow-passage phenomenon. Our results show that our model exhibits equivalent behaviors to its smooth counterpart. In particular, we identify canard solutions that are part of spike-adding transitions. Focusing on the one-spike and on the two-spike scenarios, we prove their existence in a more straightforward manner than in the smooth context. In doing so, we present several techniques that are specific to the piecewise-linear framework and with the potential to offer new tools for proving the existence of dynamical objects in a wider context.

**Keywords:** Homoclinic bifurcation, piecewise-linear systems, slow-fast dynamics, slow passage, spike adding, bursting oscillations
# 1 Introduction

Dynamical systems with multiple timescales, also often referred to as *slow-fast systems*, generate intricate dynamics while being omnipresent in applications. Whether explicit or implicit, the timescale separation allows to gather state variables in groups of slow and fast processes, which most of the time helps to break down the analysis by focusing on one group or the other. Systems with multiple timescales are ubiquitous in various application areas such as Neuroscience [1], Chemistry [2], Laser Physics [3], Population dynamics [4] or Material Science [5], to name a few.

In the application of multiple-timescale dynamical systems to brain activity, different behaviours exhibited by neurons can be mathematically understood in link with how close the dynamics of the corresponding differential system is to a bifurcation either in the full system or in the so-called *fast subsystem*; for instance, the generation of an action potential, also referred to as a spike; see [13]. In addition to action potentials, other oscillatory patterns are observed experimentally and faithfully reproduced by multiple-timescale models. In particular, one can mention *bursting oscillations* [14]).

Generally speaking, the phenomenon of bursting can be described as the result of the interplay between the fast components of the system, which in the (fast) singular limit possess coexisting stationary and periodic attractors, and the slow ones, which modulate the overall activity by driving the quasi-attractors (well defined in the fast subsystem) from quiescence to burst and vice versa. The transitions between the quasi-attractors are organised by a pair of fast-subsystem bifurcations; see [1]. These transitions from quiescence to burst and back should be seen as *dynamic bifurcations*. Dynamic bifurcations are key to understand the interesting dynamics underpinning the emergence of bursting oscillations (for the classification of bursting patterns, see [15–20]), which can only be fully deciphered by analysing both fast and slow subsystems. Dynamic bifurcation are often referred to as *slow passages through bifurcations*.

Such slow passages are characteristic of systems that exhibit slow-fast dynamics [21–25]. In fact, the slow-passage phenomenon is deeply related to connections between attracting and repelling slow manifold and associated fast-subsystem bifurcations points where these connections form. Furthermore, slow passages are related to *canard solutions*, which are generally known as objects that also stay close to repelling invariant manifolds for a considerable amount of time and are organised in parameter space along explosive branches [26–28].

As shown by Rinzel [1], the minimal configuration to obtain bursting oscillations, in the context of smooth deterministic systems, is to consider a three-dimensional system with 1 slow and 2 fast variables. In this context, in addition to the two slow-passages, a hysteresis loop is required in the full system for the bursting to occur.

From the viewpoint of bursting oscillations [1], the main profile where the burst terminates by a slow passage through a homoclinic bifurcation is
# 1 Introduction

Dynamical systems with multiple timescales, also often referred to as *slow-fast systems*, generate intricate dynamics while being omnipresent in applications. Whether explicit or implicit, the timescale separation allows to gather state variables in groups of slow and fast processes, which most of the time helps to break down the analysis by focusing on one group or the other. Systems with multiple timescales are ubiquitous in various application areas such as Neuroscience [1], Chemistry [2], Laser Physics [3], Population dynamics [4] or Material Science [5], to name a few.

In the application of multiple-timescale dynamical systems to brain activity, different behaviours exhibited by neurons can be mathematically understood in link with how close the dynamics of the corresponding differential system is to a bifurcation either in the full system or in the so-called *fast subsystem*; for instance, the generation of an action potential, also referred to as a spike; see [13]. In addition to action potentials, other oscillatory patterns are observed experimentally and faithfully reproduced by multiple-timescale models. In particular, one can mention *bursting oscillations* [14]).

Generally speaking, the phenomenon of bursting can be described as the result of the interplay between the fast components of the system, which in the (fast) singular limit possess coexisting stationary and periodic attractors, and the slow ones, which modulate the overall activity by driving the quasi-attractors (well defined in the fast subsystem) from quiescence to burst and vice versa. The transitions between the quasi-attractors are organised by a pair of fast-subsystem bifurcations; see [1]. These transitions from quiescence to burst and back should be seen as *dynamic bifurcations*. Dynamic bifurcations are key to understand the interesting dynamics underpinning the emergence of bursting oscillations (for the classification of bursting patterns, see [15–20]), which can only be fully deciphered by analysing both fast and slow subsystems. Dynamic bifurcation are often referred to as *slow passages through bifurcations*.

Such slow passages are characteristic of systems that exhibit slow-fast dynamics [21–25]. In fact, the slow-passage phenomenon is deeply related to connections between attracting and repelling slow manifold and associated fast-subsystem bifurcations points where these connections form. Furthermore, slow passages are related to *canard solutions*, which are generally known as objects that also stay close to repelling invariant manifolds for a considerable amount of time and are organised in parameter space along explosive branches [26–28].

As shown by Rinzel [1], the minimal configuration to obtain bursting oscillations, in the context of smooth deterministic systems, is to consider a three-dimensional system with 1 slow and 2 fast variables. In this context, in addition to the two slow-passages, a hysteresis loop is required in the full system for the bursting to occur.

From the viewpoint of bursting oscillations [1], the main profile where the burst terminates by a slow passage through a homoclinic bifurcation is
the so-called square-wave bursting, where the two fast-subsystem bifurcations involved are a saddle-node of equilibria, and then a saddle homoclinic bifurcation. The slow-passage through a homoclinic bifurcation is a known slow-fast scenario, as well as its relation to canard solutions [14, 31, 32].

Consider a slow-fast system of the form

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

where $\mathbf{x} \in \mathbb{R}^q$, $z \in \mathbb{R}$, $F$ is a sufficient smooth function and $0 \le \varepsilon \ll 1$ is the timescale ratio parameter; the overdot denotes differentiation with respect to the fast time $t$. The presence of $\varepsilon$ in (1) induces a timescale separation between the fast variables $\mathbf{x}$ and the slow variable $z$. Taking the $\varepsilon = 0$ limit in (1) yields the so-called fast subsystem [21], in which $z$ becomes a parameter for the $q$-dimensional differential equation $\dot{\mathbf{x}} = F(\mathbf{x}, z)$. The fast subsystem provides the first encounter with the critical manifold $S_0 = \{(\mathbf{x}, z) : F(\mathbf{x}, z) = 0\}$, which corresponds to its $z$-dependent family of equilibria.

Each subpart of the critical manifold that is formed by saddle equilibria (of the fast subsystem) admits a stable and an unstable manifold $W^s$ and $W^u$, respectively. These manifolds have a fibered structure and each fibers correspond to stable and unstable manifolds of a given saddle equilibrium of the critical manifold, respectively; see Figure 1. In the homoclinic scenario, for a specific value $z_0$ of $z$, both manifolds $W^{s,u}$ intersect along the homoclinic orbit and no other intersection exists in a neighbourhood of $z_0$; see Figure 2.

<figure><img src="image_1.png" /><figcaption>Fig. 1 Sketch of the homoclinic connection for $\varepsilon = 0$. There, without loss of generality, we can suppose that the parameter $z$ is moving: (a) from $z < z_0$; (b) passing to $z = z_0$; (c) to $z > z_0$. $SP$ corresponds to the saddle-point, while $W^s(S_0^{sd})$ and $W^u(S_0^{sd})$ correspond to the stable and unstable critical manifolds on the plane $z = z_0$.</figcaption></figure>

Following Fenichel's theorems [33], for $\varepsilon > 0$ small enough, each attracting subset $S_0^a$ (resp. repelling subset $S_0^r$ or saddle-type subset $S_0^{sd}$) of the critical manifold $S_0$ perturbs into an attracting slow manifold $S_\varepsilon^a$ (resp. a repelling slow manifold $S_\varepsilon^r$ or a saddle-type slow manifold $S_\varepsilon^{sd}$). In the saddle-type slow manifold case, its corresponding stable $W^s(S_0^{sd})$ and unstable $W^u(S_0^{sd})$ manifolds perturb into the stable $W^s(S_\varepsilon^{sd})$ and unstable $W^u(S_\varepsilon^{sd})$ manifolds of $S_\varepsilon^{sd}$. Moreover, from the unique intersection between $W^s(S_0^{sd})$ and $W^u(S_0^{sd})$ given
the so-called square-wave bursting, where the two fast-subsystem bifurcations involved are a saddle-node of equilibria, and then a saddle homoclinic bifurcation. The slow-passage through a homoclinic bifurcation is a known slow-fast scenario, as well as its relation to canard solutions [14, 31, 32].

Consider a slow-fast system of the form

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

where $\mathbf{x} \in \mathbb{R}^q$, $z \in \mathbb{R}$, $F$ is a sufficient smooth function and $0 \le \varepsilon \ll 1$ is the timescale ratio parameter; the overdot denotes differentiation with respect to the fast time $t$. The presence of $\varepsilon$ in (1) induces a timescale separation between the fast variables $\mathbf{x}$ and the slow variable $z$. Taking the $\varepsilon = 0$ limit in (1) yields the so-called fast subsystem [21], in which $z$ becomes a parameter for the $q$-dimensional differential equation $\dot{\mathbf{x}} = F(\mathbf{x}, z)$. The fast subsystem provides the first encounter with the critical manifold $S_0 = \{(\mathbf{x}, z) : F(\mathbf{x}, z) = 0\}$, which corresponds to its $z$-dependent family of equilibria.

Each subpart of the critical manifold that is formed by saddle equilibria (of the fast subsystem) admits a stable and an unstable manifold $W^s$ and $W^u$, respectively. These manifolds have a fibered structure and each fibers correspond to stable and unstable manifolds of a given saddle equilibrium of the critical manifold, respectively; see Figure 1. In the homoclinic scenario, for a specific value $z_0$ of $z$, both manifolds $W^{s,u}$ intersect along the homoclinic orbit and no other intersection exists in a neighbourhood of $z_0$; see Figure 2.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Sketch of the homoclinic connection for $\varepsilon = 0$. There, without loss of generality, we can suppose that the parameter $z$ is moving: (a) from $z < z_0$; (b) passing to $z = z_0$; (c) to $z > z_0$. SP corresponds to the saddle-point, while $W^s(S_0^{sd})$ and $W^u(S_0^{sd})$ correspond to the stable and unstable critical manifolds on the plane $z = z_0$.</figcaption></figure>

Following Fenichel's theorems [33], for $\varepsilon > 0$ small enough, each attracting subset $S_0^a$ (resp. repelling subset $S_0^r$ or saddle-type subset $S_0^{sd}$) of the critical manifold $S_0$ perturbs into an attracting slow manifold $S_\varepsilon^a$ (resp. a repelling slow manifold $S_\varepsilon^r$ or a saddle-type slow manifold $S_\varepsilon^{sd}$). In the saddle-type slow manifold case, its corresponding stable $W^s(S_0^{sd})$ and unstable $W^u(S_0^{sd})$ manifolds perturb into the stable $W^s(S_\varepsilon^{sd})$ and unstable $W^u(S_\varepsilon^{sd})$ manifolds of $S_\varepsilon^{sd}$. Moreover, from the unique intersection between $W^s(S_0^{sd})$ and $W^u(S_0^{sd})$ given
by the homoclinic connection at $\varepsilon = 0$, perturb finitely many intersections between $W^s(S_{\varepsilon}^{sd})$ and $W^u(S_{\varepsilon}^{sd})$ for $0 < \varepsilon \ll 1$.

<figure><img src="image_2.png" /><figcaption>**Fig. 2 3-dimensional representation of the homoclinic connection.** The homoclinic connection at $z = z_0$ corresponds to the intersection between the two-dimensional stable and unstable manifolds $W^{s,u}(S_0^{sd})$ of the saddle-type subset $S_0^{sd}$ of the critical manifold.</figcaption></figure>

Understanding the slow-passage through the homoclinic connection may help to describe the phenomenon of spike adding, which consists in the addition of new spikes during the bursting phase as a system parameter varies; see [14, 34, 35] for more details. Such phenomena have been widely studied in the smooth framework through phenomenological neural models, such as the Hindmarsh-Rose (HR) model [14, 36, 37], and biophysical ones, such as the Morris-Lecar (ML) model [32, 38, 39]. However, the slow-passage through the homoclinic connection has not been formally studied in piecewise linear (PWL) systems.

PWL models are quite known to reproduce qualitative and quantitative aspects of smooth systems, being more amenable to general analysis [40–42]. In recent years, they have been extended in order to incorporate slow-fast dynamics, and it has been shown in several studies that they can reproduce the most salient features of multiple-timescale systems, providing a new point of view to understand better these complex dynamics [10, 14, 43–46]. Canard-induced spike-adding transitions associated with square-wave bursting have been studied in [44] in a PWL version of the HR model. More recently, a slow passage through a Hopf-like bifurcation has been rigorously studied in [47], with application to elliptic bursting.

In the present work, we focus on the slow passage through a homoclinic connection in link with neuronal excitable dynamics in a PWL version of the classical Morris-Lecar (ML) model, that we will refer to as the PWL-ML model. The ML system is a two-dimensional excitable biophysical model characterized mainly by a cubic shaped nullcline for the membrane potential variable and a
by the homoclinic connection at $\varepsilon = 0$, perturb finitely many intersections between $W^s(S_\varepsilon^{sd})$ and $W^u(S_\varepsilon^{sd})$ for $0 < \varepsilon \ll 1$.

<figure><img src="image_3.png" /><figcaption>**Fig. 2 3-dimensional representation of the homoclinic connection.** The homoclinic connection at $z = z_0$ corresponds to the intersection between the two-dimensional stable and unstable manifolds $W^{s,u}(S_0^{sd})$ of the saddle-type subset $S_0^{sd}$ of the critical manifold.</figcaption></figure>

Understanding the slow-passage through the homoclinic connection may help to describe the phenomenon of spike adding, which consists in the addition of new spikes during the bursting phase as a system parameter varies; see [14, 34, 35] for more details. Such phenomena have been widely studied in the smooth framework through phenomenological neural models, such as the Hindmarsh-Rose (HR) model [14, 36, 37], and biophysical ones, such as the Morris-Lecar (ML) model [32, 38, 39]. However, the slow-passage through the homoclinic connection has not been formally studied in piecewise linear (PWL) systems.

PWL models are quite known to reproduce qualitative and quantitative aspects of smooth systems, being more amenable to general analysis [40–42]. In recent years, they have been extended in order to incorporate slow-fast dynamics, and it has been shown in several studies that they can reproduce the most salient features of multiple-timescale systems, providing a new point of view to understand better these complex dynamics [10, 14, 43–46]. Canard-induced spike-adding transitions associated with square-wave bursting have been studied in [44] in a PWL version of the HR model. More recently, a slow passage through a Hopf-like bifurcation has been rigorously studied in [47], with application to elliptic bursting.

In the present work, we focus on the slow passage through a homoclinic connection in link with neuronal excitable dynamics in a PWL version of the classical Morris-Lecar (ML) model, that we will refer to as the PWL-ML model. The ML system is a two-dimensional excitable biophysical model characterized mainly by a cubic shaped nullcline for the membrane potential variable and a
sigmoid shaped nullcline for the gating variable associated with the potassium ion channels. It exhibits a large spectrum of bifurcations, when considering the external current $I$ as the bifurcation parameter.

Here we define a PWL-ML differential system presenting nullclines emulating those of the original Morris-Lecar model. This allows PWL-ML to exhibit a bifurcation set similar to that of the original ML. Previous PWL versions of the ML do not present this characteristic. In particular, in [48] authors consider discontinuous nullclines what avoid some equilibrium bifurcations. Furthermore, in [49?], authors consider continuous nullclines, but without taking into account the saturated parts of the sigmoid. This makes difficult the existence of homoclinic and heteroclinic connections and the control of their bifurcations. Additionally, authors also describe the fold in the cubic shaped nullcline by a corner, what does not provide a proper canard explosion, see [50–52] for more details.

Following similar arguments than those in [12], we redefine the bifurcation parameter $I$ as a slow dynamic dependent variable, such that $I$ decreases when the membrane potential is in the oscillatory regime and $I$ increases when the membrane is resting. We show that this extended system, formed by adding the equation for the dynamics of $I$ to the PWL-ML model, is a simplified PWL burster.

The manuscript is organized as follows. In Section 2, we propose and describe the PWL-ML model, its nullclines configuration and its main bifurcation sets (Theorems 1 and 2). By adding a slow dynamics on the parameter related to the external current $I$, we build in Section 3 a 3-dimensional model exhibiting bursting oscillatory regimes, and we relate some bursting regime to the presence of a homoclinic connection in the fast subsystem, namely, the original PWL-ML model. This motivates the study of the slow-passage through a homoclinic loop, which we undertake in Section 4. Following that, we give new results related to the existence of the slow passage through a homoclinic connection in the PWL context, summarized in Theorems 3 and 4. Section 5 provides detailed proofs of all theorems. We emphasize in the advantages of the PWL context, that allows us to work with some analytical expressions, such as the slow manifolds and the local expression of the orbits, and also the use of some unusual techniques in the PWL context, in particular the contraction and expansion technique, which we introduce here.

## 2 Bifurcation structure of the PWL-ML model

We consider the following piecewise linear caricature of the Morris-Lecar (PWL-ML) model:

$$ \begin{cases} \dot{\varepsilon x} = f(x) - y + I, \\ \dot{y} = g(x) - y, \end{cases} \qquad (2) $$
sigmoid shaped nullcline for the gating variable associated with the potassium ion channels. It exhibits a large spectrum of bifurcations, when considering the external current *I* as the bifurcation parameter.

Here we define a PWL-ML differential system presenting nullclines emulating those of the original Morris-Lecar model. This allows PWL-ML to exhibit a bifurcation set similar to that of the original ML. Previous PWL versions of the ML do not present this characteristic. In particular, in [48] authors consider discontinuous nullclines what avoid some equilibrium bifurcations. Furthermore, in [49?], authors consider continuous nullclines, but without taking into account the saturated parts of the sigmoid. This makes difficult the existence of homoclinic and heteroclinic connections and the control of their bifurcations. Additionally, authors also describe the fold in the cubic shaped nullcline by a corner, what does not provide a proper canard explosion, see [50–52] for more details.

Following similar arguments than those in [12], we redefine the bifurcation parameter *I* as a slow dynamic dependent variable, such that *I* decreases when the membrane potential is in the oscillatory regime and *I* increases when the membrane is resting. We show that this extended system, formed by adding the equation for the dynamics of *I* to the PWL-ML model, is a simplified PWL burster.

The manuscript is organized as follows. In Section 2, we propose and describe the PWL-ML model, its nullclines configuration and its main bifurcation sets (Theorems 1 and 2). By adding a slow dynamics on the parameter related to the external current *I*, we build in Section 3 a 3-dimensional model exhibiting bursting oscillatory regimes, and we relate some bursting regime to the presence of a homoclinic connection in the fast subsystem, namely, the original PWL-ML model. This motivates the study of the slow-passage through a homoclinic loop, which we undertake in Section 4. Following that, we give new results related to the existence of the slow passage through a homoclinic connection in the PWL context, summarized in Theorems 3 and 4. Section 5 provides detailed proofs of all theorems. We emphasize in the advantages of the PWL context, that allows us to work with some analytical expressions, such as the slow manifolds and the local expression of the orbits, and also the use of some unusual techniques in the PWL context, in particular the contraction and expansion technique, which we introduce here.

## 2 Bifurcation structure of the PWL-ML model

We consider the following piecewise linear caricature of the Morris-Lecar (PWL-ML) model:

$$ \begin{cases} \dot{\varepsilon x} = f(x) - y + I, \\ \dot{y} = g(x) - y, \end{cases} \qquad (2) $$
where

$$f(x) = \begin{cases} -x & \text{if } x < -\sqrt{\varepsilon}, \\ \delta x + \lambda & \text{if } |x| \le \sqrt{\varepsilon}, \\ kx + \beta & \text{if } \sqrt{\varepsilon} < x < 1, \\ -x + \gamma & \text{if } x \ge 1, \end{cases} \quad g(x) = \begin{cases} 0 & \text{if } x < a, \\ lx + n & \text{if } a \le x < 1, \\ 1 & \text{if } x \ge 1. \end{cases} \quad (3)$$

Additionally, we consider {$\delta, k, \varepsilon, a, I$} to be the free parameters, being $\varepsilon > 0$, $a \le 1$, $k > 0$, and $\delta < 0$. Moreover, we consider that the rest of parameters, {$\lambda, \beta, \gamma, l, n$}, satisfy the relations

$$\lambda = \sqrt{\varepsilon}(1 + \delta), \quad \beta = \sqrt{\varepsilon}(1 - k + 2\delta), \quad \gamma = k + 1 + \sqrt{\varepsilon}(1 - k + 2\delta),$$

and for $a < 1$,

$$l = \frac{1}{1-a}, \quad n = \frac{-a}{1-a},$$

to ensure that $f$ and $g$ are continuous almost everywhere (a discontinuity arises when $a = 1$, where the central region of $g$ collapses). There exist some attempts in order to mimic the Morris-Lecar model, for instance [48, 49]. Nevertheless, these models are either discontinuous models or the minimal continuous set-up is improvable. For instance, in our model, the nullclines preserve the shape of those from the original ML model. That is, the $x$-nullcline {$y = f(x) + I$} mimics a cubic curve while the $y$-nullcline {$y = g(x)$} mimics a sigmoid. In case of [49], the author defines the $x$-nullcline with only three linearity regions, whilst we define the $x$-nullcline so that it contains a central segment in the zone {$|x| \le \sqrt{\varepsilon}$}, in order to suitably approximate the fold of the original cubic nullcline. The size of this central region is taken to be $O(\sqrt{\varepsilon})$ to be able to display a canard explosion; see [50–52]. Moreover, we describe the $y$-nullcline through three linearity regions, with the outsider ones as constants, which suits geometrically the shape of the sigmoid.

The choice of the other parameters follows from the additional dynamical considerations. First, to ensure that the bifurcation is supercritical, we take $\delta < 0$ (see [51]). Moreover, parameter $a$ will vary in the interval $[\sqrt{\varepsilon}, 1]$ so as to provide the high-threshold configuration of the ML model, which corresponds to a scenario where the system passes, as we change $a$, from a homoclinic bifurcation to a saddle point [16]. When $a = 1$, the vector field presents a discontinuity and the passage from $a = 1$ to $a < 1$ can be understood as a regulation of the discontinuous piecewise-linear vector field, which allows us to extend the existence of global dynamical objects from the discontinuous context, $a = 1$, to the continuous one, $a < 1$. Finally, we consider parameter $k$ in $(0, 1)$. This parameter can be considered as the ratio between the outer and the inner branches of the $x$-nullcline. Since $1 < l$ (given the interval that $a$ belongs to), then $k < l$. Consequently, in this particular region the slope of the $x$-nullcline is smaller than the slope of the $y$-nullcline, providing the existence of a suitable saddle-node bifurcation of equilibria at $x = a$; see Figure 3. In
where

$$f(x) = \begin{cases} -x & \text{if } x < -\sqrt{\varepsilon}, \\ \delta x + \lambda & \text{if } |x| \le \sqrt{\varepsilon}, \\ kx + \beta & \text{if } \sqrt{\varepsilon} < x < 1, \\ -x + \gamma & \text{if } x \ge 1, \end{cases} \quad g(x) = \begin{cases} 0 & \text{if } x < a, \\ lx + n & \text{if } a \le x < 1, \\ 1 & \text{if } x \ge 1. \end{cases} \quad (3)$$

Additionally, we consider {$\delta, k, \varepsilon, a, I$} to be the free parameters, being $\varepsilon > 0$, $a \le 1$, $k > 0$, and $\delta < 0$. Moreover, we consider that the rest of parameters, {$\lambda, \beta, \gamma, l, n$}, satisfy the relations

$$\lambda = \sqrt{\varepsilon}(1 + \delta), \quad \beta = \sqrt{\varepsilon}(1 - k + 2\delta), \quad \gamma = k + 1 + \sqrt{\varepsilon}(1 - k + 2\delta),$$

and for $a < 1$,

$$l = \frac{1}{1-a}, \quad n = \frac{-a}{1-a},$$

to ensure that $f$ and $g$ are continuous almost everywhere (a discontinuity arises when $a = 1$, where the central region of $g$ collapses). There exist some attempts in order to mimic the Morris-Lecar model, for instance [48, 49]. Nevertheless, these models are either discontinuous models or the minimal continuous set-up is improvable. For instance, in our model, the nullclines preserve the shape of those from the original ML model. That is, the $x$-nullcline {$y = f(x) + I$} mimics a cubic curve while the $y$-nullcline {$y = g(x)$} mimics a sigmoid. In case of [49], the author defines the $x$-nullcline with only three linearity regions, whilst we define the $x$-nullcline so that it contains a central segment in the zone {$|x| \le \sqrt{\varepsilon}$}, in order to suitably approximate the fold of the original cubic nullcline. The size of this central region is taken to be $O(\sqrt{\varepsilon})$ to be able to display a canard explosion; see [50–52]. Moreover, we describe the $y$-nullcline through three linearity regions, with the outsider ones as constants, which suits geometrically the shape of the sigmoid.

The choice of the other parameters follows from the additional dynamical considerations. First, to ensure that the bifurcation is supercritical, we take $\delta < 0$ (see [51]). Moreover, parameter $a$ will vary in the interval $[\sqrt{\varepsilon}, 1]$ so as to provide the high-threshold configuration of the ML model, which corresponds to a scenario where the system passes, as we change $a$, from a homoclinic bifurcation to a saddle point [16]. When $a = 1$, the vector field presents a discontinuity and the passage from $a = 1$ to $a < 1$ can be understood as a regulation of the discontinuous piecewise-linear vector field, which allows us to extend the existence of global dynamical objects from the discontinuous context, $a = 1$, to the continuous one, $a < 1$. Finally, we consider parameter $k$ in $(0, 1)$. This parameter can be considered as the ratio between the outer and the inner branches of the $x$-nullcline. Since $1 < l$ (given the interval that $a$ belongs to), then $k < l$. Consequently, in this particular region the slope of the $x$-nullcline is smaller than the slope of the $y$-nullcline, providing the existence of a suitable saddle-node bifurcation of equilibria at $x = a$; see Figure 3. In
this scenario, we can obtain bursting oscillations by using the PWL-ML model as a fast subsystem [16].

The piecewise linear nature of the vector field defined by equation (2) splits the phase space into different linearity regions. For clarity, we denote each of the linearity regions by $R_{i,j}$ such that, from left to right, $i \in \{1, 2, 3, 4\}$ stands for the $i$-th region of the graph of $f$, while $j \in \{1, 2, 3\}$ stands for the $j$-th region where the graph of $g$ is defined. When $a \in (\sqrt{\epsilon}, 1)$, regions $R_{1,1}$, $R_{2,1}$, $R_{3,1}$, $R_{3,2}$ and $R_{4,3}$ exist (see Figure 3). However, when $a = \sqrt{\epsilon}$, region $R_{3,1}$ vanishes while for $a = 1$, region $R_{3,2}$ is the one vanishing (see Figure 12). Local expressions of the flow used along the manuscript are given in Appendix A.

<figure><img src="image_3.png" /><figcaption>Fig. 3 Configuration of the nullclines when $\sqrt{\epsilon} < a < 1$. The green curve represents the x-nullcline $\{y = f(x) + I\}$, whose shape mimics a cubic curve; the blue curve is the y-nullcline $\{y = g(x)\}$, which mimics a sigmoid. Vertical dashed lines depict the different linearity regions.</figcaption></figure>

Consider the following 3-dimensional subset of the parameter space

$$ \mathcal{P} = \{(a, \epsilon, I) \in \mathbb{R}^3 : 0 < \epsilon < k, \sqrt{\epsilon} \le a \le 1\}. \quad (4) $$

Next, we study the different bifurcation sets that partition region $\mathcal{P}$. Those sets corresponding to saddle-node equilibria are provided in Theorem 1, whereas the bifurcation set corresponding to a global connection is provided in Theorem 2. We now define the following partition of $\mathcal{P}$:

$$ \begin{aligned} \mathcal{SN}_1 &= \{(a, \epsilon, I) \in \mathcal{P} : \sqrt{\epsilon} < a < 1, I = -\sqrt{\epsilon}(1 + 2\delta)\}, \\ \mathcal{SN}_2 &= \{(a, \epsilon, I) \in \mathcal{P} : \sqrt{\epsilon} < a < 1, I = -(ak + \beta)\}, \\ \mathcal{C} &= \{(a, \epsilon, I) \in \mathcal{P} : a = \sqrt{\epsilon}, I = -(ak + \beta)\}, \\ \mathcal{E}_3 &= \{(a, \epsilon, I) \in \mathcal{P} : \sqrt{\epsilon} < a < 1, -(ak + \beta) < I < -\sqrt{\epsilon}(1 + 2\delta)\}, \\ \mathcal{E}_1 &= \mathcal{P} \setminus (\mathcal{SN}_1 \cup \mathcal{E}_3 \cup \mathcal{SN}_2 \cup \mathcal{C}). \end{aligned} \quad (5) $$
this scenario, we can obtain bursting oscillations by using the PWL-ML model as a fast subsystem [16].

The piecewise linear nature of the vector field defined by equation (2) splits the phase space into different linearity regions. For clarity, we denote each of the linearity regions by $R_{i,j}$ such that, from left to right, $i \in \{1, 2, 3, 4\}$ stands for the $i$-th region of the graph of $f$, while $j \in \{1, 2, 3\}$ stands for the $j$-th region where the graph of $g$ is defined. When $a \in (\sqrt{\epsilon}, 1)$, regions $R_{1,1}$, $R_{2,1}$, $R_{3,1}$, $R_{3,2}$ and $R_{4,3}$ exist (see Figure 3). However, when $a = \sqrt{\epsilon}$, region $R_{3,1}$ vanishes while for $a = 1$, region $R_{3,2}$ is the one vanishing (see Figure 12). Local expressions of the flow used along the manuscript are given in Appendix A.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Configuration of the nullclines when $\sqrt{\epsilon} < a < 1$. The green curve represents the x-nullcline $\{y = f(x) + I\}$, whose shape mimics a cubic curve; the blue curve is the y-nullcline $\{y = g(x)\}$, which mimics a sigmoid. Vertical dashed lines depict the different linearity regions.</figcaption></figure>

Consider the following 3-dimensional subset of the parameter space

$$ \mathcal{P} = \{(a, \epsilon, I) \in \mathbb{R}^3 : 0 < \epsilon < k, \sqrt{\epsilon} \le a \le 1\}. \quad (4) $$

Next, we study the different bifurcation sets that partition region $\mathcal{P}$. Those sets corresponding to saddle-node equilibria are provided in Theorem 1, whereas the bifurcation set corresponding to a global connection is provided in Theorem 2. We now define the following partition of $\mathcal{P}$:

$$ \begin{aligned} \mathcal{S}\mathcal{N}_1 &= \{(a, \epsilon, I) \in \mathcal{P} : \sqrt{\epsilon} < a < 1, I = -\sqrt{\epsilon}(1 + 2\delta)\}, \\ \mathcal{S}\mathcal{N}_2 &= \{(a, \epsilon, I) \in \mathcal{P} : \sqrt{\epsilon} < a < 1, I = -(ak + \beta)\}, \\ \mathcal{C} &= \{(a, \epsilon, I) \in \mathcal{P} : a = \sqrt{\epsilon}, I = -(ak + \beta)\}, \\ \mathcal{E}_3 &= \{(a, \epsilon, I) \in \mathcal{P} : \sqrt{\epsilon} < a < 1, -(ak + \beta) < I < -\sqrt{\epsilon}(1 + 2\delta)\}, \\ \mathcal{E}_1 &= \mathcal{P} \setminus (\mathcal{S}\mathcal{N}_1 \cup \mathcal{E}_3 \cup \mathcal{S}\mathcal{N}_2 \cup \mathcal{C}). \end{aligned} \quad (5) $$
**Theorem 1** Let us consider system (2) and $(a, \varepsilon, I) \in \mathcal{P}$.

a) If $(a, \varepsilon, I) \in \mathcal{E}_3$, system (2) exhibits 3 equilibria: an attracting node, a saddle and the third one either a node or a focus.

b) If $(a, \varepsilon, I) \in \mathcal{E}_1$, system (2) exhibits one equilibrium point.

c) If $(a, \varepsilon, I) \in \mathcal{SN}_1$, system (2) exhibits two equilibria: a saddle-node and either a focus or a node. Therefore, $\mathcal{SN}_1$ is a saddle-node bifurcation manifold, being $I$ the bifurcation parameter.

d) If $(a, \varepsilon, I) \in \mathcal{SN}_2$, system (2) exhibits two equilibria: an attracting node and a saddle-node equilibrium. Therefore, $\mathcal{SN}_2$ is a saddle-node bifurcation manifold, being $I$ the bifurcation parameter.

e) If $(a, \varepsilon, I) \in \mathcal{C}$, system (2) exhibits one equilibrium point. Furthermore, $\mathcal{C}$ is a codimension 2 bifurcation curve, being $I$ and either a or $\varepsilon$ the bifurcation parameters.

The proof of Theorem 1 is given in Section 5.1.

Following Theorem 1, as we vary parameter $I$ in a neighbourhood of the curve $\mathcal{C}$, system (2) continuously change to exhibit one equilibrium point in $\mathcal{E}_1$, two equilibria on $\mathcal{SN}_1$, three equilibria in $\mathcal{E}_3$, two equilibria on $\mathcal{SN}_2$, and finally one equilibrium point when back again to $\mathcal{E}_3$. Moreover, the saddle-node surfaces $\mathcal{SN}_1$ and $\mathcal{SN}_2$ get in contact along the curve $\mathcal{C}$. Although this contact is not tangential, the scenario resembles that described by the cusp bifurcation in terms of the numbers of equilibria involved [53, 54], so we denote the curve $\mathcal{C}$ as a nonsmooth cusp bifurcation curve.

In the next result we assure the existence of a homoclinic surface rising from the parameter plane $\mathcal{P} \cap \{a = 1\}$.

**Theorem 2** There exists $\varepsilon_0 \in (0, k)$ and differentiable functions $a(\varepsilon)$ and $I(a, \varepsilon)$ defined in $(0, \varepsilon_0)$ and $U := \{(a, \varepsilon) : 0 < \varepsilon < \varepsilon_0, a(\varepsilon) \le a \le 1\}$ respectively, such that the manifold $\mathcal{H}_0 := \{(a, \varepsilon, I(a, \varepsilon)) : (a, \varepsilon) \in U\}$ is contained in $\mathcal{E}_3$ and, for parameters $(a, \varepsilon, I) \in \mathcal{H}_0$, system (2) exhibits a homoclinic connection to a saddle point. Additionally, for parameters $(a, \varepsilon, I)$ with $a < 1$, $(a, \varepsilon) \in U$ and $I > I(a, \varepsilon)$, a stable limit cycle exists and attracts the unstable manifold of the saddle.

The proof of Theorem 2 is given in Section 5.2. In Figure 4(a) we represent the bifurcation set of system (2) inside the parameter region $\mathcal{P}$. In particular, the saddle-node bifurcation surfaces $\mathcal{SN}_1$ (green surface) and $\mathcal{SN}_2$ (blue surface); the nonsmooth cusp bifurcation curve $\mathcal{C}$ at the intersection of the saddle-node bifurcations manifolds; the homoclinic connection surface $\mathcal{H}_0$, although we specifically represent the intersection of $\mathcal{H}_0$ with the plane $\{a = 1\}$ corresponding to the red curve, and the intersection of $\mathcal{H}_0$ with $\mathcal{SN}_1$ corresponding to the pink curve. Even though these curves seem to intersect, in the represented case they do not. The surface $\mathcal{H}_0$ is omitted in panel (a) for a better understanding of the image.

In Figure 4(b), we draw the intersection of the figure in panel (a) with the plane $\{\varepsilon = 1/3\}$. In particular, we draw the surface $\mathcal{H}_0$ given now by the red
**Theorem 1** Let us consider system (2) and $(a, \varepsilon, I) \in \mathcal{P}$.

a) If $(a, \varepsilon, I) \in \mathcal{E}_3$, system (2) exhibits 3 equilibria: an attracting node, a saddle and the third one either a node or a focus.

b) If $(a, \varepsilon, I) \in \mathcal{E}_1$, system (2) exhibits one equilibrium point.

c) If $(a, \varepsilon, I) \in \mathcal{SN}_1$, system (2) exhibits two equilibria: a saddle-node and either a focus or a node. Therefore, $\mathcal{SN}_1$ is a saddle-node bifurcation manifold, being $I$ the bifurcation parameter.

d) If $(a, \varepsilon, I) \in \mathcal{SN}_2$, system (2) exhibits two equilibria: an attracting node and a saddle-node equilibrium. Therefore, $\mathcal{SN}_2$ is a saddle-node bifurcation manifold, being $I$ the bifurcation parameter.

e) If $(a, \varepsilon, I) \in \mathcal{C}$, system (2) exhibits one equilibrium point. Furthermore, $\mathcal{C}$ is a codimension 2 bifurcation curve, being $I$ and either a or $\varepsilon$ the bifurcation parameters.

The proof of Theorem 1 is given in Section 5.1.

Following Theorem 1, as we vary parameter $I$ in a neighbourhood of the curve $\mathcal{C}$, system (2) continuously change to exhibit one equilibrium point in $\mathcal{E}_1$, two equilibria on $\mathcal{SN}_1$, three equilibria in $\mathcal{E}_3$, two equilibria on $\mathcal{SN}_2$, and finally one equilibrium point when back again to $\mathcal{E}_3$. Moreover, the saddle-node surfaces $\mathcal{SN}_1$ and $\mathcal{SN}_2$ get in contact along the curve $\mathcal{C}$. Although this contact is not tangential, the scenario resembles that described by the cusp bifurcation in terms of the numbers of equilibria involved [53, 54], so we denote the curve $\mathcal{C}$ as a nonsmooth cusp bifurcation curve.

In the next result we assure the existence of a homoclinic surface rising from the parameter plane $\mathcal{P} \cap \{a = 1\}$.

**Theorem 2** There exists $\varepsilon_0 \in (0, k)$ and differentiable functions $a(\varepsilon)$ and $I(a, \varepsilon)$ defined in $(0, \varepsilon_0)$ and $U := \{(a, \varepsilon) : 0 < \varepsilon < \varepsilon_0, a(\varepsilon) \le a \le 1\}$ respectively, such that the manifold $\mathcal{H}_0 := \{(a, \varepsilon, I(a, \varepsilon)) : (a, \varepsilon) \in U\}$ is contained in $\mathcal{E}_3$ and, for parameters $(a, \varepsilon, I) \in \mathcal{H}_0$, system (2) exhibits a homoclinic connection to a saddle point. Additionally, for parameters $(a, \varepsilon, I)$ with $a < 1$, $(a, \varepsilon) \in U$ and $I > I(a, \varepsilon)$, a stable limit cycle exists and attracts the unstable manifold of the saddle.

The proof of Theorem 2 is given in Section 5.2. In Figure 4(a) we represent the bifurcation set of system (2) inside the parameter region $\mathcal{P}$. In particular, the saddle-node bifurcation surfaces $\mathcal{SN}_1$ (green surface) and $\mathcal{SN}_2$ (blue surface); the nonsmooth cusp bifurcation curve $\mathcal{C}$ at the intersection of the saddle-node bifurcations manifolds; the homoclinic connection surface $\mathcal{H}_0$, although we specifically represent the intersection of $\mathcal{H}_0$ with the plane $\{a = 1\}$ corresponding to the red curve, and the intersection of $\mathcal{H}_0$ with $\mathcal{SN}_1$ corresponding to the pink curve. Even though these curves seem to intersect, in the represented case they do not. The surface $\mathcal{H}_0$ is omitted in panel (a) for a better understanding of the image.

In Figure 4(b), we draw the intersection of the figure in panel (a) with the plane $\{\varepsilon = 1/3\}$. In particular, we draw the surface $\mathcal{H}_0$ given now by the red
curve. This curve starts at the plane {$a = 1$}, corresponding with a homoclinic bifurcation in the discontinuous configuration, and ends at the bifurcation surface $SN_1$, corresponding to a homoclinic bifurcation to a saddle-node equilibrium also denoted a codimension-2 saddle-node homoclinic bifurcation. The compact region located between the red curve, the green segment and the black segment, corresponds to a bistability region where one stable node and one stable limit cycle coexist. Moreover, in region $\mathcal{E}_3$ three equilibria coexist before disappearing at two saddle-node bifurcations. Thus, by appropriately varying parameter $I$, a hysteresis loop can be obtained. Both the stability region of the PWL-ML model and the hysteresis loop are key for bursting dynamics to emerge, see [13, 15, 44].

<figure><img src="image_4.png" /><figcaption>Fig. 4 Bifurcation set of system (2) in the parameter region $\mathcal{P}$. (a) Saddle-node bifurcation surfaces $SN_1$ (green) and $SN_2$ (blue), nonsmooth cusp bifurcation curve $C$, and intersection of the homoclinic bifurcation surface $\mathcal{H}_0$ with the plane {$a = 1$} (red curve) and with the saddle-node bifurcation surface $SN_1$ (pink curve). (b) Intersection of the manifolds in (a) with the plane {$\varepsilon = 1/3$}. Black line corresponds with the plane {$a = 1$}. The red curve is $\mathcal{H}_0 \cap \{\varepsilon = 1/3\}$ corresponding with the homoclinic bifurcation. Finally, the intersection of the red curve with the green curve corresponds with a codimension-2 saddle-node homoclinic bifurcation. The pictures are obtained by considering $k = 0.52$, $\delta = -0.1$.</figcaption></figure>

## 3 Bursting in a 3D extension of the model

In this section, we extend the PWL-ML model studied so far by adding a slow dynamics on parameter $I$, so as to create precursors of bursting dynamics in the resulting 3D extended model. Namely, we consider slow passages through fast-subsystem bifurcations in order to understand key slow-fast transitions pertinent for bursting oscillations. Minimally, bursting dynamics require two fast and one slow variables [15] and this is the framework that we consider here. As a result, we define the following 3-dimensional extension of the PWL-ML
curve. This curve starts at the plane {$a = 1$}, corresponding with a homoclinic bifurcation in the discontinuous configuration, and ends at the bifurcation surface $SN_1$, corresponding to a homoclinic bifurcation to a saddle-node equilibrium also denoted a codimension-2 saddle-node homoclinic bifurcation. The compact region located between the red curve, the green segment and the black segment, corresponds to a bistability region where one stable node and one stable limit cycle coexist. Moreover, in region $\mathcal{E}_3$ three equilibria coexist before disappearing at two saddle-node bifurcations. Thus, by appropriately varying parameter $I$, a hysteresis loop can be obtained. Both the stability region of the PWL-ML model and the hysteresis loop are key for bursting dynamics to emerge, see [13, 15, 44].

<figure><img src="image_5.png" /><figcaption>Fig. 4 Bifurcation set of system (2) in the parameter region $\mathcal{P}$. (a) Saddle-node bifurcation surfaces $SN_1$ (green) and $SN_2$ (blue), nonsmooth cusp bifurcation curve $C$, and intersection of the homoclinic bifurcation surface $\mathcal{H}_0$ with the plane {$a = 1$} (red curve) and with the saddle-node bifurcation surface $SN_1$ (pink curve). (b) Intersection of the manifolds in (a) with the plane {$\varepsilon = 1/3$}. Black line corresponds with the plane {$a = 1$}. The red curve is $\mathcal{H}_0 \cap \{\varepsilon = 1/3\}$ corresponding with the homoclinic bifurcation. Finally, the intersection of the red curve with the green curve corresponds with a codimension-2 saddle-node homoclinic bifurcation. The pictures are obtained by considering $k = 0.52$, $\delta = -0.1$.</figcaption></figure>

## 3 Bursting in a 3D extension of the model

In this section, we extend the PWL-ML model studied so far by adding a slow dynamics on parameter $I$, so as to create precursors of bursting dynamics in the resulting 3D extended model. Namely, we consider slow passages through fast-subsystem bifurcations in order to understand key slow-fast transitions pertinent for bursting oscillations. Minimally, bursting dynamics require two fast and one slow variables [15] and this is the framework that we consider here. As a result, we define the following 3-dimensional extension of the PWL-ML
model

$$ egin{cases} \epsilon \dot{x} = f(x) - y + I, \\ \dot{y} = g(x) - y, \\ \dot{I} = \sigma(x^* - x). \end{cases} \tag{6} $$

In this case, the parameter $I$ becomes the slow variable, but with a feedback onto the $x$ variable, which makes it oscillate in the zone of bistability of the fast subsystem, hence allowing for robust bursting dynamics. Given the results obtained in Section 2, namely the existence of a saddle-node bifurcation of equilibria and of a saddle homoclinic connection in the PWL-ML model, for certain values of $I$, and given that this PWL-ML is effectively the fast subsystem of system (6), then we expect that the 3D extension supports solutions akin to square-wave bursting oscillations, and we showcase in Fig. 5 limit cycles that are resembling of such bursting oscillations and for which the number of spikes per burst varies brutally around very specific values of the control parameter $x$. This provides a good motivation to focus on the spike-adding process, which is closely related to the presence of a homoclinic connection in the fast subsystem, as studied in the previous section.

To illustrate this fact, we take the following parameter values: $\epsilon = 1/3$, $a = 0.99$, $\delta = -0.1$, $k = 0.52$, $\sigma = 0.01$, and we vary $x^*$. We then observe by direct simulation bursting solutions with different number of spikes per burst, sensitively depending upon the value of $x^*$. Figure 5 shows the results of such simulations. Indeed, as we increase parameter $x^*$, the number of spikes increases and the oscillatory state starts at a higher value of $I$.

<figure><img src="image_5.png" /><figcaption>Fig. 5 Spike adding phenomenon in the 3D PWL-ML continuous model (6). Labels 1-4 correspond to: (a) The transition between 1 to 2 spikes $x^* = 0.577727969500$, $x^* = 0.577727969501$; (b) The transition between 2 to 3 spikes $x^* = 0.591899523320$, $x^* = 0.591899523321$, respectively. The black lines correspond to the critical manifold, namely, the equilibria of the fast subsystem. Parameter values are: $\epsilon = 1/3$, $a = 0.99$, $\delta = -0.1$, $k = 0.52$, $\sigma = 0.01$.</figcaption></figure>

The behaviour of increasing the number of spikes as a system parameter is varied is common in many bursting oscillations and it is called *spike adding*. In the case of square-wave bursting oscillations, a key element in the spike adding phenomenon is the slow passage through the homoclinic connection exhibited by the fast subsystem of (6). In Section 4, we realize an analytical
model

$$ \begin{cases} \varepsilon \dot{x} = f(x) - y + I, \\ \dot{y} = g(x) - y, \\ \dot{I} = \sigma(x^* - x). \end{cases} \qquad (6) $$

In this case, the parameter $I$ becomes the slow variable, but with a feedback onto the $x$ variable, which makes it oscillate in the zone of bistability of the fast subsystem, hence allowing for robust bursting dynamics. Given the results obtained in Section 2, namely the existence of a saddle-node bifurcation of equilibria and of a saddle homoclinic connection in the PWL-ML model, for certain values of $I$, and given that this PWL-ML is effectively the fast subsystem of system (6), then we expect that the 3D extension supports solutions akin to square-wave bursting oscillations, and we showcase in Fig. 5 limit cycles that are resembling of such bursting oscillations and for which the number of spikes per burst varies brutally around very specific values of the control parameter $x$. This provides a good motivation to focus on the spike-adding process, which is closely related to the presence of a homoclinic connection in the fast subsystem, as studied in the previous section.

To illustrate this fact, we take the following parameter values: $\varepsilon = 1/3$, $a = 0.99$, $\delta = -0.1$, $k = 0.52$, $\sigma = 0.01$, and we vary $x^*$. We then observe by direct simulation bursting solutions with different number of spikes per burst, sensitively depending upon the value of $x^*$. Figure 5 shows the results of such simulations. Indeed, as we increase parameter $x^*$, the number of spikes increases and the oscillatory state starts at a higher value of $I$.

<figure><img src="image_6.png" /><figcaption>Fig. 5 Spike adding phenomenon in the 3D PWL-ML continuous model (6). Labels 1-4 correspond to: (a) The transition between 1 to 2 spikes $x^* = 0.577727969500$, $x^* = 0.577727969501$; (b) The transition between 2 to 3 spikes $x^* = 0.591899523320$, $x^* = 0.591899523321$, respectively. The black lines correspond to the critical manifold, namely, the equilibria of the fast subsystem. Parameter values are: $\varepsilon = 1/3$, $a = 0.99$, $\delta = -0.1$, $k = 0.52$, $\sigma = 0.01$.</figcaption></figure>

The behaviour of increasing the number of spikes as a system parameter is varied is common in many bursting oscillations and it is called *spike adding*. In the case of square-wave bursting oscillations, a key element in the spike adding phenomenon is the slow passage through the homoclinic connection exhibited by the fast subsystem of (6). In Section 4, we realize an analytical
approach of the slow passage phenomenon by considering I a slow drift; this is a typical treatment given in order to simplify the system and provide a slow-fast approach, see for instance [24].

## 4 Slow passage through the homoclinic loop

In order to study the slow passage through the homoclinic connection, we consider now system (6) by fixing $0 < \varepsilon < k$ and assuming slow constant dynamics to parameter I such that the system is now provided by

$$ \begin{cases} \varepsilon \dot{x} = f(x) - y + I, \\ \dot{y} = g(x) - y, \\ \dot{I} = \sigma. \end{cases} \qquad (7) $$

where $0 \le |\sigma| \ll 1$. Hence, I works as the slow variable and $x, y$ are the fast variables. Let $\varphi(t; \mathbf{p}, (a, \sigma))$ denote the flow of the system through a point $\mathbf{p}$ for the parameters $(a, \sigma)$. Similarly to (6), the fast subsystem (obtained when $\sigma = 0$) corresponds to the PWL-ML model given by system (2).

The critical manifold, formed by the equilibrium points of the fast subsystem, is given by the polygonal curve

$$ S_0 = \{(x, y, I) \in \mathbb{R}^3 : g(x) - f(x) = I, y = g(x), I \in \mathbb{R}\}. \qquad (8) $$

According to the character of the equilibrium points, the different branches of the critical manifold have different stability; in particular, the one defined in region $R_{3,1}$, where $f(x) = kx + \beta$ and $g(x) = 0$, and its expression is given by

$$ S_0^{sd} = \left\{ \left( -\frac{\beta + I}{k}, 0, I \right) \in \mathbb{R}^3 : I \in \left( -\beta - ka, -\beta - k\sqrt{\varepsilon} \right) \right\}, \qquad (9) $$

In this case, $S_0^{sd}$ is defined by the hyperbolic equilibrium points of saddle type of the fast subsystem of system (7). Let $W^s(S_0^{sd})$ and $W^u(S_0^{sd})$ be the stable and the unstable manifolds of $S_0^{sd}$ given by

$$ W^s(S_0^{sd}) = \left\{ \left( \frac{y}{k+\varepsilon} - \frac{I+\beta}{k}, y, I \right) \in \mathbb{R}^3 : \frac{k}{k+\varepsilon}y - (k\sqrt{\varepsilon} + \beta) < I < \frac{k}{k+\varepsilon}y - (ka + \beta) \right\}, \qquad (10) $$

$$ W^u(S_0^{sd}) = \{(x, 0, I) \in \mathbb{R}^3 : \sqrt{\varepsilon} < x < a\}. \qquad (11) $$

Following Fenichel's theorems, after perturbing the system with respect to $\sigma$, there exists a slow manifold $S_\sigma^{sd}$ at distance $O(\sigma)$ from the saddle
approach of the slow passage phenomenon by considering I a slow drift; this is a typical treatment given in order to simplify the system and provide a slow-fast approach, see for instance [24].

## 4 Slow passage through the homoclinic loop

In order to study the slow passage through the homoclinic connection, we consider now system (6) by fixing $0 < \varepsilon < k$ and assuming slow constant dynamics to parameter I such that the system is now provided by

$$ egin{cases} \varepsilon \dot{x} = f(x) - y + I, \\ \dot{y} = g(x) - y, \\ \dot{I} = \sigma. \end{cases} \tag{7} $$

where $0 \le |\sigma| \ll 1$. Hence, I works as the slow variable and $x, y$ are the fast variables. Let $\varphi(t; \mathbf{p}, (a, \sigma))$ denote the flow of the system through a point $\mathbf{p}$ for the parameters $(a, \sigma)$. Similarly to (6), the fast subsystem (obtained when $\sigma = 0$) corresponds to the PWL-ML model given by system (2).

The critical manifold, formed by the equilibrium points of the fast subsystem, is given by the polygonal curve

$$ S_0 = \{(x, y, I) \in \mathbb{R}^3 : g(x) - f(x) = I, y = g(x), I \in \mathbb{R}\}. \tag{8} $$

According to the character of the equilibrium points, the different branches of the critical manifold have different stability; in particular, the one defined in region $R_{3,1}$, where $f(x) = kx + \beta$ and $g(x) = 0$, and its expression is given by

$$ S_0^{sd} = \left\{ \left( -\frac{\beta + I}{k}, 0, I \right) \in \mathbb{R}^3 : I \in \left( -\beta - ka, -\beta - k\sqrt{\varepsilon} \right) \right\}, \tag{9} $$

In this case, $S_0^{sd}$ is defined by the hyperbolic equilibrium points of saddle type of the fast subsystem of system (7). Let $W^s(S_0^{sd})$ and $W^u(S_0^{sd})$ be the stable and the unstable manifolds of $S_0^{sd}$ given by

$$ W^s(S_0^{sd}) = \left\{ \left( \frac{y}{k+\varepsilon} - \frac{I+\beta}{k}, y, I \right) \in \mathbb{R}^3 : \frac{k}{k+\varepsilon}y - (k\sqrt{\varepsilon} + \beta) < I < \frac{k}{k+\varepsilon}y - (ka + \beta) \right\}, \tag{10} $$

$$ W^u(S_0^{sd}) = \{(x, 0, I) \in \mathbb{R}^3 : \sqrt{\varepsilon} < x < a\}. \tag{11} $$

Following Fenichel's theorems, after perturbing the system with respect to $\sigma$, there exists a slow manifold $S_\sigma^{sd}$ at distance $O(\sigma)$ from the saddle
type branch $S_0^{sd}$, having one stable manifold $W^s(S_0^{sd})$ at distance $O(\sigma)$ from $W^s(S_0^{sd})$ and one unstable manifold $W^u(S_0^{sd})$ at distance $O(\sigma)$ from $W^u(S_0^{sd})$.

From Lemma 7, we conclude that the slow manifold and both the stable and the unstable manifolds of the slow manifold are locally given, respectively, by the segment

$$S_{\sigma}^{sd} = \left\{ \left( \left( \frac{-\sigma}{k} t - \frac{\sigma\varepsilon + k(I_0 + \beta)}{k^2} \right), 0, \sigma t + I_0 \right), -\frac{k\sqrt{\varepsilon}}{\sigma} < t + C < -\frac{ka}{\sigma} \right\}, \quad (12)$$

where $C = \frac{\varepsilon}{k} + \frac{I_0 + \beta}{\sigma}$, and the pieces of plane

$$W^s(S_\sigma^{sd}) = \left\{ \left( \frac{y}{k+\varepsilon}, \frac{I}{k} - \frac{\sigma\varepsilon + k\beta}{k^2}, y, I \right) \in \mathbb{R}^3 : \begin{array}{l} \frac{k}{k+\varepsilon}y - (k\sqrt{\varepsilon} + \beta + \sigma\varepsilon/k) < I < \frac{k}{k+\varepsilon}y - (ka + \beta + \sigma\varepsilon/k) \end{array} \right\}, \quad (13)$$

$$W^u(S_\sigma^{sd}) = \{(x, 0, I) \in \mathbb{R}^3 : \sqrt{\varepsilon} < x < a\}. \quad (14)$$

As we can see from Equations (9) and (12), the critical manifold $S_0^{sd}$ and the slow manifold $S_\sigma^{sd}$ are parallel lines, separated by a factor $\sigma\varepsilon/k^2$, see Figure 6. The same applies between the stable planes arising from the critical and slow manifolds. However, the perturbation of the unstable plane of the critical manifold remains the same. Moreover, since compact subsets of planes $W^s(S_\sigma^{sd})$, $W^u(S_\sigma^{sd})$, and manifold $S_\sigma^{sd}$ are at Hausdorff distance $O(\sigma)$ from $W^s(S_0^{sd})$, $W^u(S_0^{sd})$, and $S_0^{sd}$, respectively, then the former expressions must be the canonical slow manifolds (see [55] for more details).

When $\sigma = 0$, the invariant manifolds $W^s(S_0^{sd})$ and $W^u(S_0^{sd})$ intersect exactly at the homoclinic loop, see Figure 2. Indeed, as it is proved in Proposition 1 and Lemma 5, there exists a unique orbit that, starting at the unstable manifold, achieves the stable one (see Figure 6). However, this is not the case when considering $0 < |\sigma| \ll 1$. Following [31], it can be expected that a finite number of new intersections between the manifolds $W^s(S_\sigma^{sd})$ and $W^u(S_\sigma^{sd})$ perturb from the initial intersection given by the homoclinic orbit. These intersections provide a finite number of orbits, called maximal canards [43], that starting at $W^u(S_\sigma^{sd})$, they reach $W^s(S_\sigma^{sd})$ and all these orbits perturb from the homoclinic orbit. Additionally, the number of maximal canards will depend on the singular parameter $\sigma$, let us denote this number by $N(\sigma)$. Finally, $N(\sigma)$ tipically exhibits a monotonically decreasing behavior with $\sigma$ and $\lim_{\sigma \to 0} N(\sigma) = \infty$.

In order to prove the existence of the connections between the stable and the unstable manifold, we proceed as follows. First, we study the existence of this kind of connections in the discontinuous case $a = 1$, where the situation is easy to analyze due to the fact that only two linearity regions are involved.
type branch $S_0^{sd}$, having one stable manifold $W^s(S_0^{sd})$ at distance $O(\sigma)$ from $W^s(S_0^{sd})$ and one unstable manifold $W^u(S_0^{sd})$ at distance $O(\sigma)$ from $W^u(S_0^{sd})$.

From Lemma 7, we conclude that the slow manifold and both the stable and the unstable manifolds of the slow manifold are locally given, respectively, by the segment

$$S_{\sigma}^{sd} = \left\{ \left( \left( \frac{-\sigma}{k} t - \frac{\sigma\varepsilon + k(I_0 + \beta)}{k^2} \right), 0, \sigma t + I_0 \right), -\frac{k\sqrt{\varepsilon}}{\sigma} < t + C < -\frac{ka}{\sigma} \right\}, \quad (12)$$

where $C = \frac{\varepsilon}{k} + \frac{I_0 + \beta}{\sigma}$, and the pieces of plane

$$W^s(S_{\sigma}^{sd}) = \left\{ \left( \frac{y}{k + \varepsilon}, \frac{I}{k} - \frac{\sigma\varepsilon + k\beta}{k^2}, y, I \right) \in \mathbb{R}^3 : \frac{k}{k + \varepsilon} y - (k\sqrt{\varepsilon} + \beta + \sigma\varepsilon/k) < I < \frac{k}{k + \varepsilon} y - (ka + \beta + \sigma\varepsilon/k) \right\}, \quad (13)$$

$$W^u(S_{\sigma}^{sd}) = \{(x, 0, I) \in \mathbb{R}^3 : \sqrt{\varepsilon} < x < a\}. \quad (14)$$

As we can see from Equations (9) and (12), the critical manifold $S_0^{sd}$ and the slow manifold $S_{\sigma}^{sd}$ are parallel lines, separated by a factor $\sigma\varepsilon/k^2$, see Figure 6. The same applies between the stable planes arising from the critical and slow manifolds. However, the perturbation of the unstable plane of the critical manifold remains the same. Moreover, since compact subsets of planes $W^s(S_{\sigma}^{sd})$, $W^u(S_{\sigma}^{sd})$, and manifold $S_{\sigma}^{sd}$ are at Hausdorff distance $O(\sigma)$ from $W^s(S_0^{sd})$, $W^u(S_0^{sd})$, and $S_0^{sd}$, respectively, then the former expressions must be the canonical slow manifolds (see [55] for more details).

When $\sigma = 0$, the invariant manifolds $W^s(S_0^{sd})$ and $W^u(S_0^{sd})$ intersect exactly at the homoclinic loop, see Figure 2. Indeed, as it is proved in Proposition 1 and Lemma 5, there exists a unique orbit that, starting at the unstable manifold, achieves the stable one (see Figure 6). However, this is not the case when considering $0 < |\sigma| \ll 1$. Following [31], it can be expected that a finite number of new intersections between the manifolds $W^s(S_{\sigma}^{sd})$ and $W^u(S_{\sigma}^{sd})$ perturb from the initial intersection given by the homoclinic orbit. These intersections provide a finite number of orbits, called maximal canards [43], that starting at $W^u(S_{\sigma}^{sd})$, they reach $W^s(S_{\sigma}^{sd})$ and all these orbits perturb from the homoclinic orbit. Additionally, the number of maximal canards will depend on the singular parameter $\sigma$, let us denote this number by $N(\sigma)$. Finally, $N(\sigma)$ tipically exhibits a monotonically decreasing behavior with $\sigma$ and $\lim_{\sigma \to 0} N(\sigma) = \infty$.

In order to prove the existence of the connections between the stable and the unstable manifold, we proceed as follows. First, we study the existence of this kind of connections in the discontinuous case $a = 1$, where the situation is easy to analyze due to the fact that only two linearity regions are involved.
Then, by a suitable use of the Implicit Function Theorem, we extend the
existence of these objects to the continuous context $a < 1$.

## 4.1 Slow passage in the discontinuous case $a = 1$

In this section, we indicate the set of equations that the $n$-spikes maximal canards must satisfy, where $n \le N(\sigma)$, in the discontinuous case given by $a = 1$. We also provide an analytical proof for the existence of solutions of this set of equations in the particular cases of 1-spike and 2-spikes. Finally, we find numerical evidences of the existence of maximal canards with more number of spikes.

<figure><img src="image_6.png" /><figcaption>**Fig. 6 Sketch of the invariant planes and maximal canards.** Representation of the critical manifold ($S_0^{sd}$, dotted blue line) together with its stable ($W^s(S_0^{sd})$, light blue plane) and unstable ($W^u(S_0^{sd})$, light red plane) manifolds, and the slow manifold $S_\sigma^{sd}$ together with its stable $W^s(S_\sigma^{sd})$ and unstable $W^u(S_\sigma^{sd})$ manifolds (dark blue and dark red). Locally, both unstable manifolds are coinciding planes and both stable manifolds are parallel planes. Additionally, 1-spike maximal canard (blue orbit) and 2-spikes maximal canard (green orbit) are also represented together with their intersection points with the plane $\{x=1\}$, which is the common boundary between the linearity regions $R_{3,1}$ and $R_{4,3}$.</figcaption></figure>

Along this section we deal with the case $a = 1$, hence, the region $R_{3,2}$ vanishes and the regions $R_{3,1}$ and $R_{4,3}$ share $\{x = 1\}$ as a common boundary, see Figure 3 and Figure 12. In consequence, the $n$-spikes maximal canard evolves through the linearity regions $R_{3,1}$ and $R_{4,3}$, and it intersects with the switching plane $\{x = 1\}$ at $2n$ points, that we denote by $\{\mathbf{p}_i\}_{i=1}^{2n}$, being $\mathbf{p}_1 \in W^u(S_\sigma^{sd}) \cap \{x = 1\}$ and $\mathbf{p}_{2n} \in W^s(S_\sigma^{sd}) \cap \{x = 1\}$, see Figure 6. We set the equations determining the $n$-spikes maximal canard through the flow with initial conditions at the points $\{\mathbf{p}_i\}_{i=1}^{2n}$.

In particular, for a 1-spike maximal canard orbit, we consider the existence of a point $\mathbf{p}_1 \in W^u(S_\sigma^{sd}) \cap \{x = 1\}$ and a point $\mathbf{p}_2 \in W^s(S_\sigma^{sd}) \cap \{x = 1\}$ such that the orbit of system (7) with initial condition at $\mathbf{p}_1$ achieves $\mathbf{p}_2$ at a
Then, by a suitable use of the Implicit Function Theorem, we extend the
existence of these objects to the continuous context $a < 1$.

## 4.1 Slow passage in the discontinuous case $a = 1$

In this section, we indicate the set of equations that the $n$-spikes maximal canards must satisfy, where $n \le N(\sigma)$, in the discontinuous case given by $a = 1$. We also provide an analytical proof for the existence of solutions of this set of equations in the particular cases of 1-spike and 2-spikes. Finally, we find numerical evidences of the existence of maximal canards with more number of spikes.

<figure><img src="image_7.png" /><figcaption>**Fig. 6 Sketch of the invariant planes and maximal canards.** Representation of the critical manifold ($S_0^{sd}$, dotted blue line) together with its stable ($W^s(S_0^{sd})$, light blue plane) and unstable ($W^u(S_0^{sd})$, light red plane) manifolds, and the slow manifold $S_\sigma^{sd}$ together with its stable $W^s(S_\sigma^{sd})$ and unstable $W^u(S_\sigma^{sd})$ manifolds (dark blue and dark red). Locally, both unstable manifolds are coinciding planes and both stable manifolds are parallel planes. Additionally, 1-spike maximal canard (blue orbit) and 2-spikes maximal canard (green orbit) are also represented together with their intersection points with the plane $\{x=1\}$, which is the common boundary between the linearity regions $R_{3,1}$ and $R_{4,3}$.</figcaption></figure>

Along this section we deal with the case $a = 1$, hence, the region $R_{3,2}$ vanishes and the regions $R_{3,1}$ and $R_{4,3}$ share $\{x = 1\}$ as a common boundary, see Figure 3 and Figure 12. In consequence, the $n$-spikes maximal canard evolves through the linearity regions $R_{3,1}$ and $R_{4,3}$, and it intersects with the switching plane $\{x = 1\}$ at $2n$ points, that we denote by $\{\mathbf{p}_i\}_{i=1}^{2n}$, being $\mathbf{p}_1 \in W^u(S_\sigma^{sd}) \cap \{x = 1\}$ and $\mathbf{p}_{2n} \in W^s(S_\sigma^{sd}) \cap \{x = 1\}$, see Figure 6. We set the equations determining the $n$-spikes maximal canard through the flow with initial conditions at the points $\{\mathbf{p}_i\}_{i=1}^{2n}$.

In particular, for a 1-spike maximal canard orbit, we consider the existence of a point $\mathbf{p}_1 \in W^u(S_\sigma^{sd}) \cap \{x = 1\}$ and a point $\mathbf{p}_2 \in W^s(S_\sigma^{sd}) \cap \{x = 1\}$ such that the orbit of system (7) with initial condition at $\mathbf{p}_1$ achieves $\mathbf{p}_2$ at a
14 4.1 Slow passage in the discontinuous case $a = 1$

positive time $\tau_1$, see Figure 6, that is

$$ 
\begin{aligned}
\varphi(\tau_1; \mathbf{p}_1, (1, \sigma)) &= \mathbf{p}_2, \\
\varphi(t; \mathbf{p}_1, (1, \sigma)) &\subset R_{4,3}, \text{ for } t \in (0, \tau_1).
\end{aligned} 
\qquad (15) $$

First expression in system (15) is given by three equations, one for each coordinate, with four unknowns: one coordinate of point $\mathbf{p}_1$, one coordinate of point $\mathbf{p}_2$, the time of flight $\tau_1$ and the parameter $\sigma$, see Section 5.3. Second equation in (15) provides a condition in order to avoid spurious solutions and to guarantee that any solution of the first equation corresponds with a true orbit connecting both invariant manifolds $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$.

Let us now look for a 2-spikes maximal canard orbit that starts over the unstable manifold at $\mathbf{p}_1 \in W^u(S_\sigma^{sd}) \cap \{x = 1\}$, performs a loop involving region $R_{3,1}$ and intersecting the switching plane $\{x = 1\}$ at points $\mathbf{p}_2$ and $\mathbf{p}_3$, and finally connects with the stable manifold at $\mathbf{p}_4 \in W^s(S_\sigma^{sd}) \cap \{x = 1\}$, see green curve in Figure 6. Following that description, the equations determining a 2-spikes maximal canard are given by

$$ 
\begin{aligned}
\varphi(\tau_1; \mathbf{p}_1, (1, \sigma)) &= \mathbf{p}_2, \\
\varphi(\tau_2; \mathbf{p}_2, (1, \sigma)) &= \mathbf{p}_3, \\
\varphi(\tau_3; \mathbf{p}_3, (1, \sigma)) &= \mathbf{p}_4,
\end{aligned} 
\qquad (16) $$

with restrictions

$$ 
\begin{aligned}
\varphi(t; \mathbf{p}_1, (1, \sigma)) &\subset R_{4,3}, \text{ for } t \in (0, \tau_1), \\
\varphi(t; \mathbf{p}_2, (1, \sigma)) &\subset R_{3,1}, \text{ for } t \in (0, \tau_2), \\
\varphi(t; \mathbf{p}_3, (1, \sigma)) &\subset R_{4,3}, \text{ for } t \in (0, \tau_3).
\end{aligned} 
\qquad (17) $$

System (16) is formed by 9 equations, 3 equations for each expression, and 10 unknowns: one coordinate of each of the points $\mathbf{p}_1$ and $\mathbf{p}_4$; two coordinates of $\mathbf{p}_2$ and $\mathbf{p}_3$; three times of flight $\tau_1, \tau_2$ and $\tau_3$ and the value of the parameter $\sigma$. This system provides necessary conditions for the existence of the 2-spikes maximal canards. Nevertheless, to guarantee that solutions of system (16) correspond with true orbits, it must be also satisfied the conditions in (17).

In the next result we assert the existence of two curves in the parameter space along each of them the system exhibits maximal canard orbits with different number of spikes. In particular, along one curve it exhibits 1-spike maximal canard orbits and along the other one 2-spikes maximal canard orbits.

**Theorem 3** Consider system (7) with parameters $a = 1, 0 < \varepsilon < k$ and $0 \le |\sigma|$ small enough. Then, there exist two differentiable functions $I_1(\sigma)$ and $I_2(\sigma)$ such that:

a) if $I = I_1(\sigma)$, system (7) exhibits a 1-spike maximal canard as intersection of $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$; and
14 4.1 Slow passage in the discontinuous case a = 1

positive time $\tau_1$, see Figure 6, that is

$$ 
\begin{array}{l}
\varphi(\tau_1; \mathbf{p}_1, (1, \sigma)) = \mathbf{p}_2, \\
\varphi(t; \mathbf{p}_1, (1, \sigma)) \subset R_{4,3}, \text{ for } t \in (0, \tau_1).
\end{array} 
\qquad (15) $$

First expression in system (15) is given by three equations, one for each coordinate, with four unknowns: one coordinate of point $\mathbf{p}_1$, one coordinate of point $\mathbf{p}_2$, the time of flight $\tau_1$ and the parameter $\sigma$, see Section 5.3. Second equation in (15) provides a condition in order to avoid spurious solutions and to guarantee that any solution of the first equation corresponds with a true orbit connecting both invariant manifolds $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$.

Let us now look for a 2-spikes maximal canard orbit that starts over the unstable manifold at $\mathbf{p}_1 \in W^u(S_\epsilon^{sd}) \cap \{x = 1\}$, performs a loop involving region $R_{3,1}$ and intersecting the switching plane $\{x = 1\}$ at points $\mathbf{p}_2$ and $\mathbf{p}_3$, and finally connects with the stable manifold at $\mathbf{p}_4 \in W^s(S_\epsilon^{sd}) \cap \{x = 1\}$, see green curve in Figure 6. Following that description, the equations determining a 2-spikes maximal canard are given by

$$ 
\begin{array}{l}
\varphi(\tau_1; \mathbf{p}_1, (1, \sigma)) = \mathbf{p}_2, \\
\varphi(\tau_2; \mathbf{p}_2, (1, \sigma)) = \mathbf{p}_3, \\
\varphi(\tau_3; \mathbf{p}_3, (1, \sigma)) = \mathbf{p}_4,
\end{array} 
\qquad (16) $$

with restrictions

$$ 
\begin{array}{l}
\varphi(t; \mathbf{p}_1, (1, \sigma)) \subset R_{4,3}, \text{ for } t \in (0, \tau_1), \\
\varphi(t; \mathbf{p}_2, (1, \sigma)) \subset R_{3,1}, \text{ for } t \in (0, \tau_2), \\
\varphi(t; \mathbf{p}_3, (1, \sigma)) \subset R_{4,3}, \text{ for } t \in (0, \tau_3).
\end{array} 
\qquad (17) $$

System (16) is formed by 9 equations, 3 equations for each expression, and 10 unknowns: one coordinate of each of the points $\mathbf{p}_1$ and $\mathbf{p}_4$; two coordinates of $\mathbf{p}_2$ and $\mathbf{p}_3$; three times of flight $\tau_1, \tau_2$ and $\tau_3$ and the value of the parameter $\sigma$. This system provides necessary conditions for the existence of the 2-spikes maximal canards. Nevertheless, to guarantee that solutions of system (16) correspond with true orbits, it must be also satisfied the conditions in (17).

In the next result we assert the existence of two curves in the parameter space along each of them the system exhibits maximal canard orbits with different number of spikes. In particular, along one curve it exhibits 1-spike maximal canard orbits and along the other one 2-spikes maximal canard orbits.

**Theorem 3** Consider system (7) with parameters $a = 1, 0 < \varepsilon < k$ and $0 \le |\sigma|$ small enough. Then, there exist two differentiable functions $I_1(\sigma)$ and $I_2(\sigma)$ such that:

a) if $I = I_1(\sigma)$, system (7) exhibits a 1-spike maximal canard as intersection of $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$; and
b) if $I = I_2(\sigma)$, system (7) exhibits a 2-spikes maximal canard as intersection of $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$.

In Figure 7, we represent an example of 1-spike maximal canard (panel (a)) and 2-spikes maximal canard (panel (b)) when equations (15) and (16) are numerically solved considering the parameter values $\sigma = -0.001$, $k = 0.52$, $\varepsilon = 0.5$, $\delta = -0.1$.

In general, assuming the existence of a maximum number of spikes, $N(\sigma)$, for a maximal canard orbit, the set of equations to be satisfied by any n-spikes maximal canard with $n \le N(\sigma)$ that connects a point $\mathbf{p}_1 \in W^u(S_\varepsilon^{sd}) \cap \{x = 1\}$ with $\mathbf{p}_{2n} \in W^s(S_\varepsilon^{sd}) \cap \{x = 1\}$ is given by

$$ \varphi(\tau_i; \mathbf{p}_i, (1, \sigma)) = \mathbf{p}_{i+1}, \quad i = 1, 2, \dots, 2n-1, \qquad (18) $$

with restrictions

$$ \begin{aligned} \varphi(t; \mathbf{p}_{2j-1}, (1, \sigma)) &\subset R_{4,3}, \quad \text{for } t \in (0, \tau_{2j-1}), \\ \varphi(t; \mathbf{p}_{2j}, (1, \sigma)) &\subset R_{3,1}, \quad \text{for } t \in (0, \tau_{2j}), \end{aligned} \qquad (19) $$

where $j = 1, 2, \dots, n$.

System (18)-(19) contains $6n-3$ equations with $6n-2$ unknowns, which are: the 2 coordinates of each point $\mathbf{p}_i = (1, y_i, z_i)$, with $i = 2, 3, \dots, 2n-1$; 1 coordinate of $\mathbf{p}_1$ and another one of $\mathbf{p}_{2n}$, since the rest of coordinates are fixed given their belonging in the invariant manifold and the plane $\{x = 1\}$; the $2n-1$ times of flight; and, finally, the parameter $\sigma$.

<figure><img src="image_7.png" /><figcaption>Fig. 7 1- and 2-spikes maximal canards in the discontinuous system, for $a = 1$. (a) An example of a canard segment given by a numerical method to find an approximation of the connection at $I^* = -0.674320448288416$. A perturbation on $I^*$ of order $10^{-9}$ allows passing from 1-spike to 2-spikes maximal canard. (b) Another example of a canard segment when $I^* = -0.6719579131287817$. A perturbation of $I^*$ of order $10^{-9}$ allows passing from 2-spikes to 3-spikes maximal canard.</figcaption></figure>

However, to provide a proof of the existence of a solution for system (18) requires a more complex analytical treatment, which is out of the scope of this
b) if $I = I_2(\sigma)$, system (7) exhibits a 2-spikes maximal canard as intersection of $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$.

In Figure 7, we represent an example of 1-spike maximal canard (panel (a)) and 2-spikes maximal canard (panel (b)) when equations (15) and (16) are numerically solved considering the parameter values $\sigma = -0.001$, $k = 0.52$, $\varepsilon = 0.5$, $\delta = -0.1$.

In general, assuming the existence of a maximum number of spikes, $N(\sigma)$, for a maximal canard orbit, the set of equations to be satisfied by any n-spikes maximal canard with $n \le N(\sigma)$ that connects a point $\mathbf{p}_1 \in W^u(S_\varepsilon^{sd}) \cap \{x = 1\}$ with $\mathbf{p}_{2n} \in W^s(S_\varepsilon^{sd}) \cap \{x = 1\}$ is given by

$$ \varphi(\tau_i; \mathbf{p}_i, (1, \sigma)) = \mathbf{p}_{i+1}, \quad i = 1, 2, \dots, 2n-1, \qquad (18) $$

with restrictions

$$ \begin{aligned} \varphi(t; \mathbf{p}_{2j-1}, (1, \sigma)) &\subset R_{4,3}, \quad \text{for } t \in (0, \tau_{2j-1}), \\ \varphi(t; \mathbf{p}_{2j}, (1, \sigma)) &\subset R_{3,1}, \quad \text{for } t \in (0, \tau_{2j}), \end{aligned} \qquad (19) $$

where $j = 1, 2, \dots, n$.

System (18)-(19) contains $6n-3$ equations with $6n-2$ unknowns, which are: the 2 coordinates of each point $\mathbf{p}_i = (1, y_i, z_i)$, with $i = 2, 3, \dots, 2n-1$; 1 coordinate of $\mathbf{p}_1$ and another one of $\mathbf{p}_{2n}$, since the rest of coordinates are fixed given their belonging in the invariant manifold and the plane $\{x = 1\}$; the $2n-1$ times of flight; and, finally, the parameter $\sigma$.

<figure><img src="image_8.png" /><figcaption>Fig. 7 1- and 2-spikes maximal canards in the discontinuous system, for $a = 1$. (a) An example of a canard segment given by a numerical method to find an approximation of the connection at $I^* = -0.674320448288416$. A perturbation on $I^*$ of order $10^{-9}$ allows passing from 1-spike to 2-spikes maximal canard. (b) Another example of a canard segment when $I^* = -0.6719579131287817$. A perturbation of $I^*$ of order $10^{-9}$ allows passing from 2-spikes to 3-spikes maximal canard.</figcaption></figure>

However, to provide a proof of the existence of a solution for system (18) requires a more complex analytical treatment, which is out of the scope of this
16 4.1 Slow passage in the discontinuous case $a = 1$

manuscript. Nevertheless, some numerical approach is given, in order to show the existence of maximal canard connections. For instance, in Figure 8, we provide a numerical approximation of the maximal canard orbit for 3-spikes, showing again a canard segment.

<figure><img src="image_8.png" /><figcaption>Fig. 8 A 3-spikes maximal canard in the discontinuous system, for $a = 1$. An example of a canard segment, given by a numerical method, to find an approximation of the connection at $I^* = -0.6701477318060869$. A perturbation of orbits of order $10^{-9}$ allow them to realize one more spike before leaving the slow manifold.</figcaption></figure>

The presence of a maximal canard connection with an arbitrarily number of spikes should be allowed for sufficiently small $\sigma$, as long as the determinant of the Jacobian matrix of the connection equations (18) does not nullify. However, the extension of the existence of these connections for different values of $\sigma$ may vary depending on the number of loops. In Figure 9, we find a numerical approach of these curves for 1-, 2- and 3-spikes.

<figure><img src="image_9.png" /><figcaption>Fig. 9 Maximal canard curves for $a = 1$. Numerical curves of maximal canard connections for 1-, 2- and 3-spikes, which are $I_1(\sigma)$, $I_2(\sigma)$, and $I_3(\sigma)$, respectively. As can be seen, the domain of each curve decreases as $n$ increases. Since $\sigma = -0.25$ is not contained in the domain of the function $I_3$, it follows that $N(-0.25) = 2$.</figcaption></figure>
16 4.1 Slow passage in the discontinuous case $a = 1$

manuscript. Nevertheless, some numerical approach is given, in order to show the existence of maximal canard connections. For instance, in Figure 8, we provide a numerical approximation of the maximal canard orbit for 3-spikes, showing again a canard segment.

<figure><img src="image_9.png" /><figcaption>Fig. 8 A 3-spikes maximal canard in the discontinuous system, for $a = 1$. An example of a canard segment, given by a numerical method, to find an approximation of the connection at $I^* = -0.6701477318060869$. A perturbation of orbits of order $10^{-9}$ allow them to realize one more spike before leaving the slow manifold.</figcaption></figure>

The presence of a maximal canard connection with an arbitrarily number of spikes should be allowed for sufficiently small $\sigma$, as long as the determinant of the Jacobian matrix of the connection equations (18) does not nullify. However, the extension of the existence of these connections for different values of $\sigma$ may vary depending on the number of loops. In Figure 9, we find a numerical approach of these curves for 1-, 2- and 3-spikes.

<figure><img src="image_10.png" /><figcaption>Fig. 9 Maximal canard curves for $a = 1$. Numerical curves of maximal canard connections for 1-, 2- and 3-spikes, which are $I_1(\sigma)$, $I_2(\sigma)$, and $I_3(\sigma)$, respectively. As can be seen, the domain of each curve decreases as $n$ increases. Since $\sigma = -0.25$ is not contained in the domain of the function $I_3$, it follows that $N(-0.25) = 2$.</figcaption></figure>
As it can be observed in Figure 9, each curve $I = I_n(\sigma)$ characterizing the existence of the n-spikes maximal canard seems to be defined over an interval which size decreases as $n$ increases, as it happens in [31]. Following that and fixing $\sigma$, there exists a finite number of intervals containing the specific $\sigma$. This number, $N(\sigma)$, represents the maximal number of spikes of the maximal canards that can be expected.

## 4.2 Slow passage in the continuous case $a < 1$

To finalize the study of the slow passage, we consider the continuous case focusing only on the existence of 1-spike maximal canard. To do so, we extend the result given by the discontinuous case in Theorem 3(a).

Let us consider the set of equations allowing the 1-spike maximal canard in the continuous case, that is

$$ 
\begin{aligned}
\varphi(\tau_1; \mathbf{p}_1, (a, \sigma)) &= \mathbf{p}_2, \\
\varphi(\tau_2; \mathbf{p}_2, (a, \sigma)) &= \mathbf{p}_3, \\
\varphi(\tau_3; \mathbf{p}_3, (a, \sigma)) &= \mathbf{p}_4,
\end{aligned} 
\qquad (20) $$

with restrictions

$$ 
\begin{aligned}
\varphi(t; \mathbf{p}_1, (a, \sigma)) &\subset R_{3,2}, \text{ for } t \in (0, \tau_1), \\
\varphi(t; \mathbf{p}_2, (a, \sigma)) &\subset R_{4,3}, \text{ for } t \in (0, \tau_2), \\
\varphi(t; \mathbf{p}_3, (a, \sigma)) &\subset R_{3,2}, \text{ for } t \in (0, \tau_3),
\end{aligned} 
\qquad (21) $$

Here, each equation of (20) consists of a set of three equations. Hence, we have 9 equations with 11 unknowns: 1 coordinate of points $\mathbf{p}_1$ and $\mathbf{p}_4$, 2 coordinates of points $\mathbf{p}_2$ and $\mathbf{p}_3$, the 3 times of flight $\tau_1$, $\tau_2$, and $\tau_3$, and the parameters $a$ and $\sigma$. In the next result we provide the existence of a surface in the parameter space $(a, \sigma)$ where the differential system exhibits a 1-spike maximal canard orbit.

**Theorem 4** Consider system (7) with $0 < \varepsilon < k$, and parameters $0 \le 1 - a$ and $0 \le |\sigma|$ being small enough. Let $z_1$ be the initial condition for the $I$ variable. Then there exists one differentiable function depending on $\sigma$, $I_1(\sigma; a, \varepsilon)$, such that, if $z_1 = I_1(\sigma; a, \varepsilon)$, system (7) exhibits a 1-spike maximal canard as intersection of $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$.

Similarly, as we extend the discontinuous case to the continuous one (that is, from equation (15) to equations (20)-(21)), it should be feasible to extend the equations (16)-(17) to provide the existence of a 2-spikes maximal canard also in the continuous case, and so to write the following remark.

**Remark 1** Consider system (7) with $0 < \varepsilon < k$, and parameters $0 \le 1 - a$ and $0 \le |\sigma|$ being small enough. Let $z_1$ be the initial condition for the $I$ variable. There
As it can be observed in Figure 9, each curve $I = I_n(\sigma)$ characterizing the existence of the $n$-spikes maximal canard seems to be defined over an interval which size decreases as $n$ increases, as it happens in [31]. Following that and fixing $\sigma$, there exists a finite number of intervals containing the specific $\sigma$. This number, $N(\sigma)$, represents the maximal number of spikes of the maximal canards that can be expected.

## 4.2 Slow passage in the continuous case $a < 1$

To finalize the study of the slow passage, we consider the continuous case focusing only on the existence of 1-spike maximal canard. To do so, we extend the result given by the discontinuous case in Theorem 3(a).

Let us consider the set of equations allowing the 1-spike maximal canard in the continuous case, that is

$$ 
\begin{aligned}
\varphi(\tau_1; \mathbf{p}_1, (a, \sigma)) &= \mathbf{p}_2, \\
\varphi(\tau_2; \mathbf{p}_2, (a, \sigma)) &= \mathbf{p}_3, \\
\varphi(\tau_3; \mathbf{p}_3, (a, \sigma)) &= \mathbf{p}_4,
\end{aligned} 
\qquad (20) $$

with restrictions

$$ 
\begin{aligned}
\varphi(t; \mathbf{p}_1, (a, \sigma)) &\subset R_{3,2}, \text{ for } t \in (0, \tau_1), \\
\varphi(t; \mathbf{p}_2, (a, \sigma)) &\subset R_{4,3}, \text{ for } t \in (0, \tau_2), \\
\varphi(t; \mathbf{p}_3, (a, \sigma)) &\subset R_{3,2}, \text{ for } t \in (0, \tau_3),
\end{aligned} 
\qquad (21) $$

Here, each equation of (20) consists of a set of three equations. Hence, we have 9 equations with 11 unknowns: 1 coordinate of points $\mathbf{p}_1$ and $\mathbf{p}_4$, 2 coordinates of points $\mathbf{p}_2$ and $\mathbf{p}_3$, the 3 times of flight $\tau_1$, $\tau_2$, and $\tau_3$, and the parameters $a$ and $\sigma$. In the next result we provide the existence of a surface in the parameter space $(a, \sigma)$ where the differential system exhibits a 1-spike maximal canard orbit.

**Theorem 4** Consider system (7) with $0 < \varepsilon < k$, and parameters $0 \le 1 - a$ and $0 \le |\sigma|$ being small enough. Let $z_1$ be the initial condition for the $I$ variable. Then there exists one differentiable function depending on $\sigma$, $I_1(\sigma; a, \varepsilon)$, such that, if $z_1 = I_1(\sigma; a, \varepsilon)$, system (7) exhibits a 1-spike maximal canard as intersection of $W^u(S_{\sigma}^{sd})$ and $W^s(S_{\sigma}^{sd})$.

Similarly, as we extend the discontinuous case to the continuous one (that is, from equation (15) to equations (20)-(21)), it should be feasible to extend the equations (16)-(17) to provide the existence of a 2-spikes maximal canard also in the continuous case, and so to write the following remark.

**Remark 1** Consider system (7) with $0 < \varepsilon < k$, and parameters $0 \le 1 - a$ and $0 \le |\sigma|$ being small enough. Let $z_1$ be the initial condition for the $I$ variable. There
exists one differentiable function depending on $\sigma$, $I_2(\sigma; a, \varepsilon)$, such that, if $z_1 = I_2(\sigma; a, \varepsilon)$, system (7) exhibits a 2-spike maximal canard as intersection of $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$.

# 5 Proof of main results

## 5.1 Proof of Theorem 1

In order to locate the equilibrium points of the system, we describe the shape of the x-nulcline given in (3). Hence, since $\delta < 0$ and $k > 0$, function $f(x)$ is decreasing when $x < \sqrt{\varepsilon}$ and increasing when $\sqrt{\varepsilon} < x < 1$. Therefore, it has a relative minimum at $x = \sqrt{\varepsilon}$, being $f(\sqrt{\varepsilon}) = \sqrt{\varepsilon}(1 + 2\delta)$.

Assuming $(a, \varepsilon, I) \in \mathcal{E}_3$, see equations (5), the x-nullcline intersects the y-nullcline at exactly three equilibrium points: $\mathbf{p}_{12,1} \in R_{1,1} \cup R_{2,1}$, $\mathbf{p}_{3,1} \in R_{3,1}$ and $\mathbf{p}_{3,2} \in R_{3,2}$ (see Figure 10). This proves statement (a).

<figure><img src="image_10.png" /><figcaption>**Fig. 10** Configuration of three equilibria given in $\mathcal{E}_3$. The existence of the equilibrium point $\mathbf{p}_{3,2}$ and the minimum at $y < 0$ assures the existence of a configuration of three equilibrium points.</figcaption></figure>

Assuming $(a, \varepsilon, I) \in \mathcal{E}_1$, see equations (5), and depending on the parameter $I$ we can have different scenarios. When $I < -(ak+\beta) < -\sqrt{\varepsilon}(1+2\delta)$, we have that $f(\sqrt{\varepsilon}) + I < 0$ and $f(a) + I < 0$, implying that the x-nullcline intersects the y-nullcline at exactly one point $\mathbf{p}_{12,1} \in R_{1,1} \cup R_{2,1}$, see Figure 11(a). When $-(ak + \beta) < -\sqrt{\varepsilon}(1 + 2\delta) < I$, then $f(\sqrt{\varepsilon}) + I > 0$ and $f(a) + I > 0$ implying that x-nullcline intersects y-nullcline at exactly one point $\mathbf{p}_{3,2} \in R_{3,2}$, see Figure 11(b). This proves statement (b).

The limiting case given by $I = -\sqrt{\varepsilon}(1 + 2\delta)$ (respectively, $I = -(ka + \beta)$) corresponds with the collapse of the equilibrium points $\mathbf{p}_{12,1}$ and $\mathbf{p}_{3,1}$ (respectively, $\mathbf{p}_{3,1}$ and $\mathbf{p}_{3,2}$) in a saddle-node equilibrium point whereas the equilibrium point $\mathbf{p}_{3,2}$ (respectively, $\mathbf{p}_{12,1}$) persists. As a result, as we perturb parameter $I$, we change to either $\mathcal{E}_1$ or $\mathcal{E}_3$. This proves statements (c) and (d), respectively.

Finally, the common boundary between the sets $\mathcal{S}\mathcal{N}_1$ and $\mathcal{S}\mathcal{N}_2$ is $C$, which gives us the last statement, since the collapse of two saddle-node points is a point of codimension two. There, not only a perturbation of parameter $I$
exists one differentiable function depending on $\sigma$, $I_2(\sigma; a, \varepsilon)$, such that, if $z_1 = I_2(\sigma; a, \varepsilon)$, system (7) exhibits a 2-spike maximal canard as intersection of $W^u(S_{\sigma}^{sd})$ and $W^s(S_{\sigma}^{sd})$.

# 5 Proof of main results

## 5.1 Proof of Theorem 1

In order to locate the equilibrium points of the system, we describe the shape of the x-nulcline given in (3). Hence, since $\delta < 0$ and $k > 0$, function $f(x)$ is decreasing when $x < \sqrt{\varepsilon}$ and increasing when $\sqrt{\varepsilon} < x < 1$. Therefore, it has a relative minimum at $x = \sqrt{\varepsilon}$, being $f(\sqrt{\varepsilon}) = \sqrt{\varepsilon}(1 + 2\delta)$.

Assuming $(a, \varepsilon, I) \in \mathcal{E}_3$, see equations (5), the x-nullcline intersects the y-nullcline at exactly three equilibrium points: $\mathbf{p}_{12,1} \in R_{1,1} \cup R_{2,1}$, $\mathbf{p}_{3,1} \in R_{3,1}$ and $\mathbf{p}_{3,2} \in R_{3,2}$ (see Figure 10). This proves statement (a).

<figure><img src="image_11.png" /><figcaption>**Fig. 10 Configuration of three equilibria given in $\mathcal{E}_3$.** The existence of the equilibrium point $\mathbf{p}_{3,2}$ and the minimum at $y < 0$ assures the existence of a configuration of three equilibrium points.</figcaption></figure>

Assuming $(a, \varepsilon, I) \in \mathcal{E}_1$, see equations (5), and depending on the parameter $I$ we can have different scenarios. When $I < -(ak+\beta) < -\sqrt{\varepsilon}(1+2\delta)$, we have that $f(\sqrt{\varepsilon}) + I < 0$ and $f(a) + I < 0$, implying that the x-nullcline intersects the y-nullcline at exactly one point $\mathbf{p}_{12,1} \in R_{1,1} \cup R_{2,1}$, see Figure 11(a). When $-(ak + \beta) < -\sqrt{\varepsilon}(1 + 2\delta) < I$, then $f(\sqrt{\varepsilon}) + I > 0$ and $f(a) + I > 0$ implying that x-nullcline intersects y-nullcline at exactly one point $\mathbf{p}_{3,2} \in R_{3,2}$, see Figure 11(b). This proves statement (b).

The limiting case given by $I = -\sqrt{\varepsilon}(1 + 2\delta)$ (respectively, $I = -(ka + \beta)$) corresponds with the collapse of the equilibrium points $\mathbf{p}_{12,1}$ and $\mathbf{p}_{3,1}$ (respectively, $\mathbf{p}_{3,1}$ and $\mathbf{p}_{3,2}$) in a saddle-node equilibrium point whereas the equilibrium point $\mathbf{p}_{3,2}$ (respectively, $\mathbf{p}_{12,1}$) persists. As a result, as we perturb parameter $I$, we change to either $\mathcal{E}_1$ or $\mathcal{E}_3$. This proves statements (c) and (d), respectively.

Finally, the common boundary between the sets $\mathcal{S}\mathcal{N}_1$ and $\mathcal{S}\mathcal{N}_2$ is $C$, which gives us the last statement, since the collapse of two saddle-node points is a point of codimension two. There, not only a perturbation of parameter $I$
<figure><img src="image_11.png" /><figcaption>Fig. 11 Configurations of 1 equilibrium point given in $\mathcal{E}_1$. (a) For $(a, \varepsilon, I) \in \mathcal{E}_1$ and $I < -(ak+\beta)$, at $x = \sqrt{\varepsilon}$, the y-coordinate remains negative, and thus there exists a unique equilibrium point, which is $\mathbf{p}_{12,1}$. (b) If $I > -\sqrt{\varepsilon}(1+2\delta)$, the minimum of the x-nullcline is situated at $x = \sqrt{\varepsilon}$ with $y > 0$. Therefore, it only appears the equilibrium point $\mathbf{p}_{3,2}$.</figcaption></figure>

causes the change to either $\mathcal{E}_1$ or $\mathcal{E}_3$, but also the change of either $a$ or $\varepsilon$. This ends the proof of Theorem 1.

## 5.2 Proof of Theorem 2

A singular case of system (2) is when $a = 1$, which creates a piecewise linear discontinuous vector field (see Figure 12). Discontinuous vector fields are quite common and studied, even in the piecewise linear framework (see [56, 57]). Moreover, due to the nonexistence of sliding points in the switching line $x = 1$, system (2) is of sewing type [58], which allows us to consider the continuous extension of the flow $\varphi(t; \mathbf{p}, (1, \varepsilon, I))$ to the boundary $x = 1$.

Along this section we prove the existence of a homoclinic orbit in system (2). First, we consider the discontinuous case $a = 1$. Then, by a suitable use of the Implicit Function Theorem, we extend the existence of this homoclinic orbit through $a < 1$. This is a novel technique for proving the existence of global dynamical objects. We call this technique by contraction and expansion technique, and it will be used several times along this manuscript.

Consider parameters $(a, \varepsilon, I) \in \mathcal{E}_3$ with $a = 1$, and let us establish the conditions to provide a homoclinic connection to the saddle point $\mathbf{p}_{3,1} \in R_{3,1}$. From Appendix A, the stable manifold $W^s(\mathbf{p}_{3,1})$ intersects with the switching line $\{x = 1\}$ at the point $(1, \hat{y}_1(\varepsilon))$, with

$$ \hat{y}_1(\varepsilon) = \frac{(k + \varepsilon)(k + \beta + I)}{k}, \qquad (22) $$

and the unstable manifold $W^u(\mathbf{p}_{3,1})$ intersects $\{x = 1\}$ at the point $(1,0)$, see Figure 12. Note that $\hat{y}_1(\varepsilon)$ also depends on parameters $k$ and $I$, we avoid making the dependence on $k$ and $I$ explicit to simplify the notation. Therefore, the equation

$$ \varphi(\tau; (1,0), (1, \varepsilon, I)) = (1, \hat{y}_1(\varepsilon)), \qquad (23) $$

with restriction

$$ \varphi(t; (1,0), (1, \varepsilon, I)) \subset R_{4,3}, \quad t \in (0, \tau), \qquad (24) $$

determines the existence of a homoclinic connection to the point $\mathbf{p}_{3,1}$. Therefore, according to its local expression in region $R_{4,3}$, see Appendix A,
<figure><img src="image_12.png" /><figcaption>Fig. 11 Configurations of 1 equilibrium point given in $\mathcal{E}_1$. (a) For $(a, \varepsilon, I) \in \mathcal{E}_1$ and $I < -(ak+\beta)$, at $x = \sqrt{\varepsilon}$, the y-coordinate remains negative, and thus there exists a unique equilibrium point, which is $\mathbf{p}_{12,1}$. (b) If $I > -\sqrt{\varepsilon}(1+2\delta)$, the minimum of the x-nullcline is situated at $x = \sqrt{\varepsilon}$ with $y > 0$. Therefore, it only appears the equilibrium point $\mathbf{p}_{3,2}$.</figcaption></figure>

causes the change to either $\mathcal{E}_1$ or $\mathcal{E}_3$, but also the change of either $a$ or $\varepsilon$. This ends the proof of Theorem 1.

## 5.2 Proof of Theorem 2

A singular case of system (2) is when $a = 1$, which creates a piecewise linear discontinuous vector field (see Figure 12). Discontinuous vector fields are quite common and studied, even in the piecewise linear framework (see [56, 57]). Moreover, due to the nonexistence of sliding points in the switching line $x = 1$, system (2) is of sewing type [58], which allows us to consider the continuous extension of the flow $\varphi(t; \mathbf{p}, (1, \varepsilon, I))$ to the boundary $x = 1$.

Along this section we prove the existence of a homoclinic orbit in system (2). First, we consider the discontinuous case $a = 1$. Then, by a suitable use of the Implicit Function Theorem, we extend the existence of this homoclinic orbit through $a < 1$. This is a novel technique for proving the existence of global dynamical objects. We call this technique by contraction and expansion technique, and it will be used several times along this manuscript.

Consider parameters $(a, \varepsilon, I) \in \mathcal{E}_3$ with $a = 1$, and let us establish the conditions to provide a homoclinic connection to the saddle point $\mathbf{p}_{3,1} \in R_{3,1}$. From Appendix A, the stable manifold $W^s(\mathbf{p}_{3,1})$ intersects with the switching line $\{x = 1\}$ at the point $(1, \hat{y}_1(\varepsilon))$, with

$$ \hat{y}_1(\varepsilon) = \frac{(k + \varepsilon)(k + \beta + I)}{k}, \qquad (22) $$

and the unstable manifold $W^u(\mathbf{p}_{3,1})$ intersects $\{x = 1\}$ at the point $(1,0)$, see Figure 12. Note that $\hat{y}_1(\varepsilon)$ also depends on parameters $k$ and $I$, we avoid making the dependence on $k$ and $I$ explicit to simplify the notation. Therefore, the equation

$$ \varphi(\tau; (1,0), (1, \varepsilon, I)) = (1, \hat{y}_1(\varepsilon)), \qquad (23) $$

with restriction

$$ \varphi(t; (1,0), (1, \varepsilon, I)) \subset R_{4,3}, \quad t \in (0, \tau), \qquad (24) $$

determines the existence of a homoclinic connection to the point $\mathbf{p}_{3,1}$. Therefore, according to its local expression in region $R_{4,3}$, see Appendix A,
20

5.2 Proof of Theorem 2

<figure><img src="image_12.png" /><figcaption>Fig. 12 Sketch of the nullclines configuration in the discontinuous case. The green line corresponds to the x-nullcline and the blue one to the y-nullcline, where it appears a discontinuity at x = 1. Vertical dashed lines correspond to the boundaries between regions. The black arrows define the homoclinic connection such as, the arrow above the y-nullcline denotes the unstable manifold of the saddle W^u(p_{3,1}), those above the dotted curve denotes the homoclinic path in R_{4,3}, and finally, the arrow above the black line denotes the stable manifold of the saddle W^s(p_{3,1}).</figcaption></figure>

equation (23) can be reduced to the system

$$ (1 - (k + \beta + I))e^{-\tau/\epsilon} + \frac{1}{1 - \epsilon}(e^{-\tau} - e^{-\tau/\epsilon}) + k + \beta + I = 1, \quad (25) \\ -e^{-\tau} + 1 = \hat{y}_1(\epsilon), $$

where the unknown variables are $(\epsilon, I, \tau)$. The following result assures, for every $\epsilon \in (0, k)$, the existence and uniqueness of a solution $(\epsilon, I_{\epsilon}, \tau_{\epsilon})$ of system (25), providing, consequently, the existence of a unique homoclinic connection to the saddle point $p_{3,1}$ given $\epsilon$.

**Lemma 1** Consider the differential system (2) with parameter $a = 1$. For every $\epsilon \in (0, k)$, there exists a unique value $I_{\epsilon}$ and time $\tau_{\epsilon}$ such that $(\epsilon, I_{\epsilon}, \tau_{\epsilon})$ is a solution of system (25).

*Proof* Let us consider the expression of the first integral $H$ (see equation (A2)). To guarantee the connection between points $(x, y) = (1, 0)$ and $(x, y) = (1, \hat{y}_1(\epsilon))$, the equality $H(1, v) = H(1, -u)$ must hold. This defines an implicit function $v = v_1(u)$ given by the implicit expression

$$ \frac{\epsilon - v}{(1 - v)^{1/\epsilon}} = \frac{\epsilon + u}{(1 + u)^{1/\epsilon}}, \quad (26) $$

which starts at $v_1(0) = 0$ and both equality sides are continuous functions defined for $u \ge 0$ and $0 \le v_1(u) \le 1$. We remark that $\lim_{u \to \infty} \frac{(\epsilon - v_1(u))}{(1 - v_1(u))^{1/\epsilon}} = 0$. By continuity,
20

5.2 Proof of Theorem 2

<figure><img src="image_13.png" /><figcaption>**Fig. 12** Sketch of the nullclines configuration in the discontinuous case. The green line corresponds to the x-nullcline and the blue one to the y-nullcline, where it appears a discontinuity at $x = 1$. Vertical dashed lines correspond to the boundaries between regions. The black arrows define the homoclinic connection such as, the arrow above the y-nullcline denotes the unstable manifold of the saddle $W^u(p_{3,1})$, those above the dotted curve denotes the homoclinic path in $R_{4,3}$, and finally, the arrow above the black line denotes the stable manifold of the saddle $W^s(p_{3,1})$.</figcaption></figure>

equation (23) can be reduced to the system

$$ (1 - (k + \beta + I))e^{-\tau/\epsilon} + \frac{1}{1 - \epsilon}(e^{-\tau} - e^{-\tau/\epsilon}) + k + \beta + I = 1, \quad (25) \\ -e^{-\tau} + 1 = \hat{y}_1(\epsilon), $$

where the unknown variables are $(\epsilon, I, \tau)$. The following result assures, for every $\epsilon \in (0, k)$, the existence and uniqueness of a solution $(\epsilon, I_{\epsilon}, \tau_{\epsilon})$ of system (25), providing, consequently, the existence of a unique homoclinic connection to the saddle point $p_{3,1}$ given $\epsilon$.

**Lemma 1** Consider the differential system (2) with parameter $a = 1$. For every $\epsilon \in (0, k)$, there exists a unique value $I_{\epsilon}$ and time $\tau_{\epsilon}$ such that $(\epsilon, I_{\epsilon}, \tau_{\epsilon})$ is a solution of system (25).

*Proof* Let us consider the expression of the first integral $H$ (see equation (A2)). To guarantee the connection between points $(x, y) = (1, 0)$ and $(x, y) = (1, \hat{y}_1(\epsilon))$, the equality $H(1, v) = H(1, -u)$ must hold. This defines an implicit function $v = v_1(u)$ given by the implicit expression

$$ \frac{\epsilon - v}{(1 - v)^{1/\epsilon}} = \frac{\epsilon + u}{(1 + u)^{1/\epsilon}}, \quad (26) $$

which starts at $v_1(0) = 0$ and both equality sides are continuous functions defined for $u \ge 0$ and $0 \le v_1(u) \le 1$. We remark that $\lim_{u \to \infty} \frac{(\epsilon - v_1(u))}{(1 - v_1(u))^{1/\epsilon}} = 0$. By continuity,
and since the right-hand term of (26) is always non-negative for $u > 0$, it follows that $0 \le v_1(u) < \varepsilon$ and $\lim_{u \to +\infty} v_1(u) = \varepsilon$. Further, the pair $(u, v)$ corresponding to the connection also satisfies the relation $v = v_2(u) = (\varepsilon/k)u$. Consider $V(u) = v_1(u) - v_2(u)$. Since $V(0) = 0$ and $\lim_{u \to +\infty} V(u) = -\infty$, by showing that $V'(0) = v_1'(0) - \frac{\varepsilon}{k} > 0$ we will obtain that there exists a pair $(u, v) \ne (0, 0)$ such that $v = v_1(u) = v_2(u)$. Additionally, if $V''(u) = v_1''(u) < 0$ for $u > 0$, then this intersection must be unique.

In what follows we prove that $v_1'(0) > \varepsilon/k$ and $v_1''(0) < 0$.

Hence, deriving both sides of equation (26) in respect to $u$ and simplifying, we have that

$$ \frac{v_1(u)v_1'(u)}{(1 - v_1(u))^{1/\varepsilon + 1}} = \frac{u}{(1 + u)^{1/\varepsilon + 1}}, $$

which is a tautology at $v_1(0) = 0$. Therefore, we need to derive the equation again obtaining

$$ \frac{\varepsilon v_1(u)v_1''(u)(1 - v_1(u)) + (v_1'(u))^2(\varepsilon + v_1(u))}{(1 - v_1(u))^{1/\varepsilon + 2}} = \frac{\varepsilon - u}{(1 + u)^{1/\varepsilon + 2}}, $$

where $(v_1'(0))^2 = 1$ and, since $v_1(u) \in [0, \varepsilon)$, the derivative at the origin must be positive and therefore $v_1'(0) = 1 > \varepsilon/k$.

For simplifying notation, considering now $v_1 = v_1(u) \ne 0$, we know that the derivative of $v_1$ is given by

$$ v_1' = \frac{u}{v_1} \frac{(1 - v_1)^{1/\varepsilon + 1}}{(1 + u)^{1/\varepsilon + 1}}, $$

which is well defined due to $u \ge 0$ and $v_1 \ne \{0, 1\}$ (since it will imply that either $u = -1$ or $v_1 = 1 > \varepsilon$). At $u = 0$, the derivative presents a critical point. Isolating the second derivative, we obtain

$$ \begin{aligned} v_1'' &= \frac{\varepsilon - u}{\varepsilon v_1(1 - v_1)} \left( \frac{1 - v_1}{1 + u} \right)^{1/\varepsilon + 2} - \frac{(v_1')^2(\varepsilon + v_1)}{\varepsilon v_1(1 - v_1)} \\ &= \frac{1}{\varepsilon v_1(1 - v_1)} \left( \frac{1 - v_1}{1 + u} \right)^{1/\varepsilon + 2} \left( (\varepsilon - u) - (\varepsilon + v_1) \frac{u^2}{v_1^2} \left( \frac{1 - v_1}{1 + u} \right)^{1/\varepsilon} \right) \\ &= \frac{1}{\varepsilon v_1^3(1 - v_1)(\varepsilon + u)} \left( \frac{1 - v_1}{1 + u} \right)^{1/\varepsilon + 2} \left( (\varepsilon + u)(\varepsilon - u)v_1^2 - (\varepsilon + v_1)(\varepsilon - v_1)u^2 \right) \\ &= \frac{\varepsilon(v_1 + u)(v_1 - u)}{v_1^3(1 - v_1)(\varepsilon + u)} \left( \frac{1 - v_1}{1 + u} \right)^{1/\varepsilon + 2}. \end{aligned} $$

Since at $v_1(0) = 0$ we have $v_1'(0) = 1$, by the definition of the second derivative, we only have three possibilities: $v_1 < u$, $v_1 = u$ or $v_1 > u$ for all $v_1 \in (0, \varepsilon)$. However, since $\lim_{u \to \infty} v_1(u) = \varepsilon$, we know that there must exist a value such that $v_1 < u$. Therefore, $v_1 < u$ for all $u \in (0, \infty)$.

Finally, since $u = -1 + 1/(1 - (k + \beta + I))$, we can write $I$ in terms of $u$, that is $I(u)$. On the other hand, the connection found depends on $\varepsilon$ such that $u$ writes as $u(\varepsilon)$. Therefore, $I = I(u(\varepsilon)) = I(\varepsilon)$. Considering the orbit given the connection, and let $\tau$ the time of flight. Therefore, as it depends on the initial conditions of $u$, the same happens for the time of flight, $\tau(u(\varepsilon)) = \tau(\varepsilon)$. $\square$
and since the right-hand term of (26) is always non-negative for $u > 0$, it follows that $0 \le v_1(u) < \varepsilon$ and $\lim_{u \to +\infty} v_1(u) = \varepsilon$. Further, the pair $(u, v)$ corresponding to the connection also satisfies the relation $v = v_2(u) = (\varepsilon/k)u$. Consider $V(u) = v_1(u) - v_2(u)$. Since $V(0) = 0$ and $\lim_{u \to +\infty} V(u) = -\infty$, by showing that $V'(0) = v_1'(0) - \frac{\varepsilon}{k} > 0$ we will obtain that there exists a pair $(u, v) \ne (0, 0)$ such that $v = v_1(u) = v_2(u)$. Additionally, if $V''(u) = v_1''(u) < 0$ for $u > 0$, then this intersection must be unique.

In what follows we prove that $v_1'(0) > \varepsilon/k$ and $v_1''(0) < 0$.

Hence, deriving both sides of equation (26) in respect to $u$ and simplifying, we have that

$$ \frac{v_1(u)v_1'(u)}{(1 - v_1(u))^{1/\varepsilon + 1}} = \frac{u}{(1 + u)^{1/\varepsilon + 1}}, $$

which is a tautology at $v_1(0) = 0$. Therefore, we need to derive the equation again obtaining

$$ \frac{\varepsilon v_1(u)v_1''(u)(1 - v_1(u)) + (v_1'(u))^2(\varepsilon + v_1(u))}{(1 - v_1(u))^{1/\varepsilon + 2}} = \frac{\varepsilon - u}{(1 + u)^{1/\varepsilon + 2}}, $$

where $(v_1'(0))^2 = 1$ and, since $v_1(u) \in [0, \varepsilon)$, the derivative at the origin must be positive and therefore $v_1'(0) = 1 > \varepsilon/k$.

For simplifying notation, considering now $v_1 = v_1(u) \ne 0$, we know that the derivative of $v_1$ is given by

$$ v_1' = \frac{u}{v_1} \frac{(1 - v_1)^{1/\varepsilon + 1}}{(1 + u)^{1/\varepsilon + 1}}, $$

which is well defined due to $u \ge 0$ and $v_1 \ne \{0, 1\}$ (since it will imply that either $u = -1$ or $v_1 = 1 > \varepsilon$). At $u = 0$, the derivative presents a critical point. Isolating the second derivative, we obtain

$$ \begin{aligned} v_1'' &= \frac{\varepsilon - u}{\varepsilon v_1(1 - v_1)} \left(\frac{1 - v_1}{1 + u}\right)^{1/\varepsilon + 2} - \frac{(v_1')^2(\varepsilon + v_1)}{\varepsilon v_1(1 - v_1)} \\ &= \frac{1}{\varepsilon v_1(1 - v_1)} \left(\frac{1 - v_1}{1 + u}\right)^{1/\varepsilon + 2} \left( (\varepsilon - u) - (\varepsilon + v_1) \frac{u^2}{v_1^2} \left(\frac{1 - v_1}{1 + u}\right)^{1/\varepsilon} \right) \\ &= \frac{1}{\varepsilon v_1^3(1 - v_1)(\varepsilon + u)} \left(\frac{1 - v_1}{1 + u}\right)^{1/\varepsilon + 2} \left( (\varepsilon + u)(\varepsilon - u)v_1^2 - (\varepsilon + v_1)(\varepsilon - v_1)u^2 \right) \\ &= \frac{\varepsilon(v_1 + u)(v_1 - u)}{v_1^3(1 - v_1)(\varepsilon + u)} \left(\frac{1 - v_1}{1 + u}\right)^{1/\varepsilon + 2}. \end{aligned} $$

Since at $v_1(0) = 0$ we have $v_1'(0) = 1$, by the definition of the second derivative, we only have three possibilities: $v_1 < u$, $v_1 = u$ or $v_1 > u$ for all $v_1 \in (0, \varepsilon)$. However, since $\lim_{u \to \infty} v_1(u) = \varepsilon$, we know that there must exist a value such that $v_1 < u$. Therefore, $v_1 < u$ for all $u \in (0, \infty)$.

Finally, since $u = -1 + 1/(1 - (k + \beta + I))$, we can write $I$ in terms of $u$, that is $I(u)$. On the other hand, the connection found depends on $\varepsilon$ such that $u$ writes as $u(\varepsilon)$. Therefore, $I = I(u(\varepsilon)) = I(\varepsilon)$. Considering the orbit given the connection, and let $\tau$ the time of flight. Therefore, as it depends on the initial conditions of $u$, the same happens for the time of flight, $\tau(u(\varepsilon)) = \tau(\varepsilon)$. $\square$
Next, we proof that the values $I_{\varepsilon}$ and $\tau_{\varepsilon}$ obtained in Lemma 1 define differentiable functions $I(\varepsilon)$ and $\tau(\varepsilon)$ over $(0, k)$, providing a differentiable curve $(1, \varepsilon, I(\varepsilon))$ in the parameter space $(a, \varepsilon, I)$ (red curve in Figure 4(a)). Hence, system (2) with parameters on this curve exhibits a homoclinic orbit to the saddle $p_{3,1}$. Let us first introduce the following technical lemma.

**Lemma 2** Let us consider $0 < \varepsilon < k < 1$. Therefore, for $\tau > 0$, expressions

$$ (1 - \frac{k}{k + \varepsilon}(1 - e^{-\tau}))e^{-\tau/\varepsilon} + \frac{1}{1 - \varepsilon}(e^{-\tau} - e^{-\tau/\varepsilon}) - (1 - \frac{k}{k + \varepsilon}(1 - e^{-\tau})) \quad (27) $$

and

$$ \frac{1}{k}(1 - e^{-\tau}) - e^{-\tau}(1 - e^{-\tau/\varepsilon}) \quad (28) $$

do not vanish at the same time.

*Proof* By equalling equation (28) to 0, we find that $e^{-\tau/\varepsilon} = 1 + 1/k(1 - e^{\tau})$. Substituting in equation (27) and also equalling that to 0, we find that either $e^{\tau} = 1$ or $e^{\tau} = (k^2 + k)/(\varepsilon(k-1) - 2k)$. The first case cannot happen, since $\tau > 0$. The second one neither can be, since it is a negative number for $0 < k < 1$ equalling an exponential term. Therefore, it does not exist a time $\tau > 0$ nullifying both expressions, which follows the statement. $\square$

**Lemma 3** Consider the differential system (2) with parameter $a = 1$. Given $0 < \varepsilon_0 < k$, let $I_{\varepsilon_0}$ and $\tau_{\varepsilon_0}$ be the values obtained in Lemma 1. There exist neighbourhoods $V$ of $\varepsilon_0$ and $W$ of $(\varepsilon_0, I_{\varepsilon_0}, \tau_{\varepsilon_0})$, and differentiable functions $I(\varepsilon)$ and $\tau(\varepsilon)$ defined over $V$ such that $(\varepsilon, I(\varepsilon), \tau(\varepsilon))$ is contained in $W$ and solves equation (25). Therefore, there exists a homoclinic orbit through the points $(1, 0)$ and $(1, \hat{y}_1(\varepsilon))$.

*Proof* Given $0 < \varepsilon < k$ and following Lemma 1, there exist $I_{\varepsilon}$ and $\tau_{\varepsilon}$ such that $\varphi(\tau_{\varepsilon}; (1, 0), (1, \varepsilon, I_{\varepsilon})) = (1, \hat{y}_1(\varepsilon))$, that is, $(\varepsilon, I_{\varepsilon}, \tau_{\varepsilon})$ is a solution of system (25). From there, one can deduce equation (27) equaled to 0. Moreover, the determinant of the Jacobian matrix of (25) respect to the variables $(\tau, I)$ evaluated at $(\varepsilon, I_{\varepsilon}, \tau_{\varepsilon})$ is

$$ \left| \frac{1}{1-\varepsilon} e^{-\tau_{\varepsilon}} + \frac{k}{\varepsilon(k+\varepsilon)} (1 - e^{-\tau_{\varepsilon}/\varepsilon}) \left( 1 - \frac{e^{-\tau_{\varepsilon}/\varepsilon}}{e^{-\tau_{\varepsilon}}} \right) - \frac{k+\varepsilon}{k} \right| = \frac{1}{k} (1 - e^{-\tau_{\varepsilon}}) - e^{-\tau_{\varepsilon}} (1 - e^{-\tau_{\varepsilon}/\varepsilon}) $$

which coincides with equation (25), where we have written $k + \beta + I_{\varepsilon}$ and $e^{-\tau_{\varepsilon}/\varepsilon}$ in terms of $e^{-\tau_{\varepsilon}}$ through system (25) and used equation (27) to get $\frac{1}{1-\varepsilon}e^{-\tau_{\varepsilon}} + \frac{k}{\varepsilon(k+\varepsilon)}(1 - e^{-\tau_{\varepsilon}/\varepsilon}) = -(1 - e^{-\tau_{\varepsilon}})/(k+\varepsilon)$. By Lemma 2, neither system (25) nor its Jacobian nullify at the same time, and we can use the Implicit Function Theorem, which follows the statement of the lemma.

Notice that, since the behavior in the region $R_{4,3}$ is nodal, there are not possibilities for three consecutive intersections at the line $\{x = 1\}$, therefore restriction (24) is satisfied. $\square$

Finally, we end the case $a = 1$ by showing that the homoclinic orbit satisfies a non-degeneracy condition. In particular, let $(1, y_u)$ be the second
Next, we proof that the values $I_{\varepsilon}$ and $\tau_{\varepsilon}$ obtained in Lemma 1 define differentiable functions $I(\varepsilon)$ and $\tau(\varepsilon)$ over $(0, k)$, providing a differentiable curve $(1, \varepsilon, I(\varepsilon))$ in the parameter space $(a, \varepsilon, I)$ (red curve in Figure 4(a)). Hence, system (2) with parameters on this curve exhibits a homoclinic orbit to the saddle $p_{3,1}$. Let us first introduce the following technical lemma.

**Lemma 2** Let us consider $0 < \varepsilon < k < 1$. Therefore, for $\tau > 0$, expressions

$$ (1 - \frac{k}{k + \varepsilon}(1 - e^{-\tau}))e^{-\tau/\varepsilon} + \frac{1}{1 - \varepsilon}(e^{-\tau} - e^{-\tau/\varepsilon}) - (1 - \frac{k}{k + \varepsilon}(1 - e^{-\tau})) \quad (27) $$

and

$$ \frac{1}{k}(1 - e^{-\tau}) - e^{-\tau}(1 - e^{-\tau/\varepsilon}) \quad (28) $$

do not vanish at the same time.

*Proof* By equalling equation (28) to 0, we find that $e^{-\tau/\varepsilon} = 1 + 1/k(1 - e^{\tau})$. Substituting in equation (27) and also equalling that to 0, we find that either $e^{\tau} = 1$ or $e^{\tau} = (k^2 + k)/(\varepsilon(k-1) - 2k)$. The first case cannot happen, since $\tau > 0$. The second one neither can be, since it is a negative number for $0 < k < 1$ equalling an exponential term. Therefore, it does not exist a time $\tau > 0$ nullifying both expressions, which follows the statement. $\square$

**Lemma 3** Consider the differential system (2) with parameter $a = 1$. Given $0 < \varepsilon_0 < k$, let $I_{\varepsilon_0}$ and $\tau_{\varepsilon_0}$ be the values obtained in Lemma 1. There exist neighbourhoods $V$ of $\varepsilon_0$ and $W$ of $(\varepsilon_0, I_{\varepsilon_0}, \tau_{\varepsilon_0})$, and differentiable functions $I(\varepsilon)$ and $\tau(\varepsilon)$ defined over $V$ such that $(\varepsilon, I(\varepsilon), \tau(\varepsilon))$ is contained in $W$ and solves equation (25). Therefore, there exists a homoclinic orbit through the points $(1, 0)$ and $(1, \hat{y}_1(\varepsilon))$.

*Proof* Given $0 < \varepsilon < k$ and following Lemma 1, there exist $I_{\varepsilon}$ and $\tau_{\varepsilon}$ such that $\varphi(\tau_{\varepsilon}; (1, 0), (1, \varepsilon, I_{\varepsilon})) = (1, \hat{y}_1(\varepsilon))$, that is, $(\varepsilon, I_{\varepsilon}, \tau_{\varepsilon})$ is a solution of system (25). From there, one can deduce equation (27) equaled to 0. Moreover, the determinant of the Jacobian matrix of (25) respect to the variables $(\tau, I)$ evaluated at $(\varepsilon, I_{\varepsilon}, \tau_{\varepsilon})$ is

$$ \left| \frac{1}{1-\varepsilon} e^{-\tau_{\varepsilon}} + \frac{k}{\varepsilon(k+\varepsilon)} (1 - e^{-\tau_{\varepsilon}/\varepsilon}) \left( 1 - \frac{e^{-\tau_{\varepsilon}/\varepsilon}}{e^{-\tau_{\varepsilon}}} \right) - \frac{k+\varepsilon}{k} \right| = \frac{1}{k} (1 - e^{-\tau_{\varepsilon}}) - e^{-\tau_{\varepsilon}} (1 - e^{-\tau_{\varepsilon}/\varepsilon}) $$

which coincides with equation (25), where we have written $k + \beta + I_{\varepsilon}$ and $e^{-\tau_{\varepsilon}/\varepsilon}$ in terms of $e^{-\tau_{\varepsilon}}$ through system (25) and used equation (27) to get $\frac{1}{1-\varepsilon}e^{-\tau_{\varepsilon}} + \frac{k}{\varepsilon(k+\varepsilon)}(1 - e^{-\tau_{\varepsilon}/\varepsilon}) = -(1 - e^{-\tau_{\varepsilon}})/(k+\varepsilon)$. By Lemma 2, neither system (25) nor its Jacobian nullify at the same time, and we can use the Implicit Function Theorem, which follows the statement of the lemma.

Notice that, since the behavior in the region $R_{4,3}$ is nodal, there are not possibilities for three consecutive intersections at the line $\{x = 1\}$, therefore restriction (24) is satisfied. $\square$

Finally, we end the case $a = 1$ by showing that the homoclinic orbit satisfies a non-degeneracy condition. In particular, let $(1, y_u)$ be the second
intersection point of the unstable manifold of $\mathbf{p}_{3,1}$ with the switching line $\{x=1\}$, that is, the one after passing through point $(1,0)$. Since $(1, \hat{y}_1(\varepsilon))$ is the first intersection point of the stable manifold of $\mathbf{p}_{3,1}$ with $\{x=1\}$, function $\Lambda(I) = \hat{y}_1(\varepsilon) - y_u$ is the difference between the stable and the unstable manifolds of $\mathbf{p}_{3,1}$ at the switching line $x=1$. From Lemma 3, $I(\varepsilon)$ is a zero of $\Lambda(I)$, next we prove that, in fact, it is a simple zero. This means that the homoclinic connection $\varphi(\tau; (1,0), (1, \varepsilon, I(\varepsilon))) = (1, \hat{y}_1(\varepsilon))$ appears transversely as parameter $I$ varies. In particular, we obtain $\Lambda'(I(\varepsilon)) > 0$, which means that the region limited by $W^s(\mathbf{p}_{3,1})$ and $W^u(\mathbf{p}_{3,1})$ passes, as $I$ increases, from being a backward invariant set to be a forward invariant set.

**Lemma 4** Consider the differential system (2), with parameter $a=1$,

and the function $\Lambda(I) = \hat{y}_1(\varepsilon) - y_u$, where $y_u$ is the returning y-coordinate of the orbit starting at $(1,0)$, which is given by $y_u = (k + \beta + I) + u\varepsilon(1 - (k + \beta + I))/k$, where $u = -(k + \beta + I)/(k + \beta + I - 1)$. Then, $\Lambda(I)$ has a simple zero at $I = I(\varepsilon)$, where $I(\varepsilon)$ is given in Lemma 3 and $\Lambda'(I(\varepsilon)) > 0$.

*Proof* Let us consider $\Lambda(I) = \hat{y}_1(\varepsilon) - y_u$. By Lemma 3, there exists $I(\varepsilon)$ such that $\Lambda(I(\varepsilon)) = 0$. We want to show that $\Lambda'(I) > 0$ at $I = I(\varepsilon)$, that is, for the value of the connection. Effectively, it is a direct consequence of the fact that $\Lambda'(I) = 1 + u + 1/(k + \beta + I - 1)^2$ and, in the connection, $u > 0$ as $0 < k + \beta + I(\varepsilon) < 1$. $\square$

As part of the contraction and expansion technique, the next step consists of going from the discontinuous case, $a=1$, to the continuous one, $a < 1$, by proving the persistence of the homoclinic orbit obtained in Lemma 3, as long as $1-a$ is positive and sufficiently small. That is, we first extend the system of equations that determines the homoclinic orbit for $a=1$ to a new set of equations that involves region $R_{3,2}$ and that determines the homoclinic orbit for $a < 1$. From the set of equations defined for $a=1$, we prove the existence of a family of solutions for the new set of equations. Therefore, for $a < 1$, we consider equations

$$ 
\begin{aligned}
\varphi(\tau_1; (a,0), (a, \varepsilon, I)) &= (1, y_1), \\
\varphi(\tau_2; (1, y_1), (a, \varepsilon, I)) &= (1, y_2), \\
\varphi(\tau_3; (1, y_2), (a, \varepsilon, I)) &= (a, y(a, \varepsilon)),
\end{aligned} 
\qquad (29) $$

together with the restrictions

$$ 
\begin{aligned}
\varphi(t; (a,0), (a, \varepsilon, I)) &\subset R_{3,2}, \text{ for } t \in (0, \tau_1), \\
\varphi(t; (1, y_1), (a, \varepsilon, I)) &\subset R_{4,3}, \text{ for } t \in (0, \tau_2), \\
\varphi(t; (1, y_2), (a, \varepsilon, I)) &\subset R_{3,2}, \text{ for } t \in (0, \tau_3),
\end{aligned} 
\qquad (30) $$
intersection point of the unstable manifold of $\mathbf{p}_{3,1}$ with the switching line $\{x = 1\}$, that is, the one after passing through point $(1, 0)$. Since $(1, \hat{y}_1(\varepsilon))$ is the first intersection point of the stable manifold of $\mathbf{p}_{3,1}$ with $\{x = 1\}$, function $\Lambda(I) = \hat{y}_1(\varepsilon) - y_u$ is the difference between the stable and the unstable manifolds of $\mathbf{p}_{3,1}$ at the switching line $x = 1$. From Lemma 3, $I(\varepsilon)$ is a zero of $\Lambda(I)$, next we prove that, in fact, it is a simple zero. This means that the homoclinic connection $\varphi(\tau; (1, 0), (1, \varepsilon, I(\varepsilon))) = (1, \hat{y}_1(\varepsilon))$ appears transversely as parameter $I$ varies. In particular, we obtain $\Lambda'(I(\varepsilon)) > 0$, which means that the region limited by $W^s(\mathbf{p}_{3,1})$ and $W^u(\mathbf{p}_{3,1})$ passes, as $I$ increases, from being a backward invariant set to be a forward invariant set.

**Lemma 4** Consider the differential system (2), with parameter $a = 1$,

and the function $\Lambda(I) = \hat{y}_1(\varepsilon) - y_u$, where $y_u$ is the returning y-coordinate of the orbit starting at $(1, 0)$, which is given by $y_u = (k + \beta + I) + u\varepsilon(1 - (k + \beta + I))/k$, where $u = -(k + \beta + I)/(k + \beta + I - 1)$. Then, $\Lambda(I)$ has a simple zero at $I = I(\varepsilon)$, where $I(\varepsilon)$ is given in Lemma 3 and $\Lambda'(I(\varepsilon)) > 0$.

*Proof* Let us consider $\Lambda(I) = \hat{y}_1(\varepsilon) - y_u$. By Lemma 3, there exists $I(\varepsilon)$ such that $\Lambda(I(\varepsilon)) = 0$. We want to show that $\Lambda'(I) > 0$ at $I = I(\varepsilon)$, that is, for the value of the connection. Effectively, it is a direct consequence of the fact that $\Lambda'(I) = 1 + u + 1/(k + \beta + I - 1)^2$ and, in the connection, $u > 0$ as $0 < k + \beta + I(\varepsilon) < 1$. $\square$

As part of the contraction and expansion technique, the next step consists of going from the discontinuous case, $a = 1$, to the continuous one, $a < 1$, by proving the persistence of the homoclinic orbit obtained in Lemma 3, as long as $1 - a$ is positive and sufficiently small. That is, we first extend the system of equations that determines the homoclinic orbit for $a = 1$ to a new set of equations that involves region $R_{3,2}$ and that determines the homoclinic orbit for $a < 1$. From the set of equations defined for $a = 1$, we prove the existence of a family of solutions for the new set of equations. Therefore, for $a < 1$, we consider equations

$$ 
\begin{aligned}
\varphi(\tau_1; (a, 0), (a, \varepsilon, I)) &= (1, y_1), \\
\varphi(\tau_2; (1, y_1), (a, \varepsilon, I)) &= (1, y_2), \\
\varphi(\tau_3; (1, y_2), (a, \varepsilon, I)) &= (a, y(a, \varepsilon)),
\end{aligned} 
\qquad (29) $$

together with the restrictions

$$ 
\begin{aligned}
\varphi(t; (a, 0), (a, \varepsilon, I)) &\subset R_{3,2}, \text{ for } t \in (0, \tau_1), \\
\varphi(t; (1, y_1), (a, \varepsilon, I)) &\subset R_{4,3}, \text{ for } t \in (0, \tau_2), \\
\varphi(t; (1, y_2), (a, \varepsilon, I)) &\subset R_{3,2}, \text{ for } t \in (0, \tau_3),
\end{aligned} 
\qquad (30) $$
24 5.2 Proof of Theorem 2

where $y(a, \varepsilon) = (k + \varepsilon)(ak + \beta + I)/k$. Notice that $y(1, \varepsilon) = \hat{y}_1(\varepsilon)$, see (22), and consequently, $y(a, \varepsilon)$ is the ordinate of the intersection point of the stable manifold of $\mathbf{p}_{3,1}$ with the line $\{x = a\}$ for $a \le 1$.

Then, we get a system of 6 equations with 8 unknowns grouped as $(\tau_1, \tau_2, \tau_3, y_1, y_2, I, a, \varepsilon)$. Since $(\varepsilon, I(\varepsilon), \tau(\varepsilon))$ is the solution of system (25) obtained in Lemma 3, we note that, $(0, \tau(\varepsilon), 0, 0, \hat{y}_1(\varepsilon), I(\varepsilon), 1, \varepsilon)$, is a solution of system (29). The following result asserts that this solution can be extended for $1 - a > 0$ small enough, providing the existence of the homoclinic loop.

**Proposition 1** For each $0 < \varepsilon_0 < k$, let $\tau(\varepsilon_0)$ and $I(\varepsilon_0)$ be the functions given in Lemma 3, and $\hat{y}_1(\varepsilon_0)$ the function defined in (22). Therefore, there exist neighbourhoods $U$ of $(1, \varepsilon_0)$ and $V$ of $(0, \tau(\varepsilon_0), 0, 0, \hat{y}_1(\varepsilon_0), I(\varepsilon_0))$ and differentiable functions $\tau_1, \tau_2, \tau_3, Y_1, Y_2$ and $I$ defined from $U$ to $V$ such that $\tau_1(1, \varepsilon_0) = \tau_3(1, \varepsilon_0) = 0$, $\tau_2(1, \varepsilon_0) = \tau(\varepsilon_0)$, $Y_1(1, \varepsilon_0) = 0$, $Y_2(1, \varepsilon_0) = \hat{y}_1(\varepsilon_0)$, $I(1, \varepsilon_0) = I(\varepsilon_0)$. Moreover, for every $(a, \varepsilon) \in U$, the vector $(\tau_1(a, \varepsilon), \tau_2(a, \varepsilon), \tau_3(a, \varepsilon), Y_1(a, \varepsilon), Y_2(a, \varepsilon), I(a, \varepsilon), a, \varepsilon)$ is a solution of system (29).

*Proof* Let us label by $\phi_i = 0, i = 1, \dots, 6$ the $i$th equation of system (29). Then, from the local flow in region $R_{3,2}$ and $R_{4,3}$ (see Appendix A), it follows that

$$ \phi_1 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\rho\tau_1} \left( \xi_{1,1} \frac{\sin(\theta\tau_1)}{\sqrt{1-a}} - \cos(\theta\tau_1) \right) - 1 + (k + \beta + I) \right], $$

$$ \phi_2 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I) \left( \xi_{2,1} \frac{\sin(\theta\tau_1)}{\sqrt{1-a}} e^{\rho\tau_1} + (1 - e^{\rho\tau_1}) \frac{\cos(\theta\tau_1)}{1-a} \right) - (1 - (1-a)k) \frac{y_1}{1-a} \right], $$

$$ \phi_3 = (2 - \gamma - I) \left( e^{-\frac{\tau_2}{\varepsilon}} - 1 \right) + \frac{1 - y_1}{1 - \varepsilon} \left( e^{-\tau_2} - e^{-\frac{\tau_2}{\varepsilon}} \right), $$

$$ \phi_4 = -(1 - y_1)e^{-\tau_2} + 1 - y_2, $$

$$ \phi_5 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\rho\tau_3} \left( \xi_{5,1} \frac{\sin(\theta\tau_3)}{\sqrt{1-a}} - \cos(\theta\tau_3) \right) - 1 + (k + \beta + I) \right], $$

$$ \phi_6 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\rho\tau_3} \left( \xi_{6,1} \sin(\theta\tau_3) + \xi_{6,2} \cos(\theta\tau_3) \right) + \frac{1}{1-a} \left( (ak + \beta + I) - (1-(1-a)k)(y(a, \varepsilon) + y_2) \right) \right], $$

where $\rho = \frac{k-\varepsilon}{2\varepsilon}$, $\theta = \frac{\sqrt{4\varepsilon-(k+\varepsilon)^2}}{2\varepsilon}$ and

$$ \xi_{1,1} = \frac{2 - (1 - a)(k + \varepsilon)}{\sqrt{4\varepsilon - (1 - a)(k + \varepsilon)^2}} $$

$$ \xi_{2,1} = \frac{k - \varepsilon}{\sqrt{4\varepsilon - (1 - a)(k + \varepsilon)^2}} $$
24 5.2 Proof of Theorem 2

where $y(a, \varepsilon) = (k + \varepsilon)(ak + \beta + I)/k$. Notice that $y(1, \varepsilon) = \hat{y}_1(\varepsilon)$, see (22), and consequently, $y(a, \varepsilon)$ is the ordinate of the intersection point of the stable manifold of $\mathbf{p}_{3,1}$ with the line $\{x = a\}$ for $a \le 1$.

Then, we get a system of 6 equations with 8 unknowns grouped as $(\tau_1, \tau_2, \tau_3, y_1, y_2, I, a, \varepsilon)$. Since $(\varepsilon, I(\varepsilon), \tau(\varepsilon))$ is the solution of system (25) obtained in Lemma 3, we note that, $(0, \tau(\varepsilon), 0, 0, \hat{y}_1(\varepsilon), I(\varepsilon), 1, \varepsilon)$, is a solution of system (29). The following result asserts that this solution can be extended for $1 - a > 0$ small enough, providing the existence of the homoclinic loop.

**Proposition 1** For each $0 < \varepsilon_0 < k$, let $\tau(\varepsilon_0)$ and $I(\varepsilon_0)$ be the functions given in Lemma 3, and $\hat{y}_1(\varepsilon_0)$ the function defined in (22). Therefore, there exist neighbourhoods $U$ of $(1, \varepsilon_0)$ and $V$ of $(0, \tau(\varepsilon_0), 0, 0, \hat{y}_1(\varepsilon_0), I(\varepsilon_0))$ and differentiable functions $\tau_1, \tau_2, \tau_3, Y_1, Y_2$ and $I$ defined from $U$ to $V$ such that $\tau_1(1, \varepsilon_0) = \tau_3(1, \varepsilon_0) = 0$, $\tau_2(1, \varepsilon_0) = \tau(\varepsilon_0)$, $Y_1(1, \varepsilon_0) = 0$, $Y_2(1, \varepsilon_0) = \hat{y}_1(\varepsilon_0)$, $I(1, \varepsilon_0) = I(\varepsilon_0)$. Moreover, for every $(a, \varepsilon) \in U$, the vector $(\tau_1(a, \varepsilon), \tau_2(a, \varepsilon), \tau_3(a, \varepsilon), Y_1(a, \varepsilon), Y_2(a, \varepsilon), I(a, \varepsilon), a, \varepsilon)$ is a solution of system (29).

*Proof* Let us label by $\phi_i = 0, i = 1, \dots, 6$ the $i$th equation of system (29). Then, from the local flow in region $R_{3,2}$ and $R_{4,3}$ (see Appendix A), it follows that

$$ \phi_1 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\rho\tau_1} \left( \xi_{1,1} \frac{\sin(\theta\tau_1)}{\sqrt{1-a}} - \cos(\theta\tau_1) \right) - 1 + (k + \beta + I) \right], $$

$$ \phi_2 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I) \left( \xi_{2,1} \frac{\sin(\theta\tau_1)}{\sqrt{1-a}} e^{\rho\tau_1} + (1 - e^{\rho\tau_1}) \frac{\cos(\theta\tau_1)}{1-a} \right) - (1 - (1-a)k) \frac{y_1}{1-a} \right], $$

$$ \phi_3 = (2 - \gamma - I) \left( e^{-\frac{\tau_2}{\varepsilon}} - 1 \right) + \frac{1 - y_1}{1 - \varepsilon} \left( e^{-\tau_2} - e^{-\frac{\tau_2}{\varepsilon}} \right), $$

$$ \phi_4 = -(1 - y_1)e^{-\tau_2} + 1 - y_2, $$

$$ \phi_5 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\rho\tau_3} \left( \xi_{5,1} \frac{\sin(\theta\tau_3)}{\sqrt{1-a}} - \cos(\theta\tau_3) \right) - 1 + (k + \beta + I) \right], $$

$$ \phi_6 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\rho\tau_3} \left( \xi_{6,1} \sin(\theta\tau_3) + \xi_{6,2} \cos(\theta\tau_3) \right) + \frac{1}{1-a} \left( (ak + \beta + I) - (1-(1-a)k)(y(a, \varepsilon) + y_2) \right) \right], $$

where $\rho = \frac{k-\varepsilon}{2\varepsilon}$, $\theta = \frac{\sqrt{4\varepsilon-(k+\varepsilon)^2}}{2\varepsilon}$ and

$$ \xi_{1,1} = \frac{2 - (1 - a)(k + \varepsilon)}{\sqrt{4\varepsilon - (1 - a)(k + \varepsilon)^2}} $$

$$ \xi_{2,1} = \frac{k - \varepsilon}{\sqrt{4\varepsilon - (1 - a)(k + \varepsilon)^2}} $$
$$
\begin{aligned}
\xi_{5,1} &= \frac{k(1-a)(k+\varepsilon) - 2\varepsilon}{k\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}} \\
\xi_{6,1} &= \frac{-2\varepsilon k + (k^2(1-a) - (1-(1-a)k)\varepsilon)(k+\varepsilon)}{k\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}} \\
\xi_{6,2} &= -\frac{k(1-a)(k+\varepsilon) - \varepsilon}{k(1-a)}.
\end{aligned}
$$

We remark that, in order to simplify the calculus, we have changed the orientation of the solution in $\phi_5$ and $\phi_6$ in such a way that it starts at $(a, y(a, \varepsilon))$ and finishes at $(1, y_2)$. As a result, the variable $\tau_3$ is now a non-positive value.

We notice the functions $\phi_i$ for $i = 1, 2, 5, 6$ are not defined at $a = 1$. Nevertheless, the change of variables

$$ (\tau_1, \tau_3, y_1, y_2) = \left( \frac{\sqrt{1-a}}{\theta} \tilde{\tau}_1, \frac{\sqrt{1-a}}{\theta} \tilde{\tau}_3, (1-a)g_1, y(a, \varepsilon) + (1-a)g_2 \right), \quad (31) $$

allows to express $\phi_i = (1-a)\hat{\phi}_i$ for $i = 1, 2, 5, 6$ and $\phi_i = \hat{\phi}_i$ for $i = 3, 4$, see Appendix C. Then, the system of equations

$$ \hat{\phi}_i(\tilde{\tau}_1, \tau_2, \tilde{\tau}_3, g_1, g_2, I, a, \varepsilon) = 0, \quad i = 1, 2, \dots, 6, \quad (32) $$

is equivalent to system (29) when $a < 1$. Additionally, system (32) is defined at $a = 1$ and it can be checked that the vector $v_0 = (0, \tau(\varepsilon_0), 0, 0, \hat{y}_1(\varepsilon_0), I(\varepsilon_0), 1, \varepsilon_0)$, where $0 < \varepsilon_0 < k < 1$ and $\tau(\varepsilon_0), I(\varepsilon_0)$ and $\hat{y}_1(\varepsilon_0)$ are the ones provided in Lemma 3, is a solution of (32).

We now consider the Jacobian matrix of system (32) evaluated at the solution $v_0$. In Appendix C we show that if $0 < \varepsilon_0 < k < 1$, the minor of the Jacobian matrix, obtained by omitting the columns corresponding with the derivative with respect to $a$ and to $\varepsilon$, has maximum range. By applying the Implicit Function theorem, it follows that there exist an open set $U \subset \mathbb{R}^2$ containing $(1, \varepsilon_0)$, an open set $V \subset \mathbb{R}^6$ containing $p_0 = (0, \tau(\varepsilon_0), 0, 0, \hat{y}_1(\varepsilon_0), I(\varepsilon_0))$ and a function $h: U \to V$ such that $h(1, \varepsilon_0) = p_0$ and the set of equalities $\hat{\phi}_i(h(a, \varepsilon), a, \varepsilon) = 0$ with $i = 1, 2, \dots, 6$ holds for every $(a, \varepsilon) \in U$. $\square$

From Proposition 1 we conclude for $0 < \varepsilon < k$, $0 < 1 - a$ small enough and $I = I(a, \varepsilon)$, the existence of a homoclinic connection to the saddle point $p_{3,1}$ which crosses the line $\{x = 1\}$ through the points $(1, Y_1(a, \varepsilon))$ and $(1, Y_2(a, \varepsilon))$. We note that, Proposition 1 also assures about the existence of the algebraic connection for $a > 1$; however, following the definition of the vector field, this scenario can only provide spurious solutions.

We end this case by showing that this homoclinic connection also satisfies a non degeneracy condition similar to that proved for the homoclinic connection when $a = 1$ (Lemma 4). Let $\tilde{y}_u(a, \varepsilon, I)$ and $\tilde{y}_s(a, \varepsilon, I)$ be the second coordinates of the first intersection points of the unstable manifold $W^u(p_{3,1})$ and the stable manifold $W^s(p_{3,1})$, respectively, with the line $\{x = 1\}$, and consider the function $\Lambda(a, \varepsilon, I) = \tilde{y}_s(a, \varepsilon, I) - \tilde{y}_u(a, \varepsilon, I)$.

**Lemma 5** Let us consider system (2) for $0 < \varepsilon < k$ and $0 < 1 - a$ sufficiently small, and let $\Lambda(a, \varepsilon, I) = \tilde{y}_s(a, \varepsilon, I) - \tilde{y}_u(a, \varepsilon, I)$. Then $\Lambda(a, \varepsilon, I(a, \varepsilon)) = 0$ and $\frac{\partial \Lambda}{\partial I}\bigg|_{(a, \varepsilon, I(a, \varepsilon))} > 0$. where $I(a, \varepsilon)$ is the function obtained in Proposition 1.
$$
\begin{align*}
\xi_{5,1} &= \frac{k(1-a)(k+\varepsilon) - 2\varepsilon}{k\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}} \\
\xi_{6,1} &= \frac{-2\varepsilon k + (k^2(1-a) - (1-(1-a)k)\varepsilon)(k+\varepsilon)}{k\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}} \\
\xi_{6,2} &= -\frac{k(1-a)(k+\varepsilon) - \varepsilon}{k(1-a)}.
\end{align*}
$$

We remark that, in order to simplify the calculus, we have changed the orientation of the solution in $\phi_5$ and $\phi_6$ in such a way that it starts at $(a, y(a, \varepsilon))$ and finishes at $(1, y_2)$. As a result, the variable $\tau_3$ is now a non-positive value.

We notice the functions $\phi_i$ for $i = 1, 2, 5, 6$ are not defined at $a = 1$. Nevertheless, the change of variables

$$ (\tau_1, \tau_3, y_1, y_2) = \left( \frac{\sqrt{1-a}}{\theta} \tilde{\tau}_1, \frac{\sqrt{1-a}}{\theta} \tilde{\tau}_3, (1-a)g_1, y(a, \varepsilon) + (1-a)g_2 \right), \quad (31) $$

allows to express $\phi_i = (1-a)\hat{\phi}_i$ for $i = 1, 2, 5, 6$ and $\phi_i = \hat{\phi}_i$ for $i = 3, 4$, see Appendix C. Then, the system of equations

$$ \hat{\phi}_i(\tilde{\tau}_1, \tau_2, \tilde{\tau}_3, g_1, g_2, I, a, \varepsilon) = 0, \quad i = 1, 2, \dots, 6, \quad (32) $$

is equivalent to system (29) when $a < 1$. Additionally, system (32) is defined at $a = 1$ and it can be checked that the vector $v_0 = (0, \tau(\varepsilon_0), 0, 0, \hat{y}_1(\varepsilon_0), I(\varepsilon_0), 1, \varepsilon_0)$, where $0 < \varepsilon_0 < k < 1$ and $\tau(\varepsilon_0), I(\varepsilon_0)$ and $\hat{y}_1(\varepsilon_0)$ are the ones provided in Lemma 3, is a solution of (32).

We now consider the Jacobian matrix of system (32) evaluated at the solution $v_0$. In Appendix C we show that if $0 < \varepsilon_0 < k < 1$, the minor of the Jacobian matrix, obtained by omitting the columns corresponding with the derivative with respect to $a$ and to $\varepsilon$, has maximum range. By applying the Implicit Function theorem, it follows that there exist an open set $U \subset \mathbb{R}^2$ containing $(1, \varepsilon_0)$, an open set $V \subset \mathbb{R}^6$ containing $p_0 = (0, \tau(\varepsilon_0), 0, 0, \hat{y}_1(\varepsilon_0), I(\varepsilon_0))$ and a function $h: U \to V$ such that $h(1, \varepsilon_0) = p_0$ and the set of equalities $\hat{\phi}_i(h(a, \varepsilon), a, \varepsilon) = 0$ with $i = 1, 2, \dots, 6$ holds for every $(a, \varepsilon) \in U$. $\square$

From Proposition 1 we conclude for $0 < \varepsilon < k$, $0 < 1 - a$ small enough and $I = I(a, \varepsilon)$, the existence of a homoclinic connection to the saddle point $p_{3,1}$ which crosses the line $\{x = 1\}$ through the points $(1, Y_1(a, \varepsilon))$ and $(1, Y_2(a, \varepsilon))$. We note that, Proposition 1 also assures about the existence of the algebraic connection for $a > 1$; however, following the definition of the vector field, this scenario can only provide spurious solutions.

We end this case by showing that this homoclinic connection also satisfies a non degeneracy condition similar to that proved for the homoclinic connection when $a = 1$ (Lemma 4). Let $\tilde{y}_u(a, \varepsilon, I)$ and $\tilde{y}_s(a, \varepsilon, I)$ be the second coordinates of the first intersection points of the unstable manifold $W^u(p_{3,1})$ and the stable manifold $W^s(p_{3,1})$, respectively, with the line $\{x = 1\}$, and consider the function $\Lambda(a, \varepsilon, I) = \tilde{y}_s(a, \varepsilon, I) - \tilde{y}_u(a, \varepsilon, I)$.

**Lemma 5** Let us consider system (2) for $0 < \varepsilon < k$ and $0 < 1 - a$ sufficiently small, and let $\Lambda(a, \varepsilon, I) = \tilde{y}_s(a, \varepsilon, I) - \tilde{y}_u(a, \varepsilon, I)$. Then $\Lambda(a, \varepsilon, I(a, \varepsilon)) = 0$ and $\frac{\partial \Lambda}{\partial I}\bigg|_{(a, \varepsilon, I(a, \varepsilon))} > 0$. where $I(a, \varepsilon)$ is the function obtained in Proposition 1.
26

5.3 Proof of Theorem 3

*Proof* Following the expressions appearing in Appendix C, let us consider

$$
\begin{align*}
\tilde{y}_s &= \hat{y}_1(\varepsilon) - (1-a)(k + \beta + I)\tau\sqrt{\varepsilon}(\sqrt{\varepsilon\tau} + 2(k + \varepsilon))/(2k) + O(1-a)^2, \\
\tilde{y}_u &= y_u + (1-a)g_1 e^{-\tau},
\end{align*}
$$

where $y_u$ is given by Lemma 4, and $\tilde{y}_u, \tilde{y}_s$ are taken by isolating $\hat{y}_1(\varepsilon)+(1-a)g_2$ from $\hat{\phi}_4$ and $\hat{\phi}_6$, respectively (see Appendix C). Therefore, as we calculate the derivative with respect $I$ and then we evaluate at the connection point, it follows that $\tilde{y}'_s = \hat{y}'_1(\varepsilon) = y'(1, \varepsilon)$, $\tilde{y}'_u = y'_u$ and hence, $\frac{\partial \Lambda}{\partial I}\bigg|_{(a, \varepsilon, I(a, \varepsilon))} = \hat{y}'_1(\varepsilon) - y'_u + O(1-a)$. From Lemma 4, we know that $\hat{y}'_1(\varepsilon) - y'_u > 0$. Hence, in a neighbour of $a = 1$, this derivative remains positive. $\square$

**Lemma 6** Let us consider system (2) with $0 < \varepsilon < k$ and $0 < 1 - a$ sufficiently small. For $0 < I - I(a, \varepsilon)$ small enough, there exists a stable limit cycle located in the interior of the region limited by $W^u(\mathbf{p}_{3,1})$ and $W^s(\mathbf{p}_{3,1})$.

*Proof* The statement follows as the equilibrium point on the region $R_{3,2}$ is a repelling point and, considering Lemma 5, there must be some limit object near that point. $\square$

Theorem 2 follows then as a direct consequence of the Proposition 1 and Lemma 6.

5.3 Proof of Theorem 3

Theorem 3 deals with the existence of orbits of system (7) connecting the unstable and the stable manifold of the slow manifold $S_{\sigma}^{sd}$, i.e. $W^u(S_{\sigma}^{sd})$ and $W^s(S_{\sigma}^{sd})$, respectively. Before dealing with the proof of this theorem, we provide local explicit expressions for the invariant sets $W^u(S_{\sigma}^{sd})$ and $W^s(S_{\sigma}^{sd})$.

**Lemma 7** Let us consider system (7) restricted to the region $R_{3,1}$. The segment $S_{\sigma}^{sd}$, given by equation (12), is a slow manifold of saddle type. The piece of plane $W^s(S_{\sigma}^{sd})$, given by equation (13), is a repelling invariant manifold locally containing the stable manifold of $S_{\sigma}^{sd}$, and the piece of plane $W^u(S_{\sigma}^{sd})$, given by equation (14), is an attracting invariant manifold locally containing the unstable manifold of $S_{\sigma}^{sd}$.

*Proof* The local expression of the orbits of system (7) in region $R_{3,1}$ is given by

$$
\begin{align*}
x(t) &= -\frac{\beta + k\sigma\varepsilon + k(\sigma t + C_3)}{k^2} + \left(C_1 - \frac{C_2}{k + \varepsilon} + \frac{I_0}{k}\right)e^{kt/\varepsilon} + \frac{C_2}{k + \varepsilon}e^{-t} \\
y(t) &= C_2e^{-t} \\
I(t) &= \sigma t + C_3
\end{align*}
$$

By taking $C_1 - \frac{C_2}{k+\varepsilon} + \frac{I_0}{k} = 0$, we can get the expression of the stable invariant plane $W^s(S_{\sigma}^{sd})$ by substituting $y(t)$ and $I(t)$ in the expression of $x(t)$. As the motion is given in linear terms and by a negative exponential function, it follows to be the
*Proof* Following the expressions appearing in Appendix C, let us consider

$$ \tilde{y}_s = \hat{y}_1(\varepsilon) - (1-a)(k + \beta + I)\tau\sqrt{\varepsilon}(\sqrt{\varepsilon\tau + 2(k + \varepsilon)})/(2k) + O(1-a)^2, \\ \tilde{y}_u = y_u + (1-a)g_1e^{-\tau_2}, $$

where $y_u$ is given by Lemma 4, and $\tilde{y}_u, \tilde{y}_s$ are taken by isolating $\hat{y}_1(\varepsilon)+(1-a)g_2$ from $\hat{\phi}_4$ and $\hat{\phi}_6$, respectively (see Appendix C). Therefore, as we calculate the derivative with respect $I$ and then we evaluate at the connection point, it follows that $\tilde{y}'_s = \hat{y}'_1(\varepsilon) = y'(1, \varepsilon), \tilde{y}'_u = y'_u$ and hence, $\frac{\partial \Lambda}{\partial I}|_{(a,\varepsilon, I(a,\varepsilon))} = \hat{y}'_1(\varepsilon) - y'_u + O(1-a)$. From Lemma 4, we know that $\hat{y}'_1(\varepsilon) - y'_u > 0$. Hence, in a neighbour of $a = 1$, this derivative remains positive. $\square$

**Lemma 6** Let us consider system (2) with $0 < \varepsilon < k$ and $0 < 1 - a$ sufficiently small. For $0 < I - I(a, \varepsilon)$ small enough, there exists a stable limit cycle located in the interior of the region limited by $W^u(p_{3,1})$ and $W^s(p_{3,1})$.

*Proof* The statement follows as the equilibrium point on the region $R_{3,2}$ is a repelling point and, considering Lemma 5, there must be some limit object near that point. $\square$

Theorem 2 follows then as a direct consequence of the Proposition 1 and Lemma 6.

## 5.3 Proof of Theorem 3

Theorem 3 deals with the existence of orbits of system (7) connecting the unstable and the stable manifold of the slow manifold $S_{\sigma}^{sd}$, i.e. $W^u(S_{\sigma}^{sd})$ and $W^s(S_{\sigma}^{sd})$, respectively. Before dealing with the proof of this theorem, we provide local explicit expressions for the invariant sets $W^u(S_{\sigma}^{sd})$ and $W^s(S_{\sigma}^{sd})$.

**Lemma 7** Let us consider system (7) restricted to the region $R_{3,1}$. The segment $S_{\sigma}^{sd}$, given by equation (12), is a slow manifold of saddle type. The piece of plane $W^s(S_{\sigma}^{sd})$, given by equation (13), is a repelling invariant manifold locally containing the stable manifold of $S_{\sigma}^{sd}$, and the piece of plane $W^u(S_{\sigma}^{sd})$, given by equation (14), is an attracting invariant manifold locally containing the unstable manifold of $S_{\sigma}^{sd}$.

*Proof* The local expression of the orbits of system (7) in region $R_{3,1}$ is given by

$$ \begin{aligned} x(t) &= -\frac{\beta + k\sigma\varepsilon + k(\sigma t + C_3)}{k^2} + (C_1 - \frac{C_2}{k + \varepsilon} + \frac{I_0}{k})e^{kt/\varepsilon} + \frac{C_2}{k + \varepsilon}e^{-t} \\ y(t) &= C_2e^{-t} \\ I(t) &= \sigma t + C_3 \end{aligned} $$

By taking $C_1 - \frac{C_2}{k+\varepsilon} + \frac{I_0}{k} = 0$, we can get the expression of the stable invariant plane $W^s(S_{\sigma}^{sd})$ by substituting $y(t)$ and $I(t)$ in the expression of $x(t)$. As the motion is given in linear terms and by a negative exponential function, it follows to be the
attracting one. Similarly, by taking $C_2 = 0$, $y(t) = 0$ and therefore $W^u(S_\sigma^{sd})$ become an invariant plane. This plane has a repelling trajectory given by the exponential term of $x(t)$. Finally, by intersecting both planes, we obtain the expression of $S_\sigma^{sd}$.

□

Systems of equations (15) and (16) provide the necessary and sufficient conditions to assure the existence of a maximal canard with one or two spikes, respectively, in the discontinuous scenario $a = 1$. Starting at a common solution of both systems given by the homoclinic connection when $\sigma = 0$, in the following lemmas we obtain respective $\sigma$-families of solutions by a suitable application of the Implicit Function Theorem. Let us start with the 1-spike maximal canard family.

**Lemma 8** For $a = 1$ and $0 < \varepsilon < k$, consider the values $I(\varepsilon)$, $\tau(\varepsilon)$ obtained in Lemma 3. Then, the vector $(\sigma, y_1, z_1, \tau_1) = (0, \hat{y}_1(\varepsilon), I(\varepsilon), \tau(\varepsilon))$ is a solution of the first equation in (15). Moreover, there exist a neighbourhood $U$ of the origin, a neighbourhood $V$ of $(\hat{y}_1(\varepsilon), I(\varepsilon), \tau(\varepsilon))$, and a differentiable function $M_1: U \to V$ such that $(\sigma, M_1(\sigma))$ is a solution of the first equation in (15) for every $\sigma \in U$.

*Proof* Let us consider the equations concerning system (15), which can be written as

$$ (2 - \gamma + \sigma\varepsilon - z_1)(e^{-\tau_1/\varepsilon} - 1) - \frac{1}{1-\varepsilon}(e^{-\tau_1/\varepsilon} - e^{-\tau_1}) + \sigma\tau_1 = 0, \\ 1 - e^{-\tau_1} - y_1 = 0, \quad (33) \\ \sigma\tau_1 + z_1 + k - \frac{k}{k+\varepsilon}y_1 + \frac{\sigma\varepsilon + k\beta}{k} = 0. $$

Notice that, for $\sigma = 0$ we have that the 1st and 2nd equations of system (33) coincide with the equations of system (25), considering $z_1 = I$ and $y_1 = \frac{k+\varepsilon}{k}(k + \beta + I)$. In addition, the 3rd equation is also satisfied. As a result, the vector $\mathbf{p}_0 = (0, \hat{y}_1(\varepsilon), I(\varepsilon), \tau(\varepsilon))$, formed by the functions obtained in Lemma 3, provides values for the unknowns $(\sigma, y_1, z_1, \tau_1)$ satisfying system (33), for $\tau(\varepsilon) > 0$. To assure the existence of a curve of solutions, let us consider the Jacobian matrix of the system evaluated on the theoretical solution

$$ D_J(\mathbf{p}_0) = \begin{pmatrix} \varepsilon(e^{-\tau(\varepsilon)/\varepsilon} - 1) + \tau(\varepsilon) & 0 & 1 - e^{-\tau(\varepsilon)/\varepsilon} & \frac{\gamma+I(\varepsilon)-(1+\hat{y}_1(\varepsilon))}{\varepsilon} \\ 0 & -1 & 0 & 1 - \hat{y}_1(\varepsilon) \\ \tau(\varepsilon) + \frac{\varepsilon}{k} & \frac{-k}{k+\varepsilon} & 1 & 0 \end{pmatrix}. $$

Considering the 2nd, 3rd and 4th column of $D_J(\mathbf{p}_0)$, the determinant of this submatrix is given by $\frac{k}{k+\varepsilon}(\frac{1}{k}(1 - e^{-\tau(\varepsilon)}) - e^{-\tau(\varepsilon)}(1 - e^{-\tau(\varepsilon)/\varepsilon}))$, which corresponds to equation (28). Then, considering also system (33) and joining into one equation, we get equation (27). Therefore, we are under conditions of Lemma 2.

Since the range of the Jacobian matrix is maximal for the 2nd, 3rd and 4th columns, by the Implicit Function Theorem there exist an open set $U \subset \mathbb{R}$ containing $\sigma = 0$, an open set $V$ containing $(\hat{y}_1(\varepsilon), I(\varepsilon), \tau(\varepsilon))$, and a differentiable function $M_1: U \to V$ such that $(\sigma, M_1(\sigma))$ is a solution of the first equation in (15).

□
attracting one. Similarly, by taking $C_2 = 0$, $y(t) = 0$ and therefore $W^u(S_\sigma^{sd})$ become an invariant plane. This plane has a repelling trajectory given by the exponential term of $x(t)$. Finally, by intersecting both planes, we obtain the expression of $S_\sigma^{sd}$.

□

Systems of equations (15) and (16) provide the necessary and sufficient conditions to assure the existence of a maximal canard with one or two spikes, respectively, in the discontinuous scenario $a = 1$. Starting at a common solution of both systems given by the homoclinic connection when $\sigma = 0$, in the following lemmas we obtain respective $\sigma$-families of solutions by a suitable application of the Implicit Function Theorem. Let us start with the 1-spike maximal canard family.

**Lemma 8** For $a = 1$ and $0 < \varepsilon < k$, consider the values $I(\varepsilon)$, $\tau(\varepsilon)$ obtained in Lemma 3. Then, the vector $(\sigma, y_1, z_1, \tau_1) = (0, \hat{y}_1(\varepsilon), I(\varepsilon), \tau(\varepsilon))$ is a solution of the first equation in (15). Moreover, there exist a neighbourhood $U$ of the origin, a neighbourhood $V$ of $(\hat{y}_1(\varepsilon), I(\varepsilon), \tau(\varepsilon))$, and a differentiable function $M_1: U \to V$ such that $(\sigma, M_1(\sigma))$ is a solution of the first equation in (15) for every $\sigma \in U$.

*Proof* Let us consider the equations concerning system (15), which can be written as

$$ (2 - \gamma + \sigma\varepsilon - z_1)(e^{-\tau_1/\varepsilon} - 1) - \frac{1}{1-\varepsilon}(e^{-\tau_1/\varepsilon} - e^{-\tau_1}) + \sigma\tau_1 = 0, \\ 1 - e^{-\tau_1} - y_1 = 0, \quad (33) \\ \sigma\tau_1 + z_1 + k - \frac{k}{k+\varepsilon}y_1 + \frac{\sigma\varepsilon + k\beta}{k} = 0. $$

Notice that, for $\sigma = 0$ we have that the 1st and 2nd equations of system (33) coincide with the equations of system (25), considering $z_1 = I$ and $y_1 = \frac{k+\varepsilon}{k}(k + \beta + I)$. In addition, the 3rd equation is also satisfied. As a result, the vector $p_0 = (0, \hat{y}_1(\varepsilon), I(\varepsilon), \tau(\varepsilon))$, formed by the functions obtained in Lemma 3, provides values for the unknowns $(\sigma, y_1, z_1, \tau_1)$ satisfying system (33), for $\tau(\varepsilon) > 0$. To assure the existence of a curve of solutions, let us consider the Jacobian matrix of the system evaluated on the theoretical solution

$$ D_J(\mathbf{p}_0) = \begin{pmatrix} \varepsilon(e^{-\tau(\varepsilon)/\varepsilon} - 1) + \tau(\varepsilon) & 0 & 1 - e^{-\tau(\varepsilon)/\varepsilon} & \frac{\gamma+I(\varepsilon)-(1+\hat{y}_1(\varepsilon))}{\varepsilon} \\ 0 & -1 & 0 & 1 - \hat{y}_1(\varepsilon) \\ \tau(\varepsilon) + \frac{\varepsilon}{k} & \frac{-k}{k+\varepsilon} & 1 & 0 \end{pmatrix}. $$

Considering the 2nd, 3rd and 4th column of $D_J(p_0)$, the determinant of this submatrix is given by $\frac{k}{k+\varepsilon}(\frac{1}{k}(1 - e^{-\tau(\varepsilon)}) - e^{-\tau(\varepsilon)}(1 - e^{-\tau(\varepsilon)/\varepsilon}))$, which corresponds to equation (28). Then, considering also system (33) and joining into one equation, we get equation (27). Therefore, we are under conditions of Lemma 2.

Since the range of the Jacobian matrix is maximal for the 2nd, 3rd and 4th columns, by the Implicit Function Theorem there exist an open set $U \subset \mathbb{R}$ containing $\sigma = 0$, an open set $V$ containing $(\hat{y}_1(\varepsilon), I(\varepsilon), \tau(\varepsilon))$, and a differentiable function $M_1: U \to V$ such that $(\sigma, M_1(\sigma))$ is a solution of the first equation in (15).

□
28 5.3 Proof of Theorem 3

**Lemma 9** For $a = 1$ and $0 < \varepsilon < k$, consider the solution $(\sigma, M_1(\sigma))$ of the first equation in (15) obtained in Lemma 8. Then, $(\sigma, M_1(\sigma))$ also satisfies the restriction in (15).

*Proof* Let us consider the vector field (7) with $a = 1$ and $\alpha(0)$. If we reconsider the vector field as the extension of the vector field given in the linearity region $R_{4,3}$, we note that: first, $\{y = 1\}$ is an invariant line. Thus, any orbit starting bellow $\{y = 1\}$ stays behind it. Second, $y' > 0$ for any orbit starting with $y < 1$, since $y' = 1 - y$. Finally, since $y$ is increasing, all orbits that cross $\{x = 1\}$ from $x > 1$ to $x < 1$, they will not cross again the line, since at the segment $\{x = 1\}$, $x' = 1/\varepsilon(k + \beta + I - y)$, so for $y > k + \beta + I$, it must happen $x' < 0$. As a result, any extension of the theoretical connection given in the previous result must keep at $R_{4,3}$. $\square$

From Lemmas 8 and 9, and by considering function $I_1(\sigma)$ as the third coordinate of the function $M_1(\sigma)$, we conclude statement (a) of Theorem 3. Next we deal with the proof of the statement (b) in Theorem 3.

**Lemma 10** For $a = 1$ and $0 < \varepsilon < k$, consider the values $I(\varepsilon)$, $\tau(\varepsilon)$, and $\hat{y}_1(\varepsilon)$ obtained in Lemma 3. Then, the vector $(\sigma, y_1, y_2, y_3, z_1, z_2, z_3, \tau_1, \tilde{\tau}_2, \tau_3) = (0, \hat{y}_1(\varepsilon), 0, \hat{y}_1(\varepsilon), I(\varepsilon), I(\varepsilon), I(\varepsilon), \tau(\varepsilon), 0, \tau(\varepsilon))$, where $\tilde{\tau}_2 = \tau_2^{-1}$, is a solution of system (16). Moreover, there exist a neighbourhood $U$ of $\sigma = 0$, a neighbourhood $V$ of $(\hat{y}_1(\varepsilon), 0, \hat{y}_1(\varepsilon), I(\varepsilon), I(\varepsilon), I(\varepsilon), \tau(\varepsilon), 0, \tau(\varepsilon))$, and a differentiable function $M_2: U \to V$ such that $(\sigma, M_2(\sigma))$ is a solution of system (16) for every $\sigma \in U$.

*Proof* Considering the local expressions of the flow given in Appendix A, we rewrite system (16) as a system of 9 equations. Then, changing $\tau_2$ by $\tilde{\tau}_2 = (\tau_2)^{-1}$, imposing the connection for $\sigma = 0$ and to prevent some singularities in the Jacobian matrix, we multiply the fourth equation by $e^{-k/(\varepsilon\tilde{\tau}_2)}$ and the sixth equation by $\tilde{\tau}_2$, obtaining the following set of equations

$$ 
\begin{aligned}
& (2 - \gamma + \sigma\varepsilon - z_1)(e^{-\frac{\tau_1}{\varepsilon}} - 1) - \frac{1}{1-\varepsilon}(e^{-\frac{\tau_1}{\varepsilon}} - e^{-\tau_1}) + \sigma\tau_1 = 0 \\
& 1 - e^{-\tau_1} - y_1 = 0 \\
& \sigma\tau_1 + z_1 - z_2 = 0 \\
& \left(\frac{\sigma\varepsilon + k(z_2 + \beta)}{k^2} + 1\right)(1 - e^{-\frac{k}{\varepsilon\tilde{\tau}_2}}) + \frac{y_1}{k+\varepsilon}(e^{-\frac{1}{\tilde{\tau}_2}} - \frac{k}{\varepsilon\tilde{\tau}_2} - 1) - \frac{\sigma}{k}\frac{e^{-\frac{k}{\varepsilon\tilde{\tau}_2}}}{\tilde{\tau}_2} = 0 \\
& y_1e^{-\frac{1}{\tilde{\tau}_2}} - y_2 = 0 \quad (34) \\
& \sigma + \tilde{\tau}_2(z_2 - z_3) = 0 \\
& (2 - \gamma + \sigma\varepsilon - z_3)(e^{-\frac{\tau_3}{\varepsilon}} - 1) + \frac{y_2-1}{1-\varepsilon}(e^{-\frac{\tau_3}{\varepsilon}} - e^{-\tau_3}) + \sigma\tau_3 = 0 \\
& 1 + (y_2-1)e^{-\tau_3} - y_3 = 0 \\
& \sigma\tau_3 + z_3 + k - \frac{k}{k+\varepsilon}y_3 + \frac{\sigma\varepsilon + k\beta}{k} = 0.
\end{aligned}
$$

We denote the left term of each equation by $\phi_i$, being the $i$'th equation. In order to apply the Implicit Function Theorem to prove the existence of a differentiable
**Lemma 9** For $a = 1$ and $0 < \varepsilon < k$, consider the solution $(\sigma, M_1(\sigma))$ of the first equation in (15) obtained in Lemma 8. Then, $(\sigma, M_1(\sigma))$ also satisfies the restriction in (15).

*Proof* Let us consider the vector field (7) with $a = 1$ and $\alpha(0)$. If we reconsider the vector field as the extension of the vector field given in the linearity region $R_{4,3}$, we note that: first, $\{y = 1\}$ is an invariant line. Thus, any orbit starting bellow $\{y = 1\}$ stays behind it. Second, $y' > 0$ for any orbit starting with $y < 1$, since $y' = 1 - y$. Finally, since $y$ is increasing, all orbits that cross $\{x = 1\}$ from $x > 1$ to $x < 1$, they will not cross again the line, since at the segment $\{x = 1\}$, $x' = 1/\varepsilon(k + \beta + I - y)$, so for $y > k + \beta + I$, it must happen $x' < 0$. As a result, any extension of the theoretical connection given in the previous result must keep at $R_{4,3}$. $\square$

From Lemmas 8 and 9, and by considering function $I_1(\sigma)$ as the third coordinate of the function $M_1(\sigma)$, we conclude statement (a) of Theorem 3. Next we deal with the proof of the statement (b) in Theorem 3.

**Lemma 10** For $a = 1$ and $0 < \varepsilon < k$, consider the values $I(\varepsilon)$, $\tau(\varepsilon)$, and $\hat{y}_1(\varepsilon)$ obtained in Lemma 3. Then, the vector $(\sigma, y_1, y_2, y_3, z_1, z_2, z_3, \tau_1, \tilde{\tau}_2, \tau_3) = (0, \hat{y}_1(\varepsilon), 0, \hat{y}_1(\varepsilon), I(\varepsilon), I(\varepsilon), I(\varepsilon), \tau(\varepsilon), 0, \tau(\varepsilon))$, where $\tilde{\tau}_2 = \tau_2^{-1}$, is a solution of system (16). Moreover, there exist a neighbourhood $U$ of $\sigma = 0$, a neighbourhood $V$ of $(\hat{y}_1(\varepsilon), 0, \hat{y}_1(\varepsilon), I(\varepsilon), I(\varepsilon), I(\varepsilon), \tau(\varepsilon), 0, \tau(\varepsilon))$, and a differentiable function $M_2: U \to V$ such that $(\sigma, M_2(\sigma))$ is a solution of system (16) for every $\sigma \in U$.

*Proof* Considering the local expressions of the flow given in Appendix A, we rewrite system (16) as a system of 9 equations. Then, changing $\tau_2$ by $\tilde{\tau}_2 = (\tau_2)^{-1}$, imposing the connection for $\sigma = 0$ and to prevent some singularities in the Jacobian matrix, we multiply the fourth equation by $e^{-k/(\varepsilon\tilde{\tau}_2)}$ and the sixth equation by $\tilde{\tau}_2$, obtaining the following set of equations

$$ 
\begin{aligned}
& (2 - \gamma + \sigma\varepsilon - z_1)(e^{-\frac{\tau_1}{\varepsilon}} - 1) - \frac{1}{1-\varepsilon}(e^{-\frac{\tau_1}{\varepsilon}} - e^{-\tau_1}) + \sigma\tau_1 = 0 \\
& \qquad 1 - e^{-\tau_1} - y_1 = 0 \\
& \qquad \sigma\tau_1 + z_1 - z_2 = 0 \\
& \left( \frac{\sigma\varepsilon + k(z_2 + \beta)}{k^2} + 1 \right) (1 - e^{-\frac{k}{\varepsilon\tilde{\tau}_2}}) + \frac{y_1}{k+\varepsilon} \left( e^{-\frac{1}{\tilde{\tau}_2}} - \frac{k}{\varepsilon\tilde{\tau}_2} - 1 \right) - \frac{\sigma}{k} \frac{e^{-\frac{k}{\varepsilon\tilde{\tau}_2}}}{\tilde{\tau}_2} = 0 \\
& \qquad y_1 e^{-\frac{1}{\tilde{\tau}_2}} - y_2 = 0 \quad (34) \\
& \qquad \sigma + \tilde{\tau}_2(z_2 - z_3) = 0 \\
& (2 - \gamma + \sigma\varepsilon - z_3)(e^{-\frac{\tau_3}{\varepsilon}} - 1) + \frac{y_2 - 1}{1 - \varepsilon} (e^{-\frac{\tau_3}{\varepsilon}} - e^{-\tau_3}) + \sigma\tau_3 = 0 \\
& \qquad 1 + (y_2 - 1)e^{-\tau_3} - y_3 = 0 \\
& \qquad \sigma\tau_3 + z_3 + k - \frac{k}{k+\varepsilon}y_3 + \frac{\sigma\varepsilon + k\beta}{k} = 0.
\end{aligned}
$$

We denote the left term of each equation by $\phi_i$, being the $i$'th equation. In order to apply the Implicit Function Theorem to prove the existence of a differentiable
curve $M_2$ depending on $\tilde{\tau}_2$, we need to consider a subsystem with a subset of variables which Jacobian matrix presents a non-zero determinant. To do that, we proceed as follows. Computing the partial derivatives of all $\phi_i$ with respect to the different variables and evaluating them to the connection point $z_1 = z_2 = z_3 = I(\varepsilon)$, $\tau_1 = \tau_3 = \tau(\varepsilon)$, $\tilde{\tau}_2 = 0$, $y_2 = 0$, $y_1 = y_3 = \hat{y}_1(\varepsilon)$, with $\hat{y}_1(\varepsilon)$ given in (22) (see Appendix B), we obtain that all derivatives with respect to $\tilde{\tau}_2$ are 0. Therefore, $\tilde{\tau}_2$ is removed from the set of variables. Consequently, all reminding partial derivatives of the function $\phi_6$ are zero apart from those corresponding to the derivative with respect to $\sigma$. Leaving aside $\phi_6$ and $\sigma$, the same happens if we look at the partial derivatives of $\phi_5$, where the only non-zero derivative is the one with respect to $y_2$. Consequently, the Jacobian matrix with respect to the variables ($y_1, y_3, z_1, z_2, z_3, \tau_1, \tau_3$) of the resulting system is given by

$$ Df = \begin{pmatrix} 0 & 0 & 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}} & 0 & 0 & -\frac{\hat{y}_1(\varepsilon)}{k+\varepsilon} & 0 \\ -1 & 0 & 0 & 0 & 0 & e^{-\tau(\varepsilon)} & 0 \\ 0 & 0 & 1 & -1 & 0 & 0 & 0 \\ \frac{-1}{k+\varepsilon} & 0 & 0 & \frac{1}{k} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}} & 0 & -\frac{\hat{y}_1(\varepsilon)}{k+\varepsilon} \\ 0 & -1 & 0 & 0 & 0 & 0 & e^{-\tau(\varepsilon)} \\ 0 & \frac{-k}{k+\varepsilon} & 0 & 0 & 1 & 0 & 0 \end{pmatrix}. $$

To study its determinant, let us consider the two sub-matrices

$$ A_1 = \begin{pmatrix} 0 & 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}} & 0 & -\frac{\hat{y}_1(\varepsilon)}{k+\varepsilon} \\ -1 & 0 & 0 & e^{-\tau(\varepsilon)} \\ 0 & 1 & -1 & 0 \\ \frac{-1}{k+\varepsilon} & 0 & \frac{1}{k} & 0 \end{pmatrix}, \quad A_2 = \begin{pmatrix} 0 & 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}} & -\frac{\hat{y}_1(\varepsilon)}{k+\varepsilon} \\ -1 & 0 & e^{-\tau(\varepsilon)} \\ \frac{-k}{k+\varepsilon} & 1 & 0 \end{pmatrix}. $$

We remark that $\det(Df) \neq 0$ if, and only if, any of $\det(A_1)$, $\det(A_2)$ are non-zero. Further, as we compute the determinant of $A_1$, we can see that $\det(A_1) = -1/k \det(A_2)$.

As a result, it suffices to calculate $\det(A_2)$. However, this matrix is the same as in the case of the maximal canard with a direct connection (see proof of Lemma 8), which provides the 2nd equation of system (25). The first equation of this system is directly provided by the connection equations (34). Therefore, we can apply again Lemma 2 and then the statement follows.

Since all the variables can be written in terms of $\tilde{\tau}_2$, which corresponds to the inverse of the time of flight of the orbits near the saddle point, we can express them in the following way: $\sigma = S\tilde{\tau}_2e^{-\frac{1}{\tilde{\tau}_2}}$, $z_1 = I+Z_1e^{-\frac{1}{\tilde{\tau}_2}}$, $z_2 = I+Z_2e^{-\frac{1}{\tilde{\tau}_2}}$, $z_3 = I+Z_3e^{-\frac{1}{\tilde{\tau}_2}}$, $y_1 = \hat{y}_1 + Y_1e^{-\frac{1}{\tilde{\tau}_2}}$, $y_2 = 0 + Y_2e^{-\frac{1}{\tilde{\tau}_2}}$, $y_3 = \hat{y}_1 + Y_3e^{-\frac{1}{\tilde{\tau}_2}}$, $\tau_1 = \tau(\varepsilon) + T_1e^{-\frac{1}{\tilde{\tau}_2}}$ and $\tau_3 = \tau(\varepsilon) + T_3e^{-\frac{1}{\tilde{\tau}_2}}$. By applying the unknown coefficients method to system (34), we obtain that $S < 0$. Therefore, locally we can express $\tilde{\tau}_2$ as a function of $\sigma$. Consequently, we can express function $M_2$ as a function of $\sigma$, which finishes the proof of the statement. $\square$

**Lemma 11** For $a = 1$ and $0 < \varepsilon < k$, consider the solution $(\sigma, M_2(\sigma))$ of system (16) obtained in Lemma 10. Then, $(\sigma, M_2(\sigma))$ also satisfies the conditions in (17).
curve $M_2$ depending on $\tilde{\tau}_2$, we need to consider a subsystem with a subset of variables which Jacobian matrix presents a non-zero determinant. To do that, we proceed as follows. Computing the partial derivatives of all $\phi_i$ with respect to the different variables and evaluating them to the connection point $z_1 = z_2 = z_3 = I(\varepsilon)$, $\tau_1 = \tau_3 = \tau(\varepsilon)$, $\tilde{\tau}_2 = 0$, $y_2 = 0$, $y_1 = y_3 = \hat{y}_1(\varepsilon)$, with $\hat{y}_1(\varepsilon)$ given in (22) (see Appendix B), we obtain that all derivatives with respect to $\tilde{\tau}_2$ are 0. Therefore, $\tilde{\tau}_2$ is removed from the set of variables. Consequently, all reminding partial derivatives of the function $\phi_6$ are zero apart from those corresponding to the derivative with respect to $\sigma$. Leaving aside $\phi_6$ and $\sigma$, the same happens if we look at the partial derivatives of $\phi_5$, where the only non-zero derivative is the one with respect to $y_2$. Consequently, the Jacobian matrix with respect to the variables ($y_1, y_3, z_1, z_2, z_3, \tau_1, \tau_3$) of the resulting system is given by

$$ Df = \begin{pmatrix} 0 & 0 & 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}} & 0 & 0 & -\frac{\hat{y}_1(\varepsilon)}{k+\varepsilon} & 0 \\ -1 & 0 & 0 & 0 & 0 & e^{-\tau(\varepsilon)} & 0 \\ 0 & 0 & 1 & -1 & 0 & 0 & 0 \\ \frac{-1}{k+\varepsilon} & 0 & 0 & \frac{1}{k} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}} & 0 & -\frac{\hat{y}_1(\varepsilon)}{k+\varepsilon} \\ 0 & -1 & 0 & 0 & 0 & 0 & e^{-\tau(\varepsilon)} \\ 0 & \frac{-k}{k+\varepsilon} & 0 & 0 & 1 & 0 & 0 \end{pmatrix}. $$

To study its determinant, let us consider the two sub-matrices

$$ A_1 = \begin{pmatrix} 0 & 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}} & 0 & -\frac{\hat{y}_1(\varepsilon)}{k+\varepsilon} \\ -1 & 0 & 0 & e^{-\tau(\varepsilon)} \\ 0 & 1 & -1 & 0 \\ \frac{-1}{k+\varepsilon} & 0 & \frac{1}{k} & 0 \end{pmatrix}, \quad A_2 = \begin{pmatrix} 0 & 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}} & -\frac{\hat{y}_1(\varepsilon)}{k+\varepsilon} \\ -1 & 0 & e^{-\tau(\varepsilon)} \\ \frac{-k}{k+\varepsilon} & 1 & 0 \end{pmatrix}. $$

We remark that $\det(Df) \neq 0$ if, and only if, any of $\det(A_1)$, $\det(A_2)$ are non-zero. Further, as we compute the determinant of $A_1$, we can see that $\det(A_1) = -1/k \det(A_2)$.

As a result, it suffices to calculate $\det(A_2)$. However, this matrix is the same as in the case of the maximal canard with a direct connection (see proof of Lemma 8), which provides the 2nd equation of system (25). The first equation of this system is directly provided by the connection equations (34). Therefore, we can apply again Lemma 2 and then the statement follows.

Since all the variables can be written in terms of $\tilde{\tau}_2$, which corresponds to the inverse of the time of flight of the orbits near the saddle point, we can express them in the following way: $\sigma = S\tilde{\tau}_2 e^{-\frac{1}{\tilde{\tau}_2}}$, $z_1 = I+Z_1 e^{-\frac{1}{\tilde{\tau}_2}}$, $z_2 = I+Z_2 e^{-\frac{1}{\tilde{\tau}_2}}$, $z_3 = I+Z_3 e^{-\frac{1}{\tilde{\tau}_2}}$, $y_1 = \hat{y}_1 + Y_1 e^{-\frac{1}{\tilde{\tau}_2}}$, $y_2 = 0 + Y_2 e^{-\frac{1}{\tilde{\tau}_2}}$, $y_3 = \hat{y}_1 + Y_3 e^{-\frac{1}{\tilde{\tau}_2}}$, $\tau_1 = \tau(\varepsilon) + T_1 e^{-\frac{1}{\tilde{\tau}_2}}$ and $\tau_3 = \tau(\varepsilon) + T_3 e^{-\frac{1}{\tilde{\tau}_2}}$. By applying the unknown coefficients method to system (34), we obtain that $S < 0$. Therefore, locally we can express $\tilde{\tau}_2$ as a function of $\sigma$. Consequently, we can express function $M_2$ as a function of $\sigma$, which finishes the proof of the statement. $\square$

**Lemma 11** For $a = 1$ and $0 < \varepsilon < k$, consider the solution $(\sigma, M_2(\sigma))$ of system (16) obtained in Lemma 10. Then, $(\sigma, M_2(\sigma))$ also satisfies the conditions in (17).
30 5.4 Proof of Theorem 4

Proof The 1st and 3rd equations in (16) corresponding to the flow of system (7) belong to region $R_{4,3}$, whilst the 2nd equation belongs to $R_{3,1}$. For the flow lying to region $R_{4,3}$, we note that a similar argument used in Theorem 3 can be applied, since there can only be one intersection at $\{x=1\}$ apart from the starting point. Regarding the flow in $R_{3,1}$, notice that if we extend the vector field to $\mathbb{R}^3$, if the first component vanishes, that is $x' = 0$, then it must happen that, for $\sigma$ sufficiently small, $x'' > 0$ between the invariant planes for $y > 0$. Thus, any solution that crosses the plane $\{x=1\}$ must cross at most twice, satisfying the conditions in (17). $\square$

By considering function $I_2(\sigma)$ as the fourth coordinate of the function $M_2(\sigma)$ we conclude statement (b) of Theorem 3, which finishes the proof of Theorem 3.

## 5.4 Proof of Theorem 4

Equations (20) together with restrictions (21) provide necessary and sufficient conditions for the existence of 1-spike maximal canard orbits of system (7) connecting $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$ when $a \le 1$. In the following Lemma we prove, by using the contraction and expansion technique, the existence of a set of parameters that satisfy system (20). In particular, starting at one solution of (20), provided from the discontinuous case $a = 1$, we obtain a branch of solutions of (20) by a suitable use of the Implicit Function Theorem.

**Lemma 12** For $0 < \varepsilon < k$, consider the values $I(\varepsilon)$, $\tau(\varepsilon)$ and $\hat{y}_1(\varepsilon)$ obtained in Lemma 3. Then, the vector $(a, \sigma, \tau_1, \tau_2, \tau_3, y_1, y_2, y_3, z_1, z_2, z_3) = (1, 0, 0, \tau(\varepsilon), 0, 0, \hat{y}_1(\varepsilon), \hat{y}_1(\varepsilon), I(\varepsilon), I(\varepsilon), I(\varepsilon))$ is a solution of system (20). Moreover, there exist a neighbourhood $U$ of $(a, \sigma) = (1, 0)$, a neighbourhood $V$ of $(\tau_1, \tau_2, \tau_3, y_1, y_2, y_3, z_1, z_2, z_3) = (0, \tau(\varepsilon), 0, 0, \hat{y}_1(\varepsilon), \hat{y}_1(\varepsilon), I(\varepsilon), I(\varepsilon), I(\varepsilon))$, and a differentiable function $M_1: U \to V$ such that $(a, \sigma, M_1(a, \sigma))$ is a solution of system (20) for every $(a, \sigma) \in U$.

*Proof* Let us consider the equations concerning system (20), which are

$$ 
\begin{aligned}
\psi_{1,1} \sin(\theta \tau_1) + \psi_{1,2} \cos(\theta \tau_1) + r \tau_1 + s - 1 &= 0, \\
\psi_{2,1} \sin(\theta \tau_1) + \psi_{2,2} \cos(\theta \tau_1) + (kr + \sigma)\tau_1 + ks - \varepsilon r + \beta + z_0 - y_1 &= 0, \\
\sigma \tau_1 + z_0 - z_1 &= 0,
\end{aligned}
$$

$$ 
\begin{aligned}
(2 - \gamma + \sigma\varepsilon - z_1)(e^{-\tau_2/\varepsilon} - 1) + (y_1 - 1)/(1 - \varepsilon)(e^{-\tau_2/\varepsilon} - e^{-\tau_2}) + \sigma\tau_2 &= 0, \\
1 + (y_1 - 1)^{-\tau_2} - y_2 &= 0, \\
\sigma\tau_2 + z_1 - z_2 &= 0,
\end{aligned}
$$

$$ 
\begin{aligned}
\psi_{5,1} \sin(\theta \tau_3) + \psi_{5,2} \cos(\theta \tau_3) + r \tau_3 + s - 1 &= 0, \\
\psi_{6,1} \sin(\theta \tau_3) + \psi_{6,2} \cos(\theta \tau_3) + (kr + \sigma)\tau_3 + ks - \varepsilon r + \beta + z_0 - y_2 &= 0, \\
\sigma \tau_3 + z_0 - z_1 &= 0,
\end{aligned}
$$
30 5.4 Proof of Theorem 4

Proof The 1st and 3rd equations in (16) corresponding to the flow of system (7) belong to region $R_{4,3}$, whilst the 2nd equation belongs to $R_{3,1}$. For the flow lying to region $R_{4,3}$, we note that a similar argument used in Theorem 3 can be applied, since there can only be one intersection at $\{x=1\}$ apart from the starting point. Regarding the flow in $R_{3,1}$, notice that if we extend the vector field to $\mathbb{R}^3$, if the first component vanishes, that is $x' = 0$, then it must happen that, for $\sigma$ sufficiently small, $x'' > 0$ between the invariant planes for $y > 0$. Thus, any solution that crosses the plane $\{x=1\}$ must cross at most twice, satisfying the conditions in (17). $\square$

By considering function $I_2(\sigma)$ as the fourth coordinate of the function $M_2(\sigma)$ we conclude statement (b) of Theorem 3, which finishes the proof of Theorem 3.

## 5.4 Proof of Theorem 4

Equations (20) together with restrictions (21) provide necessary and sufficient conditions for the existence of 1-spike maximal canard orbits of system (7) connecting $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$ when $a \le 1$. In the following Lemma we prove, by using the contraction and expansion technique, the existence of a set of parameters that satisfy system (20). In particular, starting at one solution of (20), provided from the discontinuous case $a = 1$, we obtain a branch of solutions of (20) by a suitable use of the Implicit Function Theorem.

**Lemma 12** For $0 < \varepsilon < k$, consider the values $I(\varepsilon)$, $\tau(\varepsilon)$ and $\hat{y}_1(\varepsilon)$ obtained in Lemma 3. Then, the vector $(a, \sigma, \tau_1, \tau_2, \tau_3, y_1, y_2, y_3, z_1, z_2, z_3) = (1, 0, 0, \tau(\varepsilon), 0, 0, \hat{y}_1(\varepsilon), \hat{y}_1(\varepsilon), I(\varepsilon), I(\varepsilon), I(\varepsilon))$ is a solution of system (20). Moreover, there exist a neighbourhood $U$ of $(a, \sigma) = (1, 0)$, a neighbourhood $V$ of $(\tau_1, \tau_2, \tau_3, y_1, y_2, y_3, z_1, z_2, z_3) = (0, \tau(\varepsilon), 0, 0, \hat{y}_1(\varepsilon), \hat{y}_1(\varepsilon), I(\varepsilon), I(\varepsilon), I(\varepsilon))$, and a differentiable function $M_1: U \to V$ such that $(a, \sigma, M_1(a, \sigma))$ is a solution of system (20) for every $(a, \sigma) \in U$.

*Proof* Let us consider the equations concerning system (20), which are

$$ 
\begin{aligned}
\psi_{1,1} \sin(\theta \tau_1) + \psi_{1,2} \cos(\theta \tau_1) + r \tau_1 + s - 1 &= 0, \\
\psi_{2,1} \sin(\theta \tau_1) + \psi_{2,2} \cos(\theta \tau_1) + (kr + \sigma)\tau_1 + ks - \varepsilon r + \beta + z_0 - y_1 &= 0, \\
\sigma \tau_1 + z_0 - z_1 &= 0,
\end{aligned}
$$

$$ 
\begin{aligned}
(2 - \gamma + \sigma\varepsilon - z_1)(e^{-\tau_2/\varepsilon} - 1) + (y_1 - 1)/(1 - \varepsilon)(e^{-\tau_2/\varepsilon} - e^{-\tau_2}) + \sigma\tau_2 &= 0, \\
1 + (y_1 - 1)^{-\tau_2} - y_2 &= 0, \\
\sigma\tau_2 + z_1 - z_2 &= 0,
\end{aligned}
$$

$$ 
\begin{aligned}
\psi_{5,1} \sin(\theta \tau_3) + \psi_{5,2} \cos(\theta \tau_3) + r \tau_3 + s - 1 &= 0, \\
\psi_{6,1} \sin(\theta \tau_3) + \psi_{6,2} \cos(\theta \tau_3) + (kr + \sigma)\tau_3 + ks - \varepsilon r + \beta + z_0 - y_2 &= 0, \\
\sigma \tau_3 + z_0 - z_1 &= 0,
\end{aligned}
$$
where $\rho \pm \theta_i$ are the eigenvalues of the system in region $R_{32}$ (see Appendix A, expressions (A3) and (A4)), $r = \sigma/(l-k)$, $s = ((l-\varepsilon)\sigma + (l-k)(\beta + I - n))/(l-k)^2$ and

$$ \psi_{1,1} = \frac{e^{\rho\tau_1}\sqrt{1-a}}{\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}}((a-s)(k+\varepsilon) + 2(ks - \varepsilon r + \beta + z_0)), $$

$$ \psi_{1,2} = -e^{\rho\tau_1}, $$

$$ \psi_{2,1} = \frac{e^{\rho\tau_1}\sqrt{1-a}}{\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}}\left(\frac{2\varepsilon(a-s)}{1-a} + (ks - \varepsilon r + \beta + z_0)(k+\varepsilon)\right), $$

$$ \psi_{2,2} = e^{\rho\tau_1}(-ks + \varepsilon r - \beta - z_0), $$

$$ \psi_{5,1} = \frac{e^{\rho\tau_3}\sqrt{1-a}}{\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}}((a-s)(k+\varepsilon) + 2(-y(a,\varepsilon) + ks - \varepsilon r + \beta + z_0)), $$

$$ \psi_{5,2} = -e^{\rho\tau_3}, $$

$$ \psi_{6,1} = \frac{e^{\rho\tau_3}\sqrt{1-a}}{\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}}\left(\frac{2\varepsilon(a-s)}{1-a} + (-y(a,\varepsilon) + ks - \varepsilon r + \beta + z_0)(k+\varepsilon)\right), $$

$$ \psi_{6,2} = e^{\rho\tau_3}(y(a,\varepsilon) - ks + \varepsilon r - \beta - z_0). $$

We label the left term of each equation by $\tilde{\phi}_i$, being $i$ the $i$'th equation, for $i = 1, \dots, 9$. Then, applying the change of variables

$$ (\tau_1, \tau_3, y_1, y_2) = \left( \frac{\sqrt{1-a}}{\theta} \tilde{\tau}_1, \frac{\sqrt{1-a}}{\theta} \tilde{\tau}_3, (1-a)g_1, y(a, \varepsilon) + (1-a)g_2 \right), $$

and rewriting the nine equations $\tilde{\phi}_i = 0$ with $i = 1, 2, \dots, 9$ concerning system (20) into the following way $\tilde{\phi}_i = (1-a)\tilde{\phi}_i$ for $i = 1, 2, 7, 8$ and $\tilde{\phi}_i = \tilde{\phi}_i$ for $i = 3, 4, 5, 6, 9$. Therefore, the new equations write as

$$ \tilde{\phi}_1 = (k + \beta + z_0)\tilde{\tau}_1/\sqrt{\varepsilon} - 1 + O(1 - a), $$

$$ \tilde{\phi}_2 = (k + \beta + z_0)\tilde{\tau}_1^2/2 - g_1 + O(1 - a), $$

$$ \tilde{\phi}_7 = (k + \beta + z_2)\sqrt{\varepsilon}\tilde{\tau}_3 - 1 + O(1 - a), $$

$$ \tilde{\phi}_8 = -(k + \beta + z_2)\tilde{\tau}_3\sqrt{\varepsilon}(\sqrt{\varepsilon}\tilde{\tau}_3 + 2(k + \varepsilon))/(2k) - g_2 + O(1 - a). $$

The resulting system of equations has a Jacobian matrix, evaluated to the solution, omitting the columns of variables $a$ and $\sigma$, is given by

$$ \begin{pmatrix} (k + \beta + I)/\sqrt{\varepsilon} & 0 & 0 & 0 & 0 & z_{0,1} & 0 & 0 & 0 \\ \tilde{\tau}_{1,2} & 0 & 0 & -1 & 0 & z_{0,2} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & \tau_{2,3} & 0 & 0 & 0 & 0 & z_{1,4} & 0 & 0 \\ 0 & \tau_{2,4} & 0 & 0 & 0 & 0 & z_{1,5} & 0 & 0 \\ 0 & \sigma & 0 & 0 & 0 & 0 & 1 & -1 & 0 \\ 0 & 0 & (k + \beta + I)\sqrt{\varepsilon}/k & 0 & 0 & 0 & z_{2,7} & 0 & 0 \\ 0 & 0 & \tilde{\tau}_{3,6} & 0 & -1 & 0 & 0 & z_{2,8} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 \end{pmatrix}, $$

where $\tau_{2,3} = \frac{k+\beta+I-1+e^{-\tau_2}}{\varepsilon}$, $\tau_{2,4} = e^{-\tau_2}$, $z_{1,4} = 1 - e^{-\tau_2/\varepsilon}$ and $z_{1,5} = -\frac{k+\varepsilon}{k}$. Computing its determinant using the co-factor expansion, we obtain that this determinant is the same to those of matrix A in Appendix C, which is proved to be non-zero. Therefore, the Implicit Function Theorem can be applied to write all the parameters in terms of $a$ and $\sigma$. $\square$

<figure><img src="image_13.png" /></figure>
where $\rho \pm \theta i$ are the eigenvalues of the system in region $R_{32}$ (see Appendix A, expressions (A3) and (A4)), $r = \sigma/(l-k)$, $s = ((l-\varepsilon)\sigma + (l-k)(\beta + I - n))/(l-k)^2$ and

$$ \psi_{1,1} = \frac{e^{\rho\tau_1}\sqrt{1-a}}{\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}}((a-s)(k+\varepsilon) + 2(ks - \varepsilon r + \beta + z_0)), $$

$$ \psi_{1,2} = -e^{\rho\tau_1}, $$

$$ \psi_{2,1} = \frac{e^{\rho\tau_1}\sqrt{1-a}}{\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}}\left(\frac{2\varepsilon(a-s)}{1-a} + (ks - \varepsilon r + \beta + z_0)(k+\varepsilon)\right), $$

$$ \psi_{2,2} = e^{\rho\tau_1}(-ks + \varepsilon r - \beta - z_0), $$

$$ \psi_{5,1} = \frac{e^{\rho\tau_3}\sqrt{1-a}}{\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}}((a-s)(k+\varepsilon) + 2(-y(a,\varepsilon) + ks - \varepsilon r + \beta + z_0)), $$

$$ \psi_{5,2} = -e^{\rho\tau_3}, $$

$$ \psi_{6,1} = \frac{e^{\rho\tau_3}\sqrt{1-a}}{\sqrt{4\varepsilon - (1-a)(k+\varepsilon)^2}}\left(\frac{2\varepsilon(a-s)}{1-a} + (-y(a,\varepsilon) + ks - \varepsilon r + \beta + z_0)(k+\varepsilon)\right), $$

$$ \psi_{6,2} = e^{\rho\tau_3}(y(a,\varepsilon) - ks + \varepsilon r - \beta - z_0). $$

We label the left term of each equation by $\tilde{\phi}_i$, being $i$ the $i$'th equation, for $i = 1, \dots, 9$. Then, applying the change of variables

$$ (\tau_1, \tau_3, y_1, y_2) = \left( \frac{\sqrt{1-a}}{\theta} \tilde{\tau}_1, \frac{\sqrt{1-a}}{\theta} \tilde{\tau}_3, (1-a)g_1, y(a, \varepsilon) + (1-a)g_2 \right), $$

and rewriting the nine equations $\tilde{\phi}_i = 0$ with $i = 1, 2, \dots, 9$ concerning system (20) into the following way $\hat{\phi}_i = (1-a)\hat{\phi}_i$ for $i = 1, 2, 7, 8$ and $\tilde{\phi}_i = \hat{\phi}_i$ for $i = 3, 4, 5, 6, 9$. Therefore, the new equations write as

$$ \hat{\phi}_1 = (k + \beta + z_0)\tilde{\tau}_1/\sqrt{\varepsilon} - 1 + O(1-a), $$

$$ \hat{\phi}_2 = (k + \beta + z_0)\tilde{\tau}_1^2/2 - g_1 + O(1-a), $$

$$ \hat{\phi}_7 = (k + \beta + z_2)\sqrt{\varepsilon}\tilde{\tau}_3 - 1 + O(1-a), $$

$$ \hat{\phi}_8 = -(k + \beta + z_2)\tilde{\tau}_3\sqrt{\varepsilon}(\sqrt{\varepsilon}\tilde{\tau}_3 + 2(k + \varepsilon))/(2k) - g_2 + O(1-a). $$

The resulting system of equations has a Jacobian matrix, evaluated to the solution, omitting the columns of variables $a$ and $\sigma$, is given by

$$ \begin{pmatrix} (k + \beta + I)/\sqrt{\varepsilon} & 0 & 0 & 0 & 0 & z_{0,1} & 0 & 0 & 0 \\ \tilde{\tau}_{1,2} & 0 & 0 & -1 & 0 & z_{0,2} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & \tau_{2,3} & 0 & 0 & 0 & 0 & z_{1,4} & 0 & 0 \\ 0 & \tau_{2,4} & 0 & 0 & 0 & 0 & z_{1,5} & 0 & 0 \\ 0 & \sigma & 0 & 0 & 0 & 0 & 1 & -1 & 0 \\ 0 & 0 & (k + \beta + I)\sqrt{\varepsilon}/k & 0 & 0 & 0 & z_{2,7} & 0 & 0 \\ 0 & 0 & \tilde{\tau}_{3,6} & 0 & -1 & 0 & 0 & z_{2,8} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 \end{pmatrix}, $$

where $\tau_{2,3} = \frac{k+\beta+I-1+e^{-\tau_2}}{\varepsilon}$, $\tau_{2,4} = e^{-\tau_2}$, $z_{1,4} = 1 - e^{-\tau_2/\varepsilon}$ and $z_{1,5} = -\frac{k+\varepsilon}{k}$. Computing its determinant using the co-factor expansion, we obtain that this determinant is the same to those of matrix A in Appendix C, which is proved to be non-zero. Therefore, the Implicit Function Theorem can be applied to write all the parameters in terms of $a$ and $\sigma$.

□
It remains to prove that these expressions lie in the domain given by the restrictions (21). Similar arguments to those in Lemma 11 can be applied in order to assure the existence of such a maximal canard connection in a neighbour of the point.

Finally, for $0 < \varepsilon < k$, and for $0 < 1 - a$ and $0 < -\sigma$ sufficiently small, the differential system with parameter $I = M_1(a, \sigma)_4 = I_1(\sigma; a, \varepsilon)$, where $M_1(a, \sigma)_4$ stands for the fourth coordinate of function $M_1(a, \sigma)$, exhibits a 1-spike maximal canard connecting the manifolds $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$, which ends the proof of Theorem 4.

# 6 Conclusion

In this paper, we have introduced and analyzed a piecewise-linear version of the so-called Morris-Lecar model, which we have called PWL-ML model. In this study, we have focused on the existence of a saddle homoclinic connection and associated dynamical phenomena.

Since the PWL-ML model is strongly inspired by the nullcline configuration of the original ML model, one would expect that the behaviour of the flows of both models to be comparable. However, it is somewhat more surprising that the bifurcation structure of both models are also similar; compare for instance Figure 4 and Figure 9.28 in [13]. It remains to find the saddle-node on invariant circle bifurcation manifold, which we conjecture will also start at the saddle-node homoclinic orbit bifurcation curve.

There also exist other approaches to reproduce the Morris-Lecar model through a piecewise linearization, see for instance [48, 49]. Nevertheless, the nullclines configuration of these models are not similar to the one of the original ML model.

In addition, we have proved the existence of a saddle homoclinic orbit. The proof involves a new method to find dynamical objects that we called the contraction and expansion technique. There, we collapse one region of the PWL continuous system, which effectively turns the system into a discontinuous one, and we find the homoclinic connection in the discontinuous limit. Using the Implicit Function Theorem through the algebraic expressions of the connection, we are able to show the persistence of the homoclinic connection to the continuous case. This novel technique can be a useful tool in order to prove the existence of other local or global dynamical objects, such as periodic orbits or heteroclinic connections.

Furthermore, we have also extended results about the slow-passage phenomenon. In general, the slow-passage phenomenon is a common behavior in multiple-timescale dynamical systems, widely studied in the smooth context. However, there are not many slow-passage type results in the PWL context, in particular related to a homoclinic connection of the fast subsystem. We recently proved rigorous results related to the slow passage through a Hopf-like bifurcation, where we have successfully reproduced the phenomenon using a minimal model [47].
It remains to prove that these expressions lie in the domain given by the restrictions (21). Similar arguments to those in Lemma 11 can be applied in order to assure the existence of such a maximal canard connection in a neighbour of the point.

Finally, for $0 < \varepsilon < k$, and for $0 < 1 - a$ and $0 < -\sigma$ sufficiently small, the differential system with parameter $I = M_1(a, \sigma)_4 = I_1(\sigma; a, \varepsilon)$, where $M_1(a, \sigma)_4$ stands for the fourth coordinate of function $M_1(a, \sigma)$, exhibits a 1-spike maximal canard connecting the manifolds $W^u(S_\sigma^{sd})$ and $W^s(S_\sigma^{sd})$, which ends the proof of Theorem 4.

# 6 Conclusion

In this paper, we have introduced and analyzed a piecewise-linear version of the so-called Morris-Lecar model, which we have called PWL-ML model. In this study, we have focused on the existence of a saddle homoclinic connection and associated dynamical phenomena.

Since the PWL-ML model is strongly inspired by the nullcline configuration of the original ML model, one would expect that the behaviour of the flows of both models to be comparable. However, it is somewhat more surprising that the bifurcation structure of both models are also similar; compare for instance Figure 4 and Figure 9.28 in [13]. It remains to find the saddle-node on invariant circle bifurcation manifold, which we conjecture will also start at the saddle-node homoclinic orbit bifurcation curve.

There also exist other approaches to reproduce the Morris-Lecar model through a piecewise linearization, see for instance [48, 49]. Nevertheless, the nullclines configuration of these models are not similar to the one of the original ML model.

In addition, we have proved the existence of a saddle homoclinic orbit. The proof involves a new method to find dynamical objects that we called the contraction and expansion technique. There, we collapse one region of the PWL continuous system, which effectively turns the system into a discontinuous one, and we find the homoclinic connection in the discontinuous limit. Using the Implicit Function Theorem through the algebraic expressions of the connection, we are able to show the persistence of the homoclinic connection to the continuous case. This novel technique can be a useful tool in order to prove the existence of other local or global dynamical objects, such as periodic orbits or heteroclinic connections.

Furthermore, we have also extended results about the slow-passage phenomenon. In general, the slow-passage phenomenon is a common behavior in multiple-timescale dynamical systems, widely studied in the smooth context. However, there are not many slow-passage type results in the PWL context, in particular related to a homoclinic connection of the fast subsystem. We recently proved rigorous results related to the slow passage through a Hopf-like bifurcation, where we have successfully reproduced the phenomenon using a minimal model [47].
The slow-passage phenomenon through a homoclinic orbit is a key element in order to understand bursting phenomena [14, 31, 34]. In particular, it is strongly related to the spike-adding phenomenon for certain types of bursting orbits, for instance square-wave bursting, but also with other complex bursting behaviors [16]. In the present work, we perturb the homoclinic connection from a fast subsystem, giving rise to a finite number of new intersections between the stable and the unstable manifold of the slow manifold. In this context, we are able to prove the existence of maximal canard connections for 1 and 2 spikes in a straightforward manner. However, the existence of more than 2 spikes seems to require a more complex treatment, and we leave it for future work. Some numerical approaches show the existence of more maximal canards, thus the mechanism to generate spikes in the PWL context is likely to be the same as in the smooth context.

Overall, the present results contribute to a wider research program aiming to reproduce various types of neural activity with both qualitatively and quantitatively accurate PWL models. Here, we extend the scenarios that contain slow passages through homoclinic connections, and these are pertinent to a number of scenarios including square-wave bursting. Further work will be needed in order to fully understand the slow passage phenomenon in this context, which goes beyond the scope of the present work and will be an interesting avenue for a follow-up study.

**Acknowledgments.** Partial financial support was received from the Ministerio de Ciencia, Innovación y Universidades (MCIU) project PID2020-118726GB-I00. AET and CV are also partially supported by the Ministerio de Economía y Competitividad through the project MTM2017-83568-P (AEI/ERDF,EU). JPV is also partially supported by a grant Beques de Mobilitat EDUIB - Santander Universitats.

# Declarations

The authors have no conflicts of interest to declare that are relevant to the content of this article.
The slow-passage phenomenon through a homoclinic orbit is a key element in order to understand bursting phenomena [14, 31, 34]. In particular, it is strongly related to the spike-adding phenomenon for certain types of bursting orbits, for instance square-wave bursting, but also with other complex bursting behaviors [16]. In the present work, we perturb the homoclinic connection from a fast subsystem, giving rise to a finite number of new intersections between the stable and the unstable manifold of the slow manifold. In this context, we are able to prove the existence of maximal canard connections for 1 and 2 spikes in a straightforward manner. However, the existence of more than 2 spikes seems to require a more complex treatment, and we leave it for future work. Some numerical approaches show the existence of more maximal canards, thus the mechanism to generate spikes in the PWL context is likely to be the same as in the smooth context.

Overall, the present results contribute to a wider research program aiming to reproduce various types of neural activity with both qualitatively and quantitatively accurate PWL models. Here, we extend the scenarios that contain slow passages through homoclinic connections, and these are pertinent to a number of scenarios including square-wave bursting. Further work will be needed in order to fully understand the slow passage phenomenon in this context, which goes beyond the scope of the present work and will be an interesting avenue for a follow-up study.

**Acknowledgments.** Partial financial support was received from the Ministerio de Ciencia, Innovación y Universidades (MCIU) project PID2020-118726GB-I00. AET and CV are also partially supported by the Ministerio de Economía y Competitividad through the project MTM2017-83568-P (AEI/ERDF,EU). JPV is also partially supported by a grant Beques de Mobilitat EDUIB - Santander Universitats.

# Declarations

The authors have no conflicts of interest to declare that are relevant to the content of this article.
# Appendix A Local flow

In this appendix, we present some local properties of the flow of system (2) when the parameters are contained in the region $\mathcal{P}$, see (4).

System (2) can be written in terms of its associated linear system

$$ \dot{\mathbf{x}} = M_{i,j}\mathbf{x} + b_{i,j}, $$

where $M_{i,j}$ and $b_{i,j}$ are the matrix and the vector components associated to the linear region $R_{i,j}$. The local flow of (2) in the region $R_{i,j}$ is determined by the trace, the determinant and the discriminant of $M_{i,j}$: tr$_{i,j}$, det$_{i,j}$, and $\Delta_{i,j}$ respectively. Let us denote $\mathbf{p}_{i,j} = -M_{i,j}^{-1}b_{i,j}$. If $\mathbf{p}_{i,j} \in R_{i,j}$, this point is an equilibrium point of the system; otherwise, it is a virtual equilibrium [59].

Depending on the region, the matrices, traces, determinants and equilibrium points are as follows. Matrix

$$ M_{1,1} = \begin{pmatrix} -1/\varepsilon & -1/\varepsilon \\ 0 & -1 \end{pmatrix}, $$

with det$_{1,1} = 1/\varepsilon > 0$, tr$_{1,1} = -1/\varepsilon - 1 < 0$ and $\Delta_{1,1} = (1/\varepsilon - 1)^2 \ge 0$. When point $\mathbf{p}_{1,1} = (I, 0)$ is an equilibrium point it is an attractor node. Matrix

$$ M_{2,1} = \begin{pmatrix} \delta/\varepsilon & -1/\varepsilon \\ 0 & -1 \end{pmatrix}, $$

with det$_{2,1} = -\delta/\varepsilon > 0$, tr$_{2,1} = \delta/\varepsilon - 1 < 0$ and $\Delta_{2,1} = (\delta/\varepsilon + 1)^2 \ge 0$. Again, when point $\mathbf{p}_{2,1} = (\frac{\lambda+I}{-\delta}, 0)$ is an equilibrium point it is an attractor node. Matrix

$$ M_{3,1} = \begin{pmatrix} k/\varepsilon & -1/\varepsilon \\ 0 & -1 \end{pmatrix}, $$

with det$_{3,1} = -k/\varepsilon < 0$. Since the determinant is negative, when the point $\mathbf{p}_{3,1} = (-\frac{\beta+I}{k}, 0)$ is an equilibrium, it is a saddle. The eigenvectors associated to the eigenvalues $\lambda_1 = k/\varepsilon$, $\lambda_2 = -1$ are $v_1 = (1, 0)$, $v_2 = (k + \varepsilon, 1)$, respectively. These directions over the equilibrium point define the unstable manifold $W^u(\mathbf{p}_{3,1})$ and the stable manifold $W^s(\mathbf{p}_{3,1})$, respectively. Matrix

$$ M_{3,2} = \begin{pmatrix} k/\varepsilon & -1/\varepsilon \\ l & -1 \end{pmatrix}. $$

with det$_{3,2} = (l - k)/\varepsilon > 0$, tr$_{3,2} = k/\varepsilon - 1$ and $\Delta_{3,2} = \frac{1}{\varepsilon^2}((k + \varepsilon)^2 - 4l\varepsilon)$. When the point $\mathbf{p}_{3,2} = \frac{1}{k-l}(n - \beta - I, kn - Il - \beta l)$ is an equilibrium point it has either a focus or a node type dynamics. In particular, for sufficiently large $l$ the point $\mathbf{p}_{3,2}$ is a repelling focus with eigenvalues $\rho \pm \theta i$, where $\rho = (k - \varepsilon)/(2\varepsilon)$
# Appendix A Local flow

In this appendix, we present some local properties of the flow of system (2) when the parameters are contained in the region $\mathcal{P}$, see (4).

System (2) can be written in terms of its associated linear system

$$ \dot{\mathbf{x}} = M_{i,j}\mathbf{x} + b_{i,j}, $$

where $M_{i,j}$ and $b_{i,j}$ are the matrix and the vector components associated to the linear region $R_{i,j}$. The local flow of (2) in the region $R_{i,j}$ is determined by the trace, the determinant and the discriminant of $M_{i,j}$: $\text{tr}_{i,j}$, $\text{det}_{i,j}$, and $\Delta_{i,j}$ respectively. Let us denote $\mathbf{p}_{i,j} = -M_{i,j}^{-1}b_{i,j}$. If $\mathbf{p}_{i,j} \in R_{i,j}$, this point is an equilibrium point of the system; otherwise, it is a virtual equilibrium [59].

Depending on the region, the matrices, traces, determinants and equilibrium points are as follows. Matrix

$$ M_{1,1} = \begin{pmatrix} -1/\varepsilon & -1/\varepsilon \\ 0 & -1 \end{pmatrix}, $$

with $\text{det}_{1,1} = 1/\varepsilon > 0$, $\text{tr}_{1,1} = -1/\varepsilon - 1 < 0$ and $\Delta_{1,1} = (1/\varepsilon - 1)^2 \ge 0$. When point $\mathbf{p}_{1,1} = (I, 0)$ is an equilibrium point it is an attractor node. Matrix

$$ M_{2,1} = \begin{pmatrix} \delta/\varepsilon & -1/\varepsilon \\ 0 & -1 \end{pmatrix}, $$

with $\text{det}_{2,1} = -\delta/\varepsilon > 0$, $\text{tr}_{2,1} = \delta/\varepsilon - 1 < 0$ and $\Delta_{2,1} = (\delta/\varepsilon + 1)^2 \ge 0$. Again, when point $\mathbf{p}_{2,1} = (\frac{\lambda+I}{-\delta}, 0)$ is an equilibrium point it is an attractor node. Matrix

$$ M_{3,1} = \begin{pmatrix} k/\varepsilon & -1/\varepsilon \\ 0 & -1 \end{pmatrix}, $$

with $\text{det}_{3,1} = -k/\varepsilon < 0$. Since the determinant is negative, when the point $\mathbf{p}_{3,1} = (-\frac{\beta+I}{k}, 0)$ is an equilibrium, it is a saddle. The eigenvectors associated to the eigenvalues $\lambda_1 = k/\varepsilon$, $\lambda_2 = -1$ are $v_1 = (1, 0)$, $v_2 = (k + \varepsilon, 1)$, respectively. These directions over the equilibrium point define the unstable manifold $W^u(\mathbf{p}_{3,1})$ and the stable manifold $W^s(\mathbf{p}_{3,1})$, respectively. Matrix

$$ M_{3,2} = \begin{pmatrix} k/\varepsilon & -1/\varepsilon \\ l & -1 \end{pmatrix}. $$

with $\text{det}_{3,2} = (l - k)/\varepsilon > 0$, $\text{tr}_{3,2} = k/\varepsilon - 1$ and $\Delta_{3,2} = \frac{1}{\varepsilon^2}((k + \varepsilon)^2 - 4l\varepsilon)$. When the point $\mathbf{p}_{3,2} = \frac{1}{k-l}(n - \beta - I, kn - Il - \beta l)$ is an equilibrium point it has either a focus or a node type dynamics. In particular, for sufficiently large $l$ the point $\mathbf{p}_{3,2}$ is a repelling focus with eigenvalues $\rho \pm \theta i$, where $\rho = (k - \varepsilon)/(2\varepsilon)$