Relaxation oscillations and canards in the Jirsa–Kelso excitator model: global flow perspective

Fenichel’s geometric singular perturbation theory and the blow-up method have been very successful in describing and explaining global non-linear phenomena in systems with multiple time-scales, such as relaxation oscillations and canards. Recently, the blow-up method has been extended to systems with flat, unbounded slow manifolds that lose normal hyperbolicity at infinity. Here, we show that transition between discrete and periodic movement captured by the Jirsa–Kelso excitator is a new example of such phenomena. We, first, derive equations of the Jirsa–Kelso excitator with explicit time scale separation and demonstrate existence of canards in the systems. Then, we combine the slow-fast analysis, blow-up method and projection onto the Poincaré sphere to understand the return mechanism of the periodic orbits in the singular case, ϵ = 0.


Introduction
The Jirsa-Kelso excitator model is a class of excitable planar systems proposed as a minimal model to describe generation of rhythmic and discrete human movement [7]. The model predictions concerning the mechanism of discrete movement generation have been tested in a number of empirical studies [3,6].
The model of a single oscillator is given as a system of two ODEs (Eq. (12) in [7]): a e-mail: K.Tsaneva-Atanasova@exeter.ac.uk

592
The European Physical Journal Special Topics Model (1) has mixed time scales. To separate them we, first, rearrange the second ODE, factor out −T v(t) and substitute the time constant T = 1/ √ , obtaining, In (1) and (2), u(t) is interpreted as a position and v(t) as a velocity of the movement, a and b are intrinsic model parameters, T and are a time scale constants and I(t) is an external stimulus input (in the rest of the paper we consider the autonomous systems, i.e. I = 0). Then, the mixed time scales of the system (2) can be separated by rescaling the time t and the variable v, respectively: The resulting system, with explicitly separated time scales, has the classical fast representation, x = f (x, y, ), y = g(x, y, ), where (·) = d dτ denotes differentiation with respect to the fast time τ [9,10]. In the rest of the paper we refer to (4) as the Jirsa-Kelso excitator (JKE). The modelling approach introduced in [7] is motivated by excitable systems with time-scale separation. In fact, the transformation given in [7], allows to transform the JKE (4) into the seminal example of such a system, the FitzHugh-Nagumo model (FHN) (with I = 0, z → −z and a → −a) [4,14,15], Figure 1 illustrates the transformation (5) between the FHN (6) and the JKE (4) models. Figure 1a shows the hypersurface H = {(x, y, z) ∈ R 3 : y = z + x − x 3 3 } together with projection of the common x-nullcline of the JKE (4) and the FHN (6) (blue curve) and y-nullcline of the JKE (4) (red curve). Additionally, it shows two representative periodic orbits (dark, a ≈ 0.78978455, and light, a ≈ 0.78978433, grey curves). Figures 1b-1d demonstrate the invariance of the x variable, and the difference in time-courses of variables y and z; the time-series correspond to the two periodic orbits in Figure 1a. Although, the variable x is invariant under the transformation (5), it has different physical interpretation in the two models. In the JKE (4) x describes position and hence allows to interpret the variable y as velocity, while in the FHN (6) x has been commonly interpreted as the membrane potential of a neuronal cell, while y represents a recovery variable associated with the slow dynamics (usually attributed to ion channels' gating) that controls the generation of action potentials [4].
Since the transformation (5) is homeomorphism (for > 0) the intrinsic dynamics of the JKE (4) and FHN (6) systems are equivalent for > 0 [7]. The equivalence between the FHN (6) and the JKE (4) means that the existing results related to the dynamics of the FHN model, e.g. existence of canard orbits or relaxation oscillations [14,15], can be extended to the JKE model. However, the singularity in (5) for = 0 implies that the global mechanism responsible for the canard cycles and relaxation oscillations in the JKE is different than the one in the FHN; i.e. the two systems have different critical manifolds. In the rest of the paper we present a detailed study of the critical manifold C 0 of the uncoupled JKE model.
2 Critical manifold C 0 of the Jirsa-Kelso excitator model As in the case of the FHN model the relaxation oscillations in the JKE appear through canard explosion, a rapid growth of periodic orbits' amplitude that happens in an exponentially small range of the control parameter [10]. However, in contrast to the FHN model, in the JKE model the canard explosion is organised by a non-generic branching point at the y-nullcline [12], rather than by a fold point [5,10]. Figure 2 depicts the canard explosion in the JKE model. Figure 2a shows the bifurcation diagram in the parameter a and Figure 2b shows the (x, y)-phase portrait, in both at which the x-coordinate is given as x H = √ 1 − b; the methods for computing of a H and x H can be found for example in [1]. The bifurcation diagram and the phase portrait in Figure 2 have been computed using continuation software AUTO [2].
We, now, focus on analysing the dynamics of the JKE system (4) in the singular limit = 0, to this end we follow [1,10]. Setting = 0 in the fast system (4) defines the layer problem of the JKE system: Since x = 0, we consider (7) as one-dimensional system with flow in the y-direction parameterised by the values of x. In other words, each value of x defines flow (7) along 1D invariant manifold (layer) in the xy-plane. The critical manifold of the JKE is given by the set of equilibria of the layer problem (7), The curve C 0 has three branches separated by two vertical asymptotes located at x = ±1. Note that, the critical manifold C 0 is the → 0 limit of the y-nullcline of JKE system (4), which has asymptotes at x = ∓ √ 1 − b. The stability of the equilibria that form C 0 can be determined from the sign of the eigenvalue λ = 1 − x 2 of the layer problem (7). For x ∈ (−∞, −1) ∪ (1, ∞) the equilibria are stable, and for x ∈ (−1, 1) they are unstable, for x = ±1 λ = 0 and C 0 is non-hyperbolic. Similarly to the FHN model the fast subsystem of (7) is tangent to the C 0 at x = ±1, however, in the case of JKE the tangency occurs at infinity. At the singular Hopf bifurcation, a H0 = a H ( = 0) = 1 − 2b/3, the middle and right parts of the critical manifold connect via a non-generic branching point [12] at (x, y) = (1, −1/2) and the critical manifold takes the form: Note that, for x = ±1 the stationary solutions of (7) are given as y(t; x) = (∓(1 − 2b/3) + a)t + const., meaning that, for a = a H0 , there is no vertical flow on the x = 1 layer.
To analyse the flow on the critical manifold C 0 , we rescale (4) by changing the fast time τ to the slow time t s = τ , In this way we obtain the so-called slow or reduced system:ẋ = f (x, y, ), ẏ = g(x, y, ), where( ·) = d dts denotes differentiation with respect to the slow time t s . Setting = 0 and eliminating y from the differential equation (10), results in a one dimensional systemẋ For 0 < b < 1 and −1 < a < 1, equation (10) has only one equilibrium, which corresponds to the equilibrium of the JKE. The equilibrium is stable for a > a H0 and unstable for a < a H0 ; b ∈ (0, 1). The layer problem and flow on the critical manifold are illustrated in Figure 3. For a = a H0 , at the singular Hopf bifurcation, there is no flow along vertical branch, x = 1, of the critical manifold. On the other branch of the critical manifold the flow is from x = ∞ to x = −1. The two branches intersect at the non-generic branching point, C 0 (x = 1) = −1/2. The flow through this point can be studied formally by blowing up the y-axis to the cylinderx 2 +¯ 2 = 1 using cylindrical blow-up transformation [9].

Global flow of the Jirsa-Keslo excitator model
In this section, we investigate the global behaviour of the layer and reduced problems of the JKE model (4). To this end, we represent the JKE (4) on the Poincaré sphere using the following transformation of variables, Transformation (12) defines a one-to-one correspondence between the points (x, y) in the plane and the points (X, Y, Z) on one of the hemisphere (we consider Z ≥ 0); it follows that x = X/Z and y = Y /Z. The advantage of using the Poincaré sphere is that the critical points at infinity are spread out along its equator, Z = 0, X 2 + Y 2 = 1. Following [13] (Theorem 1 in chapter 3.10), the equilibria and flow of the projected system on the equator, Z = 0, can be determined using the following equation, where, and, P (x, y) = y, is the JKE system (4); see [13] for details of derivation. Solving equation (13) on the equator, G(X, Y ) = 0, X 2 + Y 2 = 1, yields the four critical points where the flow between these critical points is determined by the sign of (13). For a < a H0 the two points X = 0 are unstable nodes, and the two for X = 0 are saddles. In the singular case for a = a H0 and = 0 the two points at X = 0 change from unstable nodes to non-hyperbolic nodes; the two nodes (0, ±1) provide return mechanism for the unbounded segments of the singular canards. The flow of the JKE system on the Poincaré hemisphere for Z ≥ 0 (projected orthogonally on (X, Y )-plane) is shown in Figure 4.
To understand the return mechanism at the points (0, ±1, 0) for = 0, we project the flow on a plane tangent to the sphere at these points (X, ±1, Z), and then we blow-up the non-hyperbolic nodes. We first observe that in the tangent plane, the flow in the neighbourhood of the non-hyperbolic nodes (0, ±1, 0) is radial. Specifically, from [13] (Theorem 2 in chapter 3.10) we can deduce that in the neighbourhood of (0, ±1, 0) the behaviour of the JKE model is equivalent to the behaviour of the system given by, with respective sign determined by the flow on the equator. Since = 0, P * = 0, hence the system (16) is a separable equation dZ/dX = Z/X with solution Z = CX for C ∈ R, meaning that in the neighbourhood of (0, ±1, 0) the flow has only radial direction.
To further investigate the dynamics near the projected nodes (0, ±1, 0), and to understand how they can be attracting and repelling at the same time we use polar 598 The European Physical Journal Special Topics blow-up transformation [11].
Branch B 1 changes stability at the points of intersection with branches B 2 and B 3 . Orthogonal intersection with the branch B 2 , (r, θ) = (0, π/4) (red cross in Fig. 5c), is a non-generic branching point that is of the same nature as the point (x, y) = (1, −1/2). Branches B 1 and B 3 intersect and exchange stability in a transcritical bifurcation at (r, θ) = (0, 3π/4) (red square in Fig. 5c). The intersection of branches B 1 and B 2 can be analysed formally by means of the cylindrical blowup transformation [8], while the intersection of branches B 1 and B 3 can be analysed formally by the classical techniques described in [9].
To complete the analysis of the critical manifold, we determine the direction of the flow of the reduced problem along the branch B 1 . We use again Theorem 2 from chapter 3.10 of [13] and blow-up method. Since, for the singular case = 0, the flow is defined only on the critical manifold, equations (15) for the reduced problem of the JKE model simplify to P (x, y) = y and Q(x, y) = 0. Hence, system (16) for the reduced problem is given by, After polar blow-up Φ (17), and rescaling time t s → rt s , the flow along the branch B 1 can be determined from the system, r = cos(θ) sin 2 (θ)r,

600
The European Physical Journal Special Topics Therefore, the reduced problem at infinity has two non-hyperbolic nodes on the blown-up circle, the branch B 1 , (r, θ) = (0, 0) and (r, θ) = (0, π) and the flow between the two nodes is clockwise. The behaviour near the node (0, −1, 0) can be analysed in the same way as the node (0, 1, 0), compare Figures 5b and 5d. After projection on the tangent plane and polar blow-up; the layer problem is given as, while the dynamics on the blown-up circle is described by, The layer problem at infinity (24) has again three branches of equilibria. Branches B 1 and B 2 , have the same stability and interpretation as the equivalent branches near the (0, 1, 0) node. Branch B 3 again corresponds to the projected critical manifold and has the same eigenvalue λ = −σ(θ), but near the node (0, −1, 0) it is defined over θ ∈ [0, 3π/4] rather than θ ∈ [3π/4, π) and hence is radially repelling. Branches B 1 and B 3 exchange stability at point (r, θ) = (0, 3π/4) in a transcritical bifurcation (red square in Fig. 5d). The reduced problem at infinity (25) has two non-hyperbolic nodes on the blown-up circle, (r, θ) = (0, 0) and (r, θ) = (0, π). The flow of the slow subsystem between the two nodes is counterclockwise.

Conclusion
In this paper we presented a detailed analysis of the mechanism that leads to appearance of relaxation oscillations in the Jirsa-Kelso excitator. We, first, demonstrated existence of canards in the Jirsa-Kelso excitator and showed that they are organised by a non-generic branching point at the y-nullcline of the system. We, then, combined projection onto the Poincaré sphere, slow-fast analysis and blow-up method to investigate the return mechanism of the periodic orbits in the singular case, = 0. The analysis revealed that the return mechanism is organised by two non-hyperbolic nodes at infinity, for which hyperbolicity can be recovered with a polar blow-up transformation. The blown-up flow at infinity showed the singular canard cycles can indeed be represented as concatenation of orbits, respecting the direction of time, of the layer and reduced problems. Overall, the presented analysis shows that although the Jirsa-Kelso excitator can be transformed into the FitzHugh-Nagumo model via a homeomorphic transformation (for > 0), the two models have different critical manifolds (for = 0) and different mechanisms of the global transitions. From applications point of view, the presented analysis shows that the slow-fast nature of the transition between discrete and rhythmic movements described by the JKE can be analysed and understood, not only in the abstract phase space of the FHN, but also in the experimentally relevant phase-space given by the position, u(t), and velocity, v(t), of the movement. Furthermore, our analysis suggests an alternative scenario, in addition to changes in the parameters a and b, for the transition between discrete and rhythmic movements. In this scenario, the transition is controlled by the time scale parameter . For small values of 1 the relaxation oscillations have a very long period and can be viewed as two almost steady states connected by fast