Interaction of two systems with saddle-node bifurcations on invariant circles. I. Foundations and the mutualistic case

The saddle-node bifurcation on an invariant circle (SNIC) is one of the codimension-one routes to creation or destruction of a periodic orbit in a continuous-time dynamical system. It governs the transition from resting behaviour to periodic spiking in many class I neurons, for example. Here, as a first step towards theory of networks of such units the effect of weak coupling between two systems with a SNIC is analysed. Two crucial parameters of the coupling are identified, which we call \delta_1 and \delta_2. Global bifurcation diagrams are obtained here for the"mutualistic"case \delta_1 \delta_2>0. According to the parameter regime, there may coexist resting and periodic attractors, and there can be quasiperiodic attractors of torus or cantorus type, making the behaviour of even such a simple system quite non-trivial. In a second paper we will analyse the mixed case \delta_1 \delta_2<0 and summarise the conclusions of this study.


Introduction
The saddle-node on an invariant circle (SNIC) bifurcation is one of the basic scenarios for creation of a periodic orbit in smooth continuous-time dynamical systems (number three in the list of Andronov and Pontryagin [AP]). It goes under various other names, such as saddle-node on a limit cycle (SNLC), homoclinic to a saddle-node, or saddlenode loop (though we reserve the latter name for a codimension-2 case, e.g. Z points of [BGKM]).
The SNIC bifurcation is for example the scenario by which many class I neurons are believed to make the transition from resting behaviour to periodic spiking [EK, RE] (this was also proposed by RSM to the neurophysiologist H.B.Barlow in 1986). It underlies some regenerative excitable chemical systems. It occurs in mechanical systems too, like the damped pendulum with torque for sufficiently strong damping, or its Josephson junction analogue [LHM].
Mathematically, a SNIC for a C r (r ≥ 2) vector fieldẊ = V (X) on a manifold is an elementary saddle-node (an equilibrium with a simple eigenvalue 0, no other eigenvalues on the imaginary axis and non-zero quadratic coefficient in the null direction) with a trajectory that is asymptotic to the saddle-node in both directions of time ("homoclinic") along the null direction; it follows that the resulting invariant circle is C r . Thus, taking a coordinate x along the circle, with 0 representing the saddle-node equilibrium, and as x → 0 for some C > 0 and v(x) > 0 for x = 0. Without loss of generality, by scaling x we take C = 1. We denote the length of the circle (in the rescaled coordinate x) by L (an alternative would be to scale the length of the circle to 1 and keep a general value of C).
Consider C r perturbation of a SNIC by a parameter µ, with µ = 0 representing the unperturbed case. By the elementarity condition, the invariant circle is normally hyperbolic and so persists for all C r -small perturbations. A parameter-dependent extension of the coordinate x can be chosen so that the length of the circle remains L and the vector field on it is a C r -small perturbation of that for µ = 0. By a parameter-dependent shift of origin to remove the linear term (using the implicit function theorem), the perturbed vector field has the local form for some smooth functions a, b of µ, with a(0) = 0, b(0) = 1.
It is convenient for later purposes to make b(µ) precisely 1 for µ small, by a coordinate change X(x, µ) preserving the length L. This can be achieved as in Appendix A. We will suppose that a (0) = 0 and thus for small µ we can use a as parameter instead of µ, so without loss of generality we have (2)ẋ = µ + x 2 + o(x 2 ) as x → 0. For µ < 0 it has two equilibria: a sink and a source (if the normal directions to the circle are attracting, the case of most relevance, then in the full state space these are a sink and a saddle). For µ > 0 the circle is a periodic orbit whose period T (µ) is asymptotic to π/ √ µ as µ → 0, which spends all but a bounded amount of its period in any neighbourhood of 0. An explicit example of a family with a SNIC bifurcation iṡ φ = µ + 2(1 − cos φ) on a circle of length 2π. Another artificial-looking but useful example iṡ x = µ + x 2 on the real line union a point at infinity, interpreted as the projective real line (set of lines through the origin in a plane); the projective real line is diffeomorphic to a circle with the coordinate x representing the slope of the line ; one could write x = tan θ with θ considered modulo π (half a revolution brings to itself) and theṅ θ = µ cos 2 θ + sin 2 θ on a circle of length π, which is similar to the first example (put φ = 2θ), cf. [Iz]. The goal of this paper is to find what happens when two SNIC bifurcations are coupled. This is step one towards finding how a network of class I neurons or Josephson junctions behaves, as addressed for example in chapter 8 of [HI]. Although it sounds a simple problem, we have not found a rigorous treatment in the literature and our analysis and its solution are remarkably complicated. An example was treated numerically by [GK].
We treat the general case from the point of view of determining the minimal structure that the bifurcation diagrams must possess.
We will assume the uncoupled systems are C r with r large enough that every operation we will perform makes sense (r = 5 suffices). The coupling we treat is of the form of a C r -small perturbation to the vector field given by the product of two SNIC bifurcations. This could represent the effect of gap-junction coupling ("electrical synapse") between neurons (for an introduction to neuroscience see [NMW]). We will show in Section 3.2 that gap-junction coupling produces the mutually excitatory case. On the other hand, chemical synaptic coupling requires further analysis, because it may involve adding additional degrees of freedom to represent the effects of the neuro-transmitters and time delays to represent signal travel times. The first consideration does not make much difference, because near the bifurcation the timescale for the electric dynamics is longer than for relaxation of the channels in response to the neuro-transmitters so the latter can simply be added into the normally hyperbolic directions. The second makes the dynamics infinite-dimensional and although it is likely that again the effects can simply be added into the normally hyperbolic directions, we did not pursue this yet. Similarly, the coupling can model electrical coupling between Josephson junctions. For an introduction to the theory and experiments on Josephson junction arrays, see [M, U] respectively.
Our analysis uses heavily terminology and results from section 4 of [BGKM] (on bifurcations for flows on a 2-torus) and some from [BM] (on coupling of a saddle-node periodic orbit with an oscillator). In particular, we recall some key concepts right now. A Poincaré flow on T 2 is one with a global cross-section, i.e. a transverse section such that every forward and backward orbit crosses it. With respect to a choice of coordinates (x 1 , x 2 ) on the universal cover of the torus (i.e. consider T 2 = R 2 /(L 1 Z × L 2 Z), where L j are the lengths of the cycles in the coordinate directions), the homology direction of a forward orbit of a flow on T 2 is the limit of the unit vector in the direction of the vector V of (signed) numbers of revolutions in x 1 and x 2 as time goes to +∞ (or 0 if V does not go to infinity). The winding ratio is the homology direction modulo reflection through 0. For a Poincaré flow, every orbit has the same homology direction and it is non-zero. We denote Poincaré flows by P . A Cherry flow is one with a homotopically non-trivial transverse section Σ and a direction of time such that the orbits of a nonempty subset Σ return to Σ under the flow, the induced map g : Σ → Σ is continuous, and lim x→l g(x) = lim x→r g(x) for all gaps [l, r] (components of Σ\Σ ). Every unbounded orbit of a Cherry flow has the same non-zero homology direction. We denote Cherry flows by C (in [BGKM], C denotes a larger class). If the homology direction of a flow is that of an integer vector (p 1 , p 2 ) with no common factors, we say the flow is partially mode-locked of type (p 1 , p 2 ). A flow is fully mode-locked if every orbit has homology direction 0 (equivalently if every orbit is bounded on the universal cover). We denote fully mode-locked flows by F (or F M L). A non-contractible closed curve on the torus is called rotational.

Product system
Denote the coordinates of the two systems by x j , j = 1, 2, with lengths L j for one revolution, their parameters by µ j and their (uncoupled) vector fields bẏ It will be convenient to suppose that ∂v j ∂µ j ≥ c > 0 for all x j , not just for x j near 0 (say c = 1 2 ). This can be achieved by parameter-dependent coordinate changes X j (x j , µ j ) as in Appendix B.
The product of two circles is a 2-torus. We obtain the bifurcation diagram of Figure 1 for the product system in the plane of µ = (µ 1 , µ 2 ), with global phase portraits as indicated. In the positive quadrant, the flow is smoothly conjugate to that of a constant vector field on the unit torus which varies smoothly with parameters. For the conjugacy one can use the fractions τ j (x j ) of the period T j traversed from a reference point (say x j = 0). The constant vector field is (1/T 1 , 1/T 2 ), and it has asymptotic expression ( √ µ 1 , √ µ 2 )/π as µ → 0. In particular, the flows in the positive quadrant are Poincaré flows with homology direction (T 2 , T 1 )/ T 2 1 + T 2 2 which varies smoothly and at nonzero rate with the ratio µ 1 : µ 2 . In the negative quadrant, the flow is fully mode-locked, with four invariant circles (x j ≈ ± √ −µ j ) intersecting pairwise at four equilibria. We call phase portraits topologically equivalent to this, basic tartan. The boundaries of the negative quadrant constitute two simultaneous saddle-node bifurcations of equilibria (sne), which merge into a degenerate one at µ = (0, 0). In the bottom right quadrant, the flow is a Poincaré flow of type (1, 0), with repelling and attracting periodic orbits at x 2 ≈ ± √ −µ j . Similarly, in the top left quadrant, we have Poincaré flow of type (0, 1).
The boundaries of the positive quadrant (minus the vertex 0) correspond to (elementary) saddle-node periodic orbits (snp) (periodic orbits with a Floquet multiplier +1).
3. Effect of weak coupling 3.1. First steps. The invariant 2-torus of the uncoupled system is normally hyperbolic for µ small enough, so persists under small smooth perturbation [F], and the vector field on it is a small smooth perturbation of the uncoupled case, in general µ-dependent. We restrict attention from now on to a neighbourhood of µ = 0 where the above holds. Let us denote the perturbation size in C r , for some r ≥ 2, by δ. In particular this implies that the changes to v 1 and v 2 and to their first and second derivatives are at most δ (actually in Appendix E we will require a bound on the third derivative of the perturbation, which does not need to be as strong as δ but for convenience we will assume that too).
The perturbed system iṡ x 2 ≤ 0}, and resulting homotopically non-trivial repelling and attracting sets C ± 2 (drawn as C 1 invariant circles here but could take other forms as in Figure 3).
A ± 2 be the annuli as shown in Figure 2. Then, defining φ to be the flow of the above differential system, C − 2 := ∩ t>0 φ t A − 2 and C + 2 := ∩ t<0 φ t A + 2 are homotopically nontrivial attracting, respectively repelling sets.
If we neglect the perturbation, C ± 2 are just the circles x 2 ≈ ± √ −µ 2 given by the zeroes of v 2 . Under perturbation, we will find regions in which they persist to C 1 invariant circles (either periodic orbits or chains of connecting orbits between equilibria), regions in which they are C 0 invariant circles connecting equilibria but not necessarily C 1 , and 5 in Paper II regions in which they are not even C 0 circles. Figure 3 shows some possible forms for C ± 2 . Figure 3. Some possible forms for the maximal invariant sets in the annuli. Only (d) and (g) are C 1 . Case (f) is not a C 0 circle.
In particular, in the region where both µ j ≤ −2δ, the nullclines intersect in precisely four points. As the derivative ofṽ is close to diag (2x 1 , 2x 2 ), they are two saddles, a source and a sink, arranged just as for the unperturbed case. We leave out the detailed justification. The invariant manifolds of the saddles leave close to horizontal and vertical and because of the signs of the components ofṽ between the nullclines, they are obliged to fall into the source or sink in topologically the same way as for the unperturbed case. Thus when both µ j ≤ −2δ, we continue to obtain a basic tartan phase portrait at the C 0 level, though at the C 1 level it can take forms like those in Figure 4. In particular, C ± 2 , C ± 1 are all at least C 0 circles , but not necessarily C 1 because of the ways the saddle manifolds may meet at the source or sink. . Some possible realisations of basic tartan in C 1 . The third one has C ± 2 and C − 1 non-C 1 circles.
Next we analyse µ 1 ≥ −2δ, µ 2 ≥ −Cδ for some C greater then 2 and prove that C ± 2 are C 1 invariant circles there, using normal hyperbolicity theory. The normal linearised 6 dynamics for the unperturbed casė is hyperbolic (repelling for +, attracting for −). The character of the tangential linearised dynamicsδ x 1 = v 1 (x 1 ) δx 1 depends on the sign of µ 1 .
For µ 1 ≥ 0, although tangent orbits grow a lot for 0 < x 1 L 1 /2, the growth is all cancelled out by contraction for 0 < L 1 − x 1 L 1 /2, producing zero Lyapunov exponent. Thus taking Fenichel's approach [F] to normal hyperbolicity theory, timeaveraged tangential contraction or expansion rates are less than normal ones and so the circles persist to nearby C 1 invariant circles on adding C 1 small enough perturbation.
For µ 1 < 0, there are two equilibria x 1 ≈ ± √ −µ 1 on the invariant circles of the unperturbed system, with Lyapunov exponents approximately ±2 √ −µ 1 . All other orbits are heteroclinic to these so have forward Lyapunov exponent −2 √ −µ 1 and backward Lyapunov exponent +2 √ −µ 1 . Thus tangential contraction or expansion rates are weaker than the normal ones if µ 2 < µ 1 < 0 and so under this condition the circles persist to nearby C 1 invariant circles on adding C 1 small enough perturbation for parameters in this region.
To quantify what counts as small enough perturbation, however, is not so easy. In Appendix D we show there is C > 2 such that the C 1 invariant circles persist in at least the region µ 2 ≤ −Cδ, µ 1 ≥ µ 2 + (C − 2)δ, sketched in Figure 5. -2 x 2 ≈ ± √ −µ 2 C circles at 1 C circles at 1 x 1 ≈ ± √ −µ 1 Figure 5. There are C 1 rotational invariant circles in the shaded areas.
A C 1 invariant circle may either be a periodic orbit or a C 1 string of connections between equilibria. The case of periodic orbit happens for µ 1 > Cδ (with µ 2 < −Cδ) because thenẋ 1 > 0. The case of a C 1 string of connections between equilibria happens for µ 1 ≤ −2δ because then there are two equilibria in each annulus.
Outside the region where normal hyperbolicity theory applies, however, it could happen that the attracting or repelling sets have the forms of Figure 3(a), (b) respectively, or more complicated, e.g. Figure 3(e),(f). Transition from C 1 circle (Figure 3(d)) to other phase portraits of Figure 3 would proceed via use of the fast direction to the sink (respectively source) (Figure 3(c)) or via node to focus transition leading to Figure 3(e) and possible subsequent Hopf bifurcation leading to Figure 3(f). For an example of such transitions see [BM].
Similar remarks hold for µ 1 < −Cδ and vertical invariant circles (see Figure 5). For µ 1 + µ 2 > Cδ there is a global cross-section, e.g. x 1 L 1 + x 2 L 2 = 0, becauseṽ 1 +ṽ 2 > 0 everywhere in that region of parameter space. Thus we obtain Poincaré flows. The homology direction goes continuously from (1, 0) to (0, 1) from the lower right to the upper left, but generically locks to rational values. In fact it changes monotonically with a parameter λ along the lines µ 1 = K 2 − λ, µ 2 = K 2 + λ, K constant (K > Cδ), because for the unperturbed system the component of the derivative ∂v ∂λ in the direction v ⊥ = (−v 2 , v 1 )/|v| (using Euclidean norm |v| = v 2 1 + v 2 2 ) is at least c(v 1 + v 2 )/|v| (where c is as in Section 2), which is positive in this region, and small perturbation can not change its sign, so increasing λ turns the vector field in the positive (anticlockwise) direction. On the boundaries of the regions of Poincaré flow of rational type (p, q) there is a saddle-node periodic orbit.
For µ 2 < −Cδ there are precisely two curves of saddle-node equilibria, graphs over µ 2 , separating µ 1 < −Cδ from µ 1 + µ 2 > Cδ and each creating a pair of equilibria on one of the horizontal invariant circles. This is because the equations for a saddle-node equilibrium areṽ 1 = 0,ṽ 2 = 0 and det Dṽ = 0, which can be written approximately as µ 1 ≈ −x 2 1 , x 2 2 ≈ −µ 2 , 4x 1 x 2 ≈ 0; using the implicit function theorem, the second equation determines x 2 as either of two functions of (x 1 , µ 1 , µ 2 ) near x 2 = ± √ −µ 2 for µ 2 < −Cδ; substituting these into the third determines x 1 ≈ 0 as one of two functions of (µ 1 , µ 2 ); substituting for x 2 (x 1 , µ 1 , µ 2 ) and then x 1 (µ 1 , µ 2 ) into the first provides µ 1 ≈ 0 as either of two functions of µ 2 . Without further hypotheses it is not possible to say which saddle-node bifurcation happens first (indeed the curves could cross), but in between them the flow is Cherry flow of type (1, 0). Similarly, there are two curves of saddle-node equilibrium in µ 1 < −Cδ, graphs over µ 1 , in between which we have Cherry flow of type (0, 1). The phase portrait in the region µ j ≤ −Cδ, j = 1, 2, is a basic tartan as we already discussed. It is not guaranteed to remain like this in the whole of the full mode-locked region, however. For example, if the two sne curves in µ 2 ≤ −Cδ cross then it could easily happen that to the left of this there is a heteroclinic bifurcation D → A 01 (see Fig. 14 for the notation), which would give rise to a skewed tartan (two invariant horizontal circles and two invariant circles of type (1, 1)). Nevertheless, under generic hypotheses to be formulated at the end of Sec. 3.2 (δ 1 , δ 2 = 0), we will prove in Sec. 4.3 that the phase portrait is a basic tartan in all the part of the full mode-locked region with min j∈{1,2} |µ j | δ 2 (log 1/δ) 4 for δ 1 δ 2 > 0 (the case δ 1 δ 2 < 0 will be addressed in Paper II).  Figure 6, we divide the torus into the strips |x 1 | ≤ η, |x 2 | ≤ η for some η small, their intersection B and the complement (Figure 7).
The idea is that η should be small enough so we can accurately use second order Taylor expansion of v j about x j = 0 inside the strip |x j | ≤ η, yet considerably larger than √ Cδ so that the effects of µ = (µ 1 , µ 2 ) and δ outside the strips are relatively small, in particular so thatṽ j > 0 outside |x j | ≤ η and there are no equilibria outside the intersection box. We will start with η = kδ 1/3 for a small constant k, which is the largest that suffices for neglect of the cubic and higher terms compared with µ when µ is inside the triangle, but we will at later stages reduce or increase it, so we leave η explicitly in most formulae.
The quadratic terms q(x) ofṽ can be transformed to X 2 1 +ε 1 X 2 2 , X 2 2 +ε 2 X 2 1 by putting x = M X with M some near-identity matrix, depending on µ, determined by eliminating the X 1 X 2 terms in the quadratic part M −1 q(x) ofẊ and making its "diagonal" coefficients be 1. The idea is that the derivative of the map from the matrix elements  M 22 to the coefficients of the X 2 1 and X 1 X 2 terms inẊ 1 and the X 1 X 2 and X 2 2 terms inẊ 2 is near the matrix diag(1, 2, 2, 1) which is invertible. Thus M comes out within δ of the identity (since δ bounds the size of the perturbation to the quadratic terms inẋ) and the coefficients ε j above are also of order δ. It is not strictly necessary to have reduced the quadratic terms in this way, since the quadratic part of the perturbation would have turned out to be negligible anyway, but it is tidier to eliminate terms when one can rather than have to bound their effects.
The above reduction of the quadratic terms can be equally well achieved with a global coordinate change on the torus by using instead with k j = 2π/L j , which is within order δ of the identity. Conjugating by this changes the vector field by order δ, which is the same size as the perturbation, so from now on we suppose such a change of coordinates to have been made. Now eliminate the diagonal linear terms of v by a shift of origin. This is straightforward by completing the squares and produces a shift of order δ. Finally the effect of all the above on the constant terms is to make constant terms (μ 1 ,μ 2 ) that are a near-identity diffeomorphism from (µ 1 , µ 2 ).
Thus we are left withẊ 2 + δ 2 X 1 + ε 2 X 2 1 + HOT 2 for some coefficients δ j and ε j of order δ, where HOT j denotes higher order terms of the form . We make the assumption that δ j = 0, in fact of size similar to δ (though for some purposes, δ j significantly larger than δ 4/3 would suffice). Then the ε j terms are negligible relative to the δ j terms. The tidy way to deal with this is to push the ε j terms into the HOT j , so we shall consider this done.
The signs of δ j will play a crucial role. We say the coupling coefficient δ j is excitatory if δ j > 0, inhibitory if δ j < 0, by loose analogy with neuroscience. If both δ j have the same sign we say the coupling is mutualistic (mutually excitatory if both are positive, mutually inhibitory if both are negative). If the δ j have opposite signs we say the coupling is mixed.
Gap junction coupling gives the mutually excitatory case δ 1 , δ 2 > 0, because its effect is to add current I = V 1 −V 2 R from neuron 1 to neuron 2, where V j are their electric potentials and R is the resistance of the junction (taking a linear model). This adds perturbation where C j is the capacitance of neuron j. Passage through the near resting state of a neuron corresponds to increasing V , so the angle coordinate x j for a SNIC neuron is oriented the same way as V j near the saddlenode. Up to a shift of origin of V j , we have x j ∼ V j K j for a positive scale factor K j to make the quadratic coefficient equal 1. So we read off that By reversing time and X 1 if necessary, we can always take δ 1 > 0, but one must remember to reverse time at the end of the analysis, which interchanges attractors and repellors for example (we will do this in Section 5.2 for example).
We shall switch notation back fromμ j and X j to µ j and x j , but it should be remembered that these are related to the original parameters and coordinates by a near-identity diffeomorphism.
3.3. Reduced system. In the triangle in parameter space ( Figure 6) and the box B (|x j | ≤ η) in state space we study the approximate vector fieldv Although the neglected higher order terms are small compared to µ and x 2 , they are not necessarily small compared to δ 1 x 2 and δ 2 x 1 , so one might ask why we retain the latter. The idea is that all results inside the sub-box |x| √ δ will be accurate (because the higher order terms are dominated by the δ 1 x 2 and δ 2 x 1 terms there) and the results obtained outside this sub-box will turn out insensitive to the higher order terms anyway, because the dominant ones are of the form α j x 3 j inẋ j (rather than general cubics in both variables).
Its equilibria form the graph of a function from state space to parameter space: The type of the equilibrium is determined by the determinant and trace of the derivative Dv of the vector field: In particular, the equilibrium is a saddle for det < 0, a sink for det > 0, tr < 0, and a source for det > 0, tr > 0. The two branches of hyperbola in (x 1 , x 2 ) where det = 0 correspond to saddle-node equilibria, and the line where tr = 0 to neutral saddle or Hopf bifurcation (according as det < 0, > 0). Figure 8 illustrates the mutualistic case δ 1 δ 2 > 0. Note that tr = 0, det > 0 is impossible in this case so the tr = 0 curve is all neutral saddle. . Curves where det = 0 or tr = 0 on the manifold of equilibria, considered as a graph over (x 1 , x 2 ) for the mutualistic case δ 1 δ 2 > 0.
The approximations for these curves are good for |x| √ δ. In fact they are good for all |x| ≤ η 1. To see this, the dominant correction to the equation for a saddle-node equilibrium is . Because of the factor x 1 x 2 on the right, the correction has relative size O(η) and the saddle-node curves perturb to 4x 1 x 2 = δ 1 δ 2 (1 + O(η)). The tr = 0 curve deforms to approximately x 1 + x 2 = − 3 2 (α 1 x 2 1 + α 2 x 2 2 ), which is a shift of at most O(η 2 ). The projection of the manifold (5) of equilibria to parameter space has fold curves where the determinant is zero. To find their images in parameter space it is convenient to parametrise the two branches of hyperbola. Here we make a separation of the analysis into the mutualistic and mixed cases.
In this paper we study the mutualistic case. We will treat the mixed case in Paper II.

The mutualistic case
The mutualistic case is δ 1 δ 2 > 0. By reversing time and the orientation of x 1 and x 2 if necessary, we can take δ 1 and δ 2 > 0. 4.1. Analysis of equilibria. In the mutualistic case, we parametrise the saddle-node curves by with σ = +1 for the positive branch of the hyperbola of Figure 8 and −1 for the negative branch. It follows from (5) that the saddle-node curves project to which are drawn in Figure 9. . Curves of saddle-node equilibrium (full) and of neutral saddle (dashed) in parameter space for δ 1 δ 2 > 0 (drawn for δ 1 = 0.5, δ 2 = 0.3), also indicating the number of equilibria in the regions they separate.
Note that dµ j /dθ, j = 1, 2, are non-zero along the σ = −1 curve, which we call the outer saddle-node curve, µ moving from lower right to upper left as θ increases. They have a common zero along the σ = +1 curve, however, at θ = θ c where e 3θc = δ 1 /δ 2 , which causes it to have a cusp at so we call it the cusped saddle-node curve. The position of the corresponding degenerate equilibrium in state space then is 2 . In the case δ 1 , δ 2 > 0 we are considering here, it is topologically a saddle but with zero exponent in the contracting direction. It can be useful to rewrite the position of the cusp as The approximations for the saddle-node curves are good for all |µ| ≤ Cδ, because the effect of the higher order terms is negligible anyway in |µ| δ, and for the pieces of saddle-node curve in −Cδ ≤ µ 1 ≤ −Kδ 2 , K large positive, where µ 1 ≈ −x 2 1 , µ 2 ≈ −δ 2 x 1 , we obtain The analogous result holds for the pieces in −Cδ ≤ µ 2 ≤ −Kδ 2 . In particular, the cusp in the saddle-node curve is stable to small perturbation of a two-parameter family, so survives. Readers versed in singularity theory will recognise the saddle-node curves in Figure 9 as a slice through the unfolding of a hyperbolic umbilic singularity [GG]. Indeed, define a mapping ϕ from (x 1 , x 2 , δ 1 , δ 2 ) to (µ 1 , µ 2 , δ 1 , δ 2 ) by (5). Then its singularities (points where the rank of Dϕ is less than 4) correspond to saddle-node equilibria and the origin to a hyperbolic umbilic point. Because ϕ is a stable mapping, all small smooth enough perturbations of the 4-parameter family of vector fields (4) have set of equilibria and singularity set smoothly equivalent to that of (4). Any smooth 2-parameter family of vector fields near the case δ 1 = δ 2 = 0 has bifurcations of the set of equilibria given by a slice through the unfolding of the hyperbolic umbilic.
Also shown in Figure 9 is the projection of the curve of neutral saddles to parameter space (Hopf bifurcation does not occur in the case δ 1 δ 2 > 0), which is easily computed to be a parabola (µ 1 − µ 2 ) 2 + (δ 1 + δ 2 )(δ 2 µ 1 + δ 1 µ 2 ) = 0. This is accurate, however, only for |µ| δ; the effects of higher order terms can shift the curve by order |µ| 3/2 , which becomes the same size as the distance (δ 1 + δ 2 ) |µ| between the two sides of the parabola when |µ| approaches Cδ, so even allowing the possibility that the two sides cross.
It is useful to calculate the position of the two other equilibria for parameter values on the cusped saddle-node curve. By factoring out the known double root, they are found to be at 2 e θ/2 for ε = ±1, with ε = +1 giving a saddle, −1 a sink.
Furthermore, one can calculate the eigenvalues and eigenvectors for parameters on the cusped saddle-node curve (and on the outer saddle-node curve, so we give both cases).
14 The second eigenvalue of the saddle-node is tr = 2σ √ δ 1 δ 2 cosh θ which is positive on the cusped curve, and the null eigenvector has slope −σ δ 2 /δ 1 e θ (by the slope of a vector v = (v 1 , v 2 ) we mean v 2 /v 1 ). A useful trick to save work is to note that the derivative Dv is symmetric with respect to the inner product ξ, ζ = ξ 1 ζ 1 /δ 1 + ξ 2 ζ 2 /δ 2 , so the eigenvectors are perpendicular in this inner product. This implies that the product of their slopes is −δ 2 /δ 1 . So for example, the slope of the second eigenvector of the saddle-node is σ δ 2 /δ 1 e −θ .
The eigenvalues λ and slopes s of the eigenvectors of any equilibrium are given by the following expressions: For the saddle which coexists along the cusped saddle-node curve, substitute the following expressions into (10) and (11): The ways the equilibria connect within the box B are deduced by studying the nullclines (e.g. Figure 10 for the region inside the cusped saddle-node of equilibria curve), and by computing the signs of the slopes of the eigenvectors at the equilibria from (10) and (11). In particular, for δ 1 δ 2 > 0, the source and sink are always nodes (not foci) and for δ 1 > 0, δ 2 > 0 their fast and slow directions are as indicated on Figure 11. Global connections will be analysed in subsections 4.3 and 4.4.
Thus the phase portraits in the box B for the various parameter regimes, in particular as the parameters move along the cusped saddle-node curve, are as indicated on Figure 11. 4.2. Transit map. Next, to aid in understanding the global dynamics, we consider the transit map from x 2 = η to x 2 = L 2 − η, restricting attention to those trajectories that remain within the strip |x 1 | ≤ η. This analysis is independent of the sign of δ 2 , so logically could be done in the previous section, but it would have interrupted the presentation.
By the choice of η √ Cδ and the restriction to the triangle in parameter space where in particular |µ 2 | ≤ Cδ, we haveṽ 2 > 0 for x 2 ∈ [η, L 2 − η]. The transit map x 1 → x 1 is given by integrating  Figure 10. Nullclines for parameters inside the cusped saddle-node curve (δ 1 , δ 2 > 0). sne sne cusp Figure 11. Phase portraits in the box for the indicated parameter regimes (δ 1 , δ 2 > 0). The winding ratio of the Poincaré and Cherry flows varies continuously between the extremes illustrated. Shading near a saddle-node indicates its repelling half-plane. from x 2 = η to x 2 = L 2 − η. As a first approximation, take v 1 = x 2 1 in |x 1 | ≤ η. To this level of approximation the result is that where t 2 is the time taken, which may depend on x 1 , but is dominated by the ends of the trajectory. At the lower endṽ 2 ≈ x 2 2 (µ 2 is negligible here because η ≤ x 2 ≤ L 2 − η), and at the upper endṽ 2 ≈ (L 2 − x 2 ) 2 and so t 2 ≈ 2 η . As an explicit example, if (which ∼ x 2 2 for x 2 near 0 and ∼ (L − x 2 ) 2 for x 2 near L 2 ) then t 2 = 2π L 2 cot πη L 2 ∼ 2 η for η L 2 . Thus inserting t 2 ≈ 2/η into (12) the transit map is approximately (14) x The effects on the transit map of corrections to v 1 and v 2 at µ 1 = µ 2 = 0 are shown in Appendix E to be O(x 2 1 log 1/η) and the effects of parameters µ j and the perturbation δ for |µ j | ≤ Cδ are shown to be O(δ/η).
For the approximate transit map (14), It follows that the transit map moves all points to the right by at least Kδ/η (some K > 0) except when |x 1 | = O(δ 1/2 ). Within |x 1 | = O(δ 1/2 ) it moves points possibly left or right but by at most K δ/η, some K . The analogous result holds for the transit map from x 1 = η to x 1 = L 1 − η in the strip |x 2 | ≤ η.

Extension of global dynamics into parts of the triangle.
We can now describe the global dynamics in parts of the triangle in parameter space (Figure 13 Firstly, there is K large such that for µ 1 ≤ −Kδ 4/3 between the outer sne and lower cusped sne curves, the two vertical circles C ± 1 persist (though not necessarily C 1 ). A sketch of the proof is given in Appendix F.
The invariant circles are periodic orbits above the outer saddle-node curve; the attracting one gains a saddle and sink on crossing the outer saddle-node of equilibrium curve and the repelling one gains a saddle and source on crossing the cusped saddlenode curve. Below the cusped saddle-node curve, the study of the nullclines at the end of Subsection 3.3 shows that the equilibria connect in the box as in Figure 11.
Similarly, we obtain the analogous results in a strip along the horizontal boundary of the triangle. In the intersection of the union of these boundary strips with the region inside the cusped saddle-node curve, the dynamics is fully mode-locked.
Everywhere above the outer sne curve the dynamics consists of Poincaré flows, because there are no equilibria and there is a transverse section (e.g. take x 1 = L 1 2 where µ 1 ≥ µ 2 or x 2 = L 2 2 where µ 2 > µ 1 ). The boundaries of the partial mode-locking strips are saddle-nodes of periodic orbits (snp). In particular there is a curve of snp crossing the outer sne curve and delimiting the lower right zone of type (1, 0) Poincaré flow, and an analogous one for type (0, 1). By the analysis at the beginning of this subsection, the type (0, 1) snp curve leaves the outer sne curve with |µ 1 | ≤ Kδ 4/3 ; but the most we can say for µ 2 > −µ 1 is that it remains in µ 1 ≥ −Cδ. Similarly for the type (1, 0) snp curve.
It might be possible to obtain greater control over these snp curves as follows. Consider type (0, 1), and take µ 2 sufficiently positive. The Lyapunov exponent of a periodic orbit in 2D is λ = div v dt. In the regime considered, the orbit is a graph x 1 (x 2 ), so this integral can be transformed to λ = L 2 0 div v v 2 dx 2 . The change in x 1 for an orbit segment making one revolution in x 2 is ∆x 1 = The conditions for an snp are λ = 0, ∆x 1 = 0. In the uncoupled case this happens at x 1 = 0, µ 1 = 0. The derivative of λ with respect to initial condition x 1 on x 2 = 0 is 2T 2 (µ 2 ) which is large (T 2 (µ 2 ) = L 2 0 dx 2 v 2 is the period), so under perturbation λ = 0 determines a nearby x 1 . The derivative of ∆x 1 with respect to µ 1 is T 2 , so ∆x 1 = 0 determines a nearby µ 1 . Closer analysis, however, is required to determine more precise bounds on the snp curve. Secondly, we are now able to give the phase portraits for the Cherry flows in between the outer and cusped sne curves near µ 1 = −Cδ and µ 2 = −Cδ. They are sketched in Figure 13 for δ 1 , δ 2 > 0, by analysis of the locations of the equilibria (cf. Figure 11).
Lastly, we prove that, writing ν j = −µ j , the phase portrait is a basic tartan in the part of the region bounded by the cusped sne curve with both ν j δ 2 (log 1/δ) 4 . To begin, in the region bounded by the cusped saddle-node curve, we name the saddles s and t and the branches of the local invariant manifolds of the saddles A, B, C, D as shown in Figure 14(a). On the part of the boundary going to the left from the cusp (θ > θ c ), s becomes a saddle-node, with A being in the null direction; on the part of the boundary going to the bottom (θ < θ c ), t becomes a saddle-node, with C being in the null direction. At the cusp, B merges with D. Away from the cusp (where there is ambiguity), we may continue to use the labels for the branches of invariant manifold from the remaining saddle in the region between the two saddle-node curves (viz. C, D to the left, and A, B to the bottom). We denote the sink (which is in the negative quadrant) by p. On the universal cover of the 2-torus we shall refer to the translation of p by a vector (mL 1 , nL 2 ) with m, n integer, as p mn , and similarly for s, t, A, B, C, D.
For ν 2 δ 2 (log 1/δ) 4 , the branches A and D will be shown to be close to vertical for |x 2 | ≤ η = (log 1/δ) −1 (note the new choice of η), the transit map from x 2 = η to L 2 − η has a similarly small effect, but D starts significantly to the left of A, so that D reaches To prove the near-verticality of A and D, note that they can be defined by starting from the appropriate equilibrium point which we denote in generality by (x e 1 , x e 2 ), and integrating in the appropriate direction of x 2 (increasing for D, decreasing for A). Let us consider the case of D. Then x 1 considered as a function of x 2 (which is true for η not too large) is a fixed point of the map T defined by (16) T The map T is a contraction in supremum norm on Lipschitz functions from [x e 2 , η] to R satisfying x 1 (x e 2 ) = x e 1 with Lipschitz constant K, say, if η is sufficiently less than L 2 /2, so achieved for δ small enough. Call its contraction constant λ < 1; it can be made as small as we want by suitable choice of K and δ. To bound the distance of the fixed point of T from the constant function x 1 (x 2 ) = x e 1 it suffices to estimate the distance of the image of the constant function from the constant function and divide it by (1 − λ).
Now v 1 (x e 1 , x 2 ) and v 2 (x e 1 , x 2 ) have a common factor (x 2 − x e 2 ), and for x 2 small the ratio is approximately δ 1 /(x 2 + x e 2 ). This approximation is not particularly accurate for larger x 2 but is good enough. Then (17) x 2 It follows that with k a little larger than (1 − λ) −1 . A similar bound is obtained for A, with x e 2 replaced by |x e 2 |. The effect of the transit of D from x 2 = η to L 2 − η is of order δ/η as in subsection 4.2. Now to complete the analysis, the horizontal distance ∆x 1 between the two equilibria for given ν 2 is at least that for the worst case of ν 1 , namely on the lower branch of the cusped sne curve, because as ν 1 increases from this, theẋ 1 nullcline rises, thereby separating the two equilibria further. For this worst case, and taking ν 2 significantly larger than δ 2 , the formulae for the equilibria give Thus (21) ∆x 1 ∼ 2δ 1 ν 1/4 2 . On the sne curve, we also have the approximation for both equilibria (and this does not change much if µ 1 is moved to the left). So the shifts in D and A on reaching x 2 = ±η respectively are from (18) at most approximately Choosing η = (log 1/δ) −1 (for which the value of the left hand side is of the same order as its minimum, yet η → 0 as δ → 0), we obtain the sufficient condition (24) ν 2 δ 2 (log 1/δ) 4 , for D to pass to the left of A 01 . Similarly, B passes under C 10 if ν 1 δ 2 (log 1/δ) 4 . Hence we obtain the global connections D → p 01 and B → p 10 in the regions indicated in Fig. 14(b).

Heteroclinic connections and consequences.
To determine the structure of the bifurcation diagram in the remaining central region of Fig. 13, we consider first the special case of systems which are symmetric with respect to simultaneous interchange of x 1 with x 2 and µ 1 with µ 2 , thus in particular δ 1 = δ 2 and L 1 = L 2 . Then by symmetry at the cusp we have Cherry flow of type (1, 1), but for ν 2 δ 2 log 1 δ 4 , D passes to the left of A 01 as we just showed, so a curve of heteroclinic connection D → A 01 must occur, as indicated in Fig. 15(a), separating the cusp from the region ν 2 δ 2 log 1 δ 4 . The curve of heteroclinic bifurcation leaves the lower branch of cusped sne curve tangentially; this is a bifurcation not seen in [BGKM, BM] so we give it a new name, E point, and analyse its generic unfolding in Appendix G. Where the curve of heteroclinic bifurcation hits the leftward branch of the cusped sne curve, we obtain the phase portrait of Figure 15(b). This is called an S point in the terminology of [BGKM], and it produces a "saddlenode fan" of Cherry flows. The range of winding ratios produced in the saddle-node fan depends on the copy of the sink p to which B connects there. If B → p 1n then the saddle-node fan generates all winding ratios from (0, 1) to (1, n + 1). See Figure 15(a) for the case n = 1. By reflection symmetry there is an analogous curve of heteroclinic connection B → A 10 , generating a saddle-node fan of all winding ratios from (1, 0) to (m + 1, 1) if D → p m1 . If there is no other bifurcation along the lower cusped sne curve between the D → A 01 heteroclinic and this saddle-node fan, then n = 1, and we shall suppose this case for the next five paragraphs. The curves of heteroclinic bifurcation D → A 01 and B → C 10 cross on the diagonal at a T point (in the terminology of [BGKM]). Two ways this might happen are shown in Fig. 16 (others with more crossings can be envisaged). We suspect only case (a) occurs but despite much effort to prove that the curves of heteroclinic connection are roughly vertical and horizontal we did not succeed. Nevertheless, let us concentrate attention on case (a), and we will argue at the end that case (b) can not occur.
The phase portrait is shown in Figure 17(a), and its unfolding produces a fan of curves of heteroclinic and rotational homoclinic bifurcation and a partial mode-locking tongue, as shown in Fig. 17(b). Note that in this symmetric case, x 1 + x 2 = δ 1 for the two saddles, so both are repelling; this makes the exponents of the saddles less than 1 (for a saddle with eigenvalues −λ < 0 < µ, the exponent α = λ/µ), so the bifurcation diagram is the time-reverse of the first case of Fig. 4.23(c) of [BGKM]. Now if we suppose the bifurcation diagram near the T point extends to the cusped sne curve, then we obtain a sequence of saddle-node fans from each of the indicated curves of heteroclinic connection (of which the first is the one we already found). They accumulate at what we christen an F point as shown in Fig. 18; it is an analogue of the M point of [BGKM] but with its saddle-node of periodic orbits replaced by a rotational homoclinic connection. The resulting bifurcation diagram is sketched in Figure 19. Extending further in parameter space, where the resulting partial mode-locked tongues cross the curve of neutral saddle (K points of [BGKM]), their boundaries are replaced by curves of saddle-node periodic orbits, with a curve of rotational homoclinic connection between the K point and a Z point (saddle-node loop) on the outer sne curve. We recall from [BGKM] in Figure 20, how K and Z points can produce a horn of coexistence of attracting periodic orbit with attracting equilibrium. Calculation of the generic unfolding s t 0,1 s s 0,1 s s 1,2 s t n ,n + 1 s t 1,2 s t 2, 3 s s 1 ,1 T neutral saddle snp types (m+2n, 2m+3 n) (1,2) (0,1) (2,3) (1,1) Figure 18. A sequence of S points accumulating onto an F point.  Figure 19. Simplest bifurcation diagram around the cusp for δ 1 , δ 2 > 0. The lines of heteroclinic connections that emanate from the T point cross the sne curve at saddle-node fan points (S points) and continue as lines of homoclinic bifurcations, forming boundaries of regions of Cherry flow of type indicated. The dashed lines correspond to homoclinic bifurcations that use the other branches of the manifolds of the remaining saddle.
of K points was given by Dumortier, Roussarie and Sotomayor [DRS] and of Z points by Schechter [Sch] The resulting bifurcation diagram is sketched in Fig. 21. (1,2) (0,1) (2,3) (1,1) Figure 21. Extending further in parameter space for δ 1 , δ 2 > 0, where the resulting partial mode-locked tongues cross the curve of neutral saddle (K points), their boundaries are replaced by curves of saddle-node periodic orbits, with a curve of rotational homoclinic connection between the K point and a Z point (saddle-node loop) on the outer sne curve.
It might be that some of the curves of heteroclinic connection cross the curve of neutral saddle before reaching the cusped sne curve, which would make a slightly different picture.
For families that are not symmetric, but close to symmetric, the above conclusions continue to hold. Simply the T point is not necessarily on the diagonal. Further away from symmetry, however, various changes can occur. In particular, the T point could cross the neutral saddle curve of Figure 9, changing the exponent of one saddle relative to 1, producing other cases of unfolding of the T point (see [BGKM]). Or it could happen that the first bifurcation on sliding up the lower cusped sne is the saddle-node fan generated by B → C 10 bifurcation, with a D → A 01 connection happening higher up. This would make the first saddle-node fan cover all winding ratios from (1, 0) to (1, 1) and necessitate other changes, but eventually there would be a T point (with different homotopy type of heteroclinic cycle), generating a picture similar to the above (or one of its variants depending on the exponents of the saddles).
To complete this section, we conjecture that case (b) of Figure 16 can not occur. Our reason is that the E point at the end of the D → A 01 curve produces flow of type (1, 1) just outside the cusped region, whereas the (1, 1) tongue produced by the T point is contained in the quadrant containing the cusp point, and it seems unlikely to us to see non-monotonicity of the winding ratio on turning round the cusp.

Consequences for attractors
For applications the most important feature is attractors. We can deduce what bifurcations happen to attractors in the mutually excitatory case by studying the bifurcation diagrams and phase portraits of Section 4. As mentioned, gap junction coupling gives the mutually excitatory case. To deduce those for the mutually inhibitory case, we must time-reverse the flows and the orientations of x 1 and x 2 . 5.1. Mutually excitatory case. Corresponding to the bifurcation diagram of Figure 21 for δ 1 , δ 2 > 0, we obtain the types of attractors indicated in Figure 22.
In the whole of the region to the South-West of the outer sne curve there is an attracting equilibrium. The region outside the outer sne curve is divided into strips where there are one (or more) attracting homotopically non-trivial periodic orbits. This implies a periodic spiking pattern of the pair of oscillators. There is one strip for each ratio of firing frequencies between 0 and ∞. In between the strips are curves on which the whole torus is an attracting quasiperiodic flow. Each rational strip extends across the outer sne curve in two horns terminating on the neutral saddle curve (except just one horn for the (0, 1) and (1, 0) strips). Thus in these horns two attractors coexist: an equilibrium and a homotopically non-trivial periodic orbit. This produces hysteresis effects. Deformations of Figure 21 as mentioned towards the end of Subsection 4.4 will produce deformation of Figure 22 but no qualitative change in this description of the attractors.

5.2.
Mutually inhibitory case. Time-reversing the results of Section 4, we obtain Figure 23 for the attractors corresponding to the bifurcation diagram of Figure 21 in the case δ 1 , δ 2 < 0. The situation is simpler than the mutually excitatory one. Inside the cusped sne curve there is an attracting equilibrium and it attracts almost everything.  Figure 22. Attractors in the mutually excitatory case. There is an attracting equilibrium everywhere inside the outer saddle-node curve. In the indicated strips there is an attracting periodic orbit of the indicated homotopy type. The strips extend across the outer sne curve into two horns where attracting equilibrium and periodic orbit coexist. Outside the outer sne curve the strips of periodic attractors are interleaved by curves on which the whole torus has quasi periodic flow of given homotopy type.
There is no region of coexistence of attractors, except that there may be more than one periodic orbit of given type.
The diagram may change significantly if Figure 21 deforms in some of the ways described towards the end of Subsection 4.4. Which rational gets the cusp and which are strips rather than tongues depends at least on the type of T point. But the overall result is roughly the same.

Conclusion
We have analysed the generic interaction of two dynamical systems with saddle-node on an invariant circle (SNIC) bifurcation. We identified two key coupling parameters δ 1 and δ 2 . We studied the "mutualistic" case δ 1 δ 2 > 0 in detail. We found that even the simplest sub-cases of this have fairly complicated bifurcation diagrams. We found ns sne sne (  Figure 23. Attractors in the mutually inhibitory case. Inside the cusped saddle-node curve the only attractor is an equilibrium. The region outside the cusped sne curve is divided into strips where there is an attracting homotopically non-trivial periodic orbit, interleaved by curves with a quasi-periodic attractor.
that the attractors can be equilibrium, homotopically non-trivial periodic orbit, quasiperiodic torus or quasiperiodic cantorus and that in the mutually excitatory case, which is expected to be the relevant case for gap junction coupling, there are parameter "horns" where attracting equilibrium and periodic orbit coexist. The results are expected to be a useful guide to the study of the dynamics of networks of class I neurons, Josephson junctions and other situations where SNIC bifurcations occur.
Results for the mixed case δ 1 δ 2 < 0, which turns out to be even more complicated, will be presented in a separate paper. Let X(x, µ) = x − α(µ) sin kx, with k = 2π L and |α| < 1 k for invertibility. Theṅ Now we choose α(µ) to make the coefficient of X 2 equal 1. This is a quadratic equation for α with a simple root near 0 because b is near 1. Thus we obtainẊ = A(µ) + X 2 + o(X 2 ) for some function A, with A(0) = 0.
Appendix B: Coordinate change to make ∂v ∂µ ≥ c We prove here that for a SNICẋ = v(x, µ) with form (2) and for c ∈ (0, 1), by a parameter-dependent coordinate change X(x, µ) preserving the length L of the circle, we can make ∂v ∂µ ≥ c for all x for µ small. The vector field in the new coordinate X is given byẊ = X x (x, µ)v(x, µ) (subscript denoting partial derivative), with x(X, µ) determined by X = X(x, µ).
Using x µ = −X µ /X x , we obtain Take X(x, 0) = x at µ = 0, and for all small µ keep X(x, µ) = x in a small interval |x| ≤ δ so as to preserve (2). First make ∂Ẋ ∂µ ≥ c for some c > c at µ = 0. At µ = 0 Choose a smooth function g(x) ≥ c , and g(x) = 1 in |x| ≤ δ. Let To make sure that ∂Ẋ ∂µ ≥ c for small µ, we make the additional constraint on g that L 0 g−vµ v 2 dx = 0.
Using the implicit function theorem, for eachx 1 , x 2 (x 1 ) extends to a C 1 solution of (25) for x 1 in a neighbourhood ofx 1 , thus, by uniqueness the function x 2 (x 1 ) is C 1 . Applying the chain rule to (25), The case x 2 < 0 is similar.
Appendix D: A region with C 1 invariant circles As the tangential dynamics expands a lot for 0 < x 1 L 1 /2 and then contracts a lot for 0 < L 1 − x 1 L 1 /2, the first thing we do is introduce a new horizontal coordinate in which the expansion and contraction are small.
Letv be a positive C 1 function of x 1 , close to v 1 , and write ∆v = v 1 −v. Let y 1 = x 1 0 dx v(x) be the new horizontal coordinate. It has the interpretation of time from x 1 = 0 using the vector fieldv. From the horizontal dynamicsẋ 1 = v 1 (x 1 ) + O(δ) we obtainẏ To evaluate horizontal expansion or contraction in y 1 , The vertical dynamicsẋ 2 = v 2 (x 2 ) + O(δ) has linearisation Combining (26) and (27), the slope s = δx 2 δy 1 of a tangent vector evolves by the Ricatti equatioṅ Each O(δ) term in equations (26-28) is bounded by δ, because they represent sizes of components of the perturbation or its first partial derivatives.
Let us specialise to the annulus A + 2 where v 2 > 0. To deduce that C + 2 is a C 1 circle, it suffices to find conditions on (µ 1 , µ 2 , δ) such that a cone |s| ≥ s 0 is forward invariant and vectors in it are expanded exponentially. The value s 0 = |µ 2 | will work. We show the expansion property first.
To make the cone |s| ≥ s 0 forward invariant we have to choosev. First treat the case µ 1 ≥ |µ 2 |. Then we can takev = v 1 . Since ∆v = 0 in this case, and without loss of generality considering s > 0, from equation (28) s outside a larger neighbourhood of 0,v v is bounded independently of µ 1 . Thus on s = s 0 Using δ ≤ |µ 2 |/C and µ 1 ≥ |µ 2 |, we see that at s = s 0 , Taking C a little larger than 1 + (1 + B/2) 2 , this is positive, indeed by increasing C we can make it as close as we like to 2|µ 2 |.
Next we do the case µ 1 < |µ 2 |, where we will find a constraint of the form µ 1 greater than some function of µ 2 . This time, we take ∆v to be the (negative) constant µ 1 − |µ 2 |. Then ∆v = 0 andv ≥ |µ 2 |. So (for s > 0) ; again it is bounded uniformly outside a larger neighbourhood of 0. So at s = |µ 2 |, This is non-negative if and only if This region is sketched in Figure 24. Hence the cone |s| ≥ s 0 is forward invariant for those µ 1 in (−|µ 2 |, |µ 2 |) satisfying (29). The region (29) in which we have obtained a C 1 circle has too complicated a formula to carry around with us, but it contains a region of the form as shown in Figure 24, where C * is determined by the intersection of (29) with µ 1 = −2δ.

Appendix E: Corrections to the transit map
To estimate the effects on the transit map of corrections to v 1 and v 2 and of the parameters µ j and the perturbation δ, we start from the explicit choice (13) of v 2 (x 2 ) = L 2 π 2 sin 2 πx 2 L 2 , and take v 1 = µ 1 + x 2 1 in |x 1 | ≤ η. Then Since t 2 ∼ 2/η and using √ µ 1 η and x 1 (0) ∈ [−η, η 3 ], we see that the effect of µ 1 on The dominant change to t 2 with µ 2 comes from the Taylor expansions of tan −1 x ∼ x− x 3 3 , tanh −1 x ∼ x + x 3 3 rather than the L 2 π + πµ 2 L 2 terms. So for instance for positive µ 2 η 2 , The effect of a change ∆t 2 in t 2 on x 1 , using the case µ 1 = 0 and x 1 ∈ [−η, η/3], is So the effect of the change − 2 3 µ 2 Next we consider the change induced by deforming v 2 to an arbitrary C 3 vector field with the same second order Taylor expansion and positive away from zero. We can write 33 such a deformation as α(x 2 ) sin 3 πx 2 L 2 for a bounded function |α(x 2 )| ≤ α 0 . Then t 2 is given by The leading correction to t 2 is Then using (30) the correction ∆x 1 = O(x 2 1 log 1 η ). The effect of a perturbation δ onẋ 1 is at most like changing µ 1 by δ. By our previous analysis this produces ∆x 1 = O δ η . The effect of a perturbation δ onẋ 2 can be split into two parts. Firstly, there is a part bounded by sin 3 πx 2 L 2 which can be absorbed into α so, assuming the C 3 norm of the perturbation is O(δ), it changes α by O(δ) and the resulting change to ∆x 1 is O(δx 2 1 log 1 η ). Secondly, there is a part localised near the box which fits with expansion (3), thus δ 2 x 1 + ε 2 x 2 1 + O(x 3 2 ) + O(δ|x| 3 ) , but whose O(x 3 2 ) terms can be absorbed into α, so this has size at most δ 2 |x 1 | + δη 3 . They have an effect like µ 2 so from (30) they produce ∆x 1 = O(|x 1 | 2 (δ 2 |x 1 | η 3 + δ)). Lastly the O(x 3 1 ) terms inẋ 1 produce ∆x 1 = O( x 3 1 η ), because x 1 does not change much during the transit.
The effects of the different types of change are close to additive, so putting all this together, and using |µ j | ≤ Cδ we obtain that the error in the transit map for µ 1 = µ 2 = δ = 0 is O(x 2 1 log 1 η ) and the effects of µ j and δ on the transit map are O(δ/η). Using the implicit function theorem to determine t 2 as a function of x 1 (0), one can also deduce that the corrections are small in C 1 .
Appendix F: Persistence of invariant circles for µ 1 ≤ −Kδ 4/3 In this appendix we sketch the proof that there is K large such that for µ 1 ≤ −Kδ 4/3 between the outer sne and lower cusped sne curves, the two vertical circles C ± 1 persist. Let us treat the region inside the cusp first. If δη |µ 1 |, the flow is inward, respectively outward, across the boundaries of the strips |µ 1 | − δη − ≤ ∓x 1 ≤ |µ 1 | + δη + in the box B (see Figure 25), where is a small amount to take care of the remainder terms inẋ 1 . The equilibria all lie in these strips. The upwards unstable manifold D of the top left saddle and the downward stable manifold A of the lower right saddle stay in the corresponding strips until they exit the box (the notation for branches of saddle manifold is introduced in generality in Figure 14). The transit map moves D by O(δ/η) D A (0,0) B Figure 25. Box B with strips |µ 1 | − δη − ≤ ∓x 1 ≤ |µ 1 | + δη + for parameters inside the cusp with |µ 1 | δη, µ 1 < 0.
so if this is much less than |µ 1 | then D arrives at x 2 = L 2 − η to the left of A. Because of the way the equilibria connect in the box ( Figure 11) this implies that D goes to the sink and A to the source, thus forming our two vertical invariant circles. The two conditions δη |µ 1 |, δ η |µ 1 | can be achieved for the greatest region of |µ 1 | if we choose η ∼ δ 1/3 (a choice we have already found useful) and then the conditions hold for |µ 1 | ≥ Kδ 4/3 for large enough K as claimed.
Next we treat the region between the outer sne and the upper cusped sne curves. The lefthand strip contains a saddle and a sink and the saddle manifold D reaches x 2 = L 2 −η to the left of − |µ 1 | − δη − + O(δ/η). We now need to bound where the rightward branch C of stable manifold of the saddle goes. We claim that for |µ 1 | δ 2 it crosses x 1 = 0 inside the box and hence exits x 2 = −η with x 1 > 0 as in Figure 26. . Box B with strip − |µ 1 | + δη + ≤ x 1 ≤ − |µ 1 | − δη − and the backwards orbit of (2 |µ 1 |, η), for parameters between the outer sne curve and the upper cusped sne curve with |µ 1 | δη, µ 1 < 0.
Thus for δ/η |µ 1 |, D arrives to the left of C and hence goes into the sink, making an invariant circle C − 1 . To complete the analysis of this case we construct a periodic orbit for C + 1 . The interval J on x 2 = +η between D and 2 |µ 1 | under the backward flow exits x 2 = −η by an interval J between C and something a little to the right of + |µ 1 |. This is becauseẋ 1 = µ 1 + x 2 1 in backwards time contracts x 1 > |µ 1 | towards |µ 1 | roughly like e √ |µ 1 |t . To estimate the time t it takes x 2 to flow backwards from +η to −η, useẋ 2 ≈ µ 2 + x 2 2 + δ 2 x 1 and note that |µ 2 | δ 2 |µ 1 | between the sne curves and x 1 ≤ 2 |µ 1 |, so |ẋ 2 | ≤ 3δ 2 |µ 1 | + x 2 2 . The solution ofẋ = a 2 + x 2 is x(t) = a tan at so takes time ∼ π 2a to cross |x| a. Hence x 2 takes at least time 1 3δ √ |µ 1 | . This time greatly exceeds 1/ |µ 1 |. The backwards transit map from x 2 = −η to x 2 = −L 2 + η moves points of J by O(δ/η) so if this is less than |µ 1 | the interval J is mapped by the backwards flow strictly inside itself and so has at least one attracting fixed point, making a periodic orbit. With some estimates of its derivative we could establish uniqueness.

Appendix G: Unfolding of an E point
An E point is the codimension-two situation with an elementary saddle-node equilibrium whose strong unstable manifold connects to the stable manifold of a saddle whose other branch of stable manifold lies in the repelling half-plane of the saddle-node, forming a homotopically non-trivial cycle (or the time-reverse of this situation). See Fig. 27(a).
It can be unfolded by a parameter µ to unfold the saddle-node equilibrium and a parameter ε to displace the intersection of the stable manifold of the saddle with a transverse section Σ. See Fig. 27(b). The resulting bifurcation diagram is shown in Fig. 27(c).