HOMOCLINIC CLUSTERS AND CHAOS ASSOCIATED WITH A FOLDED NODE IN A STELLATE CELL MODEL

. Acker et al ( J. Comp. Neurosci. , 15 , pp. 71–90, 2003) developed a model of stellate cells which reproduces qualitative oscillatory patterns known as mixed mode oscillations observed in experiments. This model includes dif-ferent time scales and can therefore be viewed as a singularly perturbed system of diﬀerential equations. The bifurcation structure of this model is very rich, and includes a novel class of homoclinic bifurcation points. The key to the bifurcation analysis is a folded node singularity that allows trajectories known as canards to cross from a stable slow manifold to an unstable slow manifold as well as a node equilibrium of the slow ﬂow on the unstable slow manifold. In this work we focus on the novel homoclinic orbits within the bifurcation diagram and show that the return of canards from the unstable slow manifold to the funnel of the folded node on the stable slow manifold results in a horseshoe map, and therefore gives rise to chaotic invariant sets. We also use a one-dimensional map to explain why many homoclinic orbits occur in “clusters” at exponentially close parameter values.


1.
Introduction. Stellate cells, the predominant cell type in layer II of the medial entorhinal cortex, control information flow to the hippocampus and they produce robust limbic rhythms in the theta frequency range (4)(5)(6)(7)(8)(9)(10)(11)(12). Electrophysiological studies [7] show that, when depolarized, stellate cells develop subthreshold oscillatory (STO) patterns and MMO patterns, small amplitude STOs interspersed by large amplitude oscillations. These observed oscillatory patterns are an intrinsic cell phenomenon and result mainly from the interaction between a fast persistent sodium current I Nap and a slow hyperpolarization activated current I h [7]. The stellate cell model by Acker et al [1] incorporates these currents; it consists of the current balance equation where V is the membrane potential (mV ), C is the membrane capacitance (µF/cm 2 ), I app is an applied external current (µA/cm 2 ), I Na is a sodium current, I K is the delayed rectifier current and I ℓ is the leak current. The persistent sodium current I Nap and the hyperpolarization activated current I h are specific for the STOs of stellate cells. The ionic currents are defined via Ohm's law by where the parameters G y and E y , y ∈ {Na, K, l, h, Nap}, are the maximal conductances (mS/cm 2 ) and reversal potentials (mV ) respectively (all parameter values are given in [1]). The conductances of the ionic currents are voltage dependent. This is modelled via dimensionless gating variables (m, h) for I Na , n for I K , (r f , r s ) for I h and p for I Nap . For x ∈ {m, h, n, p, r f , r s }, the dynamics of these gates where x ∞ (V ), the voltage dependent activation/inactivation curves, and τ x (V ), the voltage dependent time constants (see [1] for details). Equations (1) and (3) define the full seven dimensional model of a stellate cell (7D SC model).
In recent work [33], we presented a comprehensive bifurcation analysis of the 7D SC model. Figure 1 shows a bifurcation diagram obtained by using the software package AUTO [9]. There is a curve of equilibria, with a subcritical Hopf bifurcation in the upper branch of the curve near I app = −2.7, and a saddle-node bifurcation near I app = −2.05. In the upper branch of equilibria, those to the left of the Hopf bifurcation are stable, and those to the right have a two-dimensional unstable manifold. The equilibria in the lower branch are saddles with one-dimensional unstable manifolds.
Also shown in Figure 1 are families of periodic orbits. One family, which we denote by 0 1 , arises from the subcritical Hopf bifurcation near I app = −2.7. Along this branch of periodic orbits we also observe period-doubling bifurcations. Another family, labeled 1 0 , comes from another Hopf bifurcation (near I app = 60.7) that is not shown in the figure. The figure also includes several families of 1 s MMO patterns 1 . Families of MMOs were found by first using an initial value problem solver to find a stable periodic orbit. This orbit was used as a starting point for the AUTO continuation software. Several families of STOs also arise at SNPO bifurcations near the turning point in the 0 1 curve just below the Hopf bifurcation near I app = −2.7. For all these families of periodic orbits, AUTO was used to continue the orbits until the period exceeded 10 7 time units; we considered this to be strong numerical evidence that the family converged to a homoclinic orbit. See [33] for additional details.
A key observation in Figure 1 is that all branches of periodic orbits including STOs and MMOs terminate in novel homoclinic bifurcation points, i.e. they terminate at the lower branch of equilibria. In this work we will focus on these homoclinic bifurcation points, especially on clustering and associated chaotic dynamics. In [33], we explained all the other observed bifurcations within Figure 1. Therefore, this work extends our previous analysis.  Figure 1. Bifurcation diagram showing the equilibria, the families of periodic orbits that arise from Hopf bifurcations, and several families of 1 s mixed-mode oscillations. In all cases, the families terminate as homoclinic orbits; these are labeled using the notation described in section 2. Where more than one family appears to converge to the same point, the label is of the simplest homoclinic orbit. The numbered intervals in the lines at the bottom indicate the sectors in which φ + (p ε ) and φ − (p ε ) lie; see section 3 for details.
In Section 2 we give a summary of the bifurcation analysis presented in [33] and show examples of the homoclinic orbits observed. In Section 3 we introduce a reduced 3D SC model derived in [25] which describes the dynamics of the full 7D SC model in the subthreshold regime. Using this model, Rotstein, Wechselberger and Kopell [26] identified the canard phenomenon [4,11,21,22,31,32] to be the mechanism that generates STO and MMO patterns. In [33], we extended their work and explained the bifurcations of all MMO patterns, stable and unstable, in detail. In particular, we uncovered novel homoclinic bifurcation points as part of the canard phenomenon. This bifurcation structure in the 3D SC model is related to canards of folded node and folded saddle-node type. These type of canards are also studied in [3,5,6,14,15,16,19,23,24] and have been identified in other neuronal models as well [10,12,20,27,28]. In this paper, we review the relevant canard theory and define a global return map for the 3D SC model presented in [33] which allows us to analyse the observed bifurcations in the 7D SC model.
In section 4, we introduce a suitable Poincaré map which covers the dynamics of the 3D SC model and enables us to describe the general mechanism for bifurcations and chaos.
In Figure 1, there are several cases where it appears that multiple branches converge to a single point (e.g. two branches appear to converge to H − 0 ; multiple branches appear to converge to H − 1 , etc.), suggesting that these families might arise from a single homoclinic bifurcation. In fact, a closer look at the solutions shows that each branch appears to converge to a distinct homoclinic orbit. Specific examples are discussed in Section 2, and in Section 5, we give an argument that predicts these "clusters" of homoclinic orbits. We approximate the Poincaré map introduced in Section 4 by a one-dimensional (discontinuous) map. This is justified since the Poincaré map strongly contracts towards a one-dimensional subset due to the underlying slow/fast structure. Within that setting we show that homoclinic clusters exist. Finally, we conclude in Section 6.
2. Canard induced homoclinic orbits in the 7D SC model. Here we review the results of the bifurcation analysis performed in [33]. The 7D SC model has chaotic invariant sets, so a complete bifurcation diagram of the periodic orbits is not possible. Several of the families of the simplest MMOs are shown in Figure 1. All families of periodic orbits ultimately converge to orbits that are homoclinic to a saddle in the lower branch of equilibria. These homoclinic orbits are the main focus of this paper.
In [33], a shorthand notation was introduced to describe periodic orbits (STOs and MMOs). Each pattern is a string of symbols of the form X σ s , where X ∈ {J, C, H} describes a slow segment of a periodic orbit followed by a fast transition, σ ∈ {+, −} indicates whether, at the end of the slow segment, V quickly increases (+) or decreases (−), and the nonnegative integer s indicates the number of small O(ε) oscillations due to canards that occur near the fold. The meanings of the letters are: J The slow segment jumps at the fold. Because V always increases in a jump from the fold, the superscript σ is redundant. C The slow segment forms a canard, but the canard does not reach the saddle p. H The slow segment forms a canard that reaches the saddle; in other words, the orbit is a homoclinic connection. In Figure 1 we use this shorthand notation to denote basic homoclinic orbits observed in the model. A few of these homoclinic orbits are shown in (r s , V ) phase space in Figure 2. Clearly, slow canard segments are visible which bend around for large r s values to connect to the saddle equilibrium (shown as a red cross).  Figure 3. Three (approximate) homoclinic orbits from three different branches that appear to converge to the point labeled H − 1 in the bifurcation diagram, Figure 1. In phase space (right column) these homoclinics cannot be distinguished, but time traces (left column) reveal the different patterns (compare also with the H − 1 homoclinic in Figure 2). We show orbits with periods less than 4000, but the families were continued until the period exceeded 10 7 time units; plots of the (r s , V ) projections remain indistinguishable from those shown here.
Another important observation is that certain types of homoclinic orbits occur in "clusters". Figure 3 shows several homoclinic orbits that exist very close to the H − 1 homoclinic orbit in parameter space. These homoclinic orbits consist of more than one canard segment and they can only be distinguished by looking at the time traces. We will explain these clusters in section 5.
3. The reduced 3D SC model. It was shown by Rotstein, Wechselberger and Kopell [26] that the subthreshold dynamics (for V < −40mV ) of the 7D SC model are captured by a three dimensional model reduction (3D SC model) developed by Rotstein, et al [25]. In [33], we showed that this 3D SC model can also explain the bifurcations observed in Figure 1. The equations for the 3D SC model are Except for the applied current J app , we have maintained the same names for the variables and parameters as in the full 7D SC model. (see [25,33] for details). In [33], we provided convincing numerical evidence that this 3D SC model captures the qualitative behavior of the subthreshold dynamics of the full 7D SC model. In particular, we showed that the corresponding bifurcation diagrams (Figures 11 and  16 in [33]) are qualitatively the same.
Identifying different time scales within the 3D SC model by nondimensionalization is the key step in the analysis of the stellate cell model given in [26]. The corresponding singularly perturbed system is given by where v is the dimensionless membrane voltage and τ is the dimensionless time (see [26] for definition of all dimensionless parameters and functions).
The dimensionless 3D SC model (7)- (9) is a singularly perturbed system with small perturbation parameter ε. It has one fast variable, v, and two slow variables, (r f , r s ). Geometric singular perturbation theory [13,18] splits the analysis into two lower dimensional limiting problems ε → 0 of system (7)-(9), the two-dimensional reduced problem on the slow time scale τ and the one-dimensional layer problem on the fast time scale τ 1 = τ /ε, and pieces together the information obtained from these two limiting problems to provide a unified global description of the observed behavior in the 3D SC model.
Switching to the fast time scale τ 1 = τ /ε in system (7)-(9) and taking the limit ε → 0 gives the one-dimensional layer problem dv The two-dimensional manifold of equilibria of the layer problem, defined bȳ is called the critical manifold S. An example is shown in Figure 4. The lower branch S a is attracting (f v < 0) while the upper branch S r is repelling (f v > 0). The fold-curve L (f v = 0) corresponds to saddle-node bifurcations of equilibria in the layer problem. Therefore, equation (12) can be solved explicitly for one of the slow variables, e.g. r f = r f (r s , v), but not for the fast variable v.
The reduced problem is given by setting ε = 0 in (7)- (9): It describes the evolution of the slow variables (r f , r s ) on the critical manifold S.
To understand the reduced flow on the critical manifold we project the reduced Figure 4. The critical manifold S and the cross-section Σ 1 . Trajectories of the reduced problem include the black solid curves, which are the invariant manifolds of the saddle q, and the green curve, which is the strong stable manifold of n. The black dashed line is the fold-curve L. The magenta dashed line is the fast projection of the strong stable manifold of n from S r to S a . The blue curves are trajectories of the 3D SC model, and the red curve is a homoclinic orbit computed by AUTO, with ε = 0.01.
system onto the two-dimensional subspace (r s , v), since the critical manifold S is a graph r f = r f (r s , v) over that subspace. This projection is obtained by implicitly differentiating the algebraic constraint f = 0 which giveṡ with r f = r f (r s , v). Rescaling time by the factor −f v leads to the desingularized system which enables us to analyze the reduced flow in a neighborhood of the fold-curve f v = 0 where (16)-(17) is singular. Note that system (18)- (19) has the same phase portrait as system (16)-(17) but the orientation of the flow is reversed on S r (f v > 0) due to the rescaling of time.
The desingularized system (18)- (19) has two types of singularities: ordinary and folded. Ordinary singularities are defined by g = h = 0 and are therefore equilibria of the reduced system (16)-(17) as well. Folded singularities are defined by f v = 0 (the fold condition) and f r f g + f rs h = 0. They are classified as folded saddles, folded nodes, folded foci or folded saddle-nodes based on their classification as singularities of (18)- (19). Folded singularities are (generically) not equilibria of the  Figure 5. A sketch of the slow phase portrait with an illustration of φ + (γ 0,r ). The inset shows the intersection of φ + (γ 0,r ) and γ 1,a . This is the case shown in Figure 7.
reduced system (16)-(17) but they enable the reduced flow to cross the fold-curve in finite time (except folded foci). Such solutions within the reduced flow are called singular canards. All other points on the fold-curve L are called jump points since solutions of the reduced flow which reach the fold-curve blow-up in finite time and cease to exist. Solutions are then expected to follow close to the fast dynamics of the singularly perturbed system. A detailed bifurcation analysis of the singularities of the reduced flow can be found in [33]. A typical phase portrait of the reduced flow is shown in Figure 5 where n denotes a folded node, q ∈ S r a saddle and p ∈ S r a stable node. We denote by γ 0 the strong stable manifold of the folded node n, and define the pieces of γ 0 in S a and S r as γ 0,a = γ 0 ∩ S a , and γ 0,r = γ 0 ∩ S r . The region in S a that is bounded by γ 0,a and the fold-curve (which contains the weak eigendirection of n as well) is called the funnel. In the reduced system, all trajectories in the funnel converge to the folded node n in finite time where they cross the fold-curve L to reach the repelling branch S r of the critical manifold. Such trajectories are singular canards. Of particular interest are those singular canards which correspond to the funnel region between γ 0,a and the weak eigendirection (dashed line in Figure 5) which connects to the saddle q on S r . This sector within the funnel gets attracted to the node p on S r and is therefore (part of) the basin of attraction of the node p.
3.1. Summary of relevant canard theory. Fenichel theory [13,18] shows that those parts of the critical manifold S a and S r which are O(1) away from the foldcurve L perturb to invariant slow manifolds S a,ε and S r,ε of system (7)-(9) which are O(ε) close to their singular limit manifolds. Furthermore, the flow on these perturbed manifolds S a,ε and S r,ε is an O(ε) perturbation of the corresponding reduced flow on S a and S r . Since the node p is O(1) away from the fold-curve L it will persist as a node p ε of the flow on S r,ε . This requirement, p ε ∈ S r,ε , fixes the manifold S r,ε which is usually non-unique.
A maximal canard corresponds to a transverse intersection of the manifolds S a,ε and S r,ε extended beyond the fold-curve L. Since the equilibrium p ε is an attractor on S r,ε , a maximal canard will converge to p ε . 2 Wechselberger [32] analysed these extensions and showed that there exists generically [(1 − µ)/2µ] + 2 transverse intersections 3 of S a,ε and S r,ε near the folded node n. These (maximal) canards have distinct rotational properties near the fold-curve L. The primary weak canard corresponds to the weak eigendirection of n and serves as the axis of rotation for the other canards as well as for the invariant manifolds S a,ε and S r,ε . The primary strong canard γ 0 corresponds to the (unique) strong invariant manifold of n and makes one twist (180 • ) around the primary weak canard. For µ < 1/3, let γ k (k = 1, . . . , m, where m = [(1 − µ)/2µ]) denote the remaining secondary canards. These m secondary canards are positioned between the two primary canards (within the funnel region), and the k-th secondary canard makes (2k + 1) twists around the primary weak canard. Thus the primary strong canard and the m secondary canards are consecutively separated by a full rotation.
It was also shown in [4] that O(1) away from the fold-curve L, the secondary canards γ k are O(ε (1−µ)/2 ) close to the primary strong canard γ 0 . Consequently, there exists m + 1 sectors within the funnel bounded by the primary and secondary canards with distinctive rotational properties near the fold-curve L. Any trajectory within sector I k rotates k-times around the primary weak canard before it gets repelled. Note that sectors I k with larger k (larger rotation number) will be further from γ 0 . We refer to [32,4,16] for details and especially to [5] for a numerical study of the geometry of the invariant manifolds near a folded node singularity.
As ε → 0, each γ k converges to γ 0 . However, for ε > 0 sufficiently small, the flow in S a,ε and S r,ε is approximated by the reduced problem away from the fold. Therefore, we include the curves γ k in a phase portrait such as Figure 5, with the understanding that the figure does not capture the correct behavior near the folded node n since the rotations happen in a O( √ ε) neighborhood which vanishes in the singular limit ε → 0. These approximate trajectories γ k in the phase space of the reduced problem (16)- (17) should be viewed as the seeds of the true secondary canards of the full problem (7)- (9). This extra information is essential to understand the creation of all the oscillatory patterns (STOs and MMOs) in the SC model.

3.2.
Incorporating fast jumps and spikes. To explain STOs, MMOs and their bifurcations in the 3D SC model (as well as in the 7D SC model) we concatenate the flow of the reduced problem with the flow of the layer problem to obtain singular periodic orbits. These singular periodic orbits will persist as the observed STOs and MMOs in the full system based on the canard theory [32,4,33].
Any trajectory of the reduced flow which is not a canard will be concatenated with a fast fiber attached to the fold-curve L at the corresponding termination point of the reduced flow. This leads to fast jumps away from the fold towards increased v values.
Maximal canard trajectories γ k , will cross the fold-curve L and continue along S r towards the node p. Therefore, we concatenate γ k with the fast fiber attached to the node p. Here we have a choice of either jumping away towards increased v values or jumping back towards decreased v values and therefore back to S a . By continuation, each maximal canard γ k comes with a whole family of canards. In the singular limit ε → 0 they are all represented by the single trajectory γ k , but they can be distinguished by concatenating a family of fast fibers along γ k,r on S r . Again, there is the choice of either jumping away or jumping back along these fibers.
We focus first on the jump back mechanism associated with canards. Let φ − (x) be the base point on S a which corresponds to the (lower) fast fiber with base point x ∈ S r , i.e. the subset of a fast fiber with base point x ∈ S r along which solutions of the layer problem will converge to φ − (x) ∈ S a . Note, φ − is a projection map from S r to S a . The set φ − (γ k,r ) ⊂ S a is the set to which (one side of) the fast fibration of γ k,r projects. We call such a curve the umbra of γ k,r (see Bold, et al [2]). These are curves in S a that connect φ − (p) to n independent of the number k. In Figures 4, and 5, φ − (γ 0,r ) is shown as a dashed magenta curve.
When a solution of the 3D SC model makes a jump away in which v increases, the orbit leaves the region in which the 3D SC model provides a qualitative approximation to the full 7D SC model. Based on numerical calculations in the full 7D SC model, we formulated a list of qualitative properties of the corresponding global return map such that spiking behaviour is captured by the 3D SC model as well [33].
Definition 3.1. The global return map φ + : S r → S a has the following properties: 1. φ + preserves the orientation, in particular with respect to r s . 2. φ + contracts length scales, in particular in r s direction. 3. φ + shifts r s coordinate to a smaller value (shifts to the left). 4. The r s coordinate of φ + (n) is to the left of the r s coordinate of φ − (p).
In the 7D SC model, we generally found that the umbral curves φ + (γ k,r ) crossed at most one rotational sector boundary. Figure 5 shows a typical example, in which φ + (p) is in I 2 while φ + (n) is in I 1 . Furthermore, φ + (n) is to the left of φ − (p) as described in Definition 3.1. In the 7D SC model, as I app is decreased, we found that both φ + (p) and φ − (p) moved further into the funnel. The positions of φ + (p) and φ − (p) relative to the rotational sectors I k under variation of I app is incorporated in Figure 1.
The definitions of the global return maps φ − and φ + applied to the 3D SC model gives us sufficient information to explain the existence of chaotic invariant sets (section 4) and homoclinic clusters (section 5) in the full 7D SC model.

A return map and chaotic invariant sets.
In this section, we show that the return of canards to the funnel implies the existence of chaotic invariant sets. A (local) cross-section Σ 1 that is transverse to the flow is given by a plane r s = const near the folded node; see  Figure 6. Sketches of the set B (gray strip) and the sets B k (blue) near Γ k in the cross-section Σ 1 when φ − (p ε ) ∈ I 2 . The "curves" labeled Π(B k ) are actually exponentially thin strips.
canard must pass close through the folded node, which explains the choice of the section. By using the blow-up technique, Wechselberger [32] showed that manifolds S a,ε and S r,ε extended by the flow near a folded node intersect in Σ 1 as shown in Figure 6. The size of this rotational funnel is O( √ ε). Let Γ k denote the intersection points of S a,ε and S r,ε in Σ 1 . These points correspond to the maximal canards γ k that converge to the saddle p ε ∈ S r,ε . Γ 0 is the strong canard, Γ 1 the first secondary canard which makes one rotation around the primary weak canard, etc. The rotational sectors I k ⊂ S a,ε with rotation number k intersect Σ 1 in curves that are therefore bounded by Γ k−1 and Γ k . Recall that S r,ε contains the saddle p ε , so S r,ε is the boundary between the trajectories (exponentially close to S r,ε in Σ 1 ) which jump either forward or backward after remaining close to S r,ε . We denote by Φ + (p ε ) and Φ − (p ε ) the position of φ + (p ε ) and φ − (p ε ) within section Σ 1 . 4 For any choice of S a,ε , Φ − (p ε ) and Φ + (p ε ) are exponentially close to S a,ε . For the arguments we will make, no special choice of a stable manifold S a,ε is necessary. In Figures 6 and 7, we have sketched a S a,ε that contains Φ − (p ε ), which suggests that we chose S a,ε to be a slow manifold that contains the lower (i.e. the "jump back") one-dimensional unstable manifold of p ε .  Figure 7. A sketch of the intersections of S a,ε and S r,ε in the cross-section Σ 1 near the folded node. In this sketch, Φ − (p ε ) is in sector I 1 , between Γ 0 and Γ 1 , and Φ + (p ε ) is in sector I 2 . Each neighborhood N k of a canard point Γ k is divided into two open subsets A k (red) and B k (blue) by S r,ε . The return map Π of these sets is indicated by curves which are actually (exponentially small) strips. This is not essential, since any other choice of a slow manifold would be exponentially close to this one; any solution which returns to S a,ε an O(1) distance from the fold-curve L will be exponentially close to S a,ε in Σ 1 .
Let B be the set of points in Σ 1 sufficiently close to S r,ε so that trajectories starting in B have slow segments near S r,ε that are long enough to jump back to S a,ε to the left of Σ 1 . This is the gray shaded region in Figure 6, one boundary of which is S r,ε . Note that the "width" of this region B is exponentially small since trajectories have to follow the slow manifold S r,ε for an O(1) distance.
(B 0 and B 1 are sketched as the blue regions shown in Figure 6.) By our choice of {B k }, we can now define the return map Π to be the flow map from m k=0 B k to Σ 1 . Consider a set of points {x j } ∈ B k for which lim j→∞ x j ∈ S r,ε . By our choice of S r,ε , lim j→∞ Π(x j ) = Φ − (p ε ). The sets Π(B k ) are thin strips exponentially close to S a,ε , stretching from Φ − (p ε ) to somewhere in I 0 (past Γ 0 ) 5 .
Consider the case where φ − (p ε ) ∈ I 2 ; this is the case illustrated in Figure 6. Because the extent of each B k along S r,ε is O(ε), each strip Π(B k ), k ∈ {0, 1} necessarily covers B j , j ∈ {0, 1}. Thus the map Π has a horseshoe, and therefore an invariant set with dynamics topologically conjugate to a full shift on two symbols [29,30,17]. The transverse intersections of S a,ε with S r,ε , and the transverse intersection of the set φ − (γ k,r ) with the reduced flow in S a ensure that, for sufficiently small ε, the chaotic invariant sets are hyperbolic.
More generally, if φ − (p ε ) ∈ I k , k ≥ 2, then the map has chaotic invariant sets with dynamics topologically conjugate to a full shift on k symbols. Suppose there are m secondary canards. Because the widths of the sectors I 1 , . . . , I m are O(ε (1−µ)/2 ), the "typical" cases (for sufficiently small ε) are either φ − (p ε ) is outside the funnel, or φ − (p ε ) is in the innermost sector I m . The above argument shows that, in the latter case, there will be chaotic invariant sets with dynamics conjugate to a full shift on m symbols.
The chaotic sets considered so far involve only STOs through jump back canards. By also considering that φ + (p ε ) may be in the funnel, we can predict the existence of additional chaotic sets for the 7D SC model. Unlike φ − (γ k,r ), the umbral curve φ + (γ k,r ) generally does not extend from φ + (p ε ) to outside the funnel. It crosses at most one rotational sector boundary as shown in Figure 5.
Analogous to the definition of B, let the set A ⊂ Σ 1 be an exponentially small strip near S r,ε , with S r,ε as a border, but on the opposite side. Trajectories of this set A include jump-away canards. Here we choose the width of the exponentially small strip A such that all trajectories pass the (set of) turning points 6 in S r,ε where the fast projection map φ + (the global return) of this set is tangent to the slow flow on S a,ε . We denote by τ + k,ε the corresponding turning point of the canard γ k,r and Φ + (τ + k,ε ) denotes the corresponding point in Σ 1 which is exponentially close to S a,ε . Let A k = N k ∩ A, where N k is some O(ε) neighborhood of Γ k . Then Π(A k ) is a thin strip exponentially close to S a,ε , stretching from Φ + (p ε ) to Φ + (τ + k,ε ) and crossing near a canard point Γ j .
Consider, for example, a case where φ + (p ε ) ∈ I 2 , φ + (τ + k,ε ) ∈ I 1 for k ∈ {0, 1}, and φ − (p ε ) ∈ I 1 (see Figure 5). A sketch is shown in Figure 7 where the sets A k for k ∈ {0, 1} are shown in red. The map Π satisfies the following covering relations: Π(A 0 ) covers A 1 and B 1 ; Π(B 0 ) covers A 0 and B 0 ; Π(A 1 ) covers A 1 and B 1 ; and Π(B 1 ) covers A 0 and B 0 . Therefore the return map Π contains an invariant set on which the dynamics are topologically conjugate to a subshift of finite type whose transition diagram is shown in Figure 8(a).
Of course, more complicated cases are possible. We do not attempt to give a complete classification; instead we describe another (slightly more complicated) case to illustrate the idea. If Φ − (p ε ) is in sector I 2 , Φ + (p ε ) is in sector I 3 , and Φ + (τ + k,ε ) is in sector I 2 , then each A k (k = 0, 1, 2) covers A 2 and B 2 , and each B k (k = 0, 1, 2) covers A 0 , B 0 , A 1 and B 1 , and the return map contains an invariant set on which the dynamics are topologically conjugate to a subshift of finite type whose transition diagram is shown in Figure 8(b). Figure 8. Transition diagrams discussed in section 4. Note that some connections are bidirectional.

MARTIN WECHSELBERGER AND WARREN WECKESSER
5. Homoclinic clusters. An interesting feature of the bifurcation diagrams of the full 7D SC model (Figure 1) is the convergence of multiple families of periodic orbits to homoclinic orbits that occur within an extremely small parameter range. In this section we take advantage of the strong contraction onto S a,ε to reduce the map Π to a one-dimensional map, and we use this map to explain the mechanism that gives rise to these "homoclinic clusters". As in the previous section, we first limit the discussion to jump-back canards. Let L k = S a,ε ∩ B k , where B k was defined in Section 4. Then Π(L k ) is a curve exponentially close to S a,ε , and if {x j } is a sequence of points in L k with lim j→∞ x j = Γ k , then lim j→∞ Π(x j ) = Φ − (p); see e.g. Figure 6.
We define a one-dimensional map f (x, λ) on m k=0 L k × R as the composition of Π and a projection π from Σ 1 to the curve S a,ε ∩ Σ 1 . We abuse notation and treat Φ ± (p ε ) and {Γ k } as, say, arclength coordinates in the curve S a,ε ∩ Σ 1 , with Γ k < Γ k+1 . In this context, the projection of the curves Π(L k ) onto S a,ε is very simple, in particular orientation is preserved. The parameter λ is some system parameter which can be interpreted as a function of the main bifurcation parameter J app . 7 We have lim for 0 ≤ k ≤ m.
In the remainder of this section, we discuss specific clusters that we have observed in our computations. For each of these clusters, we construct a one-dimensional map f with the appropriate qualitative properties. Our analysis generally leads to the prediction of infinite families of homoclinic orbits. In our computations, we have observed only a small number (i.e. 2-4) homoclinic orbits in each cluster, but in all cases, they are consistent with the predictions of our qualitative arguments. λ), and the intervals exponentially close to Γ 0 and Γ 1 are where canards form, the graph of f has the form shown in Figure 9. Note that Figure 9 shows only those branches of f which with positive constants r and b, where r > (Γ 1 − Γ 0 ). 8 The reason for choosing this function is that near Γ 0 and Γ 1 , the behavior of the map will be dominated by the flow past the saddle p ε . The ratio of the eigenvalues of the saddle p ε is O(ε), so f is a plausible asymptotic form for the map. For this discussion, we assume which can be interpreted as a leading order term in an expansion in λ. With this assumption, there is a homoclinic orbit with pattern H − 1 when λ = 0. As λ is increased from 0, the homoclinic orbit is broken and Φ − (p ε , λ) increases. Therefore, an increase in λ reflects a decrease inJ app in the 3D SC model (see Section 3.2).
We have assumed r > (Γ 1 − Γ 0 ), so λ * is exponentially small. This shows that this set of homoclinic orbits occur at parameter values exponentially close to zero.
The homoclinic orbit of the SC model associated with a solution to the defining equation The above argument predicts that if there is a homoclinic orbit of type H − 1 , then for any k ≥ 0, a homoclinic orbit with the pattern C − 1 (C − 0 ) k H − 0 occurs at a parameter value exponentially close to that of H − 1 . This agrees with the numerical results shown in Figure 3. The pattern of the orbits that have convincingly converged to homoclinic orbits are The computation of one side of the 0 5 family failed to converge before reaching a very large period; however, this family also appears to be heading into the same cluster. In terms of the one-dimensional map f , this orbit has period 3, so this orbit is connected to the birth of the chaotic dynamics that exist when Φ − (p ε , λ) is in the second sector. The fate of this family as it approaches the cluster is still under investigation.
More complicated homoclinic patterns, in which the orbit makes multiple transitions between neighborhoods of Γ 0 and Γ 1 (e.g.
, may also be possible.

5.2.
The H − 0 cluster. The right branch of the 1 2 family ends at a homoclinic orbit with pattern C + 0 H − 1 , at a parameter value very close to the value at which the H − 0 homoclinic orbit occurs at the end of the 0 1 branch (see Figure 1). Here we show that the one-dimensional map predicts the closeness of these parameter values, and predicts the existence of other homoclinic orbits in this cluster. Figure 11. The map f , near the occurrence of the H − 0 homoclinic orbit. (a) An orbit of the one-dimensional map that corresponds to the orbit pattern C + 0 H − 1 in the 7D SC model; Φ − (p ε ) Γ 0 . This is, in fact, the pattern of the homoclinic orbit at the right end of the 1 2 family. It occurs at a value of I app very close to that of the For this argument, we include the "jump-away" canards. We extend the intervals on which the map is defined to be L k = S a,ε ∩ (A k ∪ B k ), where A k and B k were defined in section 4. This extends f to include small intervals to the left of and adjacent to Γ k . Note that, as shown in Figure 1, Φ + (p ε ) ∈ I 2 at the parameter value where this cluster occurs. We assume that the umbral curve of the jumpaway canard crosses γ 1,a . 9 A function with the correct qualitative properties is, for example, Because we assume that the jump-away umbral curve crosses γ 1,a , we assume here that s > Φ + (p ε , λ) − Γ 1 . Analogously to the previous section, we assume that Φ − (p ε , λ) = Γ 0 + λ so that the H − 0 homoclinic orbit occurs when λ = 0. We also assume that where Γ 1 < D < Γ 2 , so the dependence of Φ + (p ε , λ) on λ is not pathological. We look for additional homoclinic orbits when λ is small but nonzero. Figure 11 shows sketches of the hypothetical map f close to the occurrence of the H − 0 homoclinic orbit.
We first consider λ < 0, so Φ − (p ε , λ) < Γ 0 . The simplest candidate for a homoclinic orbit in this case is depicted in Figure 11(a). The defining equation is or, using the assumptions from above, Solving for λ gives which is exponentially small. The pattern of this homoclinic orbit is C + 0 H − 1 , which is the pattern of the homoclinic orbit that we computed at the end of the right branch of the 1 2 family.
The assumptions on f outlined above imply the existence of an unstable fixed point less than but close to Γ 1 , and the sequence {f −k 1 (Γ 1 , λ)} ∞ k=1 converges to this fixed point. An argument analogous to that of the previous section shows that for each k ≥ 0, the equation , λ), λ) = Γ 1 has a solution λ that is exponentially small. The corresponding homoclinic orbits have the pattern C + 0 C + 1 k H − 1 . Now consider λ > 0. Then Φ − (p ε , λ) > Γ 0 , and if λ is sufficiently small, the first iteration f (Φ − (p ε , λ), λ) is less than but close to Γ 0 . The first homoclinic orbit we expect to find as λ increases occurs when An example is sketched in Figure 11(b). Because of the unstable fixed point near Γ 1 , arguments similar to those given above lead to the prediction of orbits of the form C − 0 C + 0 (C + 1 ) k H − 1 for k ≥ 0.

5.3.
The H − s clusters, s ≥ 2. A similar analysis using one-dimensional maps as presented for the H − 0 and H − 1 clusters can be carried out for any H − s cluster, s > 1. The analysis of the one-dimensional maps predicts homoclinic orbits with patterns such as C − s (C − 0 ) k H − 0 and more complicated patterns involving arbitrary patterns of C − j before a final H − k , 0 ≤ j, k < s.

5.4.
The H + s clusters, 2 ≤ s ≤ 7. In Figure 1, we see that for 2 ≤ s ≤ 7, the left branch of the 1 s family and the right branch of the 1 s+1 family both converge to the H + s cluster. The homoclinic orbits of the two branches are H + s and C − s H + 0 , respectively. Here we use the one-dimensional map to explain this feature of the bifurcation diagram.
First we note from Figure 1 that for s ≥ 2, φ − (p ε ) is in the funnel when a homoclinic orbit of the form H + s occurs. 10 Thus Φ − (p ε ) > Γ 0 and the one-dimensional map has the form shown in Figure 12. Arguments similar to those given above show that exponentially close to the parameter value at which the H + s orbit occurs, there 10 This explains why no clusters are expected near H + 0 and H + 1 . exist homoclinic orbits of the form C − s (C − 0 ) k H + 0 ; the simplest of these is C − s H + 0 , shown in Figure 12(b). 11 6. Conclusions. Multiple time-scales govern the dynamics of the stellate cell model by Acker et al [1]. This fast/slow structure leads to a singularly perturbed system with a critical manifold S that is a folded surface, a characteristic feature of many physiological systems. The folded surface enables the system to switch between two states such as the subthreshold regime and the spiking regime.

MARTIN WECHSELBERGER AND WARREN WECKESSER
In [33], we showed that the folded node singularity n and the critical point p on the unstable sheet of the critical manifold explained the complex bifurcation diagram of mixed-mode oscillations and other periodic orbits of the 7D stellate cell model. In this paper, we have explored further consequences of the fast/slow structure of the phase space; specifically we have described how the return of canards to the funnel implies the existence of chaotic invariant sets, and we have used a onedimensional reduction of a Poincaré map to give a qualitative explanation of the homoclinic clusters in the bifurcation diagram.
The essential ingredient that leads to chaotic invariant sets is the return to the funnel of canards that passed through the folded node. In the models that we considered, the reduced flow in S r was controlled by the node p; if φ ± (p) is in the funnel, then at least parts of φ ± (γ k,r ) are sure to be in the funnel, and if these curve transversely intersect the trajectories γ k,a , then chaotic invariant sets are expected. The existence of a stable node in the reduced flow in S r is not a requirement for chaotic invariant sets, but in this model it provides control over the behavior of the canards in S r away from the fold.
Guckenheimer [14] studied a variant of the forced van der Pol oscillator which can be tuned such that the global return mechanism injects trajectories into the funnel of a folded node singularity. His numerical calculations suggest that the system has also chaotic invariant sets, therefore giving another example where canards create complex dynamics.
Doedel, Krauskopf and Osinga [8] studied a global bifurcation structure in the Lorenz system. They focused on explaining the codimension one homoclinic bifurcation also known as homoclinic explosion point by calculating global manifolds and their intersections. Their argument uses symbolic dynamics but is not based on singular perturbation analysis since the Lorenz system has no apparent multiple time-scale structure which makes their approach conceptually different to ours. The singular perturbation analysis of the stellate cell model has the advantage that we are able to relate the symbolic dynamics directly to a one-dimensional discontinuous return map. Although we describe this return map f only qualitatively in section 5, this map captures all the generic features of the stellate cell model and can therefore be viewed as its "canonical form". Numerical simulations of the Poincaré map Π (not shown) indeed verify our generic assumptions on the map f . We expect this work will motivate further theoretical study of the local and global bifurcations that occur in fast/slow systems with folded singularities.