Soledad Fernández-García

Dpto Edan&imus

Alexandre Vidal

Symmetric coupling of multiple timescale systems with Mixed-Mode Oscillations and synchronization

Keywords: Slow-Fast Dynamics, Coupled Oscillators, Mixed-Mode Oscillations, Synchronization, Nonlinear Neuron Model. 2010 MSC: Primary: 34C15, 34C23, 34D06, Secondary: 70K70, 92B25

We analyze a six-dimensional slow-fast system consisting of two symmetrically coupled identical oscillators. Each oscillator is a three-dimensional slow-fast system formed by a FitzHugh-Nagumo model with an additional slow variable. Individually, each three-dimensional subsystem admits an attractive Mixed-Mode Oscillations limit cycle, featuring small oscillations due to the presence of a folded singularity. We consider a linear symmetric coupling of the fast variable of each oscillator upon one of the slow equations of the other one and study the synchronization properties between oscillators according to the value of the coupling gain parameter.

We examine theoretically and numerically both excitatory and inhibitory cases associated with positive and negative values of the coupling gain. Apart from in-phase locking synchronization and antiphase synchronization patterns, the coupled system generates phase-locking almost-in-phase synchronization, relaxation loss of one of the oscillators and total oscillation death, intertwined with complex transitions involving period doubling cascades, period adding phenomena and chaos. We present a theoretical proof of the behavior repartition ac-

Introduction

Neural synchronization is a basic functioning principle of the nervous system resulting from the coupling between the activities of different neurons (microscopic scale), neuronal ensembles (mesoscopic scale), and brain areas (macroscopic scale). At the microscopic scale, recent advances in optical imaging, have provided evidence that calcium activities measured in-vivo (e.g. in the zebrafish embryo [START_REF] Fallani | Hierarchy of neural organization in the embryonic spinal cord: Granger-causality graph analysis of in vivo calcium imaging data[END_REF]), exhibit characteristic synchronized behaviors in the brain, as well as in the spinal cord. The data3 reveal a very strong synchronization (along the whole 4 minutes recording) of the oscillations in intracellular calcium levels between motoneurons located on the left and right side of the spinal cord respectively, while the activity profiles of these two clusters feature "antiphasic" properties: no calcium peak occur simultaneously in both clusters. A natural interpretation of these data points out the existence of an excitatory intra-cluster coupling and an inhibitory inter-cluster coupling.

To study theoretically the different features of synchronization in oscillatory patterns of intracellular calcium concentrations (ICC) between neurons, we propose a compact system composed by two coupled three-dimensional (3D) slow-fast systems, based on classical FitzHugh-Nagumo dynamics [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF]. This model has been originally designed for reproducing the changes in ICC of a single gonadotropin-releasing hormone (GnRH) expressing neuron [START_REF] Krupa | A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons[END_REF], taking particular advantage of the Mixed-Mode Oscillations feature for reaching the quiescence/active phase duration ratio. While considering two coupled oscillators for representing the activity of clusters of strongly synchronized neurons is an oversimplification of the real network phenomenon, its simplicity allows us to perform singular perturbation analysis and study, at an aggregated level, the fundamental properties of the system in a formal way, and therefore to analyze the underlying dynamical mechanisms of synchronization between neurons and clusters in case of inhibitory or excitatory coupling.

Synchronization properties of systems with multiple timescales, as relaxation oscillators, are known to be different from harmonic ones [START_REF] Bélair | On linearly coupled relaxation oscillations[END_REF][START_REF] Campbell | Multistability in coupled FitzHugh-Nagumo oscillators[END_REF][START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF][START_REF] Lee | Stable antiphase oscillations in a network of electrically coupled model neurons[END_REF][START_REF] Lee | Stability of antiphase oscillations in a network of inhibitory neurons[END_REF][START_REF] Storti | Dynamics of two strongly coupled relaxation oscillators[END_REF][START_REF] Storti | A simplified model of coupled 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].

In particular, through the Fast Threshold Modulation coupling [START_REF] Kopell | Mechanisms of phase-locking and frequency control in pairs of coupled neural oscillators[END_REF][START_REF] Kopell | Anti-phase solutions in relaxation oscillators coupled through excitatory interactions[END_REF][START_REF] Somers | Rapid synchronization through fast threshold modulation[END_REF][START_REF] Somers | Waves and synchrony in networks of oscillators of relaxation and non-relaxation type[END_REF], the synchronization is very rapid and relatively independent from the coupling strength. The impact of canard phenomena in weakly coupled slow-fast systems has been studied in several works that analyze different questions, such as synchronization and/or desynchronization [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] Ersöz | Canard-mediated (de) synchronization in coupled phantom bursters[END_REF][START_REF] Ersöz | Synchronization of weakly coupled canard oscillators[END_REF][START_REF] Rotstein | Localized and asynchronous patterns via canards in coupled calcium oscillators[END_REF], localized oscillatory patterns and formation of clusters [START_REF] Ermentrout | Canards, clusters, and synchronization in a weakly coupled interneuron model[END_REF][START_REF] Zhabotinsky | 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]. Mechanisms underlying synchronization of systems presenting bursting activity (alternation of slow quiescent phases and fast oscillatory phases) are analyzed in [START_REF] Best | The dynamic range of bursting in a model respiratory pacemaker network[END_REF][START_REF] Izhikevich | Synchronization of elliptic bursters[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] Sherman | Anti-phase, asymmetric and aperiodic oscillations in excitable cells -I. coupled bursters[END_REF].

The usual coupling between computational neuron models considered in the literature is either electrical (gap junction) or synaptic, since they are generally designed for studying the neuron-neuron communication. In both cases, the coupling between slow-fast neuron models impact the fast variable. In our study, we focus on the underlying dynamics of ICC, on which the transmission through synapses or gap-junctions between neurons has a marginal effect. For instance, it has been shown in [START_REF] Jahn | Temporospatial coupling of networked synaptic activation of AMPA-type glutamate receptor channels and calcium transients in cultured motoneurons[END_REF] that spontaneous postsynaptic currents of cultured motoneurons from the rat did not elicit significant ICC transients.

Significant elevations of ICC result from release of calcium from the endoplasmic reticulum through calcium-induced calcium release (CICR) mechanism. The dominant channel of CICR involves the non neuronal glial cells, especially the astrocytes, and the elevation of Inositol trisphosphate (IP3) from neighbor to neighbor through gap junction (see for instance [START_REF] Haydon | GLIA: Listening and talking to the synapse[END_REF][START_REF] Vandenberghe | Subcellular localization of calcium-permeable ampa receptors in spinal motoneurons[END_REF][START_REF] De Pittà | A tale of two stories: astrocyte regulation of synaptic depression and facilitation[END_REF]). This regulation mechanism is slower than neuronal transmission [START_REF] Wang | The astrocyte odyssey[END_REF], yet it is central in the functional neuron-glia system and dysfunctions of the signaling pathway have been identified as playing a major role in motoneuron diseases [START_REF] Blackburn | Astrocyte function and role in motor neuron disease: a future therapeutic target?[END_REF][START_REF] Duchen | Mitochondria, calcium-dependent neuronal death and neurodegenerative disease[END_REF]. Following a phenomenological approach, which consists in selecting the dynamics according to the consistency between the model outputs and experimental observations, we have chosen the coupling via the slow recovery variable, which represents the appropriate timescale of the above mechanism. In this type of coupling, the critical manifold of each oscillator (set of singular points of the fast dynamics) is invariant by the coupling parameter. Consequently, the coupling strength does not impact the critical manifold of the whole slow-fast coupled system either. Such kind of coupling has already been considered. In [START_REF] Clément | Mathematical modeling of the GnRH pulse and surge generator[END_REF][START_REF] Clément | Foliation-based parameter tuning in a model of the GnRH pulse and surge generator[END_REF][START_REF] Vidal | A dynamical model for the control of the gonadotrophin-releasing hormone neurosecretory system[END_REF][START_REF] Fernández-García | A multiple time scale coupling of piecewise linear oscillators. application to a neuroendocrine system[END_REF], the coupling is one-way between two planar FitzHugh-Nagumo systems to mimick the slow regulation of the GnRH secretory system by the interneuron network that conveys the steroid feedback. In [START_REF] Krupa | Mixed-mode oscillations in a multiple time scale phantom bursting system[END_REF], a rigorous study of signature transitions in Mixed-Mode Oscillations (MMOs) generated by such system has been performed. More recently, in [START_REF] Ersöz | Synchronization of weakly coupled canard oscillators[END_REF], a bidirectional coupling of this type between two planar systems was considered. In [START_REF] Ersöz | Canard-mediated (de) synchronization in coupled phantom bursters[END_REF][START_REF] Ersöz | Coupled multiple timescale dynamics in populations of endocrine neurons: Pulsatile and surge patterns of GnRH secretion[END_REF], a system is built using such type of coupling between three planar oscillators, two connected bidirectionally and a third one that acts as a regulator, connected one-way to the other two.

In this work, we consider another extension of coupled slow-fast systems using a bidirectional coupling between two 3D systems.

The new variable added to the classical FitzHugh-Nagumo system enables the emergence of MMOs in the orbits of the 3D oscillators, which makes the dynamics of the coupled six-dimensional (6D) system richer than the coupling of relaxation oscillators and harder to analyze. In this article, we propose an insight on the synchronization of coupled Mixed-Mode Oscillators, that is known to be a very difficult and yet fascinating question, from a perspective that we want to be global in two senses. First, we consider the global asymptotic behavior of the 6D system orbits and not only the local dynamics in the vicinity of the folded singularity. Second, we study the synchronization features for any coupling parameter values (both positive and negative). For medium positive values of the coupling parameter, instead of a quick in-phase synchronization pattern, we observe phase-locking almost-in-phase synchronization pattern, depending strongly on the strength of the coupling parameter. For the sake of brevity, we refer to this pattern as almost-in-phase synchronization from now on. Furthermore, for negative values, apart from antiphase synchronization patterns, the system also presents nonsymmetric behaviors: the relaxation loss of one oscillator and the total oscillation death. In this study, we describe and prove the distribution of these asymptotic behaviors according to the value of the coupling gain parameter. We also illustrate numerically and discuss the dynamical mechanisms involved in the transitions between behaviors.

The article is outlined as follows. In Section 2, we first recall the properties of the 3D slow-fast oscillator introduced in [START_REF] Krupa | A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons[END_REF], we then introduce the symmetric coupling between two such oscillators considered in this work and describe the main asymptotic patterns generated by the system. In Section 3, we describe the main asymptotic behaviors of the 6D system and prove their repartition according to the value of the coupling gain parameter. We then prove that when the coupling parameter c is negative, the oscillation frequency increases as c increases, that is, as |c| decreases. In Section 4, we show the complexity of the changes in the (quasi)-periodic behaviors for weak coupling. Basing ourselves on numerical simulations, we discuss the role of MMOs in organizing the dynamics for such coupling. Finally, in Section 5, we discuss possible extensions of the present work. The entire technical proofs of each result are available in Appendices.

Model properties

This section is divided into five subsections. We first recall the main properties of the 3D slow-fast oscillator introduced in [START_REF] Krupa | A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons[END_REF]. Then, we focus on canardinduced MMOs of the system. In the third subsection, we introduce the symmetric coupling between two identical 3D oscillators. We introduce the main asymptotic patterns generated by the 6D system. Using simulations, we introduce in the last subsection the repartition of the main behaviors that will be theoretically analyzed in the sequel.

Three-dimensional slow-fast oscillator

The 3D slow-fast system introduced in [START_REF] Krupa | A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons[END_REF] for modeling the dynamics of the ICC in a single neuron reads

ẋ = τ (-y + f (x) -φ f (z)), (1a) ẏ
= τ ε(x + a 1 y + a 2 ), (1b) ż
= τ ε φ r (x) - z -z b τ z , (1c)
where

φ f (z) = µz z + z 0 , (2)
φ r (x) = λ 1 + exp(-ρ(x -x on )) , (3)
where (x, y, z) ∈ R 3 and the overdot denotes the derivative with respect to time variable t. The timescale separation parameter fulfills 0 < ε 1. Function f is an odd cubic function such that y = f (x) is an S-shaped curve: f admits a local minimum at x = -x f < 0 and a local maximum at x = x f > 0.

Without loss of generality, we can (and, in the following, we will often) consider f (x) = -x 3 + 3x 2 f x. Parameter τ > 0 has been introduced in [START_REF] Krupa | A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons[END_REF] for the outputs to meet a given physical timescale and does not impact the phase portrait.

Moreover, we assume a 1 < 0, |a 1 | 1 (in a sense specified in the following) and parameters a 2 , z 0 , λ, τ z , z b to be strictly positive. We also consider a large enough ρ value such that the sigmoid φ r is steep at its inflection point and, therefore, represents a sharp activation function.

Subsystem (1a)-( 1b) is built on a classical FitzHugh-Nagumo system [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF],

where the fast variable x represents the electrical activity of the cell and the slow variable y is a recovery variable. Originally, the additional slow variable z represents the ICC in a single neuron. The dynamics of variable z is mainly driven by the increasing sigmoidal function φ r , and it acts as a feedback onto the x dynamics through the Hill function φ f bounded by µ. When φ r (x) is inactive (close to 0), variable z decreases to a quasi steady state close to z b which represents the baseline of the intracellular calcium level. The exponential decay rate of this motion is given by τ ε/τ z . The effect of this coupling is to reduce the electrical activity of the neuronal population in response to the rise of intracellular calcium concentration. An analogous term is used in models of single neurons to represent the hyperpolarization of the cell membrane stimulated by calcium [START_REF] Ermentrout | The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators[END_REF]. More generally, in this article, we consider system (1) fulfilling the above assumptions as a typical class of slow-fast system with 1 fast and 2 slow variables able to generate canard-induced MMOs as explained in the following.

In order to describe the behavior of the system, we specify the nullclines and possible stationary states.

• The x-nullcline is the critical manifold

S ≡ {y = f (x) -φ f (z)} . (4)
It is an S-shaped surface with two fold lines:

F -≡ {x = -x f , y = f (-x f ) -φ f (z)} , (5)
and

F + ≡ {x = x f , y = f (x f ) -φ f (z)} . (6)
The fold lines split the critical manifold into three sheets. The left and right sheets (contained in x < -x f and x > x f , respectively) are attracting for the fast dynamics, and we will note them S l and S r , respectively. The middle sheet S m (contained in |x| < x f ) is repulsive. From geometric singular perturbation theory [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], the part of these sheets away from the folds perturbs, for small enough positive values of ε, into slow manifolds S l,ε , S r,ε and S m,ε , invariant for system [START_REF] Fallani | Hierarchy of neural organization in the embryonic spinal cord: Granger-causality graph analysis of in vivo calcium imaging data[END_REF], and lying in a O(ε)-neighborhood of S l , S r and S m , respectively. Note that the O(ε) proximity between the sheets of the critical manifolds and the respective slow manifolds relies on the Hausdorff distance applied to the compact submanifolds of these sheets. The left and right slow manifolds are exponentially attractive w.r.t.

ε while the middle one is exponentially repulsive. Moreover, considering the extensions of S l,ε (resp. S r,ε ) by the positive flow, its first intersection

with x = -x f (resp. x = x f ) is a curve O(ε 2/3
)-close to the lower fold F -(resp., the upper fold F -). For more details on these results (that are now considered as classical), we refer the reader to [START_REF] Krupa | Mixed-mode oscillations in a multiple time scale phantom bursting system[END_REF][START_REF] Guckenheimer | Singular hopf bifurcation in systems with two slow variables[END_REF][START_REF] Krupa | Local analysis near a folded saddle-node singularity[END_REF][START_REF] Szmolyan | Relaxation oscillations in R 3[END_REF][START_REF] Szmolyan | Canards in R 3[END_REF][START_REF] Wechselberger | Existence and bifurcation of canards in R 3 in the case of a folded node[END_REF],

for instance.

• The y-nullcline is the plane

Π ≡ y = - x + a 2 a 1 . (7)
• The z-nullcline is the sigmoidal surface

Γ ≡ φ r (x) - z -z b τ z = 0 . (8)
In order to give the approximate location of the system equilibria, we consider the Heaviside approximation of the sigmoidal function φ r obtained for ρ → ∞. This leads to the approximation:

Γ    {z = z b } , if x < x on , {z = z b + λτ z =: z h } , if x > x on .
As in the original FitzHugh-Nagumo system, we take |a 1 | small enough to ensure that the x and y nullclines intersect only once for any fixed z value.

Thus, we define the following special curves lying on the critical manifold In the following, when necessary, we will mention the use of an approximation of the equilibrium point by the intersection between S b or S h and plane Π. In particular, such approximations allow us to greatly simplify the illustrations of the nullclines positions in a 2D view (see Fig. 1). Moreover, note that the For the simulations shown in this article, we have used f (x) = -x 3 + 4x, and the following parameter values:

S b ≡ {y = f (x) -φ f (z b ), x < x on , z = z b } , (9)
S h ≡ {y = f (x) -φ f (z h ), x > x on , z = z h } . ( 10
) (a) (b) S ⇧ -2 -1 1 2 x -6 -4 -2 2 4 6 y ⇧ x f x f x on S b S h z = z b z = z h
projection of S l ∩ {z b ≤ z ≤ z h } is
a 1 = -0.1, a 2 = 0.8, ε = 0.06, τ = 37, z 0 = 5, z b = 1, τ z = 2, x on = -0.45, λ = 1.75, ρ = 30. (11)

Mixed-Mode Oscillations of the 3D system

In Appendix A, we prove that for µ = 2.51, system (1) admits a folded saddle-node type II (FSNII) singularity in the singular limit, which implies that 175 a singular Hopf bifurcation occurs within a O(ε)-distance in the parameter space and that for µ < 2.51, system (1) admits a folded node (FN) singularity. It is well-known that in the unfolding of the FN and the FSNII singularities, systems generate canard-induced MMOs [START_REF] Guckenheimer | Singular hopf bifurcation in systems with two slow variables[END_REF][START_REF] Krupa | Local analysis near a folded saddle-node singularity[END_REF][START_REF] Szmolyan | Canards in R 3[END_REF][START_REF] Brøns | Mixed mode oscillations due to the generalized canard phenomenon[END_REF][START_REF] Guckenheimer | Unfoldings of singular hopf bifurcation[END_REF][START_REF] Wechselberger | Existence and bifurcation of canards in R 3 in the case of a folded node[END_REF]. In particular, in [START_REF] Wechselberger | Existence and bifurcation of canards in R 3 in the case of a folded node[END_REF], the author have proven the existence of a "funnel", i.e. a partition of the local phase space by "sectors of rotation", that organizes the local dynamics near FN type singularities. The sectors are separated by secondary canards corresponding to intersections between the attractive left and repulsive middle slow manifolds (S l,ε and S m,ε , respectively). Then, the number of small oscillations featured by a local orbit depends on the specific sector of rotation of the associated funnel that it enters. Together with the global return mechanism provided by the cubic term in the x-nullcline, the existence of such folded singularity allows system

(1) to generate MMOs. We refer the reader to the seminal article [START_REF] Wechselberger | Existence and bifurcation of canards in R 3 in the case of a folded node[END_REF] and to the review [START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF] and references therein for more details about canard-induced MMOs.

We recall the three qualitatively different behaviors depending on the position of the equilibrium point, since we will refer to their features in the subsequent analysis of the coupled system. We recall that the following description results from a slow-fast analysis and thus we consider ε > 0 small enough, which means, rigorously, that there exists a fixed value ε 0 > 0 such that the statements are fulfilled for any ε ∈]0, ε 0 ] (see [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF]).

The three cases described below refer to the panels displayed in Fig. 2.

(a) The equilibrium lies on the middle sheet of the critical manifold S and far from the folds: the equilibrium is unstable and the system admits a globally attractive limit cycle of relaxation type.

(b) The equilibrium lies on the middle sheet S m of S and close to a fold: the equilibrium is unstable and the system generates MMOs.

(c) The equilibrium lies on S l : the equilibrium is stable and the system asymptotically reaches a steady state. A similar situation occurs if the equilibrium lies on S r (not shown in the figure).

It is worth describing the path of x along a prototypical MMO limit cycle (case (b)) of system [START_REF] Fallani | Hierarchy of neural organization in the embryonic spinal cord: Granger-causality graph analysis of in vivo calcium imaging data[END_REF], in order to better understand the behavior of the

(a) (b) (c) 2 4 6 -4 -2 0 2 4 -2 -1 0 1 2
3D model, the periodic exploration by the current point of the regions corresponding to oscillatory regime and excitable regime of subsystem (3a)-(3b) may produce Mixed-Mode Oscillations. We will use this feature to reproduce the quiescent phase in the generated calcium pattern. System (3) is a slow-fast system with one fast and two slow variables. To describe the dynamical mechanisms of this system, we introduce the following notations. The critical manifold S 0 (or xnullcline) given by :

y = 4x x 3 IAHP (Ca) (5)
is an S-shaped surface displaying two fold lines F (lying in x < 0) and F + (lying in x > 0) that split S 0 into three parts : the left and right ones lying entirely in x < 0 and x > 0 respectively and the middle one. The y-nullcline defined by :

a0x + a1y + a2 (6)
is a plane that crosses F for a given value Ca f of Ca. Finally, the Ca-nullcline is defined by a relation between x and Ca : Ca = ⇥CaICa(x) + Ca bas [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF] that is, like ICa, a smooth sigmoidal function. This surface is attractive for the Ca dynamics.

Let us describe the typical interactions between the state variables, starting from a low level of Ca (i.e. close to Ca bas ) and (x, y) such that (x, y, Ca) lies just below F . Under the influence of the fast dynamics, the current point (x, y, Ca) quickly reaches the right branch of S 0 , x and ⇥CaICa(x) quickly increases. Consequently Ca increases while the current point goes up along the right branch of S 0 towards F + . Then, once the current point has arrived above F + , under the influence of the fast dynamics, it quickly comes back near the left branch of S 0 , x quickly decreases as well as ⇥CaICa(x) (which becomes almost null). The current point, driven by the slow dynamics, goes down along the left branch of S 0 and Ca decreases eventually down to Ca bas . Then several situation may occur depending mainly on the value of µ and related to the regime of system (3a)-(3b):

1. For small values of µ, when the current point reaches the vicinity of F , system (3a)-(3b) is in oscillatory regime. As a consequence, the current point directly and quickly reaches the right branch of S 0 , and the behavior described above repeats immediately. 9

2. For an interval of values of µ, system (3a)-(3b) is in excitable regime when Ca approaches Ca bas . Then (x, y) reaches the vicinity of the singular point of (3a)-(3b) close to the left knee.

As Ca keeps decreasing until the current point is very close to the attractive surface given by [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF]. Consequently system (3a)-(3b) passes into oscillatory regime. During this passage, the current point makes small oscillations around the fold F before the fast transition to the right branch and the whole motion repeats.

3. For great value of µ, system (3a)-(3b) remains stuck in the steady regime. Hence, after one pulse, the level of Ca remains close to the baseline.

In the following, we are specially interested in the second case since the generated pattern of Ca fits quite well the qualitative specifications of the experimental time series obtained by Terasawa in [START_REF] Bélair | On linearly coupled relaxation oscillations[END_REF], especially the quiescent phase at the baseline level between two successive pulses. Hence, we will chose the value of µ in the interval where the orbit of system (3) display MMO. It is worth noticing that this interval of values depends on the other parameters, particularly on the time scale parameters ⇧ and µ. To fix the idea, we consider the following set of parameter: a0 = 1, a1 = 0.1, a2 = 0.8, ⇧ = 0.06, µ = 2.4, Ca0 = 500, Ca bas = 100, ⌅Ca = 2, = 175, ⇤Ca = 4.5,

x on = 0.45.

(

) 8
With this values, the generated calcium pattern fulfills the typical qualitative behavior described in 3.1 and the average quantitative specifications: 10 min of IPI, peak level at 342 nM (red signal in all panels of Figure ??. Now, we will show how to mimic the variability of the quantitative features of calcium patterns between di erent cells by choosing di erent values for parameters of special importance: µ and k. Let us note that, from the above explanation, one can already see that the precise value of µ prescribes the number of small oscillations of the current point near the left fold F and, consequently, the duration of the quiescent phase. Since variations in µ does not impact much the duration of the pulses, one can consider that this parameter controls marginally the InterPulse Interval (IPI). Panel A of Figure ?? shows the results of change in µ: the red signal is obtained with parameter values (8) (IPI=10 min), the blue one with µ = 2.32 (IPI=7 min) and the yellow one with µ = 2.42 (IPI=14 min). It is worth noticing that the interval of values of µ for which a quiescent phase appear between two successive pulses in the calcium pattern depends on the other parameters.

Parameter k essentially tunes the time scale separation between y and Ca (x being much more faster). Hence, an increase in y implies a shorter time for subsystem (3a)-(3b) to complete a relaxation oscillation and, consequently, a shorter time for Ca to increase and decrease back to the baseline. In conclusion, one can increase or decrease the calcium peak level by tuning the value of parameter k. Of course, a variation in parameter k value also implies a change in the quiescence phase duration and, thus, the IPI. Panel B of Figure ?? shows, together with the initial red signal (peak level at 342 nM, IPI=10min), a blue signal obtained with k = 1.2 (peak level at 320 nM, IPI=4 min) and a yellow one obtained with k = 0.8 (peak level at 365 nM, IPI=16 min). Of course, one can choose both the peak level and the IPI by tuning first the value of k and afterwards the value of µ. Panel C of Figure Figure ?? shows calcium patterns obtained with the same set of k values as in panel B but the di erences in the IPI as been corrected by choosing appropriate values of µ : 2.448 for the blue signal, 2.238 for the yellow one. of µ, system (3a)-(3b) is in excitable regime when Ca approaches s the vicinity of the singular point of (3a)-(3b) close to the left knee. ntil the current point is very close to the attractive surface given by (3a)-(3b) passes into oscillatory regime. During this passage, the ll oscillations around the fold F before the fast transition to the le motion repeats. em (3a)-(3b) remains stuck in the steady regime. Hence, after one ains close to the baseline. cially interested in the second case since the generated pattern of Ca pecifications of the experimental time series obtained by Terasawa hase at the baseline level between two successive pulses. Hence, we e interval where the orbit of system (3) display MMO. It is worth lues depends on the other parameters, particularly on the time scale idea, we consider the following set of parameter:

1, a1 = 0. l, the periodic exploration by the current point of the regions corresponding to oscillatory nd excitable regime of subsystem (3a)-(3b) may produce Mixed-Mode Oscillations. We his feature to reproduce the quiescent phase in the generated calcium pattern. m (3) is a slow-fast system with one fast and two slow variables. To describe the dynamical ms of this system, we introduce the following notations. The critical manifold S 0 (or xgiven by : y = 4x x 3 IAHP (Ca) [START_REF] Bélair | On linearly coupled relaxation oscillations[END_REF] aped surface displaying two fold lines F (lying in x < 0) and F + (lying in x > 0) that nto three parts : the left and right ones lying entirely in x < 0 and x > 0 respectively and le one. The y-nullcline defined by :

a0x + a1y + a2 (6)
e that crosses F for a given value Ca f of Ca. Finally, the Ca-nullcline is defined by a etween x and Ca : Ca = ⇥CaICa(x) + Ca bas [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF] ke ICa, a smooth sigmoidal function. This surface is attractive for the Ca dynamics. s describe the typical interactions between the state variables, starting from a low level . close to Ca bas ) and (x, y) such that (x, y, Ca) lies just below F . Under the influence st dynamics, the current point (x, y, Ca) quickly reaches the right branch of S 0 , x and ) quickly increases. Consequently Ca increases while the current point goes up along branch of S 0 towards F + . Then, once the current point has arrived above F + , under nce of the fast dynamics, it quickly comes back near the left branch of S 0 , x quickly as well as ⇥CaICa(x) (which becomes almost null). The current point, driven by the slow , goes down along the left branch of S 0 and Ca decreases eventually down to Ca bas . Then tuation may occur depending mainly on the value of µ and related to the regime of system : small values of µ, when the current point reaches the vicinity of F , system (3a)-(3b) is scillatory regime. As a consequence, the current point directly and quickly reaches the t branch of S 0 , and the behavior described above repeats immediately. 2. For an interval of values of µ, system (3a)-(3b) is in excitable regime when Ca approaches Ca bas . Then (x, y) reaches the vicinity of the singular point of (3a)-(3b) close to the left knee.

As Ca keeps decreasing until the current point is very close to the attractive surface given by [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF]. Consequently system (3a)-(3b) passes into oscillatory regime. During this passage, the current point makes small oscillations around the fold F before the fast transition to the right branch and the whole motion repeats.

3. For great value of µ, system (3a)-(3b) remains stuck in the steady regime. Hence, after one pulse, the level of Ca remains close to the baseline.

In the following, we are specially interested in the second case since the generated pattern of Ca fits quite well the qualitative specifications of the experimental time series obtained by Terasawa in [START_REF] Bélair | On linearly coupled relaxation oscillations[END_REF], especially the quiescent phase at the baseline level between two successive pulses. Hence, we will chose the value of µ in the interval where the orbit of system (3) display MMO. It is worth noticing that this interval of values depends on the other parameters, particularly on the time scale parameters ⇧ and µ. To fix the idea, we consider the following set of parameter:

a0 = 1, a1 = 0.1, a2 = 0.8, ⇧ = 0.06, µ = 2.4, Ca0 = 500, Ca bas = 100, ⌅Ca = 2, = 175, ⇤Ca = 4.5,
x on = 0.45.

(

) 8
With this values, the generated calcium pattern fulfills the typical qualitative behavior described in 3.1 and the average quantitative specifications: 10 min of IPI, peak level at 342 nM (red signal in all panels of Figure ??. Now, we will show how to mimic the variability of the quantitative features of calcium patterns between di erent cells by choosing di erent values for parameters of special importance: µ and k. Let us note that, from the above explanation, one can already see that the precise value of µ prescribes the number of small oscillations of the current point near the left fold F and, consequently, the duration of the quiescent phase. Since variations in µ does not impact much the duration of the pulses, one can consider that this parameter controls marginally the InterPulse Interval (IPI). Panel A of Figure ?? shows the results of change in µ: the red signal is obtained with parameter values (8) (IPI=10 min), the blue one with µ = 2.32 (IPI=7 min) and the yellow one with µ = 2.42 (IPI=14 min). It is worth noticing that the interval of values of µ for which a quiescent phase appear between two successive pulses in the calcium pattern depends on the other parameters.

Parameter k essentially tunes the time scale separation between y and Ca (x being much more faster). Hence, an increase in y implies a shorter time for subsystem (3a)-(3b) to complete a relaxation oscillation and, consequently, a shorter time for Ca to increase and decrease back to the baseline. In conclusion, one can increase or decrease the calcium peak level by tuning the value of parameter k. Of course, a variation in parameter k value also implies a change in the quiescence phase duration and, thus, the IPI. Panel B of Figure ?? shows, together with the initial red signal (peak level at 342 nM, IPI=10min), a blue signal obtained with k = 1.2 (peak level at 320 nM, IPI=4 min) and a yellow one obtained with k = 0.8 (peak level at 365 nM, IPI=16 min). Of course, one can choose both the peak level and the IPI by tuning first the value of k and afterwards the value of µ. Panel C of Figure Figure ?? shows calcium patterns obtained with the same set of k values as in panel B but the di erences in the IPI as been corrected by choosing appropriate values of µ : 2.448 for the blue signal, 2.238 for the yellow one. 10 3D model, the periodic exploration by the current point of the regions corresponding to oscillatory regime and excitable regime of subsystem (3a)-(3b) may produce Mixed-Mode Oscillations. We will use this feature to reproduce the quiescent phase in the generated calcium pattern.

System ( 3) is a slow-fast system with one fast and two slow variables. To describe the dynamical mechanisms of this system, we introduce the following notations. The critical manifold S 0 (or xnullcline) given by :

y = 4x x 3 IAHP (Ca) (5)
is an S-shaped surface displaying two fold lines F (lying in x < 0) and F + (lying in x > 0) that split S 0 into three parts : the left and right ones lying entirely in x < 0 and x > 0 respectively and the middle one. The y-nullcline defined by :

a0x + a1y + a2 (6)
is a plane that crosses F for a given value Ca f of Ca. Finally, the Ca-nullcline is defined by a relation between x and Ca : Ca = ⇥CaICa(x) + Ca bas [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF] that is, like ICa, a smooth sigmoidal function. This surface is attractive for the Ca dynamics.

Let us describe the typical interactions between the state variables, starting from a low level of Ca (i.e. close to Ca bas ) and (x, y) such that (x, y, Ca) lies just below F . Under the influence of the fast dynamics, the current point (x, y, Ca) quickly reaches the right branch of S 0 , x and ⇥CaICa(x) quickly increases. Consequently Ca increases while the current point goes up along the right branch of S 0 towards F + . Then, once the current point has arrived above F + , under the influence of the fast dynamics, it quickly comes back near the left branch of S 0 , x quickly decreases as well as ⇥CaICa(x) (which becomes almost null). The current point, driven by the slow dynamics, goes down along the left branch of S 0 and Ca decreases eventually down to Ca bas . Then several situation may occur depending mainly on the value of µ and related to the regime of system (3a)-(3b):

1. For small values of µ, when the current point reaches the vicinity of F , system (3a)-(3b) is in oscillatory regime. As a consequence, the current point directly and quickly reaches the right branch of S 0 , and the behavior described above repeats immediately. coupled system later. We note p = (x, ŷ, ẑ) with x close to -x f the unique equilibrium point of system (1), which is unstable; we denote by x min (resp.

x max ) the minimal (resp. maximal) value of x along the limit cycle and P = (x, y, z) the current point. Recall that the folds F -and F + of the critical 210 manifold S given by eq. ( 4) are located at x = -x f and x = x f (see eqs.

(5)-( 6)). Also, we consider the special curves S b and S h given by eqs. ( 9)- [START_REF] Storti | Dynamics of two strongly coupled relaxation oscillators[END_REF] as the approximation of the intersections between the z-nullcline and the right and left sheets of the critical manifold, respectively.

Consider the projection F - proj of F -along the x-direction (fast fibers) onto 215 the opposite branch S r of the critical manifold (see Fig. 3). The x-component of the intersection between F - proj and the special curve

S h is the unique positive solution of f (x) -φ f (z h ) = f (-x f ) -φ f (z b
), that we denote x M . Analogously, consider the projection F + proj of F + , along the x-direction (fast fibers) onto the opposite branch S l of the critical manifold. The x-component of the intersection between F + proj and the special curve S b is the unique negative solution of

f (x) - φ f (z b ) = f (x f ) -φ f (z h ), that we denote x m .
Now along the limit cycle, starting close to S r,ε , which is attractive, and with z b < z < z h , ẏ > 0 and ż > 0, it holds that: the current point P climbs up along S r,ε until it reaches F + vicinity and z remains below z h . From classical

225 2 4 6 -4 -2 0 2 4 -2 -1 0 1 2
into three parts : the left and right ones lying entirely in x < 0 and x > 0 respectively and dle one. The y-nullcline defined by :

a 0 x + a 1 y + a 2 (6)
ne that crosses F for a given value Ca f of Ca. Finally, the Ca-nullcline is defined by a between x and Ca : Ca = ⇥ Ca I Ca (x) + Ca bas [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF] like I Ca , a smooth sigmoidal function. This surface is attractive for the Ca dynamics. us describe the typical interactions between the state variables, starting from a low level i.e. close to Ca bas ) and (x, y) such that (x, y, Ca) lies just below F . Under the influence ast dynamics, the current point (x, y, Ca) quickly reaches the right branch of S 0 , x and x) quickly increases. Consequently Ca increases while the current point goes up along t branch of S 0 towards F + . Then, once the current point has arrived above F + , under uence of the fast dynamics, it quickly comes back near the left branch of S 0 , x quickly es as well as ⇥ Ca I Ca (x) (which becomes almost null). The current point, driven by the slow cs, goes down along the left branch of S 0 and Ca decreases eventually down to Ca bas . Then situation may occur depending mainly on the value of µ and related to the regime of system ):

r small values of µ, when the current point reaches the vicinity of F , system (3a)-(3b) is oscillatory regime. As a consequence, the current point directly and quickly reaches the ght branch of S 0 , and the behavior described above repeats immediately. 9 all oscillations around the fold F before the fast transition to the ole motion repeats. tem (3a)-(3b) remains stuck in the steady regime. Hence, after one mains close to the baseline. cially interested in the second case since the generated pattern of Ca specifications of the experimental time series obtained by Terasawa phase at the baseline level between two successive pulses. Hence, we he interval where the orbit of system (3) display MMO. It is worth alues depends on the other parameters, particularly on the time scale e idea, we consider the following set of parameter:

1, a 1 = 0.1, a 2 = 0.8, ⇧ = 0.06, 2.4, Ca 0 = 500, Ca bas = 100, ⌅ Ca = 2, 175, ⇤ Ca = 4.5,
x on = 0.45.

(

) 8
ed calcium pattern fulfills the typical qualitative behavior described ative specifications: 10 min of IPI, peak level at 342 nM (red signal mimic the variability of the quantitative features of calcium patterns osing di erent values for parameters of special importance: µ and he above explanation, one can already see that the precise value f small oscillations of the current point near the left fold F and, the quiescent phase. Since variations in µ does not impact much can consider that this parameter controls marginally the InterPulse ure ?? shows the results of change in µ: the red signal is obtained I=10 min), the blue one with µ = 2.32 (IPI=7 min) and the yellow in). It is worth noticing that the interval of values of µ for which a en two successive pulses in the calcium pattern depends on the other

S h = y = f (x) -φ f (z b ), x < xon, z = z b (resp. F + proj and the special curve S b = y = f (x) -φ f (z h ), x > xon, z = z h .
results on fold-induced slow-fast transitions [START_REF] Szmolyan | Relaxation oscillations in R 3[END_REF], S r,ε crosses at O(ε 2/3 ) vicinity of F + and, under the impact of the fast dynamics, reaches S l O(ε 2/3 )-close to F + proj . Once P crosses the critical manifold, ẋ > 0 and consequently x min ≈ x m . Then, ẏ < 0 and ż < 0: P follows S l,ε towards the lower fold. The equilibrium point p is both unstable and ε-close to the lower fold. Yet, recall that we 230 consider the limit cycle to enter the funnel of the folded singularity away from any secondary canard. Hence, despite the fact that small oscillations occur near F -, a similar fast motion occurs from the fold vicinity towards the opposite branch: P reaches S r ε-close to F - proj and, since ẋ < 0 once P has crossed S r , we obtain x max ≈ x M .

235

Hence, one period of an MMO limit cycle can be divided into five phases.

In the following description, we focus on the path of variable x: I. Silent phase: P follows S l,ε exponentially close (w.r.t. ε) while x increases from x min to -x fδ ε with δ ε is of order O(ε) and positive. II. Passage through the funnel: P undergoes small oscillations during the 240 passage through the funnel near the fold F -and x ∈ (-x f -δ ε , -x f +δ ε ). III. Jump-up: P follows an almost horizontal (y and z almost constant) fast fiber to reach S r,ε vicinity and x ∈ (-x f + δ ε , x max ), IV. Active phase: P climbs up exponentially close (w.r.t. ε) to S r,ε up to the upper fold F + while x decreases from x max to x f . V. Jump-down: P follows a fast fiber to reach S l,ε vicinity and x decreases from x f to x min .

We will refer to Phases I to V of a MMO limit cycle in the whole article.

Recall that, even if all orbits (except orbits reduced to a singular point of the system) eventually adopt the above behavior, chaotic dynamics may arise for certain range of parameter values. Such a situation is characterized as follows (see [START_REF] Brøns | Mixed mode oscillations due to the generalized canard phenomenon[END_REF][START_REF] Szmolyan | Canards in R 3[END_REF][START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF] for instance). Consider the global return map π (Poincaré map), defined from a section Σ transverse to S l away from the folded singularity into itself: the image of P ∈ Σ by this map is defined as the first intersection point between Σ and the positive trajectory starting from P . Note that the trajectory between P and π(P ) necessarily features a relaxation oscillation or at least a sufficiently large excursion in the phase space. If π admits a single fixed point, it is associated with an attractive limit cycle, featuring a single relaxation and n small oscillations along one period (signature 1 n ): such cycle stays sufficiently far away from the secondary canards and spirals around the pseudo-singularity following one sector of rotation associated with n small oscillations. When the unicity of the fixed point is lost by varying a parameter value, several periodic orbits appear. In particular, the existence of multiple fixed points for π close to the intersection with a secondary canard reveals the non existence of limit cycle and the existence of quasi-periodic orbits or non isolated periodic orbits with complex signature 1 s1 1 s2 ...1 sp , where s i ∈ {n, n + 1}). In such case, along any orbit, the current point passes recurrently close to a secondary canard, on one side or the other at each passage, and therefore generates n or n + 1 small oscillations.

For the sake of simplicity, in the following, we consider ranges of parameter values such that the 3D oscillator (1) admits an attractive MMO limit cycle featuring small oscillations near the lower fold F -, which we state as hypothesis in the whole article. We consider two identical 3D systems of type [START_REF] Fallani | Hierarchy of neural organization in the embryonic spinal cord: Granger-causality graph analysis of in vivo calcium imaging data[END_REF], and aim at understanding how they interact when they are symmetrically coupled.

Symmetric coupling of two identical systems

We consider a symmetric linear coupling through the fast variable in the slow dynamics of the recovery variable between two identical oscillators:

O 1            ẋ = τ (-y + f (x) -φ f (z)), ẏ = τ ε(x + a 1 y + a 2 + c(x -X)), ż = τ ε φ r (x) - z -z b τ z , O 2            Ẋ = τ (-Y + f (X) -φ f (Z)), Ẏ = τ ε(X + a 1 Y + a 2 + c(X -x)), Ż = τ ε φ r (X) - Z -z b τ z ,                                (12)
with c either positive or negative. For the sake of simplicity, in the following, we will refer as O 1 and O 2 to each oscillator undergoing the impact of the fast variable of the other one. Hence O 1 is the X-dependent subsystem of variables (x, y, z) and O 2 is the x-dependent subsystem of variables (X, Y, Z). Remark that the coupling is not involved in the fast dynamics and, consequently, the critical manifolds of O 1 and O 2 are defined by the same equations in (x, y, z) and (X, Y, Z) spaces, respectively. Hence, we keep in mind that these manifolds do not belong to the same space, but that they are the exact same geometric object S as parameterized in eq. ( 4). In the following, for simplifying the notations as long as the context does not lead to ambiguous interpretation, we will equally refer to the critical manifold of O 1 and O 2 as S and to the left and right sheets as S l and S r .

First, remark that if we consider the same initial conditions for both oscillators (x 0 , y 0 , z 0 ) = (X 0 , Y 0 , Z 0 ), the coupling terms do not affect the dynamics:

x = X for any time and consequently y = Y and z = Z for any time as well.

Along such orbit, both oscillators O 1 and O 2 behave exactly in the same way.

Moreover, both admit the same unstable equilibrium point close to the lower fold of their critical manifold and the same MMO limit cycle. Consequently, the two oscillators are in-phase synchronized and generate MMOs. The 3D subspace {x = X, y = Y, z = Z} is therefore an invariant space where there always exist an equilibrium point of the 6D system (12) and a MMO limit cycle that is attractive inside this subspace.

If we consider different initial conditions, the coupling terms play different roles depending on the sign and strength of parameter c, resulting in different dynamical behaviors that we analyze in the sequel. Hence, for the sake of simplicity, we will only refer to the asymptotic behavior (ω-limit) of the orbits of system [START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF] with any initial data that is not a singular point and outside {x = X, y = Y, z = Z}, keeping in mind that, in this invariant subspace, the sole asymptotic behavior is obtained as the same MMO limit cycle for both oscillators with no phase-shift between them.

Main asymptotic behaviors of the coupled system

We describe the main asymptotic behaviors, i.e., those occurring for large intervals of c values. In Panel (0) of Fig. 4, we represent the time traces of (x, X) and (z, Z) variables along an orbit of system [START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF] of two identical uncoupled (c = 0) oscillators in MMO regime with different initial conditions.

If c = 0, system [START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF] with different initial conditions for O 1 and O 2 can generate different patterns depending on the value of the coupling parameter c. Note that we have produced the simulations presented in this section with µ = 2.4, but one can consider equally neighboring value of µ close to 2.51 (value obtained in the Appendix A) for which system (12) admits an attractive MMO limit cycle, and obtain the same dynamical cases yet for different ranges of the coupling parameter c. We represent an instance of each case in Fig. 4 and we explain each asymptotic behavior below.

(a) Total oscillation death: both oscillators reach a stable equilibrium.

(b) Relaxation loss: one of the oscillators produces relaxation oscillations between S l and S r , while the other oscillates confined near S l .

(0) Uncoupled case, c = 0.

(a)

(x c (t), y c (t), z c (t), X c (t), Y c (t), Z c (t)
) the solution of the coupled system for a fixed value of c starting from an initial condition that does not belong to the subspace {x = X, y = Y, z = Z}. We have simulated the solutions for each value of c on a 10 -3 step grid of [-1, 1] as long as necessary:

• either for the current point to enter the 10 -3 -radius ball defined by the euclidian distance centered at the stable equilibrium point (case of total oscillation death, panel (a) in Fig. 4),

• or for the orbit to cross {x = 0, y < 0} or {X = 0, Y < 0} a hundred times (oscillatory cases, panel (b) to (e) in Fig. 4).

For each value of c in the grid corresponding to the second case above, we have kept the 26 last crossing times, denoted (t x=0 i (c)) 25 i=0 and (t X=0 i (c)) 25 i=0 in the following. When applicable, we have calculated the successive return times:

T x=0 i (c) = t x=0 i (c) -t x=0 i-1 (c), (13)
T X=0 i (c) = t X=0 i (c) -t X=0 i-1 (c), (14)
with i ∈ [1, 25]∩N, and the non-oriented phase shifts between the two oscillators

φ i (c) = |t x=0 i (c) -t X=0 i (c)|. ( 15
)
Illustrative construction of crossing values, return times and shifts on an instance of x(t) and X(t) signals along an orbit of system ( 12) is shown in Fig. 5.

... ... ... ... ...

X(t)

x(t) i=1 (red points). An animated view of the signals generated by x and X for various c values together with the picture of Fig. 6 is provided as supplemental material (see animated Fig. 7).

t x=0 i t x=0 i+1 t x=0 i 1 t X=0 i+1 t X=0 i t X=0 i 1 T X=0 i T X=0 i+1 T x=0 i+1 T x=0 i i 1 i+1 i t - 2
From this simulation, we identify five different intervals of c values labeled from A to E and corresponding to different general behaviors. Before specifying these behaviors, we emphasize an important point: while intervals D and E are separated by a single value c p , it is not the case between A and B, B and C, C and D. Precisely, complex changes in the phase portrait occur for c values in small intervals around the specific value c d , c 0 and 0, materialized in Figure 6.

In this subsection, we mainly focus on the asymptotic behaviors displayed by the system outside small neighborhoods of these values. The changes occurring around c d and c 0 are partly explained in the analysis presented in subsection limit cycle). Since interval B is very narrow, in order to emphasize the case of relaxation loss of one oscillator, we give as supplemental material an animated view of the signals generated by x and X for c values around this interval (see animated Fig. 8).

405

(Loading Video...) In intervals C and E of c values, the sequences of return times converge for both oscillators: the 25 green and 25 blue points materializing the successive return times for each oscillators cannot be distinguished. Moreover, the limit is the same for both oscillators. Hence, the sets of blue and green points coincide and the blue points are hidden behind the green ones. These cases correspond to cases (c) and (e) of Fig. 4, respectively. It is also the case for subintervals of D for which the oscillators are perfectly synchronized (case (e) in Fig. 4) and generate a MMO orbit. A discrepancy between the return times occur for the other c values in interval D corresponding to almost-in-phase synchronization of MMO patterns (case (d) in Fig. 4) either periodic with a constant phase shift between two successive relaxation or more complex with numerous values of the phase-shifts (red dots). Above interval C, we have also plotted the half of the return time (dashed black line), which coincides precisely with the time shift between oscillators. This fact shows that the oscillators are in perfect antiphase for c values in interval C.

Note that for c > 0 the activity phases of the oscillators tend to facilitate each other, while for c < 0 the activation of one oscillator blocks the possibility of activation for the other one. Thus, we can naturally interpret c > 0 as mutual excitatory case, while c < 0 can be interpreted as mutual inhibitory case.

In the sequel, we formalize theoretically the repartition shown in this simulation and present dynamical explanations of the occurrence and properties of each behavior according to c value.

Theoretical analysis of the main behaviors repartition

The section is divided into three subsections. In subsection 3.1, we present a preliminary result about the impact of one oscillator onto the behavior of the other one. In subsection 3.2, we establish the main result of this work, where we obtain macroscopic intervals where the asymptotic behavior of system [START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF] is asymptotically stable. Finally, in subsection 3.3, we focus on the dependency of an asymptotic orbit period on the coupling gain parameter when the system displays antiphase synchronization. In particular, we prove that the frequency increases with the (negative) parameter c value.

Analysis of the impact of one 3D oscillator state upon the other one

In the following proposition, which is a preliminary result in order to prove the main result of this work in the next subsection, we discuss the effect of X on subsystem O 1 , i.e., the impact of oscillator O 2 state upon oscillator O 1 transitory dynamics. Due to the symmetry between O 1 and O 2 in the 6D dynamics, the reverse situation is exactly the same one. When c = 0, the new y-nullcline is given by

Π c ≡ y = -(1 + c)x -a 2 + cX a 1 =: r(x, X; c) , (16)
and then, parameter c changes the y-nullcline in two different ways:

1. it changes the slope of plane Π c ;

2. it changes the position of plane Π c , depending on the value of X.

In the sequel, we will refer to an X-dependent stable singular point of O 1 as a "stable quasi-equilibrium", since it is moving with the value of X. We use this notion to describe the motion of O 1 current point "following a stable quasi-equilibrium point" during a phase where it cannot experience a slow-fast transition due to the position of X-dependent Π c . We use such notion since it enables a much easier description than discussing all the possible changes in the sign of the x, y and z dynamics that are not required for depicting O 1 qualitative behaviors.

Subsequently, we include a sketch of Proposition 1 proof, where we depict the special curves of the critical manifold S b and S h given in eqs. ( 9) and ( 10), respectively, and the X-dependent position of plane Π c while X travels along phases I and IV of the limit cycle γ µ,ε (t) of the uncoupled case. We emphasize that, in the figures, we represent these objects in the same (x, y)-plane, having in mind that they are located in the plane z = z b for x < x on and in the plane z = z h for x > x on . For the sake of brevity, the remainder of the proof has been relegated to Appendix B. We recall that Phases I-IV have been defined at the end of subsection 2.1.

Proposition 1. Consider system [START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF]. Assume that, for c = 0, the values of the remainder parameters of the system are fixed so that each oscillator O 1

and O 2 admits the attractive MMO limit cycle γ µ,ε (t). There exists 0 < δ 1

(precisely δ of the same order as |a 1 |) so that, for |c| > δ, while O 2 current point

P 2 = (X, Y, Z) lies close to S l ∩ {X ∈ [x min , -x f ]} or S r ∩ {X ∈ [x f , x max ]},
the motion of O 1 current point P 1 = (x, y, z) under the impact of the coupling term c(x -X) in the y-dynamics is the following.

1. Consider c > δ.

(a) While X increases from x min to -x f (O 2 in phase I, P 2 is ε-close to S l ), P 1 follows a stable quasi-equilibrium point lying on S l and moving rightwards. When this equilibrium passes through the lower fold, P 1 generates small oscillations and then reaches quickly S r vicinity.

(b) While X decreases from x max to x f (O 2 in phase IV, P 2 is ε-close to S r ), P 1 generates relaxation oscillations.

2. For c < -δ.

(a) While X increases from x min to -x f , P 1 either oscillates or follows a stable quasi-equilibrium point lying on S r .

(b) While X decreases from x max to x f , P 1 follows a stable equilibrium point lying on S l .

Sketch of the proof

For each oscillator, the x-nullcline is the critical manifold S given in eq. ( 4)

and the z-nullcline is the sigmoidal surface Γ given in eq. ( 8). The y-nullcline is Π given in eq. ( 7) in the uncoupled case, and Π c given in eq. ( 16) in the coupled case with coupling gain c. First, note that

Π ∩ Π c = X, - X + a 2 a 1 , Z X ∈ R, Z ∈ R and in particular, for X = -x f , Π ∩ Π c ∩ {X = -x f } = -x f , x f -a 2 a 1 , Z Z ∈ R . (17)
There are four different possibilities:

1. If c > 0, (a) X ∈ (x min , -x f ). (b) X ∈ (x f , x max ). 2. If c < 0, (a) X ∈ (x min , -x f ). (b) X ∈ (x f , x max ).
We represent in Fig. 9 the dynamical configurations of the four cases, that give rise to the statements of the proposition. In particular, we depict the special curves of the critical manifold S b and S h given in eqs. ( 9) and ( 10), respectively, and the X-dependent position of plane Π c while X travels along phases I and IV of the limit cycle γ µ,ε (t) of the uncoupled case. The technical details that allow to finish the proof are included in Appendix B.

At this point, we are ready to study the different synchronization patterns generated by the interaction between both oscillators.

Repartition of the main behaviors according to the coupling gain.

We take advantage of the analysis in the preceding subsection for studying the synchronization patterns generated by the interaction between the coupled oscillators. The following theorem formalizes the macroscopic structure of the behavior repartition w.r.t. parameter c shown in Fig. 6 and the observations made in subsection 2.5.

Subsequently, we include a sketch of the proof of the Theorem. The most technical part of the proof has been relegated to Appendix C.

Theorem 3.1. Consider system [START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF]. Assume that when c = 0 the remainder parameters of the system are fixed so that O 1 and O 2 admit an attractive MMO limit cycle γ µ,ε (t). There exist 0 < δ 0 , δ 1 , . . . , δ 5 1, such that, for different initial conditions on each subsystem, the coupled system (12) displays the following synchronization patterns, depending on parameter c.

1. Case c > δ 0 (regions D and E in Fig. 6). The two oscillators are inphase locking or almost-in-phase synchronized. The MMOs persist under the coupling, i.e. each oscillator generates small oscillations close to the lower fold of their respective critical manifold.

X ∈ [x min , -x f ] X ∈ [x f , x max ] (a) Stable equilibrium plus MMOs (b) Oscillatory configuration c > 0 (c = 1) -2 -1 1 2 x -4 -2 2 4 y S b S h ⇧ X = xmin X = xf xf xf X = xf S b S h ⇧ -2 -1 1 2 x -4 -2 2 4 y X = xmax xf xf (c) Oscillatory configuration. (d) Stable equilibrium c < 0 (c = -0.25) -2 -1 1 2 x -4 -2 2 4 X = xmin X = xf y ⇧ S b S h xf xf -2 -1 1 2 x -4 -2 2 4 y X = xmax X = xf S b S h xf xf ⇧
2. Case c < -δ 1 . We can distinguish three different subcases separated by the 510 limit values c d < c 0 < 0 :

(i) For c ∈ (c 0 + δ 2 , -δ 1 ), (region C in Fig. 6), the oscillators O 1 and O 2 are in antiphase synchronization.

(ii) For c ∈ (c d + δ 4 , c 0δ 3 ) occurs the relaxation loss of one oscillator, (region B in Fig. 6).

(iii) For c < c d -δ 5 , the 6D system ( 12) admits a stable equilibrium point:

there is a total oscillation death, (region A in Fig. 6).

The MMOs disappear in the three cases, (see panels (a), (b) and (c) of Fig. 4). Moreover, c 0 can be approximated by the only solution w.r.t.

c ∈ [-1, 0] of equation a 1 f (x) + (1 + c)x + a 2 -x f c -a 1 φ f (z b ) = 0, ( 18
)
where

x(c, µ) := 1 1 + c (-a 1 f (x f ) -a 2 + cx f + a 1 φ f (z h )) , (19)
and c d can be approximated as the only real solution w.r.t. c of equation

a 1 f (x) + (1 + c)x + a 2 -x f c -a 1 φ f (z b ) = 0, ( 20
)
where

x(c, µ) := 1 c [(1 + c)x f + a 1 (f (x f ) -φ f (z h )) + a 2 ] . ( 21
)
Sketch of the proof.

The proof is based on joint slow-fast dissections taking into account that both oscillators are slow-fast systems and, therefore, any kind of oscillatory patterns of the 6D system is constituted by slow motions of P 1 and P 2 close to S l and S r , fast motions from the folds to the opposite stable branch, and occasionally small oscillations around the lower fold. Hence, except in the case of total oscillation death, during the slow motions along the 6D system orbits, x (resp., X) either increases between values ε-close to x min and -x f or decreases between values ε-close to x max and x f . Depending on the value of c, using the results of Proposition 1, we verify for each case if the bidirectional conditions on O 1 and O 2 are fulfilled along the slow motions to produce the sequence of active phase, silent phase and small oscillations of O 1 and O 2 associated with each synchronization feature (cases (a) to (e) of Fig. 4). We have relegated to 530 Appendix C the detailed proof of statements 1 and 2 of the Theorem, and we include here the approximations of c 0 and c d .

Recall that we denote P 1 = (x, y, z) and P 2 = (X, Y, Z) the current points of O 1 and O 2 , respectively, along a given orbit of system ( 12) coupling O 1 and O 2 . We know that total oscillation death occurs if both oscillators reach their respective stable equilibrium point, each of them lying on a different sheet (left or right) of the critical manifold. Without loss of generality, we focus on the case where the equilibrium point of O 1 lies on S l , bearing in mind that the other case is exactly symmetric and arises for the same value of c. Since there exists a stable equilibrium point of O 1 on S l for every c < -δ 1 when X > 0, (see Fig. 9 (d)), we search for the value of c such that there exists also a stable equilibrium point of O 2 on S r , preventing the current point from heading towards S l . Hence, we can approximate the value c d by considering the following limit case: for P 2 coinciding with the upper fold F + of the special curve S h , the equilibrium point of O 1 on the left sheet is such that the equilibrium point of O 2 is located exactly at the upper fold. Such condition is expressed by the following system of equations that the x-component of P 1 should fulfill.

   g 1 (x, x f ; c) = 0, g 2 (x, x f ; c) = 0, (22)
where

   g 1 (x, X; c) = f (x) -φ f (z b ) -r(x, X; c), g 2 (x, X; c) = f (X) -φ f (z h ) -r(X, x; c), (23)
with r(x, X; c) given by eq. ( 16). System [START_REF] Zhabotinsky | A canard mechanism for localization in systems of globally coupled oscillators[END_REF] writes

   a 1 f (x) + (1 + c)x + a 2 -cx f -a 1 φ f (z b ) = 0, a 1 f (x f ) + (1 + c)x f + a 2 -cx -a 1 φ f (z h ) = 0. ( 24
)
The second equation of ( 24) is linear in x, hence for c = 0, we can express x as a function x(c, µ) leading to statement (21) of the theorem. By substitution in the first equation of ( 24), we obtain equation [START_REF] Ersöz | Synchronization of weakly coupled canard oscillators[END_REF], which provides the roots of 535 a polynomial function of degree 3 in 1/c. For |a 1 | small enough, the equation admits only one real solution c d .

We now consider the limit case between antiphase synchronization and relaxation loss for one oscillator, and search for an approximation of c 0 at order 0 w.r.t. ε. With this aim, we consider the fast motions of the current points to be instantaneous jumps from a fold to the opposite sheet of the critical manifold.

Along an orbit in the case of relaxation loss of one oscillator (O 1 to fix the ideas), one of the current points (P 1 ) remains stuck on the left sheet, following the equilibrium point of O 1 driven by X, while the other one (P 2 ) follows a relaxation cycle. Hence, the sign of x -X should not change along the whole orbit and, consequently, neither the sign of y -Y . On the contrary, along an orbit of the antiphase synchronization case, P 1 and P 2 alternatively undergo relaxation: when P 2 arrives at the upper fold and jumps on the left sheet of the critical manifold, P 1 is close to the moving equilibrium point of O 1 on the left sheet that shall already be below (i.e. admit a lower y-component than P 1 ), so that P 1 will be ahead of P 2 and subsequently undergoes relaxation. Hence, at order 0 w.r.t. ε, c 0 corresponds to the configuration of the nullclines shown in Figure 10: we look for c such that, for X = x f (P 2 is above the upper fold), the y-nullcline of O 1 intersects the x and z-nullclines at the height y of the upper fold. Therefore, c 0 can be approximated by the solution of

   g 1 (x, x f ; c) = 0, f (x f ) -φ f (z h ) = r(x, x f ; c), (25)
where g 1 (x, X; c) is defined by eq. ( 23) and r(x, X; c) by eq. ( 16). System [START_REF] Roberts | Averaging, folded singularities, and torus canards: Explaining transitions between bursting and spiking in a coupled neuron model[END_REF] writes

   a 1 f (x) + (1 + c)x + a 2 -cx f -a 1 φ f (z b ) = 0, a 1 f (x f ) + (1 + c)x + a 2 -cx f -a 1 φ f (z h ) = 0. ( 26
)
The second equation is linear in x, so for c = -1, we can find x = x(c) given by [START_REF] Ersöz | Canard-mediated (de) synchronization in coupled phantom bursters[END_REF] and substitute in the first equation of [START_REF] Sherman | Anti-phase, asymmetric and aperiodic oscillations in excitable cells -I. coupled bursters[END_REF]. Resulting equation [START_REF] Ermentrout | Canards, clusters, and synchronization in a weakly coupled interneuron model[END_REF] gives the roots of a polynomial function in c = 1/(c + 1) of degree 3. For a 1 small , the existence of a particular periodic orbit (x(t), y(t), z(t), X(t), Y (t), Z(t)) of period T such that x(t + T /2) = X(t), y(t+T /2) = Y (t), z(t+T /2) = Z(t) for all t. In the sequel, we will refer to this particular orbit as the "perfect antiphase periodic orbit". However, either proving that this periodic orbit is a limit cycle (isolation) or that it is stable are very difficult problems to tackle theoretically for such system. Remark though that Fig. 6 shows numerically that, for a large interval of c values, the perfect antiphase orbit is actually an attractive limit cycle. Moreover, it shows that the oscillation frequency is increasing with c, i.e. the return time decreases while c increases, that is, |c| decreases.

540 -3 -2 -1 1 2 x -4 -2 2 4 6 8 y S b S h X = x max X = x f ⇧ x f x f
No matter the isolation and attractiveness of the perfect antiphase periodic orbit, in this section, we give an approximation of its period for small values of ε by neglecting the duration of the fast parts of the dynamics. We prove that this approximate period decreases while c increases and we provide a global upper bound of its partial derivative w.r.t. c.

We first give a synoptic view of the durations that will be computed in the following theorems by means of Figure 11. From the symmetry of perfect antiphase periodic orbit, the duration T 1 of the active phase during which x > 0 and X < 0 equals the duration of the active phase during which x < 0 and X > 0. Similarly, the duration T 2 of the recovery phase during which x < X < 0 equals the duration of the recovery phase during which X < x < 0. The global period is then twice the sum of the active phase and the recovery phase durations, hence 2(T 1 + T 2 ).

For the sake of simplicity, in the following, we express both approximations T 1 and T 2 of the active phase and recovery phase durations using integrals !"

x max parameterized by variable x from x max to x f and from x 0 to -x f respectively, where x 0 is the x value of the equilibrium point of O 1 at the end of the O 2 active phase. In other words, we compute the approximate durations in the cases corresponding to the situations associated with circled T 1 and T 2 in Figure 11.

x min x 0 x(t) X(t) T 1 T 2 T 1 T 2 t x f -x f
We recall that, during the active phase, the current point

P 1 = (x, y, z) of O 1
follows the right slow manifold ε-close to y = f (x)-φ f (z), x > x f , for the Hausdorff distance. Along the corresponding trajectory, x decreases from a maximal value x max to the right fold value x f , and the state variables of the second oscillator remains close to a stable equilibrium point, slowly moving with the value of x. Once x becomes less than x f , the first oscillator experiences a fast motion of negligible duration before the beginning of the subsequent recovery phase.

The following theorem takes advantage of reductions to the slow manifolds and quasi-stationary approximations to state an integral expression for the duration of P 1 slow motion (under the impact of the bidirectional coupling) when x goes from x max to x f . Theorem 3.2. Assume the hypotheses of Theorem 3.1, 2.(i), c ∈ (c 0 +δ 2 , -δ 1 ), leading to antiphase synchronization between oscillators. Consider a parameterization (x(t), y(t), z(t), X(t), Y (t), Z(t)) of the perfect antiphase periodic orbit C(c, ε). The duration of O 1 (resp. O 2 ) active phase, i.e. the time spent by P 1 (resp. P 2 ) close to S r during one cycle, can be approximated at the first order in ε by T 1 (c, ε) + O(1) with

T 1 (c, ε) = 1 ετ x f xmax f (x) + ∂ψr ∂x (x, z, ε) a 1 f (x) + (1 + c)x -c X(x, c, ε) + k r (x, z, ε) dx , (27)
where

X(x, c, ε) = X eq (x, c, ε) + O(ε), (28)
X eq (x, c, ε) denotes the X component of O 2 equilibrium point moving with x and lying on the left sheet of the critical manifold,

k r (x, z, ε) = a 2 + a 1 ψ r (x, z, ε) -a 1 φ f (z) -φ f (z) + ∂ψ r ∂z (x, z, ε) φ r (x) - z -z b τ z , (29)
with φ f (x) given by eq. ( 2), φ r (x) given by eq. ( 3), and ψ r is a differentiable function defined on

(x f , x max ) × I r z × (0, ε 0 ), with I r z the interval of z values along C(c, ε) while x ∈ (x f , x max ), such that ∃ξ r (ε) = O(ε 2/3 ), such that, ∀(x, z) ∈ (x f , x max ) × I r z , |ψ r (x, z, ε)| < ξ r (ε). ( 30
)
Sketch of the proof.

595

The detailed proof of Theorem 3.2 is available in Appendix D. After rescaling the system, we use the symmetry of the system in the antiphase regime for giving an approximation of the time that the current point P 1 of O 1 spends along S r during one cycle. We parametrize the corresponding part of the cycle by (x, z) and use an implicit expression of the slow manifold S r,ε together with 600 the differentiable property of this parameterization. From the reduction of the X-dependent O 1 flow to the slow manifold, we obtain a formal expression of T 1 corresponding to eq. ( 27).

Theorem 3.3. With the same hypotheses and notations as in Theorem 3.2, the recovery time of the perfect antiphase periodic orbit can be approximated at first order in ε by T 2 (c, ε) + O(1) where

T 2 (c, ε) = 1 ετ -x f x0(c) f (x) + ∂ψ l ∂x (x, z, ε) a 1 f (x) + (1 + c)x -c X(x, c, ε) + k l (x, z, ε) dx , (31)
x 0 (c) is the x value of the equilibrium point of O 1 at the end of the O 2 active phase, i.e. the unique solution for x < -x f of equation

a 1 f (x) + (1 + c)x + a 2 -cx f -a 1 φ f (z b ) = 0, (32)
function X(x, c, ε) is the parameterization of the X value along the part of C(c, ε)

while x ∈ (x min , -x f ), k l (x, z, ε) = a 2 + a 1 ψ l (x, z, ε) -a 1 φ f (z) -φ f (z) + ∂ψ l ∂z (x, z, ε) φ r (x) - z -z b τ z , with φ r (x) given in (2)
, and ψ l is a differentiable function defined on

(x min , -x f ) × I l z × (0, ε 0 ), with I l z the z-path while x ∈ (x min , -x f ) along C(c, ε), such that ∃ξ l (ε) = O(ε 2/3 ), such that,∀(x, z) ∈ (x min , -x f ) × I l z , |ψ l (x, z, ε)| < ξ l (ε). ( 33
)
Sketch of the proof.

The detailed proof of Theorem 3.3 is available in Appendix E. It is mainly based on the same idea as the proof Theorem 3.2: we give an approximation of the time that the current point P 1 of O 1 spends along S l during one cycle. Yet an important difference consists in the parameterization X of X path during the silent phase, which naturally depends on x. After the reduction of the Xdependent O 1 flow to the left slow manifold, we keep this dependency in the formal expression of T 2 corresponding to eq. ( 31). Nevertheless, this does not impact the differentiability properties of the implicit functions introduced in the statement.

In the subsequent Theorem 3.4, we prove that T 1 and T 2 decrease as c

increases in [c 0 + δ, -δ]. Theorem 3.4. For c ∈ [c 0 + δ 2 , -δ 1 ] and ε small enough, ∂T 1 ∂c (c, ε) < 0 and ∂T 2 ∂c (c, ε) < 0. (34)
Consequently, the first order approximation in ε of the global period 2(T 1 (c, ε) + T 2 (c, ε)) of the perfect antiphase periodic orbit, decreases while parameter c increases and the global oscillation frequency is increasing with c.

Sketch of the proof. The proof of Theorem 3.4 is available in Appendix F.

We first prove a technical Lemma stating that x 0 (c) introduced in Theorem 3.3

(the x value of O 1 equilibrium point at the end of O 2 active phase) increases with c in the antiphase case. From Leibniz formula applied to the integral expression of T 1 , we obtain easily ∂T1 ∂c (c, ε) < 0. We use the same tool for T 2 , yet the dependency on X results in the sum of two terms in the integral defining ∂T1 ∂c (c, ε). The proof is based on the fact that the integral of the term known to be negative (from the Lemma) compensates the integral of the other one, no matter its sign. This is proved by bounding this integrand by a positive function that we can explicitly integrate.

Remark that the decrease of the period can only be obtained for a subinterval of (c 0 , 0) from the above approach. This restriction results indirectly from the fact that the approximation of the trajectories durations using the reduction to the slow manifold only makes sense if we consider ε small enough. The result can be stated in two different ways depending on whether we fix ε or δ. The first one is stated in Theorem 3.4: we define a fixed value of δ from geometrical properties (independently from ε), then there exists ε * > 0 such that for any [START_REF] Clément | Mathematical modeling of the GnRH pulse and surge generator[END_REF] hold. A second way could be to consider a value of ε small but fixed, then there exists a non-empty interval [c 0 + δ 2 , -δ 1 ] of c values along which the period of the perfect antiphase periodic orbit decreases. This statement implicitly link the optimal (i.e. smallest) values of δ 1 and δ 2 with ε that would tend to 0 with ε. Yet it is important to note that the proofs given in Appendices are not designed to allows δ to tend to 0 with ε.

ε ∈ (0, ε * ], for any c ∈ [c 0 + δ 2 , -δ 1 ], inequalities
It is worth noticing that numerical simulations for a fixed value of ε show an increase of the period for small (in absolute value) negative values of c. But, more interestingly, they also show an interval [-η, 0], (with η > 0) of c values for which the behavior is chaotic. In the next section, we provide numerical simulations and discuss the transitions undergone by the system behavior induced by the passage of c value from one macroscopic interval mentioned in Theorem 3.1 to another one.

Insight into almost-in-phase synchronization of MMOs

In Theorem 3.1, we have described the different behaviors of system [START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF], depending on the value of parameter c, out of the small regions in grey in Fig. 6, where tangled bifurcations take place. This section focuses on these intervals as well as on the sequence of bifurcations occurring in [0, c p ] where the MMOs play an essential role. This section is divided into two subsections. In the first one, we focus on the complex cascade of bifurcations for medium excitatory coupling. The second one is devoted to the chaotic behaviors observed for weak inhibitory and excitatory coupling.

Complex cascade of bifurcations for medium excitatory coupling

In the first item of Theorem 3.1, we have proven the qualitative property of almost-in-phase synchronization of the two coupled oscillators O 1 and O 2 resulting in MMO patterns for c > 0. Nevertheless, the existence of an attractive periodic orbit of the whole coupled system and the precise link between the value of c and the shift between the two oscillator patterns come under the theoretical study of the synchronization of MMOs, which remains a major and highly challenging problem. As a matter of fact, a large panel of complex bifurcations appears w.r.t. c in the case of weak excitatory coupling. In this section, we illustrate the high complexity of the bifurcation cascades occurring for c ∈ [0, c p ] using numerical simulations and discuss the associated types of transitions, in particular in terms of MMO signature. i=0 out of a hundred. In addition to the return times and phase shifts, we have retrieved the value of the first variable of each oscillator at the times when the other oscillator crosses the section {x = 0} (resp. {X = 0}), i.e. for each value of c and for i = [0, 25] ∩ N,

x X=0 i (c) = x c (t X=0 i (c)), (35)
X x=0 i (c) = X c (t x=0 i (c)). (36)
The two top panels of Figure 12 show, for each value of c, the 25 return times • For a given value of c, the unicity of the X and x return points together with a single non vanishing shift value indicates that the asymptotic behavior is an almost-in-phase synchronization of a simple kind: there exists a stable limit cycle such that, along one period, each oscillator has a single activity phase.

T x=0 i (c) = t x=0
• The existence of two (and only two) different return values for x and X, with a single phase shift value, corresponds to an attractive limit cycle such that, along one period, each oscillator has two activity phases. In Consequently, the phase shifts remain extremely close to 0.

Consider a starting c value leading to in-phase synchronization associated with a limit cycle of period T . We describe below the four types of transitions numbered from 1 to 4 observed in the simulation shown in Fig. 13 and refer to (both in Zoom4).

2. Infinite period doubling cascade: it results in chaotic dynamics for a small interval of c values, and leads to a limit cycle of period close to T . Such transition occurs, for instance, for c decreasing from 0.394 to 0.388 (Zoom1) and for c increasing from 0.465 to 0.475 (Zoom2). Note that the phase shift admits a unique value that is rather small, which evidences almost-in-phase synchronization. The resulting limit cycle is characterized by a single return value for one oscillator, greater than x f , while the single return value of the other oscillator is lower than -x f : the same oscillator is always ahead of the other one, i.e. if O 1 activates before O 2 at a certain time, O 1 always activates first afterwards. The whole transition consists in the addition of one small oscillation in the MMO generated by both oscillators (signature transition from 1 s to 1 s+1 ).

3. Infinite period doubling cascade resulting and priority exchange: it results in emergence of chaos and a subsequent structural re-stabilization of the dynamics and the existence of an attractive limit cycle. Yet, along this limit cycle, the oscillators exchange their priority from one activation to the subsequent one, i.e. if O 1 activates before O 2 at a certain time, O 2 activates first during the subsequent activation and the process repeats.

Such transitions occur for instance for c increasing from 0.396 to 0.398 (in Zoom1) and for c decreasing from 0.46 to 0.45 (in Zoom2). In that case, the signature of the MMO generated by each oscillator is 1 s 1 s+1 .

4. Sudden appearance of chaos emerging from a long period limit cycle:

this phenomenon occurs when decreasing c value from 0.52 to 0.515 (in Zoom3). This transition is rather mysterious since the numerical simulations do not show evidence for a period doubling cascade. Nevertheless, it is clearly different from the three transitions mentioned above since the long period orbit includes a return value for one of the oscillator that is below -x f . In other terms, another mechanism than the interactions between oscillators during their fast dynamics is involved in the emergence of the chaotic behavior.

Zoom5 in Fig. 13 (c ∈ values bounded below by 0, the whole system admits a limit cycle along which the oscillators exchange their priority from one activation to the subsequent one.

A period doubling cascade occurs for increasing c values, and therefore, chaotic behavior only emerges from a strictly positive value of c (around 0.02 in the simulation shown in Fig. 14).

Conclusions and future work

We have analyzed the synchronization of two identical 3D slow-fast systems generating MMOs under symmetric linear coupling according to the sign and strength of the coupling parameter. We have described the different synchronization patterns occurring in both the excitatory and inhibitory cases, corresponding to positive and negative values of coupling parameter c, respectively.

Using geometric arguments and tools based on singular perturbation theory, we have proven that inhibition between cells only preserves sustained MMOs for weak coupling (small |c| value), while excitatory coupling preserves MMO for any strength. We have also derived accurate approximations of the values of the (Loading Video...) This study leads to different conclusions depending on the scale at which we consider the model. Considering that each 3D oscillator represents an individual cell following MMOs in its dynamical behavior, this study proves that strong enough coupling leads to in-phase synchronization between cells. Hence, such system can also be interpreted as an idealized model of two clusters of neurons, with excitatory coupling within each cluster and strong enough for the cells to be in-phase synchronized. From this viewpoint, this study is a first step for understanding the impact of inhibition or weak excitation on model of clustered cells in (at least) two different contexts. On the one hand, the inhibition case finds a direct application in the study of motoneurons that are often organized as two clusters (ipsi-versus contra-lateral side of the spinal chord). The role of inhibition between the two clusters, that is suggested by experimental data [START_REF] Fallani | Hierarchy of neural organization in the embryonic spinal cord: Granger-causality graph analysis of in vivo calcium imaging data[END_REF],

is central leading to different reactions of the low motor system compartment to external stimulus. On the other hand, intracellular calcium concentrations have been recorded in cultured cells and these data (that we have used originally for fitting the single cell model [START_REF] Krupa | A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons[END_REF]) evidenced the role of the synchronization of calcium oscillations in the GnRH secretion. The medium excitation case of the present study can be used for understanding the transition from synchronized to desynchronized behaviors of neurons.

X < 0), the y-component of the intersection with x = 0 is smaller (resp. larger) than in the uncoupled case.

(a) Consider X ∈ (x min , -x f ), i.e. O 2 in phase I. Taking into account that x u is located in the middle sheet and Π c has a larger x-slope than Π, together with relation ( 17), we have Π c ∩ S ∩ Γ ∩ {X = -x f } = q := (x q , y q , z q ), with x q ∈ (-x f , x u ).

Since x q is close to the lower fold but still on the middle sheet, there exists a minimal value δ (0) > 0 such that, for c > δ (0) , Π c ∩ S ∩ Γ ∩ {X = x min } lies on S l , which implies the existence of a stable quasi-equilibrium point. Its

x-component increases with X values (see Fig.

2.

If c < 0, Π c has a lower x-slope than Π, which possibly gives rise to the existence of more than one quasi-equilibrium point for O 1 . Moreover, if X > 0, the y-component of the intersection with x = 0 is larger than in the uncoupled case and if X < 0, it is smaller.

(a) Consider X ∈ (x min , -x f ), i.e. O 2 in phase I. Bearing in mind that x u is located on the middle sheet, Π c has a lower x-slope than Π and from eq. ( 17), we have

Π c ∩ S ∩ Γ ∩ {X = -x f } = v := (x v , y v , z v ), with x v > x u .
So, while X ∈ (x min , -x f ), the possible quasi-equilibrium points of O 1 are on the right of the lower fold (see Fig. For c ∈ [-1a 1 f (0), 0) and a 1 f (0) > -1, O 1 admits a single quasiequilibrium point while X ∈ (x min , -x f ), and if c < -1a 1 f (0), subsystem O 1 can have one, two, or three quasi-equilibrium points while X ∈ (x min , -x f ). Hence, either a stable quasi-equilibrium point exists, lying on S r , and P 1 follows it, or O 1 generates relaxation oscillations.

(b) Consider X ∈ (x f , x max ), i.e. O 2 in phase IV. There exists a minimal δ (1) > 0 such that, for c < -δ (1) , Π c intersects S on its left sheet S l (see Fig. 9 (d)). Then, in this case, P 1 follows the corresponding stable quasiequilibrium moving rightwards as X decreases.

Finally, we take δ = max δ (0) , δ (1) . (i) For c ∈ (c 0 + δ 2 , -δ 1 ), the nullclines of O 1 when P 2 comes back near S l are represented in panel (d) of Fig. 9. Thus, when P 2 approaches the upper fold, P 1 is approaches the black line, so y < Y and P 2 comes back ε-close to S l when P 1 is below it. Therefore, P 1 follows S l and reach S r vicinity before P 2 , since P 2 pushes the y-nullcline of O 1 further to the right. remains below P 1 , both ε-close to S l and X > x. Hence, P 2 subsequently reaches S r vicinity before P 1 . Therefore, oscillator O 2 keeps on oscillating between S l and S r , and prevents P 1 from leaving S l vicinity. Thus, P 1 oscillates yet remains confined ε-close to S l and above the projection of the upper fold upon S l along the x-direction. In this case, the asymptotic periodic behavior is characterized by the relaxation loss of one oscillator, see Fig. 4 (a).

(iii) If c < c dδ 5 , the 6D system admits an attractive equilibrium corresponding to equilibria (x eq , y eq , z eq ) ∈ S l of O 1 with X = X eq and (X eq , Y eq , Z eq ) ∈ S r of O 2 with x = x eq (see Fig. C.18). Note that the symmetry in the system implies the existence of another stable equilibrium obtained by exchanging the coordinates between O 1 and O 2 . Depending on the initial condition of the 6D system, the associated orbit reaches either one or the other, both cases corresponding to symmetric total oscillation death. In cases 2.(i), (ii) and (iii), the oscillator does not generate MMOs (see panels (a), (b) and (c) of Fig. 4), since each oscillator prevents the y-nullcline of the other one from passing through the folds of the critical manifold, where the small oscillations would occur. We now prove the decrease of T 2 while c increases. We introduce the following additional notations:

α(x, c, ε) = a 1 f (x) + (1 + c)xc X(x, c, ε) + k l (x, z, ε), which is the denominator of the integrand in expression (31) defining T 2 . Hence, eq. ( 31) can be written as

T 2 (c, ε) = 1 ετ -x f x0(c)
g(x, c, ε)dx, with g(x, c, ε) = f (x) + ∂ψ l ∂x (x, z, ε) α(x, c, ε) .

Note that function g is directly linked with the x-dynamics of the 2D reduced system (E.5) by ẋ = ε g(x, c, ε) .

Since ẋ > 0 for x < -x f and X < -x f , we have g(x, c, ε) > 0 for x ∈ [x min , -x f ] and ε small enough. Moreover, f (x)+ ∂ψ l ∂x (x, z, ε) < 0 over [x min , -x f ], since it is Note that, from g > 0 and Lemma 1, g(x 0 (c), c, ε) ∂x0 ∂c (c) > 0. Moreover, x-X > 0 along all the recovery phase, c < 0 and, as already mentioned, for x < -x f , f (x) + ∂ψ l ∂x (x, z, ε) < 0.

Hence, the sign of the integrand ∂g ∂c (x, c, ε) depends on the sign of x -Xc ∂ X ∂c , which is hard to determine: x -X > 0 during the studied phase, c < 0 and 1180 ∂ X ∂c < 0, but difficult to evaluate. In the following, we prove that the negative term g(x 0 (c), c, ε) ∂x0 ∂c (c) compensates the integral (i.e. its opposite is greater than the integral) even if x -Xc ∂ X ∂c > 0 for some x ∈ [x 0 (c), -x f ]. In order to help the track of the signs along the sequence of following equations and inequalities, we write (whenever it is possible) the terms and factors depending 1185 on parameters and variables under a positive form with an explicit minus sign ahead if needed.

We first evaluate, using calculation (F.1) in the proof of Lemma 1, g(x 0 (c), c, ε)x 0 (c) = -f (x 0 (c)) -∂ψ l ∂x (x 0 (c), z 0 , ε) -α(x 0 (c), c, ε)

x fx 0 (c) a 1 f (x 0 (c)) + 1 + c . (F.4)

We now bound the integral from above, as follows. Consider Moreover, note that ∂ X ∂x is O(1) w.r.t. ε. On the other hand, ∂k l ∂x (x, z, ε) = O(ε 2/3 ) and, therefore, its absolute value is bounded from above by 1 + c, for ε small enough. Consequently, ∂α ∂x (x, c, ε) > 0 and we can write Using this upper bound in (F.5) together with max(0, x -Xc ∂ X ∂c ) < max(x -X), we obtain Since x min -2x f for ε small enough and x 0 (c) < -x f , max(x -X) < -x fx min < x fx 0 (c).

Finally, using this inequality, β(x) > 0, as well as inequality (F.8) and equation (F.7), the derivative of the recovery time w.r.t. c given by (F.3) is bounded from above as follows:

∂T 2 ∂c (c, ε) ≤ β(x 0 (c)) max(x -X) -(x fx 0 (c)) (-α(x 0 (c), c, ε)) -max(x -X) (-α(-x f , c, ε)) ≤ -β(x 0 (c)) max(x -X) (-α(-x f , c, ε)) < 0.

Figure 1 :

 1 Figure 1: Nullclines of system (1). Panel (a): x, y and z-nullclines are the purple, blue and yellow surfaces, respectively. Panel (b): Representations in the same (x, y)-plane of the half-planes z = z b for x < xon and z = z h for x > xon. Thus, the purple lines are the approximations S b and S h given in eqs. (9) and (10). The blue straight line represents the projection of the y-nullcline Π. The grey area close to S b (resp. Sr) is the projection of the attractive sheet S l (resp. Sr) of the fast dynamics for z ∈ [z b , z h ].

  well-approximated by the part of S b on the 170 left of the lower fold, and the projection of S r is well-approximated by the part of S h ∩ {z b ≤ z ≤ z h } on the right of the upper fold.

10

 10

10

 10

Figure 2 :

 2 Figure 2: Main behaviors of the single oscillator. Panel (a): relaxation oscillation. Panel (b): Mixed-Mode Oscillation (with a magnified view on the small oscillations). Panel (c): no sustained oscillation; the red point is an attractive equilibrium.

Figure 3 :

 3 Figure 3: Construction of x M (resp. xm) as the x-component of the intersection point between F - proj and the special curveS h = y = f (x)φ f (z b ), x < xon, z = z b (resp. F + proj and the special curve S b = y = f (x)φ f (z h ), x > xon, z = z h .

Figure 4 : 2 . 5 .

 425 Figure 4: Patterns generated by system (12) according to the values of c. The other parameter values are given by 11 and µ = 2.4. We emphasize the difference between almost-in-phase and in-phase locking synchronization cases with magnified views on the shifted signals generated by O 1 and O 2 in panel (d) and on coincident signals in panel (e), despite the different initial conditions taken for the two oscillators.

2 4 0Figure 5 :

 245 Figure 5: Illustration of the definitions of O 1 and O 2 crossing times and return times sequences and sequence of shifts between them on an instance of system (12) orbit.

Figure 6

 6 Figure 6 shows, for each c value in the 10 -3 step grid of [-1, 1], the return times (T x=0 i (c)) 25 i=1 (blue points) and (T X=0 i (c)) 25 i=1 (green points), respectively and the phase shifts (φ i (c)) 25 i=1 (red points). An animated view of the signals generated by x and X for various c values together with the picture of Fig. 6 is

3. 2 Figure 6 :Figure 4 )

 264 Figure 6: Synoptic view of the different synchronization patterns described in Theorem 3.1, depending on the value of the coupling parameter c ∈ [-1, 1]. The other parameter values are given by 11. For each value of c in intervals B to E, 25 return times for each oscillator (blue and green points) and 25 phase shifts are plotted. Blue points are hidden behind green points since they almost coincide for c value in intervals C and E and most of values in interval D.

Figure 7 :

 7 Figure 7: Repartition of synchronization types between O 1 and O 2 for c from -1 to 1 with a 2 10 -3 step value. Top panel : 50 seconds signals generated by x and X after a transient of 250 seconds starting from an initial data outside {x = X, y = Y, z = Z}. Bottom panel : partition into 5 intervals A to E of c values according to the type of synchronization between oscillators, sequences of 25 O 1 and O 2 return times (blue and green points), shifts (red points), half return times in case of antiphase synchronization (dashed black line). The moving vertical red bar materializes the c value corresponding to the associated frame in the top panel.

Figure 8 :

 8 Figure 8: Transition from total death oscillation to antiphase synchronization through relaxation loss of one oscillator for values of c varying from -0.53 to -0.48 with 10 -3 step value. The animation shows signals generated by x and X along an orbit of the coupled system after a transient of 250 seconds starting from an initial data outside {x = X, y = Y, z = Z}. In the manuscript, it corresponds to transitions between interval of c values A→B→C. Interval B associated to relaxation loss of one oscillator is very narrow. The detailed transitions (not shown in this animation) between the above synchronization types involve complex changes in the phase portrait (canard explosion and cascade of bifurcations) that occur in exponentially small (w.r.t. ε) intervals of c values.

Figure 9 :

 9 Figure 9: Projections of the nullclines of oscillator O 1 with parameter values given by 11 in different cases: c > 0 (top) and c < 0 (bottom), X ∈ [x min , -x f ] (left) and X ∈ [x f , xmax](right). We represent in the same (x, y)-plane the two planes z = z b for x < xon and z = z h for x > xon. In purple, the special curves of the critical manifold S b and S h given in eq. (9); in blue, plane Π corresponding to the uncoupled case; the grey area materializes the locus of plane Πc for X increasing from x min (red line) to -x f (black line) for (a) and (c), and for X decreasing from xmax (red line) to x f (black line) for (b) and (d). The dashed green lines materialize the x-component of the folds x = -x f and x = x f . In case (a) and (d), O 1 admits a stable equilibrium moving rightwards (cyan trajectories) while X decreases from xmax to x f : in (a), this point eventually crosses the fold leading to small oscillations, while in (d) this point lies on the left sheet for any X ∈ [x f , xmax]. In case (b) and (c), O 1 is in oscillatory configuration since the X-dependent singular point always lies on the middle sheet Sm of S.

Figure 10 : 3 . 3 .

 1033 Figure 10: Projections of the nullclines of oscillator O 1 with parameter values given by 11, c = c 0 ≈ -0.482 and X decreasing from xmax to x f . We represent in the same plane the two planes z = z b for x < xon and z = z h for x > xon. In purple, the special curves of the critical manifold S b and S h given in eq. (9); in blue, plane Π corresponding to the uncoupled case; the grey area materialize the locus of plane Πc for X between xmax (red line) and X = x f (black line). Considering X decreasing between these values, O 1 current point stays close to the stable singular point, which follows the cyan trajectory. The dashed green lines materialize the x-component of the folds x = -x f and x = x f and in brown the fast fibers corresponding to the upper fold F + . The c value corresponding to the transition between antiphase synchronization and loss of relaxation of one oscillator is approximated by the c value such that O 1 singular point for X = x f has the same y-component as the upper fold (black, purple and dashed brown lines intersect at the same point).

Figure 11 :

 11 Figure 11: Decomposition of the global period of the perfect antiphase periodic orbit into two active phase durations (T 1 ) and two recovery phase durations (T 2 ).

Figure 12

 12 Figure 12 shows return times, shifts and return values associated with the signals generated by O 1 and O 2 according to the value of c ∈ [0, c p ]. Return times and shifts are defined as explained at the beginning of section 2.5 and the same protocol has been used for the numerical simulations. Recall that we note (x c (t), y c (t), z c (t), X c (t), Y c (t), Z c (t)) the solution of the coupled system for a fixed value of c starting from an initial condition that does not belong to the subspace {x = X, y = Y, z = Z}. For this simulation, we have used a 10 -4 step grid of interval [0.1, 0.8] for c values and kept the 26 last values (t x=0 i (c)) 25 i=0 and (t X=0 i (c)) 25i=0 out of a hundred. In addition to the return times and phase shifts, we have retrieved the value of the first variable of each oscillator at the times

i

  (c)t x=0 i-1 (c) (green points) and T X=0 i (c) = t X=0 i (c)t X=0 i-1 (c) (blue points) for i = [1, 25] ∩ N. For the sake of comparison between the two sequences, we have plotted both on each panel, with the green points on top of the blue ones in the first panel and the blue points on top of the green ones in the second panel. The middle panel shows, for each value of c on the grid, the 25 evaluated shifts between the oscillators (φ i (c)) 25 i=1 . Two bottom panels show the sequences of crossed return values (x X=0 i (c)) 25 i=1 (green points) and (X x=0 i (c)) 25 i=1 (blue points), respectively. Once again, from Fig. 12, we can isolate macroscopic subintervals of interval [0.1, 0.8] of c values for which the system admits an attractive asymptotic behavior. Due to the symmetry between O 1 and O 2 assumed in this study (same parameter values and symmetric coupling), we obtain direct information on the MMO signature (number of small oscillations) displayed by x(t) and X(t) in the two following cases.

Figure 12 :

 12 Figure 12: Return times, shifts and return values associated with the signals generated by O 1 and O 2 according to the value of c ∈ [0, c p ]. Two top panels: return times (t x=0 i (c)) 25 i=1 (green points) and (t X=0 i (c)) 25 i=1 (blue points), respectively, according to c. Both are plotted on each panel, with the green points on top of the blue ones in the first panel and the blue points on top of the green ones in the second panel. Middle panel: for each value of c on the grid, the 25 evaluated shifts between the oscillators (φ i (c)) 25 i=1 . Two bottom panels: sequences of crossed return values (x X=0 i (c)) 25 i=1 (green points) and (X x=0 i (c)) 25 i=1 (blue points), respectively.

  that case, if O 1 begins an activity phase just before O 2 , this activity phase is preceded by n small oscillations for O 1 and n + 1 for O 2 . Then, just before the subsequent activity phase, O 2 and O 1 display n and n + 1 small oscillations, implying that O 2 becomes active slightly ahead of O 1 .In addition, chaotic behaviors arise when a large amount of different return values and phase shifts are recorded.We can isolate four patterns of bifurcation cascades occurring in the transition between the two "simple" cases mentioned above: they are emphasized and numbered from 1 to 4 in 5 zooms on small intervals of c values included in [0, c p ] shown in Fig.13. Due to the symmetry between the oscillators, return value 0 has a particular place in the structure of the bifurcation diagram: it is associated with a perfect in-phase synchronization between the two oscillators, as shown by the associated vanishing phase shift. Yet, for neighboring values of c, cascades of period doubling bifurcation occur. Such a cascade implies a route to chaos if an infinite number of period doubling occur in a compact interval of c values. It is worth noticing that, along such a cascade, the return value of each oscillator from one quasi-period to the subsequent one is necessarily non zero. Yet the return value corresponds to a fast dynamics of the oscillator.

Figure 13 : 1 .

 131 Figure 13: Magnified views of bifurcations in the transition from the cases of existence of an attractive limit cycle to chaotic behaviors.

4 . 2 .

 42 [0.68, 0.72]) shows a transition of the third type, i.e. an infinite period doubling cascade resulting in a almost-in-phase synchronized orbit with alternation of priority between oscillators O 1 and O 2 after each active phase. We emphasize the simulation on this interval of c values around c p since it corresponds to the stabilization of the symmetric orbit of the 6D system lying in {x = X, y = Y, z = Z} into an attractive limit cycle for strengthened coupling gain (transition between intervals D and E in Fig.6). Hence, for c > c p ≈ 0.716 corresponding to interval E, in-phase locking synchronization between oscillators O 1 and O 2 is the unique asymptotic behavior of the 6D system.This relatively large variety of transitions in the dynamics highlights the complexity resulting from the interactions between the MMO generated by each oscillators. Despite almost-in phase or in-phase synchronization feature has been shown in Theorem 3.1 (case 1.), the changes in the fine structure of the phase portrait during these transitions are highly complex and their analysis necessitate more acute tools than geometric arguments based on singular perturbation theory. Chaotic behaviors for weak inhibitory and excitatory coupling Chaotic behaviors are expected for weak coupling, yet, in this subsection, we point out the clear difference in the emergence of chaos for excitatory versus inhibitory weak coupling.

Figure 14 :

 14 Figure 14: Return times and return values for c ∈ [-0.08, 0.08] and 25 return values for each value of c in the 10 -4 step grid.

Fig. 14 shows

 14 Fig. 14 shows the return times and return values for c values in [-0.08, 0.08] and, as before, 25 return values for each value of c in the 10 -4 step grid for this interval. An animated view of the signals generated by x and X for varying c values together with the picture of Fig. 14 is provided as supplemental material (see animated Fig. 15). For c increasing along [-0.08, 0], the oscillators are first in perfect antiphase synchronization, and, from a precise value c chaos of c up to 0, the model enters a chaotic behavior, evidenced in Fig. 14 by the large amount of different values of return times and return values for any fixed c ∈ (c chaos , 0). Note that the generated signals are anti-phasic in the sense that, along an orbit that does not belong to {x = X, y = Y, z = Z}, x(t) and X(t) are never positive at the same time. On the contrary, weak excitatory coupling results in almost-in-phase synchronization between oscillators O 1 and O 2 . Moreover, for an interval of c

Figure 15 :

 15 Figure 15: Repartition of synchronization type between oscillator O 1 and O 2 for values of the coupling gain parameter from -0.06 to 0.08 with a 2 10 -4 step value. The top panel shows the 100 seconds signals generated by x and X along an orbit of the coupled system after a transient of 250 seconds starting from an initial data outside {x = X, y = Y, z = Z}. The bottom panel the sequences of 25 O 1 and O 2 return times (blue and green points), shifts (red points), half return times in case of antiphase synchronization (dashed black line). We refer to Fig. 14 in the manuscript for extended explanation. The moving vertical red bar materializes the c value corresponding to the associated frame of the top panel.

Figure B. 17 :

 17 Fig. B.17 (a) for a zoom around the lower fold). Then, the current pointP 1 = (x, y, z) of O 1 isforced to follow a stable quasi-equilibrium that is moving rightwards, then to display small oscillations and finally to quickly reach S r vicinity.

  9 (c) for a global view and Fig.B.17 (b) for a zoom around the lower fold). In this case, depending on the value of c, different possibilities may arise. Depending whether the x-slope of Π c is smaller or larger than f (0) (the slope of the critical manifold at its inflection w.r.t the x direction), the system possesses either one, two or three equilibria. Moreover, function r defining Π c in eq. (16) fulfills ∂r ∂x(0, X; c) = -c + 1 a 1 .

(the nullclines of O 2 Figure C. 18 :

 218 Figure C.18: Projections of the nullclines of oscillators O 1 (panel (a)) and O 2 (panel (b))with parameter values given by 11 and c = -0.7 < c d . We represent in the same plane the two planes z = z b for x < xon and z = z h for x > xon. In purple, the special curves of the critical manifold S b and S h given in eq. (9); in blue, the plane Π (uncoupled case); in red, y-nullcline for X = xmax (panel (a)) and the Y -nullcline for x = x min (panel (b)); in orange, the y and Y -nullclines for x = xeq and X = Xeq corresponding to the stable quasi-equilibrium point of the coupled system. We also represent in green the lines x = -x f and x = x f .

  Hence, we have the same situation as during the previous oscillation, yet O 1 and O 2 have exchanged their roles and the sequence repeats, leading to antiphase synchronization regime. (ii) For c ∈ (c d + δ 4 , c 0δ 3 ), the nullclines of O 1 when P 2 comes back εclose to S l are represented in Fig. C.19. Thus, when P 2 approaches to the upper fold, P 1 approaches the black line. Consequently, Y < y, P 2

Figure C. 19 :

 19 Figure C.19: Projections of the nullclines of oscillator O 1 with parameter values given by 11, c = -0.502 and X decreasing from xmax to x f . We represent in the same plane the two planes z = z b for x < xon and z = z h for x > xon. In purple, the special curves of the critical manifold S b and S h given in eq. (9); in blue, plane Π corresponding to the uncoupled case; the grey area materialize the locus of plane Πc for X between xmax (red line) and X = x f (black line). Considering X decreasing between these values, O 1 current point stays close to the stable singular point, which follows the cyan trajectory. The dashed green lines materialize the x-component of the folds x = -x f and x = x f and in brown the fast fibers corresponding to the upper fold F + . In the case of relaxation loss of O 1 , its singular point remains above (i.e. has a greater y-component than) the upper fold: the intersection of the black line with the purple curve is above the dashed brown line.

1175

  the partial derivative w.r.t. x of the parameterization of C l lying exponentially close (w.r.t. ε) to the left slow manifold. Consequently, α(x, c, ε) < 0 above [x min , -x f ] as well.Using again the Leibniz formula, we obtain, c, ε)dxg(x 0 (c), c, ε)x 0 (c) ,(F.3) Since ∂α ∂c (x, c, ε) = x -X(x, c, ε)c ∂ X ∂c (x, c, ε),we have∂g ∂c (x, c, ε) = -(f (x) + ∂ψ l ∂x (x, z, ε))(x -X(x, c, ε)c ∂ X ∂c (x, c, ε)) α 2 (x, c, ε).

  c, ε) = a 1 f (x) + 1 + c + (-c) ∂ X ∂x (x, c, ε) + ∂k l ∂x (x, z, ε).On the one hand, since x(t) and X(t) are strictly increasing along C l (c, ε), we directly obtain, for any x ∈ [x 0 (c), -x f ],Ẋ(t) = ẋ(t) ∂ X ∂x (x(t), c, ε) > 0 =⇒ ∂ X ∂x (x, c, ε) > 0.

(F. 5 )

 5 (x) + ∂ψ l ∂x (x, z, ε))( X(x, c, ε)x + c ∂ X ∂c (x, c, ε)) ∂α ∂x (x, c, ε) (-∂α ∂x (x, c, ε)) α 2 (x, c, ε) dx.Then, for ε small enough, for anyx ∈ [x 0 (c), -x f ], (-c) ∂ X ∂x (x, c, ε) + ∂k l ∂x (x, z, ε) > 0,which allows us to obtain-(f (x) + ∂ψ l ∂x (x, z, ε)) ∂α ∂x (x, c, ε) = -(f (x) + ∂ψ l ∂x (x, z, ε)) a 1 f (x) + 1 + c + (-c) ∂ X ∂x (x, c, ε) + ∂k l ∂x (x, z, ε) ≤ -(f (x) + ∂ψ l ∂x (x, z, ε)) a 1 f (x) + 1 + c = β(x). (F.6)Remark that β(x) > 0 and that the above notation allows us to write the term evaluated in equation (F.4) asg(x 0 (c), c, ε) ∂x 0 ∂c (c) = β(x 0 (c)) x fx 0 (c) -α(x 0 (c), c; ε) . (F.7)Since f (x) is positive for x < 0 and both ∂ψ l ∂x and ∂ 2 ψ l ∂x 2 are O(ε 2/3 ), for ε small enough,β (x) = -f (x) 1 + c + (-a 1 ) ∂ψ l ∂x (x, z, ε) -(1 + c) ∂ 2 ψ l ∂x 2 (x, z, ε) (a 1 f (x) + 1 + c) 2 < 0,and, for any x ∈ [x 0 (c), -x f ],-f (x) -∂ψ l ∂x (x, z, ε) ∂α ∂x (x, c, ε) < β(x 0 (c)).

  c, ε)dx ≤ β(x 0 (c)) max(x -X) -x f x0(c) -∂α ∂x (x, c, ε) (α(x, c, ε)) 2 dx ≤ β(x 0 (c)) max(x -X) 1 α(-x f , c, ε) -1 α(x 0 (c), c, ε). (F.8)

  Ca to increase and decrease back to the n increase or decrease the calcium peak level by tuning the value of ation in parameter k value also implies a change in the quiescence PI. Panel B of Figure ?? shows, together with the initial red signal min), a blue signal obtained with k = 1.2 (peak level at 320 nM, tained with k = 0.8 (peak level at 365 nM, IPI=16 min). Of course, level and the IPI by tuning first the value of k and afterwards the Figure ?? shows calcium patterns obtained with the same set of k erences in the IPI as been corrected by choosing appropriate values , 2.238 for the yellow one.

	1, a2 = 0.8, ⇧ = 0.06,
	2.4, Ca0 = 500, Ca bas = 100, ⌅Ca = 2,	(8)
	175, ⇤Ca = 4.5,	x on = 0.45.
	d calcium pattern fulfills the typical qualitative behavior described
	tive specifications: 10 min of IPI, peak level at 342 nM (red signal
	imic the variability of the quantitative features of calcium patterns
	sing di erent values for parameters of special importance: µ and
	e above explanation, one can already see that the precise value
	small oscillations of the current point near the left fold F and, the quiescent phase. Since variations in µ does not impact much
	can consider that this parameter controls marginally the InterPulse
	re ?? shows the results of change in µ: the red signal is obtained
	=10 min), the blue one with µ = 2.32 (IPI=7 min) and the yellow
	n). It is worth noticing that the interval of values of µ for which a
	n two successive pulses in the calcium pattern depends on the other
	es the time scale separation between y and Ca (x being much more
	y implies a shorter time for subsystem (3a)-(3b) to complete a
	equently, a shorter time for

  Total oscillation death, c = -0.7.

					x(t), X(t)								x(t), X(t)
											4
											2
											0
		-2									2
		-4	10	20	30 z(t), Z(t) 40	50	60	70	0 4 4 x(t), X(t) 10 20	30 z(t), Z(t) 40	50	60	70
	4	x(t), X(t)					2
			10	20	30	40	50	60	70		0 0	10	20	30	40	50	60	70
	2	(b) Relaxation loss, c = -0.502. x(t), X(t)			x(t), X(t) (c) antiphase synchronization, c = -0.25.
											4
	0										0 2
											2
	-2		10	20	30 z(t), Z(t) 40	50	60	70		0 4	10	20	30 z(t), Z(t) 40	50	60	70
											4
	-4 20 0	10 x(t), X(t) 20 30 40 (d) Almost-in-phase synchronization, c = 0.1. (e) In-phase locking synchronization, c = 1. 50 60 70 0 10 20 30 40 50 60 0 2 70 x(t), X(t) 30 40 50 60 x(t), X(t) 10 20 30 40 50 70 2 2 z(t), Z(t) 4 x(t), X(t) 4 z(t), Z(t)	60	7
											0
											-2
	4	0	10	20	30 z(t), Z(t) 40	50	60	70		0 -4	10	20	30 z(t), Z(t) 40	50	60	70
		-2									4
											2
	2	10 0 -4	20	30	40 10	50	60	70 20	0 0 30 z(t), Z(t) 10 20 30 40	40	50	60 50	70	60	70
	20 20 0 0	30 z(t), Z(t) 40 4 2 30 40 10 20			50 50 30			60 60 40		70 50 70	60	7
		0
		0			10			20	30		40		50	60	70

Corresponding author: soledad@us.es

alexandre.vidal@univ-evry.fr

Fig.1in[START_REF] Fallani | Hierarchy of neural organization in the embryonic spinal cord: Granger-causality graph analysis of in vivo calcium imaging data[END_REF] : https://arxiv.org/ftp/arxiv/papers/1407/1407.4380.pdf#page=19

Acknowledgments

The authors would like to thank Frédérique Clément and Fabrizio de Vico Fallani (INRIA Paris Centre) for funding and helpful discussions about this study. The first author is now supported by VPPI-US and partially supported by grants P12-FQM-1658 and MTM2015-65608-P.

Several perspectives of research are open by the results of the present study:

we evoke the two most natural ones. The first question is related to the fine characterization of the transition between orbits of different periods and signatures in the almost-in-phase synchronization case. Different analyses should be performed depending on the type of transitions shown and discussed in Section 4. Note that geometric arguments are no longer sufficient to understand the fine mechanisms of period-doubling cascades (finite or infinite) and the precise properties of the MMO orbits (period, signature, asymptotic stability). The second natural extension consists in considering (at least weak) heterogeneity between cells, i.e. (slightly) different values of the parameters. Indeed, the strong assumption that the two cells are identical has allowed us to benefit from the symmetry of the 6D system for performing theoretical analysis. A natural way to extend the study is to introduce heterogeneity in parameter values among cells and consider the system as a perturbation of the identical oscillator case. Away from the transition between qualitatively different behavior, the structural stability of the unperturbed vector field ensures that a hyperbolic limit cycle persists under small perturbation, i.e. in case of sufficiently weak heterogeneity between cells. In particular, a central question is to characterize the new synchronization types obtained with asymmetric coupling.

Appendix A. Computation of folded sigularities 1010

Here we compute the folded singularities of system (1) [START_REF] Guckenheimer | Singular hopf bifurcation in systems with two slow variables[END_REF][START_REF] Krupa | Local analysis near a folded saddle-node singularity[END_REF][START_REF] Szmolyan | Canards in R 3[END_REF][START_REF] Brøns | Mixed mode oscillations due to the generalized canard phenomenon[END_REF][START_REF] Guckenheimer | Unfoldings of singular hopf bifurcation[END_REF][START_REF] Wechselberger | Existence and bifurcation of canards in R 3 in the case of a folded node[END_REF]. First, we consider the new time variable s = τ t, and obtain the rescaled system

The associated reduced problem is given by

where we have denoted

Now, the reduced problem (A.1) projected onto the (x, z)-plane is given by

where f (x, z) defines the critical manifold S (see ( 4)). A time rescaling by function -h x drives us to the desingularized system

It is known that folded singularities are fold points p = (x, f (x)φ f (z), z) of the critical manifold such that

Remember that the critical manifold has two fold lines F -and F + given in Eq. (A.3) for x = -x f is linear in µ, so for every z, there exists μ(z) such that (-x f , z) is a folded singularity for the 3D system with µ = μ(z).

The linearization of the desingularized system (A. values are given by [START_REF] Storti | A simplified model of coupled relaxation oscillators[END_REF].

for z ∈ (1.14093, 1.54561) and a folded saddle-node for z = 1.14093, which corresponds to µ = 2.51. Moreover, it is possible to check that it corresponds to a real equilibrium of the full system (1), and then it is a folded saddle-node of type II [START_REF] Guckenheimer | Singular hopf bifurcation in systems with two slow variables[END_REF][START_REF] Krupa | Local analysis near a folded saddle-node singularity[END_REF][START_REF] Szmolyan | Canards in R 3[END_REF][START_REF] Brøns | Mixed mode oscillations due to the generalized canard phenomenon[END_REF][START_REF] Guckenheimer | Unfoldings of singular hopf bifurcation[END_REF][START_REF] Wechselberger | Existence and bifurcation of canards in R 3 in the case of a folded node[END_REF].

Appendix B. Technical details of the proof of Proposition 1.

Let us denote by x u the first component of the only equilibrium point in the uncoupled case, which lies on the middle sheet of S but close to the lower fold.

We proceed now through the different cases.

1.

If c > 0, Π c has a larger x-slope than Π. Moreover, if X > 0 (resp.

Appendix C. Technical details of the proof of Theorem 3.1.

Recall that we denote P 1 = (x, y, z) and P 2 = (X, Y, Z) the current points of O 1 and O 2 , respectively, along a given orbit of system [START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF] coupling O 1 and O 2 .

We consider δ 0 > δ and δ 1 > δ, with δ defined in Proposition 1, and proceed case by case. To fix the idea in the following description, without loss of generality considering the symmetry between the oscillators, we assume that P 1 is ahead of First, we consider in system [START_REF] Terman | Stability of anti-phase and in-phase locking by electrical coupling but not fast inhibition alone[END_REF] the new time variable

We obtain the rescaled system

where, now, the dot denotes the derivative w.r.t. the new time variable s. Using the symmetry of the system in the antiphase regime, we will approximate the time that the current point P 1 of O 1 spends along S r during one cycle. Consider the following set where ẋ < 0:

with I r z the z-path while x ∈ (x f , x max ) along C(c, ε). When c ∈ (c 0 + δ 2 , -δ 1 ), the oscillators do not possess any equilibrium point on S r . Since ẋ = 0 along S r and ẏ > 0 in U r , every orbit starting in U r escapes from U r across the half-plane

From the dynamics of subsystem (D.2), we have also, along any trajectory close to the critical manifold in U r ,

It follows from eqs. (D.6)-(D.7) that such a trajectory can be represented as the graph of a differentiable function of (x, z), in particular, the trajectory corresponding to C(c, ε). Let us parameterize this trajectory (the intersection of

and define

We can directly deduce that ψ r > 0 and is differentiable w.r.t. x and z.

1130 Moreover, from classical results [START_REF] Szmolyan | Relaxation oscillations in R 3[END_REF][START_REF] Szmolyan | Canards in R 3[END_REF] on fold-induced slow-fast transitions, there exists ε 0 such that for ε ∈ (0,

, and fulfills statement [START_REF] De Pittà | A tale of two stories: astrocyte regulation of synaptic depression and facilitation[END_REF]. By setting y = f (x)φ f (z) + ψ r (x, z, ε) into eq. (D.2), we obtain the following 2D reduction of O 1 along the slow manifold

From Proposition 1, statement 2.(b), when x ∈ (x f , x max ), it forces P 2 to follow a stable quasi-equilibrium point X eq (x, c) on the left sheet of the critical manifold.

Thus, we approximate X X given in eq. [START_REF] Haydon | GLIA: Listening and talking to the synapse[END_REF]. Then, the time that variable x spends along the right sheet of the slow manifold during one cycle is given by: As previously, due to the symmetry between the oscillators, we consider, without loss of generality, the recovery phase with x > X, i.e. O 1 (variables

The following proof is based on the same technics as the preceding one. However, we point out a capital difference between the stated results and linked with the boundaries of the integral expressions and quasi-stationary approximations. We recall that the recovery phase that we consider starts just after the fast motion of P 2 . Since we can neglect the O(ε) duration of P 2 fast motion, we can approximate O 1 starting point of the recovery phase by its stable equilibrium, lying on S l , when X = x f (corresponding to the end of O 2 active phase). The corresponding x value is then the unique solution x 0 (c) ∈ (x min , -x f ) of eq. ( 32).

Similarly as in the proof of Theorem 3.2, we consider the following set where ẋ > 0:

where I l z denotes the z-path while x ∈ (x min , -x f ) along C(c, ε). From the same argument as in the previous proof, every orbit starting in U l escapes from U l across the half-plane

relations (D.6) and (D.7) are fulfilled, and such a trajectory can be represented as the graph of a differentiable function of (x, z). We denote the intersection of

and we pose

We can directly deduce that ψ l < 0 and is differentiable with respect to x and z. As previously, there exists ε 0 such that for ε ∈ (0,

, and verifies relation 1155 [START_REF] Duchen | Mitochondria, calcium-dependent neuronal death and neurodegenerative disease[END_REF].

By setting y = f (x)φ f (z) + ψ l (x, z, ε) into system (D.2), we obtain the reduced 2D system,

During the recovery phase, the current point P 2 of oscillator O 2 follows the left slow manifold and, for any x ∈ [x 0 (c), -x f ], Ẋ > 0. Hence, we can parameterize the X component of trajectory C l (c, ε) by variable x. We keep the dependency w.r.t. c and ε explicit and note this parameterization X(x, c, ε). Then, at the beginning of the recovery phase (after neglecting the O(ε) duration corresponding to the fast transition from X = x f to X = X min ), x is ε-close to x 0 (c).

It increases and remains greater than X during the recovery phase which ends when x reaches -x f . Hence, from (E.5), the time spent by P 1 along the left slow manifold from x = x 0 (c) to x = -x f is given by

where we have already taken into account rescaling (D.1). From relation (E.6)

and system (E.5), direct computations lead to expression [START_REF] Wang | The astrocyte odyssey[END_REF]. Proof. By differentiating where k r (x, z, ε) is given in [START_REF] Vandenberghe | Subcellular localization of calcium-permeable ampa receptors in spinal motoneurons[END_REF]. Since c < 0, x > x f , X < 0, and ∂ Hence, both T 1 and T 2 decrease while c increases in [c 0 + δ 2 , -δ 1 ], and the approximation 2(T 1 + T 2 ) of the global period decreases, concluding the proof.