and $\theta = \sqrt{4l\epsilon - (k + \epsilon)^2/(2\epsilon)}$. In such a case, the local expression of the orbit of system (2) with initial condition $\mathbf{q}_0 = (x_0, y_0)^T \in \mathbb{R}_{3,2}$ is given by:

$$ \begin{pmatrix} x(t) \\ y(t) \end{pmatrix} = P^{-1}e^{Jt}P(\mathbf{q}_0 + A^{-1}\mathbf{a}_0) - A^{-1}\mathbf{a}_0, $$

where,

$$ P^{-1} = \begin{pmatrix} \frac{\rho+1}{l} & \frac{\theta}{l} \\ 1 & 0 \end{pmatrix}, \quad e^{Jt} = e^{\rho t} \begin{pmatrix} \cos(\theta t) & \sin(\theta t) \\ -\sin(\theta t) & \cos(\theta t) \end{pmatrix}, \\ A = \begin{pmatrix} k/\epsilon & -1/\epsilon \\ l & -1 \end{pmatrix}, \quad \mathbf{a}_0 = \begin{pmatrix} \frac{1}{\epsilon}(\beta + I), & n \end{pmatrix}^T. $$

We finally consider matrix

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

which satisfies that $M_{4,3} = M_{1,1}$. Hence, when $\mathbf{p}_{4,3} = (\gamma - 1 + I, 1)$ is an equilibrium point then it is an attractor node. The expressions of the orbits in the half-space $\{x > 1\}$ are given by

$$ \begin{aligned} x(t) &= C_1 e^{-t/\epsilon} - \frac{C_2}{1-\epsilon} e^{-t} - 1 + \gamma + I, \\ y(t) &= C_2 e^{-t} + 1, \end{aligned} \qquad (A1) $$

where $C_1, C_2 \in \mathbb{R}$.

Let us consider the point $\mathbf{p} = (1, -1+\gamma+I)^T$ on the separation line $\{x = 1\}$ of the region $\mathbb{R}_{4,3}$. Denoting by $\dot{\mathbf{p}}$ the vector field of system (2) evaluated at $\mathbf{p}$ it follows that $\dot{\mathbf{p}} = M_{4,3}\mathbf{p} + b_{4,3} = (0, 2-\gamma-I)^T$. Hence, provided $2-\gamma-I \neq 0$, the set $\{\mathbf{p}, \dot{\mathbf{p}}\}$ is a base of $\mathbb{R}^2$ called a Krylov base, and the orbits of system (2) can be expressed in terms of this base as $\Gamma(t) = u(t)\mathbf{p} + v(t)\dot{\mathbf{p}}$. From Theorem 5 of [60], as long as $\Gamma(t) \subset \mathbb{R}_{4,3}$, function

$$ H(u, v) = (u + \lambda_1 v)^{\lambda_2} (u + \lambda_2 v)^{-\lambda_1}, \qquad (A2) $$

is constant along $\Gamma(t)$, where $\lambda_1, \lambda_2$ are the eigenvalues of matrix $M_{4,3}$. Therefore, the point transformation defined by the flow in region $\mathbb{R}_{4,3}$ on the line $\{x = 1\}$ can be expressed as $\pi(z)$, where the point $\mathbf{p} + z\dot{\mathbf{p}}$ maps onto the point $\mathbf{p} + \pi(z)\dot{\mathbf{p}}$, through the implicit relation $H(1, z) = H(1, \pi(z))$.

For a sufficiently smooth function $h(z)$, and to simplify the notation, let us consider the following abuse of notation $H(z) = H(1, h(z))$. It follows that

$$ H''(z) = D(\lambda_1 - \lambda_2)(\lambda_1 h(z) + 1)^{\lambda_2 - 2}(\lambda_2 h(z) + 1)^{-\lambda_1 - 2} \\ \qquad (h'(z))^2(D(\lambda_1 - \lambda_2 + 1)h(z)^2 - 1) $$
and $\theta = \sqrt{4l\epsilon - (k + \epsilon)^2/(2\epsilon)}$. In such a case, the local expression of the orbit of system (2) with initial condition $\mathbf{q}_0 = (x_0, y_0)^T \in \mathbb{R}_{3,2}$ is given by:

$$ \begin{pmatrix} x(t) \\ y(t) \end{pmatrix} = P^{-1}e^{Jt}P(\mathbf{q}_0 + A^{-1}\mathbf{a}_0) - A^{-1}\mathbf{a}_0, $$

where,

$$ P^{-1} = \begin{pmatrix} \frac{\rho+1}{l} & \frac{\theta}{l} \\ 1 & 0 \end{pmatrix}, \quad e^{Jt} = e^{\rho t} \begin{pmatrix} \cos(\theta t) & \sin(\theta t) \\ -\sin(\theta t) & \cos(\theta t) \end{pmatrix}, \\ A = \begin{pmatrix} k/\epsilon & -1/\epsilon \\ l & -1 \end{pmatrix}, \quad \mathbf{a}_0 = \begin{pmatrix} \frac{1}{\epsilon}(\beta + I), & n \end{pmatrix}^T. $$

We finally consider matrix

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

which satisfies that $M_{4,3} = M_{1,1}$. Hence, when $\mathbf{p}_{4,3} = (\gamma - 1 + I, 1)$ is an equilibrium point then it is an attractor node. The expressions of the orbits in the half-space $\{x > 1\}$ are given by

$$ \begin{aligned} x(t) &= C_1 e^{-t/\epsilon} - \frac{C_2}{1-\epsilon} e^{-t} - 1 + \gamma + I, \\ y(t) &= C_2 e^{-t} + 1, \end{aligned} \tag{A1} $$

where $C_1, C_2 \in \mathbb{R}$.

Let us consider the point $\mathbf{p} = (1, -1+\gamma+I)^T$ on the separation line $\{x = 1\}$ of the region $\mathbb{R}_{4,3}$. Denoting by $\dot{\mathbf{p}}$ the vector field of system (2) evaluated at $\mathbf{p}$ it follows that $\dot{\mathbf{p}} = M_{4,3}\mathbf{p} + b_{4,3} = (0, 2-\gamma-I)^T$. Hence, provided $2-\gamma-I \neq 0$, the set $\{\mathbf{p}, \dot{\mathbf{p}}\}$ is a base of $\mathbb{R}^2$ called a Krylov base, and the orbits of system (2) can be expressed in terms of this base as $\Gamma(t) = u(t)\mathbf{p} + v(t)\dot{\mathbf{p}}$. From Theorem 5 of [60], as long as $\Gamma(t) \subset \mathbb{R}_{4,3}$, function

$$ H(u, v) = (u + \lambda_1 v)^{\lambda_2} (u + \lambda_2 v)^{-\lambda_1}, \tag{A2} $$

is constant along $\Gamma(t)$, where $\lambda_1, \lambda_2$ are the eigenvalues of matrix $M_{4,3}$. Therefore, the point transformation defined by the flow in region $\mathbb{R}_{4,3}$ on the line $\{x = 1\}$ can be expressed as $\pi(z)$, where the point $\mathbf{p} + z\dot{\mathbf{p}}$ maps onto the point $\mathbf{p} + \pi(z)\dot{\mathbf{p}}$, through the implicit relation $H(1, z) = H(1, \pi(z))$.

For a sufficiently smooth function $h(z)$, and to simplify the notation, let us consider the following abuse of notation $H(z) = H(1, h(z))$. It follows that

$$ H''(z) = D(\lambda_1 - \lambda_2)(\lambda_1 h(z) + 1)^{\lambda_2 - 2}(\lambda_2 h(z) + 1)^{-\lambda_1 - 2} \\ \qquad (h'(z))^2(D(\lambda_1 - \lambda_2 + 1)h(z)^2 - 1) $$
$$ -h(z)(\lambda_1 h(z) + 1)(\lambda_2 h(z) + 1)h''(z) $$

$$ H'''(z) = D(\lambda_1 h(z) + 1)^{\lambda_2-3} (\lambda_2 h(z) + 1)^{-\lambda_1-3} \\ \left( (\lambda_2 - \lambda_1)h(z)h^{(3)}(z)(\lambda_1 h(z) + 1)^2 (\lambda_2 h(z) + 1)^2 - (\lambda_1 - \lambda_2)h'(z)^3 \right. \\ \left. + 3(\lambda_1 - \lambda_2)(\lambda_1 h(z) + 1)(\lambda_2 h(z) + 1)h'(z)h''(z) \right. \\ \left. (D(\lambda_1 - \lambda_2 + 1)h(z)^2 - 3) - 2(T) \right) \\ (D(\lambda_1 - \lambda_2 + 1)h(z)^2 - 1) $$

where $T = \lambda_1 + \lambda_2$ and $D = \lambda_1\lambda_2$. Hence, by equalling $H''(0) = H''(\pi(0))$ and $H'''(0) = H'''(\pi(0))$, we get the first terms of the time series expansion of $\pi(z)$, where

$$ \pi(z) = -z + \frac{2T}{3}z^2 + O(z^3), $$

which describes the flow in $R_{4,3}$ in a neighbourhood of $z = 0$, corresponding with the coordinate in the Krylov base of the contact point $(1, k + \beta + I)$.

Finally, by direct integration, the solution of system (7) in the half-space $\{x > 1\}$ is given by

$$ \begin{aligned} x_{4,3}(t) &= C_1 e^{-t/\epsilon} - \frac{C_2}{1-\epsilon} e^{-t} + \sigma t - 1 + \gamma - \sigma\epsilon + I_0 \\ y_{4,3}(t) &= C_2 e^{-t} + 1 \\ I(t) &= \sigma t + I_0 \end{aligned} \qquad (A3) $$

We notice that $I(t) = \sigma t + I_0$ is the same expression for any region. On the other hand, orbits in region $R_{3,2}$ for the $\mathbf{x}(t) = (x_{3,2}(t), y_{3,2})^t$ variables are given by

$$ \mathbf{x}(t) = P^{-1}e^{Jt}P\mathbf{x}_0 + \mathbf{x}_p, \qquad (A4) $$

where $P, e^{Jt}$ are the ones given previously, $\mathbf{x}_0 \in \mathbb{R}^2$ and $\mathbf{x}_p$ is a particular solution; for instance, in our case we have taken the solution $\mathbf{x}_p = (a_p t + b_p, (ka_p + \sigma)t + kb_p - \epsilon a_p + \beta + I_0)$, with $a_p = \sigma/(l-k)$ and $b_p = \frac{(k-\epsilon)\sigma/(l-k)+\sigma+\beta-n+I_0}{l-k}$.

## Appendix B 2-spikes maximal canards

In this appendix, we provide all the partial derivatives of the equations in system (34) with respect to the variables $\{\sigma, y_1, y_2, y_3, z_1, z_2, z_3, \tau_1, \tau_2, \tau_3\}$ and evaluated at the point $(0, \hat{y}_1(\epsilon), 0, \hat{y}_1(\epsilon), I(\epsilon), I(\epsilon), I(\epsilon), \tau(\epsilon), 0, \tau(\epsilon))$.

We first compute the partial derivative with respect to $y_1$. Since this variable only appears in the 2nd, 4th and 5th equation, it follows that $\frac{d\phi_1}{dy_1} = 0$ for $i \neq 2, 4, 5$. Moreover

$$ \frac{d\phi_2}{dy_1} = -1, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_4}{dy_1} = \frac{-1}{k+\epsilon}, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_5}{dy_1} = 0. $$
$$ -h(z)(\lambda_1 h(z) + 1)(\lambda_2 h(z) + 1)h''(z) ), $$

$$ H'''(z) = D(\lambda_1 h(z) + 1)^{\lambda_2 - 3} (\lambda_2 h(z) + 1)^{-\lambda_1 - 3} \\ \qquad \left( (\lambda_2 - \lambda_1)h(z)h^{(3)}(z)(\lambda_1 h(z) + 1)^2 (\lambda_2 h(z) + 1)^2 - (\lambda_1 - \lambda_2)h'(z)^3 \right. \\ \qquad \left. + 3(\lambda_1 - \lambda_2)(\lambda_1 h(z) + 1)(\lambda_2 h(z) + 1)h'(z)h''(z) \right. \\ \qquad \left. (D(\lambda_1 - \lambda_2 + 1)h(z)^2 - 3) - 2(T) \right) \\ \qquad + 3(\lambda_1 - \lambda_2)(\lambda_1 h(z) + 1)(\lambda_2 h(z) + 1)h'(z)h''(z) \\ \qquad (D(\lambda_1 - \lambda_2 + 1)h(z)^2 - 1) ), $$

where $T = \lambda_1 + \lambda_2$ and $D = \lambda_1\lambda_2$. Hence, by equalling $H''(0) = H''(\pi(0))$ and $H'''(0) = H'''(\pi(0))$, we get the first terms of the time series expansion of $\pi(z)$, where

$$ \pi(z) = -z + \frac{2T}{3}z^2 + O(z^3), $$

which describes the flow in $R_{4,3}$ in a neighbourhood of $z = 0$, corresponding with the coordinate in the Krylov base of the contact point $(1, k + \beta + I)$.

Finally, by direct integration, the solution of system (7) in the half-space $\{x > 1\}$ is given by

$$ \begin{aligned} x_{4,3}(t) &= C_1 e^{-t/\epsilon} - \frac{C_2}{1-\epsilon} e^{-t} + \sigma t - 1 + \gamma - \sigma\epsilon + I_0 \\ y_{4,3}(t) &= C_2 e^{-t} + 1 \\ I(t) &= \sigma t + I_0 \end{aligned} \qquad (A3) $$

We notice that $I(t) = \sigma t + I_0$ is the same expression for any region. On the other hand, orbits in region $R_{3,2}$ for the $\mathbf{x}(t) = (x_{3,2}(t), y_{3,2})^t$ variables are given by

$$ \mathbf{x}(t) = P^{-1}e^{Jt}P\mathbf{x}_0 + \mathbf{x}_p, \qquad (A4) $$

where $P, e^{Jt}$ are the ones given previously, $\mathbf{x}_0 \in \mathbb{R}^2$ and $\mathbf{x}_p$ is a particular solution; for instance, in our case we have taken the solution $\mathbf{x}_p = (a_p t + b_p, (ka_p + \sigma)t + kb_p - \epsilon a_p + \beta + I_0)$, with $a_p = \sigma/(l-k)$ and $b_p = \frac{(k-\epsilon)\sigma/(l-k)+\sigma+\beta-n+I_0}{l-k}$.

## Appendix B 2-spikes maximal canards

In this appendix, we provide all the partial derivatives of the equations in system (34) with respect to the variables $\{\sigma, y_1, y_2, y_3, z_1, z_2, z_3, \tau_1, \tilde{\tau}_2, \tau_3\}$ and evaluated at the point $(0, \hat{y}_1(\epsilon), 0, \hat{y}_1(\epsilon), I(\epsilon), I(\epsilon), I(\epsilon), \tau(\epsilon), 0, \tau(\epsilon))$.

We first compute the partial derivative with respect to $y_1$. Since this variable only appears in the 2nd, 4th and 5th equation, it follows that $\frac{d\phi_1}{dy_1} = 0$ for $i \neq 2, 4, 5$. Moreover

$$ \frac{d\phi_2}{dy_1} = -1, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_4}{dy_1} = \frac{-1}{k + \epsilon}, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_5}{dy_1} = 0. $$
Variable $y_2$ only appears in the 5th, 7th and 8th equations, so

$$ \frac{d\phi_5}{dy_2} = -1, \quad \frac{d\phi_7}{dy_2} = \frac{e^{\frac{\tau(\varepsilon)}{\varepsilon}} - e^{-\tau(\varepsilon)}}{1 - \varepsilon}, \quad \frac{d\phi_8}{dy_2} = e^{-\tau(\varepsilon)}. $$

Similarly, variable $y_3$ only appears in the 8th and 9th equations, so

$$ \frac{d\phi_8}{dy_3} = -1, \quad \frac{d\phi_9}{dy_3} = \frac{k}{k + \varepsilon}. $$

Variable $z_1$ only appears in 1st and 3rd equations, thus

$$ \frac{d\phi_1}{dz_1} = 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}}, \quad \frac{d\phi_3}{dz_1} = 1. $$

Variable $z_2$ appears in the 3rd, 4th and 5th equations, thus

$$ \frac{d\phi_3}{dz_2} = -1, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_4}{dz_2} = \frac{1}{k}, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_5}{dz_2} = 0. $$

Variable $z_3$ appears in the 6th, 7th and 9th equations, and so

$$ \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_6}{dz_3} = 0, \quad \frac{d\phi_7}{dz_3} = 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}}, \quad \frac{d\phi_9}{dz_3} = 1. $$

Variable $\sigma$ appears in the 1st, 3rd, 4th, 6th, 7th and 9th equations, and its derivatives are

$$ \frac{d\phi_1}{d\sigma} = \frac{d\phi_7}{d\sigma} = \varepsilon(e^{-\frac{\tau(\varepsilon)}{\varepsilon}} - 1) + \tau(\varepsilon), \quad \frac{d\phi_3}{d\sigma} = \frac{d\phi_9}{d\sigma} = \tau(\varepsilon), $$

$$ \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_4}{d\sigma} = \frac{\varepsilon}{k^2}, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_6}{d\sigma} = 1. $$

Variable $\tau_1$ appears in the 1st, 2nd, 3rd equations, and so

$$ \frac{d\phi_1}{d\tau_1} = \frac{k + \beta + I(\varepsilon)}{\varepsilon}, \quad \frac{d\phi_2}{d\tau_1} = e^{-\tau(\varepsilon)}, \quad \frac{d\phi_3}{d\tau_1} = 0. $$

Variable $\tilde{\tau}_2$ appears in the 4th, 5th and 6th equations, and its derivatives are

$$ \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_4}{d\tilde{\tau}_2} = 0, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_5}{d\tilde{\tau}_2} = 0, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_6}{d\tilde{\tau}_2} = z_2 - z_3. $$

Variable $\tau_3$ only appears in the 7th, 8th and 9th equations, so

$$ \frac{d\phi_7}{d\tau_3} = \frac{k + \beta + I(\varepsilon) - \hat{y}_1(\varepsilon)}{\varepsilon}, \quad \frac{d\phi_8}{d\tau_3} = (1 - \hat{y}_1(\varepsilon))e^{-\tau(\varepsilon)}, \quad \frac{d\phi_9}{d\tau_3} = 0. $$
Variable $y_2$ only appears in the 5th, 7th and 8th equations, so

$$ \frac{d\phi_5}{dy_2} = -1, \quad \frac{d\phi_7}{dy_2} = \frac{e^{\frac{\tau(\varepsilon)}{\varepsilon}} - e^{-\tau(\varepsilon)}}{1 - \varepsilon}, \quad \frac{d\phi_8}{dy_2} = e^{-\tau(\varepsilon)}. $$

Similarly, variable $y_3$ only appears in the 8th and 9th equations, so

$$ \frac{d\phi_8}{dy_3} = -1, \quad \frac{d\phi_9}{dy_3} = \frac{k}{k + \varepsilon}. $$

Variable $z_1$ only appears in 1st and 3rd equations, thus

$$ \frac{d\phi_1}{dz_1} = 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}}, \quad \frac{d\phi_3}{dz_1} = 1. $$

Variable $z_2$ appears in the 3rd, 4th and 5th equations, thus

$$ \frac{d\phi_3}{dz_2} = -1, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_4}{dz_2} = \frac{1}{k}, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_5}{dz_2} = 0. $$

Variable $z_3$ appears in the 6th, 7th and 9th equations, and so

$$ \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_6}{dz_3} = 0, \quad \frac{d\phi_7}{dz_3} = 1 - e^{-\frac{\tau(\varepsilon)}{\varepsilon}}, \quad \frac{d\phi_9}{dz_3} = 1. $$

Variable $\sigma$ appears in the 1st, 3rd, 4th, 6th, 7th and 9th equations, and its derivatives are

$$ \frac{d\phi_1}{d\sigma} = \frac{d\phi_7}{d\sigma} = \varepsilon\left(e^{-\frac{\tau(\varepsilon)}{\varepsilon}} - 1\right) + \tau(\varepsilon), \quad \frac{d\phi_3}{d\sigma} = \frac{d\phi_9}{d\sigma} = \tau(\varepsilon), $$

$$ \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_4}{d\sigma} = \frac{\varepsilon}{k^2}, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_6}{d\sigma} = 1. $$

Variable $\tau_1$ appears in the 1st, 2nd, 3rd equations, and so

$$ \frac{d\phi_1}{d\tau_1} = \frac{k + \beta + I(\varepsilon)}{\varepsilon}, \quad \frac{d\phi_2}{d\tau_1} = e^{-\tau(\varepsilon)}, \quad \frac{d\phi_3}{d\tau_1} = 0. $$

Variable $\tilde{\tau}_2$ appears in the 4th, 5th and 6th equations, and its derivatives are

$$ \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_4}{d\tilde{\tau}_2} = 0, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_5}{d\tilde{\tau}_2} = 0, \quad \lim_{\tilde{\tau}_2 \to 0} \frac{d\phi_6}{d\tilde{\tau}_2} = z_2 - z_3. $$

Variable $\tau_3$ only appears in the 7th, 8th and 9th equations, so

$$ \frac{d\phi_7}{d\tau_3} = \frac{k + \beta + I(\varepsilon) - \hat{y}_1(\varepsilon)}{\varepsilon}, \quad \frac{d\phi_8}{d\tau_3} = (1 - \hat{y}_1(\varepsilon))e^{-\tau(\varepsilon)}, \quad \frac{d\phi_9}{d\tau_3} = 0. $$
# Appendix C Homoclinic connection

In this appendix we obtain system (32) from system (29), and we show that both systems are equivalent for $a < 1$. Additionally, we compute the Jacobian matrix of system (32) at the solution $v_0 = (0, \tau(\varepsilon_0), 0, 0, \hat{y}_1(\varepsilon), I(\varepsilon_0), 1, \varepsilon_0)$. Finally, we show that this Jacobian has a minor of maximal range.

In order to do that, let us consider the change of variables (31). From this, we obtain the following expression for the first equation in system (29):

$$ \tilde{\phi}_1 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\frac{\beta}{\theta}\sqrt{1-a}\tilde{\tau}_1} \left( \xi_{1,1} \frac{\sin(\sqrt{1-a}\tilde{\tau}_1)}{\sqrt{1-a}} - \cos(\sqrt{1-a}\tilde{\tau}_1) \right) - 1 + (k + \beta + I) \right]. $$

We notice that functions

$$ \frac{\sin(\sqrt{1-a}\tilde{\tau}_1)}{\sqrt{1-a}}, \quad \cos(\sqrt{1-a}\tilde{\tau}_1), \quad e^{\frac{\beta}{\theta}\sqrt{1-a}\tilde{\tau}_1}, $$

are well-defined analytical functions at $a = 1$, and

$$ \begin{aligned} \frac{\sin(\sqrt{1-a}\tilde{\tau}_1)}{\sqrt{1-a}} &= \tilde{\tau}_1 + \frac{\tilde{\tau}_1^3 (a-1)}{6} + O((1-a)^2), \\ \cos(\sqrt{1-a}\tilde{\tau}_1) &= 1 + \tilde{\tau}_1^2 (a-1)/2 + O((1-a)^2), \\ e^{\frac{\beta}{\theta}\sqrt{1-a}\tilde{\tau}_1} &= 1 + \frac{(\varepsilon-k)}{2\sqrt{\varepsilon}}\tilde{\tau}_1 (a-1) + O((1-a)^2). \end{aligned} $$

Hence, we may consider $\tilde{\phi}_1 = (1-a)\hat{\phi}_1$, where $\hat{\phi}_1 = \tilde{\tau}_1(k+\beta+I)/\sqrt{\varepsilon}-1+O(1-a)$. This gives us the theoretical value for the variable $\tilde{\tau}_1 = \sqrt{\varepsilon}/(k+\beta+I)$.

The 2nd equation in system (29) is given by:

$$ \tilde{\phi}_2 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I) \left( \xi_{2,1} \frac{\sin(\sqrt{1-a}\tilde{\tau}_1)}{\sqrt{1-a}} e^{\frac{\beta}{\theta}\sqrt{1-a}\tilde{\tau}_1} + (1-e^{\frac{\beta}{\theta}\sqrt{1-a}\tilde{\tau}_1}) \frac{\cos(\sqrt{1-a}\tilde{\tau}_1)}{1-a} \right) - (1-(1-a)k)g_1 \right] $$

Again, since the expression $(1 - e^{\frac{\beta}{\theta}\sqrt{1-a}\tilde{\tau}_1}) \frac{\cos(\sqrt{1-a}\tilde{\tau}_1)}{1-a}$ is well-defined and analytical at $a = 1$, we can consider $\tilde{\phi}_2 = (1-a)\hat{\phi}_2$ where $\hat{\phi}_2 = \frac{\tilde{\tau}_1^2(k+\beta+I)}{2\varepsilon} - g_1 + O(1-a)$. This gives us the theoretical value for variable $g_1 = \frac{\varepsilon}{2(k+\beta+I)}$.

The 3rd and 4th equations in system (29) become

$$ \begin{aligned} \hat{\phi}_3 &= (2 - \gamma - I)(e^{-\frac{\tau_2}{\varepsilon}} - 1) + \frac{1 - (1 - a)g_1}{1 - \varepsilon}(e^{-\tau_2} - e^{-\frac{\tau_2}{\varepsilon}}) = 0 \\ \hat{\phi}_4 &= -(1 - (1 - a)g_1)e^{-\tau_2} + 1 - \hat{y}_1(\varepsilon) - (1 - a)g_2 = 0. \end{aligned} $$
# Appendix C Homoclinic connection

In this appendix we obtain system (32) from system (29), and we show that both systems are equivalent for $a < 1$. Additionally, we compute the Jacobian matrix of system (32) at the solution $v_0 = (0, \tau(\varepsilon_0), 0, 0, \hat{y}_1(\varepsilon), I(\varepsilon_0), 1, \varepsilon_0)$. Finally, we show that this Jacobian has a minor of maximal range.

In order to do that, let us consider the change of variables (31). From this, we obtain the following expression for the first equation in system (29):

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

We notice that functions

$$ \frac{\sin(\sqrt{1-a}\tilde{\tau}_1)}{\sqrt{1-a}}, \quad \cos(\sqrt{1-a}\tilde{\tau}_1), \quad e^{\frac{\rho}{\theta}\sqrt{1-a}\tilde{\tau}_1}, $$

are well-defined analytical functions at $a = 1$, and

$$ \begin{aligned} \frac{\sin(\sqrt{1-a}\tilde{\tau}_1)}{\sqrt{1-a}} &= \tilde{\tau}_1 + \frac{\tilde{\tau}_1^3 (a-1)}{6} + O((1-a)^2), \\ \cos(\sqrt{1-a}\tilde{\tau}_1) &= 1 + \tilde{\tau}_1^2 (a-1)/2 + O((1-a)^2), \\ e^{\frac{\rho}{\theta}\sqrt{1-a}\tilde{\tau}_1} &= 1 + \frac{(\varepsilon-k)}{2\sqrt{\varepsilon}}\tilde{\tau}_1 (a-1) + O((1-a)^2). \end{aligned} $$

Hence, we may consider $\tilde{\phi}_1 = (1-a)\hat{\phi}_1$, where $\hat{\phi}_1 = \tilde{\tau}_1(k+\beta+I)/\sqrt{\varepsilon}-1+O(1-a)$. This gives us the theoretical value for the variable $\tilde{\tau}_1 = \sqrt{\varepsilon}/(k+\beta+I)$. The 2nd equation in system (29) is given by:

$$ \tilde{\phi}_2 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I) \left( \xi_{2,1} \frac{\sin(\sqrt{1-a}\tilde{\tau}_1)}{\sqrt{1-a}} e^{\frac{\rho}{\theta}\sqrt{1-a}\tilde{\tau}_1} + (1-e^{\frac{\rho}{\theta}\sqrt{1-a}\tilde{\tau}_1}) \frac{\cos(\sqrt{1-a}\tilde{\tau}_1)}{1-a} \right) - (1-(1-a)k)g_1 \right] $$

Again, since the expression $(1 - e^{\frac{\rho}{\theta}\sqrt{1-a}\tilde{\tau}_1}) \frac{\cos(\sqrt{1-a}\tilde{\tau}_1)}{1-a}$ is well-defined and analytical at $a = 1$, we can consider $\tilde{\phi}_2 = (1-a)\hat{\phi}_2$ where $\hat{\phi}_2 = \frac{\tilde{\tau}_1^2(k+\beta+I)}{2\varepsilon} - g_1 + O(1-a)$. This gives us the theoretical value for variable $g_1 = \frac{\varepsilon}{2(k+\beta+I)}$.

The 3rd and 4th equations in system (29) become

$$ \begin{aligned} \hat{\phi}_3 &= (2 - \gamma - I)(e^{-\frac{\tau_2}{\varepsilon}} - 1) + \frac{1 - (1 - a)g_1}{1 - \varepsilon}(e^{-\tau_2} - e^{-\frac{\tau_2}{\varepsilon}}) = 0 \\ \hat{\phi}_4 &= -(1 - (1 - a)g_1)e^{-\tau_2} + 1 - \hat{y}_1(\varepsilon) - (1 - a)g_2 = 0. \end{aligned} $$
The 5th equation is then given by:

$$ \tilde{\phi}_5 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\frac{a}{b}\sqrt{1-a}\tilde{\tau}_3} \left( \xi_{5,1} \frac{\sin(\sqrt{1-a}\tilde{\tau}_3)}{\sqrt{1-a}} - \cos(\sqrt{1-a}\tilde{\tau}_3) \right) - 1 + (k + \beta + I) \right], $$

which can be simplified again as $\tilde{\phi}_5 = (1-a)\hat{\phi}_5$ where $\hat{\phi}_5 = -(k + \beta + I)\sqrt{\varepsilon}\tilde{\tau}_3/k - 1 + O(1-a)$. This gives us the value $\tilde{\tau}_3 = -k/(\sqrt{\varepsilon}(k + \beta + I))$. Finally, $\tilde{\phi}_6$ is given by

$$ \tilde{\phi}_6 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\frac{a}{b}\sqrt{1-a}\tilde{\tau}_3} \left( \xi_{6,1} \frac{\sin(\sqrt{1-a}\tilde{\tau}_3)}{\sqrt{1-a}} + \xi_{6,2} \cos(\sqrt{1-a}\tilde{\tau}_3) \right) + \frac{1}{1-a}((ak + \beta + I) - (1-(1-a)k)(\hat{y}_1(\varepsilon) + (1-a)g_2)) \right]. $$

This expression can be written as $\tilde{\phi}_6 = (1-a)\hat{\phi}_6$ where $\hat{\phi}_6 = -(k + \beta + I)\tilde{\tau}_3\sqrt{\varepsilon}(\sqrt{\varepsilon}\tilde{\tau}_3 + 2(k + \varepsilon))/(2k) - g_2 + O(1-a)$. This provides the value $g_2 = (2(k + \varepsilon) - k/(k + \beta + I))/2$.

We conclude that system (32) is equivalent to system (29) when $a < 1$. Moreover, system (32) is defined at $a = 1$ and

$$ \tilde{v}_0 = \left( \frac{\sqrt{\varepsilon}}{k + \beta + I}, \tau, -\frac{k/\sqrt{\varepsilon}}{k + \beta + I}, \frac{\varepsilon/2}{k + \beta + I}, k + \varepsilon - \frac{k/2}{k + \beta + I}, I, 1, \varepsilon \right), $$

is a solution, where $\tau = \tau(\varepsilon)$ and $I = I(\varepsilon)$ are given in Lemma 3.

Let us now compute the Jacobian matrix of system (32) at $\tilde{v}_0$ and show that there exists a minor with maximal range. In particular, the partial derivatives of equation $\hat{\phi}_i$ with $i = 1, 2, \dots, 6$ with respect to the variables $(\tilde{\tau}_1, \tau_2, \tilde{\tau}_3, g_1, g_2, I)$ produce the minor

$$ A = \begin{pmatrix} \frac{k+\beta+I}{\sqrt{\varepsilon}} & 0 & 0 & 0 & 0 & I_1 \\ \tilde{\tau}_{1,2} & 0 & 0 & -1 & 0 & I_2 \\ 0 & \frac{k+\beta+I-1+e^{-\tau}}{\varepsilon} & 0 & 0 & 0 & 1-e^{-\frac{\tau}{\varepsilon}} \\ 0 & e^{-\tau} & 0 & 0 & 0 & -\frac{k+\varepsilon}{k} \\ 0 & 0 & \frac{(k+\beta+I)\sqrt{\varepsilon}}{k} & 0 & 0 & I_5 \\ 0 & 0 & \tilde{\tau}_{3,6} & 0 & -1 & I_6 \end{pmatrix} $$
The 5th equation is then given by:

$$ \tilde{\phi}_5 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\frac{a}{b}\sqrt{1-a}\tilde{\tau}_3} \left( \xi_{5,1} \frac{\sin(\sqrt{1-a}\tilde{\tau}_3)}{\sqrt{1-a}} - \cos(\sqrt{1-a}\tilde{\tau}_3) \right) - 1 + (k + \beta + I) \right], $$

which can be simplified again as $\tilde{\phi}_5 = (1-a)\hat{\phi}_5$ where $\hat{\phi}_5 = -(k + \beta + I)\sqrt{\varepsilon}\tilde{\tau}_3/k - 1 + O(1-a)$. This gives us the value $\tilde{\tau}_3 = -k/(\sqrt{\varepsilon}(k + \beta + I))$. Finally, $\tilde{\phi}_6$ is given by

$$ \tilde{\phi}_6 = \frac{1-a}{1-(1-a)k} \left[ (ak + \beta + I)e^{\frac{a}{b}\sqrt{1-a}\tilde{\tau}_3} \left( \xi_{6,1} \frac{\sin(\sqrt{1-a}\tilde{\tau}_3)}{\sqrt{1-a}} + \xi_{6,2} \cos(\sqrt{1-a}\tilde{\tau}_3) \right) + \frac{1}{1-a}((ak + \beta + I) - (1-(1-a)k)(\hat{y}_1(\varepsilon) + (1-a)g_2)) \right]. $$

This expression can be written as $\tilde{\phi}_6 = (1-a)\hat{\phi}_6$ where $\hat{\phi}_6 = -(k + \beta + I)\tilde{\tau}_3\sqrt{\varepsilon}(\sqrt{\varepsilon}\tilde{\tau}_3 + 2(k+\varepsilon))/(2k) - g_2 + O(1-a)$. This provides the value $g_2 = (2(k+\varepsilon) - k/(k+\beta+I))/2$.

We conclude that system (32) is equivalent to system (29) when $a < 1$. Moreover, system (32) is defined at $a = 1$ and

$$ \tilde{v}_0 = \left( \frac{\sqrt{\varepsilon}}{k + \beta + I}, \tau, -\frac{k/\sqrt{\varepsilon}}{k + \beta + I}, \frac{\varepsilon/2}{k + \beta + I}, k + \varepsilon - \frac{k/2}{k + \beta + I}, I, 1, \varepsilon \right), $$

is a solution, where $\tau = \tau(\varepsilon)$ and $I = I(\varepsilon)$ are given in Lemma 3.

Let us now compute the Jacobian matrix of system (32) at $\tilde{v}_0$ and show that there exists a minor with maximal range. In particular, the partial derivatives of equation $\hat{\phi}_i$ with $i = 1, 2, \dots, 6$ with respect to the variables $(\tilde{\tau}_1, \tilde{\tau}_2, \tilde{\tau}_3, g_1, g_2, I)$ produce the minor

$$ A = \begin{pmatrix} \frac{k+\beta+I}{\sqrt{\varepsilon}} & 0 & 0 & 0 & 0 & I_1 \\ \tilde{\tau}_{1,2} & 0 & 0 & -1 & 0 & I_2 \\ 0 & \frac{k+\beta+I-1+e^{-\tau}}{\varepsilon} & 0 & 0 & 0 & 1-e^{-\frac{\tau}{\varepsilon}} \\ 0 & e^{-\tau} & 0 & 0 & 0 & -\frac{k+\varepsilon}{k} \\ 0 & 0 & \frac{(k+\beta+I)\sqrt{\varepsilon}}{k} & 0 & 0 & I_5 \\ 0 & 0 & \tilde{\tau}_{3,6} & 0 & -1 & I_6 \end{pmatrix} $$
where not all the entries has been explicitely computed. In fact, $\tilde{\tau}_{i,j}$ stands for $\frac{\partial \hat{\phi}_j}{\partial \tilde{\tau}_i}$ and $I_j$ stands for $\frac{\partial \hat{\phi}_j}{\partial I}$. Therefore,

$$ \det(A) = - \frac{(k + \beta + I)^2}{k} \left| \begin{array}{cc} \frac{k + \beta + I - 1 + e^{-\tau}}{\varepsilon} & 1 - e^{-\frac{\tau}{\varepsilon}} \\ e^{-\tau} & -\frac{k + \varepsilon}{k} \end{array} \right| \quad (C5) $$

where, $(k + \beta + I) \neq 0$ and the second term is precisely the determinant calculated in Lemma 3, which is non-zero.

## References

[1] Rinzel, J.: A formal classification of bursting mechanisms in excitable systems. In: International Congress of Mathematicians, Berkeley, California, USA, August 3-11, 1986, vol. II, pp. 1578-1593. American Mathematical Society, Providence (1987)

[2] Rotstein, H.G., Kopell, N., Zhabotinsky, A.M., Epstein, I.R.: Canard phenomenon and localization of oscillations in the belousov-zhabotinsky reaction with global feedback. The Journal of Chemical Physics 119(17), 8824-8832 (2003)

[3] Mbé, J.H.T., Talla, A.F., Chengui, G.R.G., Coillet, A., Larger, L., Woafo, P., Chembo, Y.K.: Mixed-mode oscillations in slow-fast delayed optoelectronic systems. Physical Review E 91(1), 012902 (2015)

[4] Rinaldi, S., Muratori, S.: Slow-fast limit cycles in predator-prey models. Ecological Modelling 61(3), 287-308 (1992). https://doi.org/10.1016/0304-3800(92)90023-8

[5] Rajesh, S., Ananthakrishna, G.: Relaxation oscillations and negative strain rate sensitivity in the portevin-le chatelier effect. Physical Review E 61(4), 3664 (2000)

[6] Guckenheimer, J., Hoffman, K., Weckesser, W.: The forced van der pol equation i: The slow flow and its bifurcations. SIAM Journal on Applied Dynamical Systems 2(1), 1-35 (2003). https://doi.org/10.1137/S1111111102404738

[7] Bold, K., Edwards, C., Guckenheimer, J., Guharay, S., Hoffman, K., Hubbard, J., Oliva, R., Weckesser, W.: The forced van der pol equation ii: Canards in the reduced system. SIAM Journal on Applied Dynamical Systems [electronic only] 2 (2003). https://doi.org/10.1137/S1111111102419130
where not all the entries has been explicitely computed. In fact, $\tilde{\tau}_{i,j}$ stands for $\frac{\partial \hat{\phi}_j}{\partial \tilde{\tau}_i}$ and $I_j$ stands for $\frac{\partial \hat{\phi}_j}{\partial I}$. Therefore,

$$ \det(A) = - \frac{(k + \beta + I)^2}{k} \left| \begin{array}{cc} \frac{k + \beta + I - 1 + e^{-\tau}}{\varepsilon} & 1 - e^{-\frac{\tau}{\varepsilon}} \\ e^{-\tau} & -\frac{k + \varepsilon}{k} \end{array} \right| \quad (C5) $$

where, $(k + \beta + I) \neq 0$ and the second term is precisely the determinant calculated in Lemma 3, which is non-zero.

## References

[1] Rinzel, J.: A formal classification of bursting mechanisms in excitable systems. In: International Congress of Mathematicians, Berkeley, California, USA, August 3-11, 1986, vol. II, pp. 1578-1593. American Mathematical Society, Providence (1987)

[2] Rotstein, H.G., Kopell, N., Zhabotinsky, A.M., Epstein, I.R.: Canard phenomenon and localization of oscillations in the belousov-zhabotinsky reaction with global feedback. The Journal of Chemical Physics 119(17), 8824-8832 (2003)

[3] Mbé, J.H.T., Talla, A.F., Chengui, G.R.G., Coillet, A., Larger, L., Woafo, P., Chembo, Y.K.: Mixed-mode oscillations in slow-fast delayed optoelectronic systems. Physical Review E 91(1), 012902 (2015)

[4] Rinaldi, S., Muratori, S.: Slow-fast limit cycles in predator-prey models. Ecological Modelling 61(3), 287-308 (1992). https://doi.org/10.1016/0304-3800(92)90023-8

[5] Rajesh, S., Ananthakrishna, G.: Relaxation oscillations and negative strain rate sensitivity in the portevin-le chatelier effect. Physical Review E 61(4), 3664 (2000)

[6] Guckenheimer, J., Hoffman, K., Weckesser, W.: The forced van der pol equation i: The slow flow and its bifurcations. SIAM Journal on Applied Dynamical Systems 2(1), 1-35 (2003). https://doi.org/10.1137/S1111111102404738

[7] Bold, K., Edwards, C., Guckenheimer, J., Guharay, S., Hoffman, K., Hubbard, J., Oliva, R., Weckesser, W.: The forced van der pol equation ii: Canards in the reduced system. SIAM Journal on Applied Dynamical Systems [electronic only] 2 (2003). https://doi.org/10.1137/S1111111102419130
[8] Krupa, M., Sandstede, B., Szmolyan, P.: Fast and slow waves in the fitzhugh-nagumo equation. Journal of Differential Equations 133(1), 49-97 (1997). https://doi.org/10.1006/jdeq.1996.3198

[9] Kakiuchi, N., Tchizawa, K.: On an explicit duck solution and delay in the fitzhugh-nagumo equation. Journal of Differential Equations 141, 327-339 (1997)

[10] Konishi, K., Takeuchi, M., Shimizu, T.: Design of external forces for eliminating traveling wave in a piecewise linear fitzhugh-nagumo model. Chaos: An Interdisciplinary Journal of Nonlinear Science 21(2), 023101 (2011) https://doi.org/10.1063/1.3545162. https://doi.org/10.1063/1.3545162

[11] Sriram, K., Bernard, S.: Complex dynamics in the oregonator model with linear delayed feedback. Chaos: An Interdisciplinary Journal of Nonlinear Science 18(2), 023126 (2008)

[12] Rinzel, J., Ermentrout, G.: In: Koch, C., Segev, I. (eds.) Analysis of neural excitability and oscillations, 2nd edn., pp. 251-291. MIT Press, Cambridge (MA) (1989)

[13] Izhikevich, E.M.: Dynamical Systems In Neuroscience. MIT Press, Cambridge (MA) (2007)

[14] Desroches, M., Kaper, T.J., Krupa, M.: Mixed-mode bursting oscillations: Dynamics created by a slow passage through spike-adding canard explosion in a square-wave burster. Chaos: An Interdisciplinary Journal of Nonlinear Science 23(4), 046106 (2013)

[15] Rinzel, J.: A formal classification of bursting mechanisms in excitable systems. In: Teramoto, E., Yumaguti, M. (eds.) Mathematical Topics in Population Biology, Morphogenesis and Neurosciences. Lecture Notes in Biomathematics, vol. 71, pp. 267-281. Springer, Berlin (1987)

[16] Izhikevich, E.M.: Neural excitability, spiking and bursting. International Journal of Bifurcation and Chaos 10(06), 1171-1266 (2000)

[17] Bertram, R., Butte, M.J., Kiemel, T., Sherman, A.: Topological and phenomenological classification of bursting oscillations. Bulletin of Mathematical Biology 57(3), 413-439 (1995)

[18] Golubitsky, M., Josic, K., Kaper, T.J.: An unfolding theory approach to bursting in fast-slow systems. In: Broer, H.W., Krauskopf, B., Vegter, G. (eds.) Global Analysis of Dynamical Systems, pp. 277-308. IOP Publishing, Bristol (2001)
[8] Krupa, M., Sandstede, B., Szmolyan, P.: Fast and slow waves in the fitzhugh-nagumo equation. Journal of Differential Equations 133(1), 49-97 (1997). https://doi.org/10.1006/jdeq.1996.3198

[9] Kakiuchi, N., Tchizawa, K.: On an explicit duck solution and delay in the fitzhugh-nagumo equation. Journal of Differential Equations 141, 327-339 (1997)

[10] Konishi, K., Takeuchi, M., Shimizu, T.: Design of external forces for eliminating traveling wave in a piecewise linear fitzhugh-nagumo model. Chaos: An Interdisciplinary Journal of Nonlinear Science 21(2), 023101 (2011) https://doi.org/10.1063/1.3545162. https://doi.org/10.1063/1.3545162

[11] Sriram, K., Bernard, S.: Complex dynamics in the oregonator model with linear delayed feedback. Chaos: An Interdisciplinary Journal of Nonlinear Science 18(2), 023126 (2008)

[12] Rinzel, J., Ermentrout, G.: In: Koch, C., Segev, I. (eds.) Analysis of neural excitability and oscillations, 2nd edn., pp. 251-291. MIT Press, Cambridge (MA) (1989)

[13] Izhikevich, E.M.: Dynamical Systems In Neuroscience. MIT Press, Cambridge (MA) (2007)

[14] Desroches, M., Kaper, T.J., Krupa, M.: Mixed-mode bursting oscillations: Dynamics created by a slow passage through spike-adding canard explosion in a square-wave burster. Chaos: An Interdisciplinary Journal of Nonlinear Science 23(4), 046106 (2013)

[15] Rinzel, J.: A formal classification of bursting mechanisms in excitable systems. In: Teramoto, E., Yumaguti, M. (eds.) Mathematical Topics in Population Biology, Morphogenesis and Neurosciences. Lecture Notes in Biomathematics, vol. 71, pp. 267-281. Springer, Berlin (1987)

[16] Izhikevich, E.M.: Neural excitability, spiking and bursting. International Journal of Bifurcation and Chaos 10(06), 1171-1266 (2000)

[17] Bertram, R., Butte, M.J., Kiemel, T., Sherman, A.: Topological and phenomenological classification of bursting oscillations. Bulletin of Mathematical Biology 57(3), 413-439 (1995)

[18] Golubitsky, M., Josic, K., Kaper, T.J.: An unfolding theory approach to bursting in fast-slow systems. In: Broer, H.W., Krauskopf, B., Vegter, G. (eds.) Global Analysis of Dynamical Systems, pp. 277-308. IOP Publishing, Bristol (2001)
[19] Saggio, M.L., Spiegler, A., Bernard, C., Jirsa, V.K.: Fast-slow bursters in the unfolding of a high codimension singularity and the ultra-slow transitions of classes. The Journal of Mathematical Neuroscience **7**(1), 7 (2017)

[20] Desroches, M., Rinzel, J., Rodrigues, S.: Classification of bursting patterns: A tale of two ducks. PLoS Computational Biology **18**(2), 1009752 (2022)

[21] Baer, S., Erneux, T., Rinzel, J.: The slow passage through a hopf bifurcation: Delay, memory effects, and resonance. SIAM Journal on Applied Mathematics **49**(1), 55-71 (1989)

[22] Neishtadt, A.I.: Persistence of stability loss for dynamical bifurcations I. Differential Equations **23**, 1385-1391 (1987)

[23] Neishtadt, A.I.: Persistence of stability loss for dynamical bifurcations II. Differential Equations **24**, 171-176 (1988)

[24] Neishtadt, A.I.: On stability loss delay for dynamical bifurcations. Discrete & Continuous Dynamical Systems - Series S **2**(4), 897 (2009)

[25] Diminnie, D.C., Haberman, R.: Slow passage through homoclinic orbits for the unfolding of a saddle-center bifurcation and the change in the adiabatic invariant. Physica D: Nonlinear Phenomena **162**(1-2), 34-52 (2002)

[26] Benoît, E., Callot, J.-L., Diener, F., Diener, M.: Chasse au canard. Collectanea Mathematica **32**(1-2), 37-119 (1981)

[27] Dumortier, F., Roussarie, R.: Canard cycles and center manifolds. Memoirs of the American Mathematical Society **577** (1996)

[28] Krupa, M., Szmolyan, P.: Relaxation oscillation and canard explosion. Journal of Differential Equations **174**(2), 312-368 (2001)

[29] Kuznetsov, Y.A., Muratori, S., Rinaldi, S.: Homoclinic bifurcations in slow-fast second order systems. Nonlinear Analysis: Theory, Methods & Applications **25**(7), 747-762 (1995). https://doi.org/10.1016/0362-546X(94)E0005-2

[30] Matsue, K.: Rigorous numerics for fast-slow systems with one-dimensional slow variable: Topological shadowing approach. Topological Methods in Nonlinear Analysis **49** (2015). https://doi.org/10.12775/TMNA.2016.072

[31] Carter, P.: Spike-adding canard explosion in a class of square-wave bursters. Journal of Nonlinear Science **30** (2020). https://doi.org/10.
[19] Saggio, M.L., Spiegler, A., Bernard, C., Jirsa, V.K.: Fast-slow bursters in the unfolding of a high codimension singularity and the ultra-slow transitions of classes. The Journal of Mathematical Neuroscience **7**(1), 7 (2017)

[20] Desroches, M., Rinzel, J., Rodrigues, S.: Classification of bursting patterns: A tale of two ducks. PLoS Computational Biology **18**(2), 1009752 (2022)

[21] Baer, S., Erneux, T., Rinzel, J.: The slow passage through a hopf bifurcation: Delay, memory effects, and resonance. SIAM Journal on Applied Mathematics **49**(1), 55-71 (1989)

[22] Neishtadt, A.I.: Persistence of stability loss for dynamical bifurcations I. Differential Equations **23**, 1385-1391 (1987)

[23] Neishtadt, A.I.: Persistence of stability loss for dynamical bifurcations II. Differential Equations **24**, 171-176 (1988)

[24] Neishtadt, A.I.: On stability loss delay for dynamical bifurcations. Discrete & Continuous Dynamical Systems - Series S **2**(4), 897 (2009)

[25] Diminnie, D.C., Haberman, R.: Slow passage through homoclinic orbits for the unfolding of a saddle-center bifurcation and the change in the adiabatic invariant. Physica D: Nonlinear Phenomena **162**(1-2), 34-52 (2002)

[26] Benoît, E., Callot, J.-L., Diener, F., Diener, M.: Chasse au canard. Collectanea Mathematica **32**(1-2), 37-119 (1981)

[27] Dumortier, F., Roussarie, R.: Canard cycles and center manifolds. Memoirs of the American Mathematical Society **577** (1996)

[28] Krupa, M., Szmolyan, P.: Relaxation oscillation and canard explosion. Journal of Differential Equations **174**(2), 312-368 (2001)

[29] Kuznetsov, Y.A., Muratori, S., Rinaldi, S.: Homoclinic bifurcations in slow-fast second order systems. Nonlinear Analysis: Theory, Methods & Applications **25**(7), 747-762 (1995). https://doi.org/10.1016/0362-546X(94)E0005-2

[30] Matsue, K.: Rigorous numerics for fast-slow systems with one-dimensional slow variable: Topological shadowing approach. Topological Methods in Nonlinear Analysis **49** (2015). https://doi.org/10.12775/TMNA.2016.072

[31] Carter, P.: Spike-adding canard explosion in a class of square-wave bursters. Journal of Nonlinear Science **30** (2020). https://doi.org/10.
1007/s00332-020-09631-y

[32] Guckenheimer, J., Kuehn, C.: Computing slow manifolds of saddle type. SIAM Journal on Applied Dynamical Systems **8**(3), 854–879 (2009)

[33] Fenichel, N.: Geometric singular perturbation theory for ordinary differential equations. Journal of Differential Equations **31**(1), 53–98 (1979)

[34] Terman, D.: Chaotic spikes arising from a model of bursting in excitable membranes. SIAM Journal on Applied Mathematics **51**(5), 1418–1450 (1991)

[35] Channell, P., Cymbalyuk, G., Shilnikov, A.: Origin of bursting through homoclinic spike adding in a neuron model. Physical Review Letters **98**(13), 134101 (2007)

[36] Hindmarsh, J.L., Rose, R.M.: A model of neuronal bursting using three coupled first order differential equations. Proceedings of the Royal Society of London. Series B, Biological Sciences **221**(1222), 87–102 (1984)

[37] Shilnikov, A., Kolomiets, M.: Methods of the qualitative theory for the hindmarsh-rose model: A case study-a tutorial. International Journal of Bifurcation and chaos **18**(08), 2141–2168 (2008)

[38] Lecar, H.: Morris-lecar model. Scholarpedia **2**(10), 1333 (2007)

[39] Morris, C., Lecar, H.: Voltage oscillations in the barnacle giant muscle fiber. Biophysical Journal **35**(1), 193–213 (1981)

[40] Llibre, J., Ponce, E., Teruel, A.E.: Horseshoes near homoclinic orbits for piecewise linear differential systems in $\mathbb{R}^3$. International Journal of Bifurcation and Chaos **17**, 1171–1184 (2007). https://doi.org/10.1142/S0218127407017756

[41] Carmona, V., Fernández-García, S., Fernández-Sánchez, F., García-Medina, E., Teruel, A.E.: Reversible periodic orbits in a class of 3d continuous piecewise linear systems of differential equations. Nonlinear Analysis: Theory, Methods & Applications **75**(15), 5866–5883 (2012)

[42] Carmona, V., Fernández-García, S., Teruel, A.E.: Saddle-node of limit cycles in planar piecewise linear systems and applications. Discrete & Continuous Dynamical Systems **39**(9), 5275–5299 (2019)

[43] Desroches, M., Guillamon, A., Ponce, E., Prohens, R., Rodrigues, S., Teruel, A.E.: Canards, folded nodes, and mixed-mode oscillations in piecewise-linear slow-fast systems. SIAM Review **58**(4), 653–691 (2016) https://doi.org/10.1137/15M1014528. https://doi.org/10.1137/15M1014528
1007/s00332-020-09631-y

[32] Guckenheimer, J., Kuehn, C.: Computing slow manifolds of saddle type. SIAM Journal on Applied Dynamical Systems **8**(3), 854–879 (2009)

[33] Fenichel, N.: Geometric singular perturbation theory for ordinary differential equations. Journal of Differential Equations **31**(1), 53–98 (1979)

[34] Terman, D.: Chaotic spikes arising from a model of bursting in excitable membranes. SIAM Journal on Applied Mathematics **51**(5), 1418–1450 (1991)

[35] Channell, P., Cymbalyuk, G., Shilnikov, A.: Origin of bursting through homoclinic spike adding in a neuron model. Physical Review Letters **98**(13), 134101 (2007)

[36] Hindmarsh, J.L., Rose, R.M.: A model of neuronal bursting using three coupled first order differential equations. Proceedings of the Royal Society of London. Series B, Biological Sciences **221**(1222), 87–102 (1984)

[37] Shilnikov, A., Kolomiets, M.: Methods of the qualitative theory for the hindmarsh-rose model: A case study-a tutorial. International Journal of Bifurcation and chaos **18**(08), 2141–2168 (2008)

[38] Lecar, H.: Morris-lecar model. Scholarpedia **2**(10), 1333 (2007)

[39] Morris, C., Lecar, H.: Voltage oscillations in the barnacle giant muscle fiber. Biophysical Journal **35**(1), 193–213 (1981)

[40] Llibre, J., Ponce, E., Teruel, A.E.: Horseshoes near homoclinic orbits for piecewise linear differential systems in $\mathbb{R}^3$. International Journal of Bifurcation and Chaos **17**, 1171–1184 (2007). https://doi.org/10.1142/S0218127407017756

[41] Carmona, V., Fernández-García, S., Fernández-Sánchez, F., García-Medina, E., Teruel, A.E.: Reversible periodic orbits in a class of 3d continuous piecewise linear systems of differential equations. Nonlinear Analysis: Theory, Methods & Applications **75**(15), 5866–5883 (2012)

[42] Carmona, V., Fernández-García, S., Teruel, A.E.: Saddle-node of limit cycles in planar piecewise linear systems and applications. Discrete & Continuous Dynamical Systems **39**(9), 5275–5299 (2019)

[43] Desroches, M., Guillamon, A., Ponce, E., Prohens, R., Rodrigues, S., Teruel, A.E.: Canards, folded nodes, and mixed-mode oscillations in piecewise-linear slow-fast systems. SIAM Review **58**(4), 653–691 (2016) https://doi.org/10.1137/15M1014528. https://doi.org/10.1137/15M1014528
[44] Desroches, M., Fernández-García, S., Krupa, M.: Canards in a minimal piecewise-linear square-wave burster. Chaos: An Interdisciplinary Journal of Nonlinear Science 26(7), 073111 (2016) https://doi.org/10.1063/1.4958297. https://doi.org/10.1063/1.4958297

[45] Prohens, R., Teruel, A.E., Vich, C.: Slow-fast n-dimensional piecewise linear differential systems. Journal of Differential Equations 260(2), 1865–1892 (2016). https://doi.org/10.1016/j.jde.2015.09.046

[46] Desroches, M., Fernández-García, S., Krupa, M., Prohens, R., Teruel, A.E.: Piecewise-linear (pwl) canard dynamics - simplifying singular perturbation theory in the canard regime using piecewise-linear systems. In: Carmona, V., Cuevas-Maraver, J., Fernández-Sánchez, F., García-Medina, E. (eds.) Nonlinear Systems, Vol. 1: Mathematical Theory and Computational Methods, pp. 67–86. Springer, Cham (2018)

[47] Penalva, J., Desroches, M., Teruel, A.E., Vich, C.: Slow passage through a hopf-like bifurcation in piecewise linear systems: Application to elliptic bursting. Chaos: An Interdisciplinary Journal of Nonlinear Science 32(12), 123109 (2022) https://doi.org/10.1063/5.0101778. https://doi.org/10.1063/5.0101778

[48] Tonnelier, A., Gerstner, W.: Piecewise linear differential equations and integrate-and-fire neurons: insights from two-dimensional membrane models. Phys Rev E Stat Nonlin Soft Matter Phys 67(2 Pt 1), 021908 (2003). https://doi.org/10.1103/PhysRevE.67.021908

[49] Coombes, S.: Neuronal networks with gap junctions: A study of piecewise linear planar neuron models. SIAM Journal on Applied Dynamical Systems 7(3), 1101–1129 (2008) https://doi.org/10.1137/070707579. https://doi.org/10.1137/070707579

[50] Arima, N., Okazaki, H., Nakano, H.: A generation mechanism of canards in a piecewise linear system. Trans. IEICE, A 80(3), 447–453 (1997)

[51] Fernández-García, S., Desroches, M., Krupa, M., Teruel, A.: Canard solutions in planar piecewise linear systems with three zones. Dynamical Systems An International Journal 31, 173–197 (2016). https://doi.org/10.1080/14689367.2015.1079304

[52] Carmona, V., Fernández-García, S., Teruel, A.E.: Saddle-node canard cycles in planar piecewise linear differential systems (2020) arXiv:2003.14112 [math.DS]

[53] Kuznetsov, Y.A.: Elements of Applied Bifurcation Theory, Second Edition, 2nd edn. Applied Mathematical Sciences. Springer, ??? (1998). libgen.li/file.php?md5=4ccf8c905f92abfbed3aa98edeec5dbf
[44] Desroches, M., Fernández-García, S., Krupa, M.: Canards in a minimal piecewise-linear square-wave burster. Chaos: An Interdisciplinary Journal of Nonlinear Science 26(7), 073111 (2016) https://doi.org/10.1063/1.4958297. https://doi.org/10.1063/1.4958297

[45] Prohens, R., Teruel, A.E., Vich, C.: Slow-fast n-dimensional piecewise linear differential systems. Journal of Differential Equations 260(2), 1865–1892 (2016). https://doi.org/10.1016/j.jde.2015.09.046

[46] Desroches, M., Fernández-García, S., Krupa, M., Prohens, R., Teruel, A.E.: Piecewise-linear (pwl) canard dynamics - simplifying singular perturbation theory in the canard regime using piecewise-linear systems. In: Carmona, V., Cuevas-Maraver, J., Fernández-Sánchez, F., García-Medina, E. (eds.) Nonlinear Systems, Vol. 1: Mathematical Theory and Computational Methods, pp. 67–86. Springer, Cham (2018)

[47] Penalva, J., Desroches, M., Teruel, A.E., Vich, C.: Slow passage through a hopf-like bifurcation in piecewise linear systems: Application to elliptic bursting. Chaos: An Interdisciplinary Journal of Nonlinear Science 32(12), 123109 (2022) https://doi.org/10.1063/5.0101778. https://doi.org/10.1063/5.0101778

[48] Tonnelier, A., Gerstner, W.: Piecewise linear differential equations and integrate-and-fire neurons: insights from two-dimensional membrane models. Phys Rev E Stat Nonlin Soft Matter Phys 67(2 Pt 1), 021908 (2003). https://doi.org/10.1103/PhysRevE.67.021908

[49] Coombes, S.: Neuronal networks with gap junctions: A study of piecewise linear planar neuron models. SIAM Journal on Applied Dynamical Systems 7(3), 1101–1129 (2008) https://doi.org/10.1137/070707579. https://doi.org/10.1137/070707579

[50] Arima, N., Okazaki, H., Nakano, H.: A generation mechanism of canards in a piecewise linear system. Trans. IEICE, A 80(3), 447–453 (1997)

[51] Fernández-García, S., Desroches, M., Krupa, M., Teruel, A.: Canard solutions in planar piecewise linear systems with three zones. Dynamical Systems An International Journal 31, 173–197 (2016). https://doi.org/10.1080/14689367.2015.1079304

[52] Carmona, V., Fernández-García, S., Teruel, A.E.: Saddle-node canard cycles in planar piecewise linear differential systems (2020) arXiv:2003.14112 [math.DS]

[53] Kuznetsov, Y.A.: Elements of Applied Bifurcation Theory, Second Edition, 2nd edn. Applied Mathematical Sciences. Springer, ??? (1998). libgen.li/file.php?md5=4ccf8c905f92abfbed3aa98edeec5dbf
[54] Perko, L.: Differential Equations and Dynamical Systems. Texts in Applied Mathematics. Springer, New York (1996). https://books.google.es/books?id=-j3vAAAAMAAJ

[55] Prohens, R., Teruel, A.E., Vich, C.: Slow-fast n-dimensional piecewise linear differential systems. Journal of Differential Equations **260**(2), 1865–1892 (2016)

[56] Euzébio, R.D., Llibre, J.: On the number of limit cycles in discontinuous piecewise linear differential systems with two pieces separated by a straight line. Journal of Mathematical Analysis and Applications **424**(1), 475–486 (2015). https://doi.org/10.1016/j.jmaa.2014.10.077

[57] Freire, E., Ponce, E., Torres, F.: Canonical discontinuous planar piecewise linear systems. SIAM Journal on Applied Dynamical Systems **11**, 181–211 (2012). https://doi.org/10.1137/11083928X

[58] Medrado, J.C., Torregrosa, J.: Uniqueness of limit cycles for sewing piecewise linear systems. J. Math. Anal. Appl. (431), 529–544 (2015)

[59] Di Bernardo, M., Champneys, A., Budd, C., Kowalczyk, P.: Piecewise-smooth Dynamical Systems: Theory and Applications, (2008). https://doi.org/10.1007/978-1-84628-708-4

[60] Llibre, J., Núñez, E., Teruel, A.: Limit cycles for planar piecewise linear differential systems via first integrals. Qualitative Theory of Dynamical Systems **3**, 29–50 (2002). https://doi.org/10.1007/BF02969332
[54] Perko, L.: Differential Equations and Dynamical Systems. Texts in Applied Mathematics. Springer, New York (1996). https://books.google.es/books?id=-j3vAAAAMAAJ

[55] Prohens, R., Teruel, A.E., Vich, C.: Slow-fast n-dimensional piecewise linear differential systems. Journal of Differential Equations **260**(2), 1865–1892 (2016)

[56] Euzébio, R.D., Llibre, J.: On the number of limit cycles in discontinuous piecewise linear differential systems with two pieces separated by a straight line. Journal of Mathematical Analysis and Applications **424**(1), 475–486 (2015). https://doi.org/10.1016/j.jmaa.2014.10.077

[57] Freire, E., Ponce, E., Torres, F.: Canonical discontinuous planar piecewise linear systems. SIAM Journal on Applied Dynamical Systems **11**, 181–211 (2012). https://doi.org/10.1137/11083928X

[58] Medrado, J.C., Torregrosa, J.: Uniqueness of limit cycles for sewing piecewise linear systems. J. Math. Anal. Appl. (431), 529–544 (2015)

[59] Di Bernardo, M., Champneys, A., Budd, C., Kowalczyk, P.: Piecewise-smooth Dynamical Systems: Theory and Applications, (2008). https://doi.org/10.1007/978-1-84628-708-4

[60] Llibre, J., Núñez, E., Teruel, A.: Limit cycles for planar piecewise linear differential systems via first integrals. Qualitative Theory of Dynamical Systems **3**, 29–50 (2002). https://doi.org/10.1007/BF02969332