Elif Köksal Ersöz

Mathieu Desroches

Martin Krupa

Synchronization of weakly coupled canard oscillators

Keywords: canards, phase response curves, slow-fast systems, synchronization, weak coupling 2010 MSC: 34C15, 34C26, 34D06, 34E17, 92C20

come

Introduction

Synchronization is a research topic of its own, which has produced a large body of knowledge, in particular for so-called weakly coupled oscillators [START_REF] Ermentrout | Frequency plateaus in a chain of weakly coupled oscillators[END_REF][START_REF] Van Vreeswijk | When inhibition not excitation synchronizes neural firing[END_REF][START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF][START_REF] Hoppenteadt | Weakly Connected Neuroal Networks[END_REF][START_REF]Phase Response Curves in Neuroscience: Theory, Experiment, and Analysis[END_REF][START_REF] Ashwin | Mathematical framework for oscillatory network dynamics in neuroscience[END_REF]. A classical object of interest in this context is the (infinitesimal) phase response curve or (i)PRC which encodes how a small perturbation affects the phase of an oscillator when applied all along the associated stable limit cycle solution. The derivation of the PRC relies on the linearization of the system along the unperturbed (i.e. uncoupled) cycle and corresponds to the solution of the adjoint variational equation. Solutions to the adjoint problem and PRCs give insights on the synchronization properties of coupled oscillating systems [START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF][START_REF]Phase Response Curves in Neuroscience: Theory, Experiment, and Analysis[END_REF] when the coupling strength is small enough. Such studies are gathered under the name "weakly coupled oscillator theory" [START_REF] Hoppenteadt | Weakly Connected Neuroal Networks[END_REF]. This theory has been linked with earlier studies from Malkin [START_REF] Malkin | Methods of Poincare and Liapunov in Theory of Non-Linear Oscillations[END_REF][START_REF] Malkin | Some Problems in Nonlinear Oscillation Theory[END_REF] by Izhikevich and Hoppensteadt in [START_REF] Hoppenteadt | Weakly Connected Neuroal Networks[END_REF]; an explicit proof was given in [START_REF] Buicǎ | Periodic solutions of nonlinear periodic differential systems with a small parameter[END_REF] based on the work of Roseau [START_REF] Roseau | Vibrations nonlinéaires et théorie de la stabilité[END_REF][START_REF] Françoise | Oscillations en biologie: Analyse qualitative et modèles[END_REF].

Weakly coupled oscillator theory has been used in many studies, especially to investigate the effects of slowly-varying parameters, underlying bifurcations and coupling strengths on collective dynamics. In one of the pioneering papers on this topic [START_REF] Van Vreeswijk | When inhibition not excitation synchronizes neural firing[END_REF], out-of-phase (OP) synchronization (intermediate modes between in-phase (IP) and anti-phase (AP) solutions) was shown to emerge from a pitchfork bifurcation in the phase difference as a function of the coupling parameter. A similar bifurcation structure has been found in type-I spiking neuron models, see e.g. [START_REF] Crook | Spike frequency affects the synchronization properties of networks of cortical oscillators[END_REF][START_REF] Chow | Dynamics of spiking neurons with electrical coupling[END_REF][START_REF] Lewis | Dynamics of spiking neurons connected by both inhibitory and electrical coupling[END_REF]. Another recent study related to type-I membranes [START_REF] Urban | Formation of antiwaves in gap-junction-coupled chains of neurons[END_REF] focused on the transition from IP to OP synchronous states in chains of Wang-Buszaki models coupled by gap junctions. This transition was investigated both analytically and numerically as a function of intrinsic system properties by using phase models and interaction function. In the framework of type-II neuron models, the impact of the Hopf bifurcation on the possible synchronization patterns has been studied, e.g., in [START_REF] Feingold | Phase locking, period doubling, and chaotic phenomena in externally driven excitable systems[END_REF][START_REF] Aronson | Amplitude response of coupled oscillators[END_REF][START_REF] Ermentrout | Mechanisms of phase-locking and frequency control in pairs of coupled neural oscillators[END_REF][START_REF] Ermentrout | Type I membranes, phase resetting curves and synchrony[END_REF][START_REF] Hansel | Synchrony in excitatory neural networks[END_REF]. Furthermore, variations of the PRC across a Hopf bifurcation were analyzed in cortical excitatory neuron models in [START_REF] Ermentrout | The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators[END_REF]. Qualitative changes in the behaviour of the PRC were also looked at directly from experimental data in [START_REF] Kiss | Predicting mutual entrainment of oscillators with experiment-based phase models[END_REF] and [START_REF] Wickramasinghe | Spatially organized dynamical states in chemical oscillator networks: Synchronization, dynamical differentiation, and chimera patterns[END_REF] where the interaction functions were analysed during the transition from Hopf and relaxation oscillators. The existence of different synchronization modes and of bistable regions in weakly coupled slow-fast systems interacting via gap junctions has been underlined in [START_REF] Belair | On linearly coupled relaxation oscillators[END_REF][START_REF] Bem | Short duty cycle destabilizes a half-center oscillator, but gap junctions can restabilize the anti-phase pattern[END_REF][START_REF] Chakraboty | The transitions from phase locking to drift in a system of two weakly coupled van der Pol oscillators[END_REF][START_REF] Cymbalyuk | In-phase and antiphase selfoscillations in a model of two electrically coupled pacemakers[END_REF][START_REF] Sherman | Rhythmogenic effects of weak electrotonic coupling in neuronal models[END_REF][START_REF] Chow | Dynamics of spiking neurons with electrical coupling[END_REF][START_REF] Kopell | Anti-phase solutions in relaxation oscillators coupled through excitatory interactions[END_REF][START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF][START_REF] Pfeuty | The combined effects of inhibitory and electrical synapses in synchrony[END_REF][START_REF] Mancilla | Synchronization of electrically coupled pairs of inhibitory interneurons in neocortex[END_REF]. Synchronization has also been studied in the context of coupled piecewise linear models, in particular in [START_REF] Coombes | Phase-locking in networks of pulse-coupled McKean relaxation oscillators[END_REF][START_REF] Coombes | Neuronal networks with gap junctions: A study of piecewise linear planar neuron models[END_REF].

Slow-fast oscillators are an important source of complicated dynamics, and particularly in relation to the canard phenomenon [START_REF] Benoît | Chasse au canard[END_REF][START_REF] Krupa | Relaxation oscillation and canard explosion[END_REF]. The term canard cycle refers to a class of periodic solutions of slow-fast systems which follow for a long time interval a repelling slow manifold. Canards occur in slow-fast systems near regions of the critical manifold (fast nullsurface) where the key assumption of normal hyperbolicity fails. The most common points of this kind in systems with one slow and one fast variables are generically fold points, so-called canard points. A canard solution flows from an attracting slow manifold to a repelling slow-manifold by passing close to such a canard point. In planar systems, canard cycles exist in a very narrow range of bifurcation parameters, an interval that is exponentially small in the time scale separation parameter ε (0 < ε 1). These sharp transitions upon parameter variation through the canard regime are called canard explosions [START_REF] Brøns | Bifurcations and instabilities in the Greitzer model for compressor system surge[END_REF]. Combinations of advanced theoretical techniques, such as blow-up methods [START_REF] Krupa | Local analysis near a folded saddle-node singularity[END_REF], and numerical methods [START_REF] Guckenheimer | Numerical computation of canards[END_REF] have introduced a new understanding of canard-induced complex oscillations in systems with multiple time scales (in R n , n ≥ 3), in particular mixed-mode oscillations (MMOs) [START_REF] Desroches | Mixed-mode oscillations in a multiple time scale[END_REF] and bursting oscillations [START_REF] Desroches | Mixed-mode bursting oscillations: Dynamics created by a slow passage through spike-adding canard explosion in a square-wave burster[END_REF][START_REF] Vo | Multiple geometric viewpoints of mixed mode dynamics associated with pseudo-plateau bursting[END_REF], and extended their applications to neuroscience [START_REF] Krupa | Mixed-mode oscillations in a three time-scale model for the dopaminergic neuron[END_REF][START_REF] Rotstein | The dynamic structure underlying subthreshold oscillatory activity and the onset of spikes in a model of medial entorhinal cortex stellate cells[END_REF][START_REF] Rubin | Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model[END_REF]. In (weakly) coupled slow-fast systems, the effect of canard solutions has been considered in several aspects such as the formation of clusters, synchrony, phase and amplitude dynamics [START_REF] Drover | Analysis of a canard mechanism by which excitatory synaptic coupling can synchronize neurons at low firing frequencies[END_REF][START_REF] Ermentrout | Canards, clusters and synchronization in a weakly coupled interneuron model[END_REF][START_REF] Roberts | Averaging, folded singularities, and torus canards: Explaining transitions between bursting and spiking in a coupled neuron model[END_REF][START_REF] Rotstein | A canard mechanism for localization in systems of globally coupled oscillators[END_REF][START_REF] Rotstein | Localized and asynchronous patterns via canards in coupled calcium oscillators[END_REF]. Recently, canard-mediated variability has been investigated in coupled phantom bursting systems addressing issues on synchronization and desynchronization [START_REF] Ersöz | Canard-mediated (de)synchronization in coupled phantom bursters[END_REF].

In this work we extend previous results on adjoint solutions and weakly cou-pled slow-fast oscillators to the case of canard cycles. Analytical formulations of adjoints and interaction functions were studied in [START_REF] Brown | On the phase reduction and response dynamics of neural oscillator populations[END_REF], which also provides a review of the behavior near bifurcation points. In the framework of relaxation cycles, an expression for the adjoints could be obtained in [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF] by taking the singular limit approximation, considering the attracting branches of the critical manifold in place of the slow segments of relaxation cycles, and instantaneous jumps in place of the fast flow. The consequence of using this setup is that the canard regime has not been dealt with. In the present study, we propose an alternative numerical strategy, based on numerical continuation, for the computation of solutions to the adjoint variational equation associated with planar slow-fast systems along a canard explosion.

In parameter space, canards organize the transition between the Hopf regime and the relaxation regime. Therefore we may expect to link the synchronous behavior between these two families [START_REF] Kiss | Predicting mutual entrainment of oscillators with experiment-based phase models[END_REF][START_REF] Wickramasinghe | Spatially organized dynamical states in chemical oscillator networks: Synchronization, dynamical differentiation, and chimera patterns[END_REF] by computing adjoints for canard solutions. When performing such computations, we observe a qualitative change in the sign and shape of the adjoint (or equivalently, of the iPRC) near the maximal canard (the cycle with the longest repelling segment). This phenomenon occurs in both canard-explosive systems that we consider here, namely the van der Pol (VDP) oscillator and a two-dimensional (2D) reduction of the Hodgkin-Huxley (HH) model. We propose an explanation to this qualitative change through the period function of the canard family which has a non-monotonic behavior across the explosion from the Hopf bifurcation point to the relaxation regime, namely, it increases during the headless-canard regime and it decreases during the canard-with-head regime. This remarkable property of canard explosions singles out the maximal (period) canard, for which we highlight a key role in synchronization, which to the best of our knowledge was not reported in previous studies. Similar dependence of the frequency upon a bifurcation parameter has been studied in [START_REF] Curtu | Mechanisms for frequency control in neural competition models[END_REF] in the context of "escape-release" mechanisms of central patterns generators. The authors of [START_REF] Zhang | Phase response properties of half-center oscillators[END_REF] have then linked this dependence to transitions in PRCs and phase-locking properties occurring in the low-frequency region.

In the second half of this work, we explore the dependence of the phase difference between the two weakly coupled identical VDP systems on system parameters. By investigating the effect of the main parameter displaying the canard explosion, we observe that the transition in synchronization properties occurring at the maximal canard of the coupled system manifests itself as the AP synchronization state changing its stability through a pitchfork bifurcation in the phase difference. Furthermore, we reveal the presence of 2nT -periodic synchronous states in the maximal canard neighborhood appearing via multiple period-doubling (PD) bifurcations. Finally, we consider the effect of the coupling strength on the synchronous states of the oscillators in the maximal canard regime. We give numerical evidence of the presence of PD cascades not predicted by the theory of weakly coupled oscillators (which is valid for moderate coupling strengths in various systems [START_REF] Lewis | Dynamics of spiking neurons connected by both inhibitory and electrical coupling[END_REF][START_REF] Netoff | Beyond two-cell networks: experimental measurement of neuronal responses to multiple synaptic inputs[END_REF]) but that can be justified using phase plane analysis of the single canard oscillators under scrutiny. We also propose in an analytical formula to compute adjoints associated with limit cycles of slow-fast systems in Liénard systems, which gives satisfactory yet improvable results. This paper is organized as follows. In Section 2, we introduce the main objects required to compute adjoint solutions along a limit cycle and we present our numerical strategy to do so along families of canard cycles. In Section 3, we analyze numerically the effect of the main parameter controlling the canard explosion on the synchronous states of the coupled VDP system and report a qualitative change occurring near the maximal canard solution. We then explain this change by invoking the properties of the period function associated with such a canard-explosive branch of limit cycles. In Section 4, we focus on the effect of the coupling strength parameter on the synchronous structure in the coupled VDP system near a maximal canard. After concluding and proposing a few perspectives to this work, we present an analytical formula to compute adjoint solutions for the type of systems we have investigated and test this formula numerically in Appendix A.

Computations of adjoint solutions along a family of cycles

Phase Response Curve and Adjoints

PRCs describe the phase shifts along a stable limit cycle of a dynamical system in response to a stimulus. Weakly coupled oscillator theory [START_REF] Ermentrout | Frequency plateaus in a chain of weakly coupled oscillators[END_REF][START_REF] Kuramoto | Chemical oscillators, waves, and turbulence[END_REF][START_REF] Hoppenteadt | Weakly Connected Neuroal Networks[END_REF] is used to predict the phase-locking properties of coupled oscillating system with a "small enough" coupling strength. This theory, which reduces the dynamics of oscillators to a phase variable, implies that coupling has small effects that can accumulate over time and lead to phase-locking behaviors. IPRCs correspond to PRCs in the limit of infinitesimal stimulus. One way to compute iPRCs is by means of non-trivial solutions to the adjoint variational equation associated to the stable limit cycle under consideration; there are numerous other approaches, see e.g. [START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF][START_REF]Phase Response Curves in Neuroscience: Theory, Experiment, and Analysis[END_REF].

We now recall the main elements necessary to introduce adjoint solutions associated with a limit cycles. Consider a dynamical system in

R n dX dt = F (X) (1)
that possesses a T -periodic asymptotically stable limit cycle γ. A phase variable φ ∈ [0, T ) is defined along the limit cycle γ parameterized by time and it is typically normalized to 1 or to 2π. It can be associated with points on the cycle by writing φ = Θ(x) for x ∈ γ. Then, perturbing a point x on the limit cycle with corresponding phase φ = Θ(x) (which we can also write as x = X(φ)) by a small quantity y ∈ R n leads to a delay or an advance of the phase. The new phase φ is given by

φ = φ + ∇ X Θ(x).y + O(||y|| 2 )
and the difference between the old and new phases for small perturbations are expressed as φφ = ∇ X Θ(x).y.

The vector function Z defined by Z(φ) = ∇ X Θ(X(φ))) is the gradient of the phase map describing how infinitesimal perturbations on any system variable along the limit cycle changes its phase. The function Z (which depends on φ or equivalently on t ∈ [0, T ]) is the solution of the adjoint variational equation

dZ(t) dt + A(t) T Z(t) = 0 (2)
which satisfies the normalization condition

Z(t) dX 0 (t) dt = 1, (3)
where

A(t) = D X F (X)| γ
is the linearization of system (1) around the limit cycle γ. The adjoint equation should be integrated backwards in time to eliminate all the transient components except the periodic one, which gives the solution. An algorithm to compute solutions to adjoint equations, based on backward integration, is embedded in software package xppaut [START_REF] Ermentrout | Simulating, analyzing, and animating dynamical systems: a guide to XPPAUT for researchers and students[END_REF], or can be coded in matlab [START_REF] Izhikevich | Slowly coupled oscillators: phase dynamics and synchronization[END_REF], in addition to a continuation-based approach in MatCont [START_REF] Govaerts | Computation of the phase response curve: A direct numerical approach[END_REF].

Canard explosions occur in slow-fast systems in a very narrow parameter range which is exponentially small in the timescale separation parameter 0 < ε 1. Naturally, this parameter range gets narrower as ε tends to 0, and limits the usage of classical tools to compute family of canard orbits and their adjoints. This limit has been acknowledged by Govaerts and Sautois who have introduced a direct numerical approach in the continuation package MatCont [START_REF] Govaerts | Computation of the phase response curve: A direct numerical approach[END_REF]. In addition to existing methods, we propose an alternative and simpler continuation-based strategy using the software package auto [START_REF] Doedel | auto-07p : Continuation and bifurcation software for ordinary differential equations[END_REF]. We formulate a periodic continuation problem which allows us to compute rapidly and reliably a family of cycles with associated non-trivial periodic solution of the adjoint equation. Note that a boundary-value problem (BVP) approach has been proposed in [START_REF] Ermentrout | Phase response curves to measure ion channel effects on neurons[END_REF], outside of a continuation setup given that the system was with reset. Here, for simplicity, we avoid dealing specifically with boundary conditions and opt for the most natural periodic setting of this numerical

problem. An extension of the analytic approach for solving adjoint variational equations in slow-fast systems is given in Appendix A.

Numerical continuation alternative for adjoints

The numerical continuation approach proposed in the present work allows us to compute limit cycles and associated adjoint solutions along a canard-explosive branch. One of the main advantages of the continuation is the possibility to find solutions in the limit ε → 0. We extend the continuation setting of the original system (1), solved in order to find limit cycles, by including equation (2) (once written in first-order form) to find periodic solutions of the adjoint problem along these cycles. In order to compute a limit cycle γ together with a periodic solution of the associated adjoint problem along γ, one needs to solve the following system of equations

Ẋ = F (X), Ż = -D X F (X) T |γ Z. (4)
Our numerical continuation strategy requires two steps: first, to find a nontrivial solution of the adjoint problem along the (initial) cycle γ, and second, to follow the extended system (4) (as a periodic continuation problem) in a 165 bifurcation parameter in order to find a branch of such solutions. In the following section we describe these steps by considering two examples of coupled slowfast systems in the canard regime, namely, the VDP system and and a twodimensional reduction of the HH model for action potential generation whose slow-fast structure and associated canard dynamics were analyzed in [START_REF] Moehlis | Canards for a reduction of Hodgkin-Huxley equations[END_REF]. In the case of the VDP system, the extended continuation setting (4) reads

x = y -f (x) y = ε(c -x) z 1 = f (γ 1 (t))z 1 + εz 2 z 2 = -z 1 , (5)
with f (x) = x 3 /3x, 0 < ε 1 and c is the bifurcation parameter displaying the canard explosion. The system given in (5) merges the original VDP system with the adjoint equation. As hinted at above, the continuation procedure is divided into two steps.

In the first step, we initialize system (5) with the limit cycle γ for the first two equations, and the trivial solution for the remaining two (which is trivially periodic). We need to obtain a non-trivial periodic solution of the adjoint equation which can be found by continuing system [START_REF]Phase Response Curves in Neuroscience: Theory, Experiment, and Analysis[END_REF] in an extra parameter. Indeed,

given that the trivial solution to the adjoint equation exists for all values of parameters c and ε, by continuing in any of these we can only hope to find a branch point and switch at this bifurcation to the non-trivial solution branch.

An alternative is to introduce a dummy parameter µ, such that system (5) becomes

x = y -f (x) y = ε(c -x) z 1 = f (γ 1 (t))z 1 + εz 2 z 2 = -z 1 + µ, (6)
and to continue the starting solution in µ along a very small interval, as small as possible. It turns out that we can compute a branch in µ and stop at µ = 10 -8 , which is indeed very small but sufficient to find a non-trivial solution of the extended problem [START_REF] Ashwin | Mathematical framework for oscillatory network dynamics in neuroscience[END_REF].

Given that µ is very small, we can, in the second step, impose back µ = 0 and run a simple Newton iteration so as to converge to a non-trivial solution of the original extended problem [START_REF]Phase Response Curves in Neuroscience: Theory, Experiment, and Analysis[END_REF]. The advantage of using numerical continuation to compute a non-trivial solution to the adjoint equation along a canard cycle is that we can then continue the extended problem (5) in parameter c and follow both the cycle and the associated periodic solution of the adjoint equation along the entire canard explosion.

Finally, the normalization condition ( 3) is required to close to the linear problem corresponding to the adjoint equation. Implementing this condition as part of our numerical continuation procedure can be a little delicate for small values of ε, therefore we decided to use a periodic continuation in auto and apply the scaling that corresponds to (3) as a post-processing step. Note that we refrain from computing the Floquet bundle to obtain the non-trivial solution of the adjoint equation for this numerical problem since we only need any non-trivial solution to the adjoint equation, which we can then normalize appropriately.

Starting from the Hopf bifurcation and continuing all the way to the relaxation regime, we can therefore follow the canard cycles by varying c together with their associated adjoint solutions. Figure 1 shows some of the orbits lying in the headless canard and in the canard with head regimes. We observe a qualitative change in the adjoint solution, where max (x(t)) point on the periodic orbit corresponds to the zero phase φ = 0, as the limit cycle γ passes through the maximal canard.

In order to see whether or not the transitions that we observe in coupled VDP oscillators are system dependent, we next compute adjoints solutions associated with canards in a planar reduction of the HH model.

Adjoints of canard cycles in a reduced Hodgkin-Huxley model

A reduction of the classical HH model to two variables was analyzed from the viewpoint of canard dynamics in [START_REF] Moehlis | Canards for a reduction of Hodgkin-Huxley equations[END_REF]; the planar system has the form

C V = I -ḡNa [m ∞ (V )] 3 (0.8 -n)(V -V N a ) -ḡK n 4 (V -V K ) -g L (V -V L ) ṅ = α n (V )(1 -n) -β n (V )n, (7)
where

α n (V ) = (0.01(V +55))/(1-exp[-(V +55)/10]), β n (V ) = 0.125 exp[-(V + 65)/80], m ∞ (V ) = α m /(α m + β m ) with α m = (0.1(V + 40))/(1 -exp[-(V + 40)/10]), β m = 0.4 exp[-(V + 65)/18]
). Moehlis has shown in [START_REF] Moehlis | Canards for a reduction of Hodgkin-Huxley equations[END_REF] that system (7) displays a canard explosion when parameter I is varied, for the following fixed values of the other parameters: ḡNa = 120, ḡK = 36, ḡL = 0.3,

V N a = 50, V K = -77, V L = -54.4, C = 1.
After verifying numerically that the dynamics of V are much faster than the dynamics of n and, hence, that the system effectively displays slow-fast dynamics, a formal asymptotic analysis was performed in ε which appeared in the rescaled form of the slow equation [START_REF] Moehlis | Canards for a reduction of Hodgkin-Huxley equations[END_REF]. Following the treatment of ε as a small parameter in asymptotic analysis and obtaining an ε-expansion of the I-value at which the canard explosion occurs, ε = 1 was plugged in the final formula.

n V t/T t/T t/T Z 1 Z 1 Z 1 S 0 1 2 3 4 5 6 7 1.5 ⇥10 7 ⇥10 11 -20 20 0 Z 1 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0 0 1 0 1 t/T t/T t/T Z 1 Z 1 Z 1 -5 5 10 ⇥10 6 0.2 0.4 0.6 0.8 1 0 t/T -1.5 ⇥10 7 -1.5 ⇥10 6
ṅ = ε(α n (V )(1 -n) -β n (V )n)) in
Despite the instability of part of the canard branch in system ( 7), the continuation strategy allows to find solutions to adjoint equations. Since we are interested in the shape of the adjoints of the canard cycles lying on different sides of the repelling slow manifold, we can ignore the stability issue. Following the same continuation procedure described above, we compute adjoints of the 2D reduced HH system. Canard cycles and corresponding adjoint solutions are visualized in Figure 2. As in the VDP system, the transition from headless canards to canards with head changes qualitatively the adjoint solution.

Consequences of a non-monotonic period function on the iPRC

As shown in Figure 3, the period function is non-monotonic along the canard explosion. It increases in the headless canard regime, reaches its maximum at the maximal canard and then decreases in the canard-with-head regime. The non-monotonicity of the period function along the explosive branch of canard cycles is one key aspect of the canard phenomenon in VDP-type systems, and the maximum of the period function can be used to detect numerically the maximal canard [START_REF] De Maesschalck | Numerical continuation techniques for planar slow-fast systems[END_REF]. The shape of this period function is sufficient to understand the effect of a perturbation of a canard cycle close enough to the lower fold of the critical manifold S 0 . Indeed, O(1) away from this fold point, a sufficiently small perturbation from the slow manifold takes the perturbed trajectory back to it very rapidly and therefore the effect of this perturbation is largely attenuated. This justifies that the solution to the adjoint equation along a canard cycle is close to zero for most of the cycle apart from the time interval corresponding to when the cycle is close to the lower fold (where the canard point is). On the other hand, near the lower fold of the critical manifold, the attraction to the slow manifold associated to the chosen canard cycle is weaker and the effect of the perturbation becomes large; see Figure 4(a2),(b2) for an illustration of this point. This effect can be understood by invoking the period function of the (x = 1), then becomes negative as the solution (x(t), y(t)) moves away from the fold region. Note that invoking the period function to explain a change of shape and sign of the adjoint solution has been used in [START_REF] Zhang | Phase response properties of half-center oscillators[END_REF] in the context of so-called escape-release mechanism for the synchronization of half-center oscillators. Here we show that it also applies in the context of coupled canard oscillators. 280

(c) (d) (b) (c) (d) -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -0.

Synchronization properties of weakly coupled canard oscillators

The behavior of the adjoint solutions (or equivalently, of the iPRCs) provides predictions on the collective behavior in the weak coupling regime via the interaction function which is the convolution of adjoint solutions and the coupling function [START_REF] Kuramoto | Chemical oscillators, waves, and turbulence[END_REF][START_REF] Malkin | Methods of Poincare and Liapunov in Theory of Non-Linear Oscillations[END_REF][START_REF] Malkin | Some Problems in Nonlinear Oscillation Theory[END_REF][START_REF]Phase Response Curves in Neuroscience: Theory, Experiment, and Analysis[END_REF][START_REF] Hoppenteadt | Weakly Connected Neuroal Networks[END_REF]. In coupled identical systems the interaction function of each oscillator reads:

H(φ j -φ i ) = 1 T T 0 Z(t)U j (γ(t), γ(t + φ j -φ i ))dt. ( 8
)
where φ jφ i (i = {1, 2}, j = 3i) is the phase difference between the two oscillators and U is the coupling function. The dynamics of the phase difference, φ = φ jφ i , is described by the following equation

dφ dt = α[H(-φ) -H(φ)] = αG(φ). (9)
where 1 α > 0 is the coupling strength. Equation ( 9) has a stable solution at In the case of coupled identical oscillators, that both IP and AP solutions are guaranteed to exist [START_REF] Ashwin | The dynamics of n weakly coupled identical oscillators[END_REF].

φ * if G (φ * ) < 0,
IP synchronization of two identical relaxation cycles (coming from oscillators with cubic-shaped fast nullclines) that are weakly coupled via fast to fast (FF) connections has been shown in [START_REF] Somers | Rapid transition though fast threshold modulation[END_REF][START_REF] Kopell | Anti-phase solutions in relaxation oscillators coupled through excitatory interactions[END_REF][START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF][START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF][START_REF] Lee | Stable antiphase oscillations in a network of electrically coupled neurons[END_REF][START_REF] Cymbalyuk | In-phase and antiphase selfoscillations in a model of two electrically coupled pacemakers[END_REF], outside the canard regime.

In addition to FF coupling -which is the coupling function generally considered since it acts as a prototype for the electrical interaction between neuronal systems-we consider fast to slow (FS) coupling, which is not physiologically realistic but provides insight into understanding the interactions between perturbation and canards. The FF-coupled VDP oscillators read:

ε ẋi = y i + x i - x 3 i 3 + α(x j -x i ), ẏi = (c -x i ), (10)
and the FS-coupled system is given by The effect of a small perturbation on the canard cycles in the neighborhood of the lower fold of the critical manifold S 0 , is different for canards with head than for headless canard cycles, as revealed by the corresponding adjoint solutions;

ε ẋi = y i + x i - x 3 i 3 , ẏi = c -x i + α(x i -x j ). ( 11
)
see Figure 1. This qualitative change occurs at the maximal canard. Figure 5 shows the interaction functions H associated with the cycles in the neighborhood of the maximal canard (shown in Figure 3 The function G (see Figure 6) is computed for the cycles (whose adjoints are presented in Figure 1) interacting via FF coupling. The location of the zeros φ * of G, and the sign of its derivative at such points, determine the type and stability of synchronized state of the coupled system. The IP synchronized solution which exists for the Hopf cycles (not shown on this figure) loses stability along the canard explosion (due to the high sensitivity to perturbation resulting from the passage near the fold of the critical manifold S 0 ) and a stable OP solution appears for the headless canard cycles (orbits 1-4). The phase difference of the stable OP solution increases as the cycle approaches the maximal canard (Panels 1-4). Bistability appears for the canards-with-head (orbits 5-7), where IP and AP solutions are the stable synchronous solutions and the OP is the unstable solution (Panels 5-7).

The information obtained with the function G about synchronized states of the weakly coupled VDP system with FF coupling, can be confirmed by a numerical bifurcation analysis of the coupled system in question. We have performed this analysis by continuing synchronous states of system (10) (including the ones which are not visualized in Figure 6) in parameter c. The result is presented in Figure 7 where the chosen solution measure is the difference between the x-component of each oscillator at time t = 0, x 2 (0)x 1 (0), regardless of its varying amplitude as a function of c. That measure has the same interpretation as the phase difference for these simple orbits and it is often used in the analysis of weakly coupled oscillators [START_REF] Van Vreeswijk | When inhibition not excitation synchronizes neural firing[END_REF]. Panels (b) to (d) are successive zooms of panel (a) in the region corresponding to maximal canards for each oscillator. The properties of the synchronized states of the FF-coupled cycles are tracked starting from the double Hopf bifurcation point at c = c Hopf = 1 down to the relaxation regime near c ≈ 0.615. We consider a fixed coupling strength α = 10 -5 for which the weakly coupled oscillators theory is expected to be valid; a detailed discussion on the effect of α is presented in Section 4.

One stable and two (symmetric) unstable branches, which correspond to IP and AP solutions, respectively, appear at c = c Hopf . The IP solution undergoes a pitchfork bifurcation through which it loses its stability as a stable OP solution appears (Panel (b)). The OP branches become unstable at a PD bifurcation which is followed by a PD cascade corresponding to 2nT -periodic stable synchronous solutions (Panel (d)), where the interaction function analysis is not valid. The T -periodic OP branches become stable again via a second PD bifurcation. It changes its stability two times via a couple of fold bifurcations before connecting to the second pitchfork bifurcation point on the IP branch that restabilizes the IP state.

The unstable AP branch that appears at c Hopf becomes stable at the maximal canard of the coupled system through a pitchfork bifurcation (Panel (c)). The stable AP and OP solution related to this pitchfork bifurcation coexist with stable IP solutions for a some range of c in the neighborhood of the maximal canard. For smaller values of c, IP and AP remain stable, while OP states are unstable.

The bistability regions (illustrated in Figure 8) already hinted at with the investigation of the function G, are well identified through the continuation analysis, in particular the coexisting stable IP and stable AP states (canards with head and relaxation cycles) born near the maximal canard solutions. This intricate bifurcation structure unveils a main connection between the stable IP and the AP states through the double Hopf point at c = c Hopf , which gives rise to both the IP stable state and a branch of unstable AP states. Decreasing c further, additional bifurcations occur, in particular pitchfork bifurcation points (black dots in Figure 7 (a) to (c)) which correspond to events where the synchronous state loses some symmetry. Indeed, on both the IP and the AP branches these bifurcations lead to additional solution branches along which the two canard oscillators do not follow identical cycles; in each case, the synchronous state becomes identical again through fold bifurcations. Note that these non-identical branches emanating from both the IP and the AP states come close to each other (near a second pair of fold bifurcations) forming a structure that seems to be a broken transcritical bifurcation. This perturbed bifurcation is only conjectured here; a more detailed analysis of the ε-dependence of the synchronous states goes beyond the scope of this paper and will be a question for future work. We simply remark that this structure seems to perturb from an additional connection between the stable IP and stable AP coupled canard states. Finally we note the presence of several sequences of PD bifurcations (colored dots in Figure 7 (d)) which are likely to indicate small zones of chaotic dynamics in this region of parameter space.

One striking element about the bifurcation diagram described above is the fact that most of the connecting branches between the IP and the AP synchronous states are organized near solutions that correspond to maximal canards. It is therefore natural to ask about the effect of the coupling strength α on such synchronous states containing maximal canard segments; we focus on this aspect in the next section.

Effect of the coupling strength α

The interaction function analysis reveals the existence and stability of synchronous states for weakly coupled oscillators, although how "weak" the coupling should be in order that the theory applies is questionable. For instance, it was shown in [START_REF] Lewis | Dynamics of spiking neurons connected by both inhibitory and electrical coupling[END_REF] that for leakly integrate-and-fire type of oscillators the H function analysis is valid for moderate coupling strengths, whereas other papers (see e.g. [START_REF] Oh | Loss of phase-locking in non-weakly coupled inhibitory networks of type-I model neurons[END_REF][START_REF] Maran | Using phase resetting to predict 1:1 and 2:2 locking in two neuron networks in which firing order is not always preserved[END_REF]) have mentioned a loss of 1:1 phase locking estimated by the interaction function analysis. In the case of coupled canard-explosive systems where the properties of the underlying oscillators vary brutally in parameter ranges that are exponentially small in time-scale parameter ε, the notion of weak coupling can be even more vague. For instance the region with cascades of PD bifurcations, highlighted in Figure 7 (d) and corresponding to cycles that are close to the maximal canard regime (under weak coupling of strength α = 10 -5 ), gives a good numerical evidence that canard orbits are very sensitive to perturbations and that the validity of the interaction function analysis is limited in such cases.

In order to investigate this aspect further, we next consider the phase difference dynamics of two coupled identical headless canard cycles for a c-value in the neighborhood of the maximal canard, as a function of the coupling strength α > 0. This numerical continuation study will focus both FF and FS interac- A stable OP synchronous state with a phase difference φ * = 0.34 is predicted by the interaction function analysis for the case of two headless canard cycles with FF-coupling, that is, for system [START_REF] Roseau | Vibrations nonlinéaires et théorie de la stabilité[END_REF]; see Figure 5 (a) and Figure 6 panel 4. Using the bifurcation diagram presented in Figure 9, we can conclude that this OP regime persists for α ∈ (0, 6.63371×10 -5 ]. It loses its stability at α ≈ 6.63371 × 10 -5 via a PD bifurcation where the interaction function result is violated, and consequently, not valid for greater coupling strengths.

Switching branch at this PD point reveals the presence of a PD cascade, for which we compute only a few subsequent branches. Among the stable part of these branches of period-2nT synchronous solutions (near which chaotic orbits (black), 2T -periodic (green) and 4T -periodic (blue) branches coexist with stable (solid) and unstable (dashed) solutions. surely exist too), that is, for a coupling strength α ∈ (6.63371×10 -5 , 0.0083195],

x 2 (0) x 1 (0) ↵ -1.5 -0.75 0 0.1 0.2 0 0.16 0.28 -0.08 0 (a) • • • • • • • • • (a) (b) -1.4 0 0 0.008 (b) • • • • • • • • •• •• • •
• • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • 0 0.00015 0.0003 -0.9 -1.4 -0.4 ↵ x 2 (0) x 1 (0) -0.75 -0.4 0.00028 0.000292 ↵ x 2 (0) x 1 (0) -1.03 -1.01 0.00005 0.00008 ↵ x 2 (0) x 1 (0) (a) (b) (a) (b)
y 1 , y 2 2 -2 0 x 1 , x 2 -0.5 0.5 1.5 y 1 , y 2 2 -2 0 x 1 , x 2 x1(t) x2(t) x1(t) x2(t)
there exists a family of solutions displaying what we call "spike suppression".

This scenario corresponds to when one of the oscillators spikes by following a canard with head while the other always remain in the headless canard regime.

Regarding the IP solution branch, it becomes stable at α ≈ 0.0085633416545 and coexists, for α ∈ [0.0085633416545, 0.289498], with the 2nT -periodic headless canard solution branch.

Fast-to-Slow (FS) coupling. The bifurcation structure in α for this case in presented in Figure 13; associated solution profiles are shown in Figure 14. The stable IP synchronization state predicted by the interaction function analysis for the FS-coupling (Figure 5 (b)) becomes unstable at α ≈ 0.007498445 (Figure 13 (a)) via a subcritical PD bifurcation that introduces an unstable 2T -periodic branch which becomes stable at α ≈ 8.74785268 × 10 -5 , where the interaction function analysis loses its validity. Continuing that branch leads to the detec-

S1 S2 S1 S2 S1 S2 2 -0.5 0.5 1.5 -2 0 x 1 , x 2 y 1 , y 2 S 0 2 -0.5 0.5 1.5 -2 0 x 1 , x 2 y 1 , y 2 S 0 2 -0.5 0.5 1.5 -2 0 x 1 , x 2 y 1 , y 2 S 0 x1(t) x2 (t) x1(t) x2(t) x1(t) x2(t
Figure 12: Period T (top panels, identical), period 2T (middle panels, non-identical) and period 4T (bottom panels, non-identical) stable OP synchronous states of the FF-coupled VDP system in the maximal canard regime. The phase differences for these states are coherent with the interaction function analysis. Values of the coupling strength α are 6.63371 × 10 -5

in panels (a1)-(a2), 7.04717 × 10 -5 in panels (b1)-(b2) and 7.13322 × 10 -5 in panels (c1)-(c2).

Left panels: Trajectories projected onto the (x i , y i ) planes. Right panels: Time series of the

x i coordinates.
tion of further PD bifurcations organized in a cascade, which we compute only the beginning of; see Figure 13 (b). These nT -periodic branches correspond to families of solutions displaying what we call "spike alternation", that is, a scenario for which both oscillators of the FS-coupled system follow subsequently a headless canard segment and then a canard-with-head segment, hence performing an MMO [START_REF] Desroches | Mixed-mode oscillations in a multiple time scale[END_REF]; see Figure 14 for an illustration on such MMO cycles with period T , 2T and 4T on (a), (b) and (c) panels, respectively. Depending on the value of the coupling strength α, the oscillators may follow the same or different canard trajectories.

On both FF-and FS-coupled canard systems, we have observed using a numerical bifurcation analysis the proximity of several stable solution branches 0.0032 0.00295 -0.18 -0.11 -0.04 0.0027 with complicated oscillatory patterns mixing passages along headless canards and along canards with head. In the context of neuronal systems, these solutions alternate subthreshold oscillations and spikes. These solutions are not predicted by the interaction function analysis typically employed in weakly coupled oscillator studies. However, one can justify their existence by invoking the presence in such systems of repelling (Fenichel) slow manifolds, which are known to be exponentially close to each other (in the timescale separation parameter ε).

x 2 (0) x 1 (0) ↵ (b) x 2 (0) x 1 (0) -0.9 -0.4 0.1 0 0.005 0.01 ↵ (a) • • • • • • • • • • •
Therefore, the presence of these manifolds near the middle branch of the critical manifold S 0 of each individual slow-fast oscillator can allow to explain why, for values of the coupling strength α that are larger that exponentially small quantities, synchronized states of the coupled system may follow these manifolds on one side (subthreshold regime) or the other (spiking regime) while staying very close to the boundary (well approximated by maximal canards). In other words, when two identical slow-fast systems are weakly coupled, the existence of a repelling slow manifold and of an associated maximal canard trajectory in the uncoupled system can give rise to solutions to the full (coupled) system for which each node follows this maximal canard on opposite sides, hence, separate after an O(1) time. This can happen as soon as the coupling is stronger than an exponentially small function of ε and can therefore be responsible for the presence of canard-induced states in the coupled system; see [START_REF] Ersöz | Canard-mediated (de)synchronization in coupled phantom bursters[END_REF] for an example of this phenomenon in weakly coupled folded-singularity systems.

Discussion

In this paper, we have extended previous results on weakly coupled slow-fast oscillators to the canard regime, both from theoretical and numerical perspectives. Our main finding is that the behavior of adjoint solutions (or equivalently, of iPRCs) changes qualitatively when the canard cycle under consideration is moving (as the canard parameter is varied) along the associated explosive branch. Indeed, the sign and shape of the adjoint solutions flip as the underlying canard cycle goes from the headless canard regime to the canard-with-head regime, the transition taking place at the maximal canard cycle, which in particular sheds new light onto the previously unnoticed role played by this special canard in the context of coupled slow-fast oscillators, deeply connected to the fact that it corresponds to a critical point of the period function [START_REF] De Maesschalck | Numerical continuation techniques for planar slow-fast systems[END_REF]. This change of behavior of adjoints of canard cycles upon infinitesimal perturbations can be explained by the peculiar known property of the period function of a canard-explosive branch, which can be summarized as follows: larger headless canards have greater periods, whereas larger canards-with-head have smaller periods. As explained in Section 3, this arguments is fully applicable when the perturbation is applied near the fold point of the critical manifold corresponding to the canard point, and its validity is weakened as the perturbation is applied further away from this fold point, where the contraction towards the unperturbed cycle rapidly annihilates the effect of the perturbation. This justifies that adjoints computed along canard cycles are very close to zero during most of the cycle except along a time interval corresponding to when the canard cycle passes near the fold (canard) point of the critical manifold S 0 . Nevertheless, the explanation that we provide is valid for the most informative part of the adjoint solutions and bears consequences on the synchronized solutions of coupled canard systems.

We have shown this mechanism for a prototypical canard oscillator, namely the VDP system, but it is clearly applicable to all excitable systems of this form, in particular, to slow-fast type-II neuron models such as the reduced HH model studied in [START_REF] Moehlis | Canards for a reduction of Hodgkin-Huxley equations[END_REF]. This opens the way to a renewed understanding of iPRCs in such neuron models, from the Hopf cycles (whose adjoint solutions will qualitatively look like those associated with small headless canard cycles) to the spiking cycles (whose adjoint solutions will qualitatively look like those associated with canards with head). In particular, our findings can be related to recent work on isochrones since PRC analysis originates in the study of phase models and isochrones [START_REF] Winfree | The geometry of biological time[END_REF]. Recently the isochrons of canard cycles were investigated numerically in [START_REF] Langfield | Solving Winfree's puzzle: The isochrons in the FitzHugh-Nagumo model[END_REF][START_REF] Langfield | Forward-time and backwardtime isochrons and their interactions[END_REF] where evidence was given that their properties change in the vicinity of the maximal canard neighborhood; this is likely to be closely linked with the results presented here. While a full comparison of these two aspects of canard cycles' phase properties goes beyond the scope of the present work, it is certainly an interesting topic for future work.

While studying adjoint solutions along canard cycles, we have also proposed a numerical strategy based on numerical continuation in auto to compute these objects as a system parameter is varied, that is, to reliably compute a family of limit cycles and at the same time a family periodic solutions to the associated adjoint problem. Making use of the boundary-value solver of auto, we could easily identify the flip in the solution to the adjoint problem as the cycle goes through the maximal canard. Our strategy with a one-step homotopy approach to compute a non-trivial solution to the adjoint equation associated with a given limit cycle, and then continue both this solution and the limit cycle as an extended periodic continuation problem is more elementary than the one developed by Govaerts and Sautois in [START_REF] Govaerts | Computation of the phase response curve: A direct numerical approach[END_REF] for MatCont, as it simply relies on a periodic continuation, yet giving access to the objects of interest. Moreover, the simplicity of our approach makes it easily adaptable to other systems and we believe that it is a interesting computation for the community of auto users.

In Section 3, we looked at the bifurcation structure of the synchronous states of the weakly coupled identical VDP systems when varying the main system parameter, which in this case controls the position of the slow nullcline but would likely be an applied current in the neuronal context. We found an intricate structure of solution branches of IP, AP and OP states, connected through both PD and pitchfork bifurcations, which are organized around the maximal canard solution. While the synchronization properties of relaxation cycles were already known, we believe that the bifurcation structure of the weakly coupled canard regime is by-and-large novel, in particular the role of the maximal canard as an organizing center for the IP, AP and OP families.

In Section 4, we focused on the bifurcation structure of synchronous states of coupled identical VDP systems in the maximal canard regime, depending on coupling strength α. PD bifurcations and chaotic trajectories in VDP-like systems under periodic perturbation have been studied in e.g. [START_REF] Ermentrout | Period doublings and possible chaos in neural models[END_REF][START_REF] Feingold | Phase locking, period doubling, and chaotic phenomena in externally driven excitable systems[END_REF][START_REF] Schwartz | Subharmonic hysteresis and period doubling bifurcations for a periodically driven laser[END_REF]. In the present study, we unveiled a complex web of period-2nT branches suggestive of the presence of nearby chaotic attractors, which we chose not to investigate.

Instead, we highlighted these further synchronous states, all existing close to maximal canard solutions but not all predicted by standard interaction function analysis. Being close to a maximal canard, hence to threshold, these solutions may contain both passages near headless canards and near canards-with-head, therefore an alternation between subthreshold oscillations and spikes. Even when the classical weakly coupled theory may not apply, the slow-fast phase plane structure of the underlying single canard oscillator enables one to understand why such mixed-mode oscillatory synchronous states can arise for small to moderate coupling strength, owing to the geometry and proximity between families of repelling slow manifolds. As a question for future work, we plan to investigate the relevance of these complicated synchronous states in the context of neuron models, where canards-with-head may be considered as not so rare events but rather as spikes with a slow activation.

Control aspects of canard cycles have been studied in [START_REF] Durham | Feedback control of canards[END_REF] where the authors have obtained MMOs, cascades of PD bifurcations and chaotic behavior in a FHN-type relaxation oscillator depending on the control setup. Developing control strategies for reaching desired spiking behavior in coupled systems can be a future direction of our study.

Finally, as an appendix, we also provided an analytical formula for the adjoint solutions associated with limit cycles of Liénard systems, which gives reasonable numerical results for headless canard cycles.

This work is only a first step towards extending canard studies to the realm of weakly coupled oscillators and, more generally, to weakly connected networks.

It is not rigorous yet but we have identified the main geometrical structures that play a pivotal role in shaping the main family of synchronous solutions to coupled planar slow-fast systems in the canard regime. Moreover, we have highlighted the central role of the maximal canard in organizing the synchronization properties of such systems. Maximal canards have been identified as the spiking threshold in excitable neuron models [START_REF] Desroches | Inflection, canards and excitability threshold in neuronal models[END_REF][START_REF] Mitry | Exitable neurons, firing threshold manifolds and canards[END_REF]. Our study on weakly coupled pla-nar systems, which are reduced models for excitable neurons, has underlined the relation between synchronisation, excitability and maximal canard by relating it to the properties of the period function of the canard explosion. The interaction between the slow-fast structure and the weak coupling has given quite rich dynamics in these simple planar systems having already a two slow/two fast structure under coupling. Similar possible links between synchronisation, excitability and maximal canards certainly deserves further attention in coupled slow-fast systems of higher dimensions as well as in larger networks. Beyond the effect of canard-explosive dynamics on synchronization, we plan in the near future to investigate similar effects in (at least three-dimensional) systems with canards organized by folded singularities [START_REF] Krupa | Weakly coupled two-slow-twofast systems, folded singularities and mixed mode oscillations[END_REF] as well as in systems with slowly varying quantities, such as bursting systems where spike-adding canard explosions will be likely to have a dramatic effect on the synchronization properties of coupled bursters [START_REF] Vries | Beyond synchronization: Modulatory and emergent effects of coupling in square-wave bursting[END_REF][START_REF] Desroches | Mixed-mode bursting oscillations: Dynamics created by a slow passage through spike-adding canard explosion in a square-wave burster[END_REF].

The adjoint equation associated with system (A.1) along the limit cycle γ is given by Ż = -J(γ(t)) T Z, (A.7)

where Z is a two-dimensional real vector and J(γ(t)) is the Jacobian matrix evaluated along the solution γ. Following [START_REF] Schecter | Adjoint equation and Melnikov function[END_REF], we write the solution to equation (A. where s = (s 1 , s 2 ) is a solution to the linearized equation (A.3). We apply this formula to the two solutions (γ 1 (t), γ 2 (t)) and (v(t), w(t)) of the linearized equation, which gives us two solutions of the adjoint equation. What we wish to get is a periodic solution of the adjoint; to get it, we will find a suitable linear combination of the two solutions obtained using Schecter's strategy, imposing periodicity. Namely, we will find scalars α and β such that αZ γ (T ) + βZ s (T ) = αZ γ (0) + βZ s (0), (A.9)

where α and β are reals, Z γ and Z s being obtained using formula (A.8) from the linearization of the limit cycle γ and the solution (v(t), w(t)) described above, respectively. Therefore, focusing on the second component only, the periodicity condition (A.9) becomes α exp Limits of the formula. The strategy we proposed overcomes the singularities due to the presence of folds on the critical manifold, however it has limitations.

First, the approximation of adjoint solutions of canard cycles with this formula can be considered as successful for headless canards (see Figure A.15), yet a lot of care in the numerical simulations used is required. However, even with such care we have been unable to compute adjoints associated with large canards using this formula. The reason for this can be understood by looking the expression in (A.13) which is singular when γ 1 (t) = 0, that is, at extrema of γ 1 . We can try to integrate these equations by splitting the solution into two branches excluding the extrema. With this strategy, our formula can be used to compute adjoints for all canard cycles and, hence, extend Izhikevich's approach. The second drawback of formula (A.8) is that it assumes a Liénard form for the system under consideration. Hence, it is not directly applicable to more general planar slow-fast systems, in particular, to biophysical neuron models such as the 2D reduction of the HH system that we considered in Section 2.2.2.

170 2 . 2 . 1 .

 221 Adjoint solutions of the VDP system

Figure 1 :

 1 Figure 1: Top left panel: Canard orbits of the VDP system in the phase plane, for ε = 0.1. Panels 1-7: time profile of the first component of the adjoint solution associated with each canard cycle shown in the phase plane (together with the critical manifold S 0 := {y = f (x)}), keeping the same color coding. A qualitative change in the adjoint solution occurs in between Orbit 4 and Orbit 5, corresponding to the passage through the maximal canard cycle.

Figure 2 :

 2 Figure 2: Top left panel: Canard orbits of the reduced HH system in the phase plane. Panels 1-7: time profile of the first component of the adjoint solution associated with each canard cycle shown in the phase plane (together with the critical manifold S 0 := { V = 0}), keeping the same color coding. A qualitative change in the adjoint solution occurs in between Orbit 4 and Orbit 5, corresponding to the passage through the maximal canard cycle.

(

Figure 3 :

 3 Figure 3: (a) Period of limit cycles along the canard explosion in the VDP system for ε = 0.1;the parameter that varies is c. The period is increasing along the headless canard part of the branch, it reaches its maximum at the maximal canard and then decreases along the canard-with-head cycles. (b) Three headless canard cycles and their periods marked on the period curve. Smaller cycles have smaller periods. (c) Three cycles in the neighborhood of the maximal canard, together with their periods marked on the period curve. Canards with head and headless canards have very close periods in this vicinity. (d) Three canards with head and their periods marked on the period curve. Larger cycles have smaller periods. Also shown on panels (b) to (d) is the critical manifold S 0 , on which solid (resp. dashed) parts represent stable (resp. unstable) branches.

Figure 4 :Figure 5 :

 45 Figure 4: (a1, b1) Transient effect (dashed curves) of a small perturbation of the canard cycles (red solid curves) in the positive x-direction. (a2)-(b2) time profile of the first component of the adjoint solution associated with the red canard cycles and (inset) (x(t), y(t)) during one cycle. Perturbing a headless canard (resp. a canard with head) away from the attracting slow manifold (yellow asterisk) delays (resp. advances) its phase by driving it to a larger yet slower (resp. faster) yellow dashed cycle. Perturbing a headless canard (resp. a canard with head) towards the repelling slow manifold (blue asterisk) advances (resp. delays) its phase by driving it to a smaller yet faster (resp. slower) blue dashed cycle.

  meaning that the two oscillators will synchronize with a phase difference φ * . The solution φ * = 0 corresponds to IP synchronization, φ * = π (or equivalently φ * = 0.5 if the phase is rescaled to [0,1]) to AP synchronization, and any other values of φ * corresponds to OP synchronization of coupled oscillators.

Figure 6 :

 6 Figure 6: Selection of canard cycles of the VDP oscillator in the phase plane (x, y) (top left panel) together with the corresponding G functions (panels 1 to 7; the phase φ is rescaled to [0,1]).

  (c)) interacting via FF (panel (a)) and FS (panel (b)) connections. Given that a headless canard cycle resembles more the maximal canard (the maximal canard being a maximal headless canard), the amplitude of the corresponding function H decreases while the number of zeros, that is, solutions to H(φ * ) = 0, and the sign of H (φ * ) remain the same. The sign of H (φ * ) changes when the cycle moves to the canard-with-head regime, while the number of zeros φ * is preserved.

FFigure 7 :Figure 8 :

 78 Figure 7: Bifurcation diagram of system (10) with respect to variations of c for α = 10 -5 , from the Hopf regime to the relaxation regime. The output solution measure is the difference between the first components of each oscillator at time t = 0. The region of the maximal canard is enlarged from left to right and top to bottom panels. Black dots in panels (a) to (c) denote pitchfork bifurcation points; the black star in panel (b) corresponds to the double Hopf point that initiate the periodic regime in this coupled system; colored dots in panel (d) denote PD bifurcation points.

  tions. The aim is to identify what range of the perturbation strength can give rise to interesting canard-mediated dynamics that are not predicted by the interaction function analysis but whose existence can be justified using slow-fast arguments.Fast-to-Fast (FF) coupling. The bifurcation structure in α for this case in presented in Figures9 and 10(zoomed views); associated solution profiles are shown in Figures 11 and 12.

Figure 9 :

 9 Figure 9: Continuation in α for the FF-coupled VDP system for a c value in the vicinity of the maximal canard. Inset panels (a) and (b) are zoomed viewed of different parts of the main panel. Bifurcation points (mainly PD bifurcations) are indicated by red dots. T-periodic

Figure 10 :

 10 Figure 10: Continuation in α for the FF-coupled VDP system in the maximal canard regime: zoomed view from Figure 9 in the region of PD cascades (most of the computed PD bifurcation points being highlighted by colored dots).

Figure 11 :

 11 Figure 11: Period 2T (top panels) and period 4T (bottom panels) non-identical OP synchronous states for the FF-coupled system in the maximal canard regime, illustrating the spike suppression scenario. Values of the coupling strength α are 2.86959 × 10 -4 in panels (a1)-(a2) and 2.85359 × 10 -4 in panels (b1)-(b2).

Figure 13 :

 13 Figure 13: Continuation in α for the FS-coupled VDP system. Bifurcation points (PD bifurcations) are indicated by red dots. Both stable (solid) and unstable (dashed) parts of T -periodic (black), 2T -periodic (red), 4T -periodic (green) and 8T -periodic (blue) branches are shown.

x 1 , x 2 y 1

 121 , y 2

Figure 14 :

 14 Figure 14: Period T (top panels), 2T (middle panels) and 4T (bottom panels) stable nonidentical OP synchronous states displaying spike alternation for the FS-coupled system near the maximal canard regime. Values of the coupling strength α are 2.8502655978 × 10 -3 in panels (a1)-(a2), 2.8939484985 × 10 -3 in panels (b1)-(b2) and 2.9039987077 × 10 -3 in panels (c1)-(c2). Left panels: Trajectories projected onto the (x i , y i ) planes. Right panels: Time series of the x i coordinates.

7 ) 0 f

 70 in the formZ T (t) = exp t (γ 1 (s))ds -s 2 s 1 , (A.8)

T 0 fFigure A. 15 : 0 f (γ 1

 01501 Figure A.15: Adjoint solutions for the headless canard cycles shown in Figure 1 computed analytically using formula (A.8). Blue curve: Z 1 . Red curve: Z 2 .

Appendix A. Analytical expression for adjoints of canard cycles

Here we use classical results from the theory of linear differential equations [START_REF] Coddington | Theory of ordinary differential equations[END_REF] as well as unpublished results by Schecter [START_REF] Schecter | Adjoint equation and Melnikov function[END_REF] in order to derive an expression for the periodic solution of the adjoint problem associated with a limit cycle of a Liénard system. This extends the approach taken by Izhikevich in [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF], who considered the case of relaxation cycles by taking the limit ε = 0.

Izhikevich's formulation is not applicable to canard cycles due to the presence of the folds of the critical manifold S 0 which gives rise to canard dynamics and requires to have ε = 0 in the computation of the adjoints.

We consider the following VDP type slow-fast system written in Liénard form

where f (x) = x 3 /3x is a cubic function and the prime denotes differentiation with respect to the fast time t. We consider a canard cycle solution of system (A.1), that is, a periodic solution γ(t) = (x(εt), y(εt)).

The linearized system associated with (A.1) along γ is given by

which we can recast as a second-order linear differential equation

An obvious solution of (A.3) is (γ 1 (t), γ 2 (t)). Recall that if one knows a particular solution y * of the second-order linear differential y (t) + p(t)y (t) + q(t)y(t) = 0, then one can obtain another solution y # , non-proportional to the first onehence forming a basis of the space of solutions together with the first oneusing a variation of constant type formula, that is, Therefore, knowing the solution (γ 1 (t), γ 2 (t)) of the linearized system written as a second-order equation (A.3), a non-proportional solution is given by