Jin Gyu Lee

Rodolphe Sepulchre

Rapid and robust synchronization via weak synaptic coupling Extended arXiv version

Keywords: Rapid synchronization, heterogeneous multi-agents, excitable systems, weak synaptic coupling, slow-fast models

This paper examines how weak synaptic coupling can achieve rapid synchronization in heterogeneous networks. The assumptions aim at capturing the key mathematical properties that make this possible for biophysical networks. In particular, the combination of nodal excitability and synaptic coupling are shown to be essential to the phenomenon.

Introduction

Biophysical neuronal networks exhibit rapid transitions between asynchronous and synchronous modes in spite of weak coupling and significant heterogeneity in the network. Such properties might be desirable in engineered networks, but they do not seem immediately aligned with the mathematical analysis of synchronization. A standard model of synchronization in physical models is through diffusive coupling in a network of Van der Pol nonlinear oscillators, see, e.g., (Wang & Slotine, 2005;Stan & Sepulchre, 2007;[START_REF] Angeli | A Lyapunov approach to incremental stability properties[END_REF][START_REF] Arcak | Passivity as a design tool for group coordination[END_REF]. Synchronization in those models is the result of a contraction property for the error dynamics. The analysis requires identical oscillators and the coupling strength must be strong enough to overcome the shortage of contraction in the decoupled model. The robustness of this mechanism to the heterogeneity of the network has been studied by several authors, see, e.g., [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF][START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF][START_REF] Lee | Heterogeneous Van der Pol oscillators under strong coupling[END_REF][START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF]. However, the common conclusion of such analysis is that stronger coupling is required to compensate for stronger This paper was not presented at any IFAC meeting. Corresponding author R. Sepulchre. The research leading to these results has received funding from the European Research Council under the Advanced ERC Grant Agreement Switchlet n.670645 and SpikyControl n.101054323. This work was done while Jin Gyu Lee was with University of Cambridge.

Email addresses: jin-gyu.lee@inria.fr (Jin Gyu Lee), r.sepulchre@eng.cam.ac.uk (Rodolphe Sepulchre).

heterogeneity. This conclusion is consistent with synchronization studies in more abstract phase models such as the celebrated Kuramoto model [START_REF] Kuramoto | Self-entrainment of a population of coupled non-linear oscillators[END_REF]Strogatz, 2000;[START_REF] Dörfler | Synchronization and transient stability in power networks and nonuniform Kuramoto oscillators[END_REF].

The goal of the present paper is to highlight the additional modeling properties needed to reconcile the apparent discrepancy between biophysical observations and mathematical analysis. We will analyze synchronization in a model of excitable neurons under synaptic coupling.

Those two properties are central to biophysical models of synchronization: a mechanism akin to fast threshold modulation (Somers & Kopell, 1993). The nonlinear nature of the synaptic coupling makes it possible to bound the coupling by a given constant and at the same time make its gain high in a narrow amplitude range. Synaptic coupling hence reconciles the weak coupling biophysical assumption with the localized strong coupling property required for mathematical contraction. Excitable models have a threshold, that is, a highly localized amplitude range where the spike can be triggered. Thanks to the threshold, the range of strong coupling can be narrowed down to where it is most beneficial to synchronization. The two properties combined together drastically simplify the analysis of the overall network because the neurons are nearly decoupled almost everywhere in the state-space. As a result, the convergence analysis decomposes into a sequence of distinct phenomena:

(1) Rapid convergence of each individual neuron to a neighborhood of its limit cycle attractor, deter-mined by the absolute bound on the coupling term treated as an exogenous disturbance; the weaker the coupling, the better the bound (Sections 3.1-3.2). (2) Owing to the fast-slow nature of the oscillations and to the synaptic nature of the coupling, the limit cycle trajectory of each neuron is rapidly entrained by the short coupling pulses triggered by each spike of presynaptic neurons. The rate of entrainment is independent of the strength of the coupling (Section 5.1).

(3) Phase synchronization of the spiking events occurs because the coupling pulses of the presynaptic neurons only affect the phase of the postsynaptic neuron near the threshold, the narrow range where they reduce the phase difference (Section 5.3).

Our objective is to show that those three intuitive properties can be demonstrated in a standard mathematical model using standard analysis. In particular, the simple linearization analysis of the differential equation highlights how the synaptic coupling provides strong diffusive coupling only near the threshold, making the overall network nearly decoupled anywhere else.

The specificity of pulse coupling with respect to diffusive coupling has been thoroughly studied in the dynamical systems literature, see, e.g., the seminal paper (Mirollo & Strogatz, 1990a) and the tutorial paper [START_REF] Mauroy | Kick synchronization versus diffusive synchronization[END_REF]. However, combining the slow-fast nature of spiking oscillations with a reduction to phase models has proven challenging [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF][START_REF] Sacré | Sensitivity analysis of oscillator models in the space of phase-response curves: Oscillators as open systems[END_REF][START_REF] Sacré | Singularly perturbed phase response curves[END_REF]. Those technical obstacles have made it difficult to connect the literature on diffusive coupling in ordinary differential equations with the literature on abstract phase models. The analysis in the present paper does not perform any phase reduction and proceeds straight from linearization analysis of the full differential equation model: we instead perform a network-wise singular perturbation analysis.

A strong source of inspiration for the present paper is the work by Somers andKopell (Somers &Kopell, 1993, 1995). To the best of our knowledge, they were the first to observe the rapid synchronization phenomenon in a model that combines excitability and synaptic coupling (which yields fast threshold modulation). They provided a detailed analysis of two identical relaxation oscillators (two-time scale oscillators characterized by a squarewave oscillation) under weak excitatory synaptic coupling and provided numerical evidence that a ring of relaxation oscillators synchronizes quicker than a ring of sinusoidal oscillators, under excitatory synaptic coupling. This basic but key observation is illustrated in Figure 1 (a) and is a starting point for the analysis in the present paper.

The phase portrait of Figure 1 (a) illustrates why localizing the coupling near the threshold is key to rapid synchronization. In this example, the vertical direction is the direction where coupling is exerted. Under uniform diffusive coupling, most of the coupling force is balanced by the strong attraction of individuals to their own limit cycles, hence resulting in a slower convergence rate for synchronization (Figure 1 (c) illustrates that the strong diffusive coupling significantly deforms the trajectory of the oscillation). In contrast, the synaptic coupling is active only near the threshold, where it matches the direction of pulling. This mechanism is quite robust to heterogeneity as illustrated by Figure 1 (b). An idea to design nonlinear coupling to synchronize identical agents efficiently under the same philosophy has been introduced in [START_REF] Pavlov | Nonlinear integral coupling for synchronization in networks of nonlinear systems[END_REF], where the goal is to make the error dynamics contractive everywhere.

The remainder of the paper is organized as follows. We motivate the discussions made earlier with a concrete example in Section 2. Then, we introduce a generalized model in Section 3. The latter two steps of the convergence analysis are detailed in Section 5 based on a network-wise reduction illustrated in Section 4 followed by a discussion in Section 6. The first step of the convergence is straightforward from the definition of the model (Section 3.2). A simulation result illustrating the robustness of these biophysical assumptions with respect to heterogeneity is given in Section 7. Finally, we conclude in Section 8. Additional technical details are given in the Appendix.

Notation: A (directed) graph is a pair G = (N , E) consisting of a finite nonempty set of nodes N and an edge set of ordered pairs of nodes E ⊆ N ×N , where (i, i) / ∈ E for all i ∈ N (i.e., the graph does not contain self-loops).

The set

N i = {j ∈ N : (j, i) ∈ E} denotes the neigh- bors of the node i. A tuple (i 0 , i 1 , ..., i l ) is called a path (of length l) from node i 0 to node i l if i k ∈ N i k+1 for all k = 0, . . . , l -1. Given a graph G = (N , E), let N ⊆ N and E ⊆ E| N := {(j, i) ∈ E : i, j ∈ N }. Then the pair G = (N , E ) is called a subgraph of G. If N = N , then
G is a spanning subgraph. A graph G is strongly connected if, for any pair of agents i and j, there exists a path from i to j. A spanning subgraph T of G having a root i such that there is exactly one path from i to any other node is called a spanning tree of G. An independent strongly connected component (iSCC) of G is a subgraph G = (N , E ) such that it is maximal subject to being strongly connected and satisfies that (j, i) / ∈ E for any j ∈ N \ N and i ∈ N . For vectors or matrices a and b, col(a, b) := [a T , b T ] T .

A motivating example

The simulation in Figure 1 involves the standard Morris-Lecar excitable model

τ (v) ẋ = -x + x ∞ (v), v = g L (-0.4 -v) + g Ca m ∞ (v)(1 -v) + g K x(-0.7 -v) + 0.4 + I (1)
where the small parameter = 0.02 controls the time scale separation of the model, and

τ (v) = 1 cosh v+0.1 0.29 , x ∞ (v) = 1 2 + 1 2 tanh v + 0.1 0.145 , m ∞ (v) = 1 2 + 1 2 tanh v 0.15
.

Figure 1 (a) considers two identical neurons, with (voltage) states v and v. The numerical values of the leak, Ca ++ , and K + maximal conductances are g L = 0.5, g Ca = 1, and g K = 2, respectively. The coupling term from the (presynaptic) neuron v to the (postsynaptic) neuron v is the synaptic current

I = g coup (1 -v)m ∞ (v)
called excitatory for v ≤ 1 and weak if the maximal conductance parameter g coup > 0 is small. The simulation uses g coup = 0.67. We note that the weaker g coup , the weaker the external input I in (1). We assume that the coupling is symmetric, just for simplicity.

A simple calculation provides the linearized error dynamics along the synchronized solution v(•) = v(•). The linearized error dynamics have the following expression:

δ x = - 1 τ (v(t)) δ x + b(t)δṽ, δṽ = [a(t) -k(t)] δṽ + g K (-0.7 -v(t))δ x, (2) with b(t) = τ (v(t)) [x(t) -x ∞ (v(t))] + τ (v(t))x ∞ (v(t)) τ (v(t)) 2 , a(t) = -g L -(g Ca + g coup )m ∞ (v(t)) + g Ca m ∞ (v(t))(1 -v(t)) -g K x(t), k(t) = g coup (1 -v(t))m ∞ (v(t)).
System (2) has the familiar form of error dynamics of a second-order linear system with diffusive coupling. The coupling gain k(t) is always non-negative and can be large in the local range of the sigmoidal nonlinearity, even if the synaptic coupling strength g coup is small. In other words, the synaptic coupling can provide a strong diffusive coupling gain in a localized voltage range. According to the classical theory of diffusive coupling, this large gain can make the linearized error dynamics contractive, which guarantees the exponential synchronization of solutions starting in the vicinity of the limit cycle. In the rest of the paper, we will see that this results in an open set of synchronous spiking in the limiting case in which the spikes take place infinitely fast: that is akin to fast threshold modulation, yielding rapid convergence.

It will be also apparent that owing to the spiky nature of the oscillation, the coupling will almost vanish along the slow branch of the limit cycle attractor, only contributing to the contraction during a synchronous spiking with the help of the sufficiently steep synaptic coupling and the sufficient time scale separation, which gives nodal excitability.

This will be not so different for heterogeneous systems as illustrated in Figure 1 (b). Here, we consider the heterogeneous network where the second neuron has the maximal conductance parameters g L = 0.25, g Ca = 0.5, and g K = 4. According to the blended dynamics theory [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF][START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF][START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF], a sufficiently large gain k(t) can also make the error dynamics of heterogeneous networks contractive.

For the numerical example shown in this section, the value

a(t) -k(t) = -0.5 -(1 + g coup )m ∞ (v(t)) + (1 -g coup )m ∞ (v(t))(1 -v(t)) -2x(t)
becomes uniformly negative when (since x ∈ [0, 1])

g coup = 0.67 > max v∈[-0.7,1] -0.5 -m ∞ (v) + m ∞ (v)(1 -v) m ∞ (v) + m ∞ (v)(1 -v) .
Meanwhile, further simulation with a smaller g coup suggests an example where the periodic linear system is contracting over one period, but not point-wise in time. In the general model treated in this paper, this mechanism will lead to synchronization of the spiking times, rather than arbitrary precision approximate synchronization to a common limit cycle trajectory, especially for heterogeneous networks. The spiking time refers to the timing of the fast time scale behavior. This timing shrinks to a single instant of time in the singular limit where the spikes occur infinitely fast. See Section 3 for more detail. This weaker form of synchronization (called phase synchronization) is the one considered in this paper.

3 Model assumptions

Spiking model for individual neurons

We consider N heterogeneous neurons individually modeled by the second-order ultrafast-fast-slow system

τ i, (v i ) ẋi = -x i + x i,∞ (v i ) vi = f i (x i , v i ) + (E i -v i )G i (3)
with i ∈ N := {1, . . . , N }, and G i is the external input of neuron i. Here, τ i, (v) > 0 is a state-dependent fast-slow time constant. In the absence of coupling, the external input is G i = 0. The phase portrait is then the classical Morris-Lecar phase portrait. the lower (upper) branch of the nullcline f i (x, v) = 0, hence there is only one intersection between the two, which is located in the middle. We have f i (x, v) > 0 (< 0) in the left (right) region of the nullcline f i (x, v) = 0.

Under Assumption E1, the Poincare-Bendixson theorem ensures that all solutions (except the single unstable equilibrium) converge to a limit cycle. One easily determines a positively invariant compact set D i that includes the limit cycle (the grey box in Figure 2). The external input (E i -v i )G i is called 'excitatory,' as the input gain (E i -v i ) is positive inside the positively invariant set D i .

The spiking nature of the limit cycle is ensured by the additional assumption that the time constant

τ i, (•) is O(1)
in the lower branch, while O( q ) with some q ∈ (0, 1) in the upper branch of the nullcline f i (x, v) = 0. Indeed, as > 0 becomes small, this ensures that the limit cycle exhibits a slow-fast-ultrafast behavior of spiking: the slow part of the limit cycle oscillation is a slow slide along the slow (lower) branch, whereas the ultrafast part is a switch between the lower branch and the upper branch, and finally the fast part is a fast slide along the fast (upper) branch. Such assumptions are standard in biophysical models of neuronal models: significant time scale separation between the voltage dynamics and the ion channel kinetics is only pronounced in the subthreshold region.

Assumption E2 (Biophysical nature of spiking) We assume that τi (v) := lim →0 τ i, (v) ≥ 0 and λi (v) := lim

→0 q τ i, (v) < ∞
are well-defined in the set D i , and satisfy

τi (v) = 0, if v ≥ E th i > 0, if v < E th i and λi (v) > 0, if v > E th i = 0, if v ≤ E th i
with some E th i that is located between the upper branch and the lower branch, as illustrated in Figure 2.

The spiking behavior described above for an isolated neuron (G i = 0) is robust to the perturbations induced by synaptic coupling. This is formalized in the following assumption.

Assumption E3 (Robustness margin) We assume that there exists M i > 0, such that for any m ∈ [0, M i ], the perturbed system of (3) represented as

τ i, (v i ) ẋi = -x i + x i,∞ (v i ), vi = f i (x i , v i ) + (E i -v i )m =: f m i (x i , v i ),
has the same D i as its positively invariant set, and satisfies Assumptions E1 and E2 with the same E th i . In particular, we assume max m∈

[0,Mi] v i (m) < E th i < min m∈[0,Mi] v i (m), where (x i (m), v i (m)) ((x i (m), v i (m))
) denotes the left (right) knee (that is, the extreme left (right) point of the nullcline f m i (x, v) = 0) of the perturbed system. We also assume x i (M i ) < x i (0), as illustrated in Figure 2.

Lemma 1 Under Assumptions E1 and E3, the nullcline f m i (x, v) = 0 moves to the right as m increases, i.e., for each v < E i , we have (∂x m i /∂m)(v) > 0 for all m ∈ (0, M i ), where x m i (v) denotes the unique root of f m i (x, v) = 0. In particular, the left knee, the right knee, and the intersection move to the right as m increases.

PROOF. Since Assumption E1 is satisfied for all m ∈ [0, M i ], for each v < E i , there exists a unique root

x m i (v) of f m i (x, v) = 0. Moreover, since f m i (x, v) < 0 (> 0) for x > x m i (v) (< x m i (v)), we have (∂f i /∂x)(x m i (v), v) = (∂f m i /∂x)(x m i (v), v) < 0. Now, this implies from ∂f i ∂x (x m i (v), v) ∂x m i ∂m (v)+(E i -v) = ∂ ∂m f m i (x m i (v), v) = 0 that (∂x m i /∂m)(v) > 0 for v < E i .
Assumption E3 also illustrates the nature of the excitatory coupling: a weak pulse G i creates an open set near the left knee of the isolated system (3), advancing the timing of the next spike.

Synaptic coupling

A standard model of excitatory synaptic coupling has the form

G i (t) = g i j∈Ni α ij S ij (v j (t)) (4)
where N i is the index set of presynaptic neurons, α ij > 0 is the adjacency weight, and S ij (•) : R → [0, 1] is the coupling function of sigmoidal type (S ij (•) is non-decreasing). The constant g i controls the coupling strength:

Assumption S1 (Weak coupling) The constant g i > 0 is chosen small enough to satisfy

g i j∈Ni α ij =: g i d i ≤ M i ,
where M i is the robustness margin in Assumption E3.

Our next assumption highlights the essential characteristics of excitatory synaptic coupling in the rapid and robust synchronization phenomenon: the coupling gain (determined by the derivative of the synaptic current) is high only near the threshold.

Assumption S2 (Localized coupling) For j ∈ N i , S ij (v) = 1, if v ≥ E th j , 0, if v ≤ max m∈[0,Mj ] v j (m),
where v j (m) is defined in Assumption E3.

As illustrated in the Introduction, under Assumptions E1-E3 and S1-S2, rapid convergence of each neuron to a neighborhood of its limit cycle attractor is guaranteed by treating the synaptic coupling as a weak exogenous disturbance (the first step of the convergence): every trajectory converges to D := D 1 ו • •×D N . Hence, individual neurons persistently spike regardless of the coupling. Next, we turn our attention to the network behavior of those spiking oscillators.

Recall that any directed graph G consists of multiple iS-CCs and their followers. The behavior of each iSCC is independent of the other iSCCs, due to the local interaction described in (4). Therefore, we focus our attention on a single iSCC. In particular, we assume the following. The behavior of the followers corresponds to the phenomenon of entrainment, rather than synchronization.

Assumption G1 (Connectivity) The communication graph G = (N , E) determined by the neighborhoods N i for i ∈ N is strongly connected.

A network-wise phase reduction

The analysis of the network relies on its three-time scale structure. Our biophysical assumptions (Assumptions E1-E3 and S1-S2) illustrated in the previous subsections ensure that, for the = 0 limit, a trajectory that starts from the slow region

S := {col(x 1 , v sm 1 (x 1 ), . . . , x N , v sm N (x N )) : x i ∈ [x i (0), x i (M i )], i ∈ N } , (5)
where

v sm i (x) denotes the smallest root of f i (x, v) = 0 (i.e., min{v : x 0 i (v) = x})
, returns back to the slow region S after an ultrafast-fast transition of spiking.

Note that the network is completely decoupled in the slow region S:

ẋi = -x i + x i,∞ (v sm i (x i )) τi (v sm i (x i )) =: h i (x i ), v i = v sm i (x i ), i ∈ N , (6)
each neuron independently travels to the left along the lower branch according to (6), until one of them reaches the left knee (x i (0), v i (0)). The set of all such endpoints is

P sp := {col(x 1 , v 1 , . . . , x N , v N ) ∈ S : ∃i s.t. x i = x i (0)} . ( 7
)
Any point in the set P sp is a jump point where the flow must be continued as a solution of the ultrafast-fast behavior. Assumptions E1-E3 and S1-S2 ensure that the ultrafast-fast behavior that starts from a jump point in P sp eventually returns to the region S, and hence further back to another jump point of P sp via another round of slow travel. The spiking map from P sp to S is denoted by

X + (x 1 , . . . , x N ) := col(x + 1 , . . . , x + N ) : P sp → S. (8)
Note that the mapping x + i (x 1 , . . . , x N ) only uses x j , j ∈ N as their argument for simplicity of notation, as v i is entirely determined by the relation v sm i (x i ), in the slow region S.

Singular perturbation analysis (see Appendix A for details of the argument) rests on the slow behavior in the singular limit = 0. This behavior is a pulse-coupled integrate-and-fire model of the type studied in (Mirollo & Strogatz, 1990b): the neurons follow the decoupled integrating flow (6) between instantaneous spikes that define the spiking times. This singular behavior is studied in the next section.

Integrate-and-fire model

The integrate-and-fire model is represented as

ẋi = h i (x i ) < 0, x i ∈ [x i , x i ], x + i = x + i (x 1 , . . . , x N ) ∈ [x i , x i ], ∃j s.t. x j = x j , i ∈ N := {1, . . . , N }. (9)
The state-space of this behavior is the set S:

S := {col(x 1 , . . . , x N ) : x i ∈ [x i , x i ], i ∈ N }.
The integrate-and-fire model has the following properties with

x mid i , x mid i ∈ (x i , x i ) such that x mid i < x mid i , i ∈ N .
This integrate-and-fire model corresponds to the singular limit = 0 of the network dynamics. The goal of the two following sections is to establish that this model possesses a hyperbolic fixed point that corresponds to synchronous firing.

Property 1 (A network-wise lower threshold) When neuron j ∈ N reaches a lower threshold x j , the network spikes: there exists a nonempty index set N F (x 1 , . . . , x N ) ⊆ N of synchronously firing neurons. In particular, j ∈ N F (x 1 , . . . , x N ).

Let us denote the set of all points where a network-wise spiking happens (uniformly) as P sp :

P sp := {col(x 1 , . . . , x N ) ∈ S : ∃j s.t. x j = x j }.
Then, the return map is well-defined as

X ++ (X + (p)) : P sp → P sp , where X ++ (p + ) := col(x ++ 1 (p + ), • • • , x ++ N (p + )), x ++ i (x + 1 , . . . , x + N ) := τ -1 i τ i (x + i ) -min j∈N τ j (x + j ) , τ i (x i ) := x i xi 1 h i (x) dx ≥ 0, x i ∈ [x i , x i ],
for i ∈ N , and X + (p) := col(x + 1 (p), . . . , x + N (p)). Note that τ i (x i ) represents the time of travel of neuron i from x i to the lower threshold x i .

Property 2 (Local interaction) Neuron i with x i > x i spikes only when a presynaptic neuron spikes and

x i ≤ x mid i . In particular, i ∈ N F (x 1 , . . . , x N ) only if there exists a set of neurons i 0 , . . . , i l-1 ∈ N F (x 1 , . . . , x N ) such that x i0 = x i0 and i k-1 ∈ N i k for all k = 1, . . . , l with i l = i (chain reaction), where N i ⊂ N is the index set of presynaptic neurons.
On the other hand, whenever a presynaptic neuron spikes, neuron i with x i > x i synchronously spikes provided that it is sufficiently near the threshold x i .

We define the set of network synchrony as

P syn := {p ∈ P sp : N F (p) = N }.
The connectivity of the graph determined by the neighborhoods N i for i ∈ N (Assumption G1) ensures by Property 2 that for each j ∈ N , there exists a point col(x 1 , . . . , x N ) inside the closure P syn of P syn such that x j = x j and x i > x i for all i = j.

Property 3 (Convexity) For any two points p = col(x 1 , . . . , x N ), p = col(x 1 , . . . , x N ) ∈ P sp such that

x i ≤ x i , ∀i ∈ N , we have N F (x 1 , . . . , x N ) ⊇ N F (x 1 , . . . , x N ).
In particular, for any monotonically increasing (i.e., at least one is strictly increasing) path that connects p and p , the points p where N F (p ) changes are isolated.

By Property 3, with a finite subset B of B syn , P syn can be approximated as

p=col(x1,...,x N )∈B x j(p) ≡ x j(p) × i =j(p) [x i , x i ) ⊆ p=col(x1,...,x N )∈B syn {x j(p) } × i =j(p) [x i , x i ) = P syn , ( 10
)
where B syn is the set of all points p = col(x 1 , . . . , x N ) on the boundary of P syn such that there exists an index j(p) such that x i = x i if and only if i = j(p). The equality holds except for a measure zero set. Note that for (10) to be a good approximation, for each j ∈ N , there exists at least one point p ∈ B such that j(p) = j.

Finally, regarding the spiking map x + i (•, . . . , •), i ∈ N , we have the following properties.

Property 4 (Localized coupling) If neuron i spikes, then

x + i (x 1 , . . . , x N ) ∈ [x mid i , x i ].
On the other hand, if neuron i does not spike, then (η ≥ 0)

x + i (x 1 , . . . , x N ) ∈ [max{x i -η, x i }, x i ].
If, in addition,

N i ∩ N F (x 1 , . . . , x N ) = ∅, then x + i (x 1 , . . . , x N ) = x i .
Property 4 illustrates that the network dynamics is almost decoupled away from the synchronous spike times.

Property 5 (Stable spiking) The spiking map is continuous and differentiable almost everywhere in each region where N F (x 1 , . . . , x N ) is invariant. A network-wise spiking is stable: there exist functions H i : [x i , x i ] → R >0 , i ∈ N and r ≥ 1 such that when the spiking map is differentiable, for any infinitesimal perturbation col(δx 1 , . . . , δx N ) that preserves N F (x 1 , . . . , x N ) and an index j such that x j = x j (so to remain in P sp ), we have

max i∈N F (p) 1 H i (x + i (p)) δx + i -min i∈N F (p) 1 H i (x + i (p)) δx + i ≤ r max i∈N F (p) 1 H i (x i ) δx i -min i∈N F (p) 1 H i (x i ) δx i (11) and δx + i = δx i , i / ∈ N F (p)
, where p := col(x 1 , . . . , x N ).

Properties 1-5 are the key properties enabling synchrony in earlier integrate-and-fire models of the literature such as those studied in (Mirollo & Strogatz, 1990b;[START_REF] Mauroy | Clustering behaviors in networks of integrate-and-fire oscillators[END_REF][START_REF] Mauroy | On the dichotomic collective behaviors of large populations of pulse-coupled firing oscillators[END_REF]. In the rest of the section, we will show that they also hold for the spiking map X + (•, . . . , •) defined in (8), where x i , x i , x mid i , and andx i (0), respectively. 4.1 Absorption: synchronous spiking Property 1 has already been established in Section 3.3. So, let us detail more on what happens after when the trajectory leaves the slow region S from P sp . This corresponds to the ultrafast-fast behavior as explained in Section 3.3. In detail, a trajectory that leaves from P sp follows the ultrafast behavior

x mid i are x i (0), x i (M i ), x i (M i ),
x i = 0, v i = f i (x i , v i ) + g i (E i -v i ) j∈Ni α ij S ij (v j ), i ∈ N , (12)
in the ultrafast region UF ⊂ D \ S (where the ultrafast time scale dominates), which is defined as the region where the vector field of the network ( 12) is non-zero. This corresponds to the phenomenon of absorption: firing neurons absorb neighboring neurons that are sufficiently near the threshold, resulting in synchronous firing (Mirollo & Strogatz, 1990b;[START_REF] Mauroy | Clustering behaviors in networks of integrate-and-fire oscillators[END_REF]. More precisely, a set N rt of neurons i at their left knee (x i (0), v i (0)) spike, i.e., go over the threshold

E th i (because v i > 0 for all v i ∈ (v i (0), E th i ]), resulting in S * i (v i ) = 1,
which yields a shift of the nullcline of other neurons (neuron j such that i ∈ N j ) into the right, yielding another spike, all in the same ultrafast time scale.

In particular, for each synchronous spiking, there exists a unique subgraph T = (N T , E T ) of G with the hierarchical structure, such that a spike of the 0-th layer N T 0 = N rt (i.e., set of all roots) yields a spike of the next layer N T 1 , and so on to all the neurons in N T ⊆ N , where

N T l+1 := i∈N T l {j : i ∈ N j,T } \ N T l , l ≥ 0,
N j,T is the set of neighbors associated with the subgraph T , and

N T l := ∪ l m=0 N T m .
For point col(x 1 , . . . , x N ) ∈ P sp , the corresponding subgraph can be found as follows.

1) Let l = 1 and N T 0 = {j :

x j = x j (0)}. 2) Let N T l = {j : x j < x j (g j d T j )} \ N T l-1 and E T l = ∪ k∈N T l ((N k ∩ N T l-1 ) × {k}), where d T j := k∈Nj ∩N T l-1 α jk .
3) If N T l = ∅, then let l = l+1 and repeat 2). Otherwise, let l * = l-1 and define T as (N T , E T ), where

N T := N T l * and E T := ∪ l * l=1 E T l .
Note that this means

x j ∈ X T j :=    x j g j d T j , x j g j d T j , if j ∈ N T 0 , x j g j d T j , x j g j d T j , otherwise, where d T j := k∈Nj ∩N T l-1
α jk for neuron j in the (l +1)th layer and

N T -1 = N T -2 = ∅.
Property 2 (Fast threshold modulation) Neuron i with x i > x i spikes if and only if there exists a subgraph T = (N T , E T ) of G such that i ∈ N T and

x k ∈ X T k ⊂ [x k , x mid k ], ∀k ∈ N T .
In particular, for each p ∈ P sp , there exists such a subgraph T satisfying

N T = N F (p).
This stronger property ensures that a collection B of explicit points indexed by a spanning tree T of G as

p B T := col(x 1 (g 1 d T 1 ), . . . , x N (g N d T N )) ∈ B syn ,
exactly characterizes P syn , according to (10). In particular,

T of G x j(p B T ) × i =j(p B T ) x i , x i (g i d T i ) = P syn (13)
holds except for a measure zero set.

We emphasize that the nodal excitability (as illustrated at the end of Section 3.1) and the weak localized synaptic coupling (as emphasized in Section 3.2) are essential to provide this open set of synchronous spiking: a strong contraction for synchronous spiking behavior, which provides both the rapid convergence and the robustness with respect to heterogeneity (see Sections 5 and 6). This is akin to the mechanism of fast threshold modulation: the coupling acts mainly to create an open transitional region.

Remark 1 If we linearize the ultrafast behavior (12), which corresponds to a synchronous spiking, then we observe diffusive coupling with a time-varying edge-wise coupling gain k ij (•) that is temporarily strong:

δv i = ∂f i ∂x (x i (t), v i (t))δx i + a i (t)δv i + j∈Ni k ij (t)[δv j -δv i ],
where

k ij (t) := g i (E i -v i (t))α ij S ij (v j (t)) and a i (t) := ∂fi ∂v (x i (t), v i (t))-g i j∈Ni α ij S ij (v j (t))+ j∈Ni k ij (t).
This strong pulling in the ultrafast time scale is what provides a synchronous spiking of the next layer.

Property 3 follows from the construction of the subgraph

T for point p ∈ P sp . In particular, N T l ⊇ N T l
for all l ≥ 0, where the subgraph T corresponds to the point p ∈ P sp . In particular, N F (x 1 , . . . , x N ) changes only when

x j = x j (g j d T j ) > x j (0) or x j = x j (g j d T j ) for some j / ∈ N T 0 .
Now, since the ultrafast behavior (12) preserves the value x i during its entire transition, a trajectory that starts from P sp ends up (via ultrafast behavior) in the fast region

F (that is S → U F → F), which is the set of all points col(x 1 , v 1 , . . . , x N , v N ) ∈ D \ (S ∪ UF) such
that there exists a nonempty index set N F ⊆ N satisfying

x i ∈ [x i (0), x i (m i )] ⊃ [x i (g i d T i ), x i (g i d T i )], if i ∈ N F , [x i (m i ), x i (M i )], if i / ∈ N F , v i = v lg i (m i , x i ), if i ∈ N F , v sm i (m i , x i ), if i / ∈ N F , m i = g i j∈Ni∩N F α ij , i ∈ N , where v lg i (m, x) (v sm i (m, x)) denotes the largest (small- est) root of f m i (x, v) = 0 (i.e., max{v : x m i (v) = x} (min{v : x m i (v) = x})) for x ≤ x i (m) (x ≥ x i (m)). Note that v sm i (•) = v sm i (0, •).
In the next subsection, we will detail on the fast behavior in the fast region F, where the fast time scale q dominates. In doing so, we will more detail on the exact spiking map x + i (•, . . . , •), i ∈ N defined in (8), and its properties.

Spiking map

To establish Property 4, we first note that the network is again completely decoupled in the fast region F:

x i = λi (v lg i (m i , x i )) -x i + x i,∞ (v lg i (m i , x i )) =: H mi i (x i ), (v i = v lg i (m i , x i )), i ∈ N F , x i = 0, (v i = v sm i (m i , x i )), i / ∈ N F , (14)
each neuron i ∈ N F independently travels to the right along the upper branch according to ( 14), until one of them reaches the right knee (x i (m i ), v i (m i )). Any such point is again a jump point where the flow must be continued as a solution of the ultrafast behavior (12).

Similar to the mechanism of synchronous spiking illustrated in the previous subsection, this yields that a subset N F ⊆ N F of neurons jumps to their lower branch while preserving their x i value, hence

x i ∈ [x i (0), x i (m i )].
If there are still some neurons remaining in the upper branch, i.e., N F \ N F = ∅, then it means that we again enter the fast region F, but with a different index set N F \ N F . This will change m i , hence neurons remaining in the upper branch have a small jump before it moves again along the fast upper branch according to ( 14).

After possibly a few of these iterations (F → UF → F), all the neurons in the upper branch will eventually jump to their lower branch (F → UF → S), which completes the mechanism of spiking. Note that the x i value of the neurons that previously jumped to their lower branch is preserved. The number of these iterations

(F → UF → F) is bounded by |N F | -1 ≤ N -1.
Property 4 (Localized coupling) If neuron i spikes, then

x + i (x 1 , . . . , x N ) ∈ [x i (0), x i (g i d i )] ⊂ [x mid i , x i ].
If neuron i does not spike, then x + i (x 1 , . . . , x N ) = x i .

Next, consider the propagation of an infinitesimal perturbation col(δx 1 , . . . , δx N ) (that preserves N F (x 1 , . . . , x N ) and an index j such that x j = x j ) through one iteration of F → U F → F or S. For this purpose, assume that a subset N F ⊆ N F of neurons jumps to their lower branch after this iteration. Note that this iteration is differentiable if N F is preserved in an open neighborhood, and is continuous even if N F changes.

Since the fast behavior just corresponds to the shift of time spent on the fast upper branch, our infinitesimal perturbation can be redefined as

δT i := δx i /H mi i (x i ), i ∈ N F
and accordingly

δT + i := δx + i /H mi i (x + i ), i ∈ N F ,
for simplicity of representation, where x + i and δx + i are the corresponding values after the iteration. Then, we have

δT + k -δT k = δT + i -δT i , ∀i, k ∈ N F , because x + k x k 1 H m k k (x) dx = x + i xi 1 H mi i (x) dx, ∀i, k ∈ N F .
Note that there exists an index k ∈ N F such that δx + k = 0, as

x + k = x k (m k ) (neuron k reaches the right knee), hence max i∈N F δT + i -min i∈N F δT + i ≤ max i∈N F δT i -min i∈N F δT i .
Note that between two iterations, we have to multiply

H mi i (x + i )/H m + i i (x + i )
to δT + i for neurons remaining in the upper branch, i ∈ N F \ N F , as they experience a small jump before flowing again along the fast upper branch. This establishes Property 5.

Property 5 (Stable spiking) Property 5 holds with

r := (r F →F ) N +1 ≥ 1, r F →F := max i∈N max m,m ∈[0,Mi] max x∈[x i ,xi(min{m,m })] H m i (x) H m i (x)
,

and H i (•) := H Mi i (•).
Almost any trajectory that starts from the rest of the region D \ (S ∪ UF ∪ F) moves to S in a finite time, hence never going back to the rest of the region (see Appendix B).

Remark 2 This subsection shows the importance of a network-wise reduction for nonstandard singularly perturbed systems (nonstandard as in (Wechselberger, 2020)) that represent networks of neurons, instead of a neuron-wise reduction. In particular, all the important interaction via weak localized synaptic coupling takes place in the ultrafast-fast time scale, making it difficult to characterize the network-wise spiking map X + (•) via phase response curves.

5 Main result: integrate-and-fire model

It is known that synchronous spiking is not necessarily persistent in the presence of heterogeneity and a general network topology (Mirollo & Strogatz, 1990b;[START_REF] Mauroy | Clustering behaviors in networks of integrate-and-fire oscillators[END_REF]. A synchronous spike of the entire network does not necessarily lead to a synchronous spike of the entire network after a return.

The goal of the present section is to detail the conditions that make perfect network synchrony stable and robust.

Condition for the local stability of network synchrony

Assume for a moment that there exists a fixed point p * of the return map, i.e., X ++ (X + (p * )) = p * , in the interior of P syn so that a synchronous spiking that starts from p * returns back to p * . Then, the local stability of the fixed point can be guaranteed by a contraction property, that is, the existence of c S→F > 0 and c F →S > 0 such that

|h i (x * i )| ≤ c S→F H i (x * i ), H i (x + i (p * )) ≤ c F →S |h i (x + i (p * ))|, ∀i ∈ N , c S→F c F →S < 1/r ≤ 1, ( 15
)
where col(x * 1 , . . . , x * N ) := p * ∈ P syn , and the constant r ≥ 1 and the function H i (•) are given in Property 5. Given that

x * i ≤ x mid i < x mid i ≤ x + i (p * ), this is merely equivalent to |h i (•)|/H i (•)
being monotone increasing in the mean:

|h i (x + i (p * ))| H i (x + i (p * )) - |h i (x * i )| H i (x * i ) ≥ 1 c F →S -c S→F > 0.
Theorem 1 (Local stability) Given the integrateand-fire model (9) with Properties 1-5, in addition to Assumption G1 assume that there is a fixed point p * of the return map X ++ (X + (•)) in the interior of P syn that satisfies (15) with some c S→F , c F →S . Then, it is locally exponentially stable with the convergence rate c S→F rc F →S < 1.

Remark 3 For this fixed point to have an open neighborhood of persistent synchronous spiking, there must exist a point p ∈ P syn such that x i > x i for all neurons except one say 1, i.e., x 1 = x 1 . Then, by Property 2, this implies that there exists a spanning tree having a root 1. So, as illustrated in Section 3.2, we need connectivity of G.

PROOF. Since there is an open neighborhood of the fixed point p * =: col(x * 1 , . . . , x * N ), its local exponential stability can be determined by the linearization of the return map. In particular, since by Property 5, the spiking map is differentiable almost everywhere in P syn , and for any infinitesimal perturbation col(δx 1 , . . . , δx N ) that preserves N F (p * ) = N and an index j such that x j = x j (hence δx j = 0), we have (11) with N F = N and p = p * . Now, analogously as in Section 4.2, X ++ (•) just corresponds to the shift of time spent, hence let us redefine our infinitesimal perturbation as

δτ + i := 1 h i (x + i (p * ))
δx + i and accordingly

δτ ++ i := 1 h i (x ++ i (X + (p * ))) δx ++ i and δτ i := 1 h i (x * i ) δx i ,
for simplicity of representation. Then, we have

δτ ++ k -δτ + k = δτ ++ i -δτ + i , ∀i, k ∈ N , because x ++ k x + k 1 h k (x) dx = x ++ i x + i 1 h i (x) dx, ∀i, k ∈ N .
Since p * is a fixed point, the infinitesimal perturbation col(δx + 1 , . . . , δx + N ) again preserves the index j such that

x ++ j (X + (p * )) = x * j = x j , hence δτ ++ j = δx ++ j = 0, hence max i∈N δτ ++ i -min i∈N δτ ++ i ≤ max i∈N δτ + i -min i∈N τ + i ,
and thus, by the stability of the spiking (Property 5),

max i∈N δτ ++ i -min i∈N δτ ++ i ≤ c F →S max i∈N 1 H i (x + i (p * )) δx + i -min i∈N 1 H i (x + i (p * )) δx + i ≤ c F →S r max i∈N 1 H i (x * i ) δx i -min i∈N 1 H i (x * i ) δx i ≤ c F →S rc S→F max i∈N δτ i -min i∈N δτ i . (16)
This concludes the proof.

The local stability of a fixed point p * , which is in the interior of P syn , ensures the existence of a neighborhood P pos which is stable, that is P pos is a positively invariant set in P syn : any synchronous spiking that starts from P pos returns back to P pos . In particular, there exists P pos ⊆ P syn , c S→F > 0, and c F →S > 0 that satisfies (15) for all p ∈ P pos instead of just p * . Note that by Properties 2 and 4, the region that corresponds to x i and the region that corresponds to x + i (p) are separated, hence (15) is well justified. In particular, we have

x i ≤ x i ≤ x mid i < x mid i ≤ x + i (x 1 , . . . , x N ) ≤ x i ,
for any col(x 1 , . . . , x N ) ∈ P syn . In this sense, (15) for all p ∈ P pos can be guaranteed by the following assumption.

Assumption M1 (Monotonicity in the mean)

There exists c S→F , c F →S > 0 such that

max x∈[x i ,x mid i ] |h i (x)| H i (x) ≤ c S→F , max x∈[x mid i ,xi] H i (x) |h i (x)| ≤ c F →S , ∀i ∈ N ,
and c S→F c F →S < 1/r, where the constant r ≥ 1 and the function H i (•) are given in Property 5.

A simple sufficient condition for Assumption M1 is the positivity of the derivative: that there exists a sufficiently large λ > 0 such that d dx

|h i (x)| H i (x) ≥ λ > 0, ∀x ∈ (x i , x i ), ∀i ∈ N .
Under this stability condition (Assumption M1), in particular, we can show contraction inside P syn : the contraction of the metric defined for any two points p = col(x 1 , . . . , x N ) and p = col(x 1 , . . . , x N ) in P sp :

d(p, p ) := max i∈N [τ i (x i ) -τ i (x i )] -min i∈N [τ i (x i ) -τ i (x i )].
As in (Somers & Kopell, 1993), the notion of phase difference depends on time difference along the trajectory, not differences in position in phase space.

Theorem 2 (Contraction inside P syn ) Given the integrate-and-fire model (9) with Properties 1-5, let Assumptions G1 and M1 holds. Then, we have

d(X ++ (X + (p)), X ++ (X + (p ))) ≤ c S→F rc F →S d(p, p ), (17)
for all p, p ∈ P syn .

PROOF. From the argument in the proof of Theorem 1, we can analogously conclude ( 16) for any point along the path that connects p and p inside P syn for any eligible infinitesimal perturbation col(δx 1 , . . . , δx N ). Note that such a path exists, as P syn is convex, from Property 3. Then, by integrating this along the path, we can conclude (17). In particular, we can construct a path p(s) = col(x 1 (s), . . . , x N (s)) : [0, 2] → P syn as follows.

(1) Let p(0) = p.

(2) Let N ⊆ N be the set of all the indices i such that τ i (x i ) ≤ τ i (x i ) and let τ := max i∈N [τ i (x i )τ i (x i )] ≥ 0. Then, for each s ∈ [0, 1], we define p(s)

as x i (s) =              x i , if i / ∈ N , τ -1 i (τ i (x i ) -sτ ), if i ∈ N and τ i (x i ) -τ i (x i ) > sτ , x i , if i ∈ N and τ i (x i ) -τ i (x i ) ≤ sτ .
(3) Finally, let τ := -min i / ∈N [τ i (x i )-τ i (x i )] ≥ 0, and for each s ∈ [1, 2], define p(s) as

x i (s) =              x i , if i ∈ N , τ -1 i (τ i (x i ) + (s -1)τ ), if i / ∈ N and τ i (x i ) -τ i (x i ) > (s -1)τ , x i , if i / ∈ N and τ i (x i ) -τ i (x i ) ≤ (s -1)τ .
By its construction, for s ∈ [0, 1], δτ i (s) ∈ {0, τ }, i ∈ N . Similarly, we have for s ∈ [1, 2], δτ i (s) ∈ {0, -τ }, i ∈ N . This now implies from ( 16) that

max i∈N δτ ++ i (s) -min i∈N δτ ++ i (s) ≤ c S→F rc F →S τ , if s ∈ [0, 1], c S→F rc F →S τ , if s ∈ [1, 2].
Therefore, we can conclude that

d(X ++ (X + (p)), X ++ (X + (p ))) = max i∈N [τ i (x ++ i (X + (p(0)))) -τ i (x ++ i (X + (p(2))))] -min i∈N [τ i (x ++ i (X + (p(0)))) -τ i (x ++ i (X + (p(2))))] ≤ 2 0 D + max i∈N [τ i (x ++ i (X + (p(0)))) -τ i (x ++ i (X + (p(s))))]ds - 2 0 D + min i∈N [τ i (x ++ i (X + (p(0)))) -τ i (x ++ i (X + (p(s))))]ds ≤ 2 0 max i∈N δτ ++ i (s) -min i∈N δτ ++ i (s)ds ≤ c S→F rc F →S [τ + τ ] = c S→F rc F →S d(p, p ),
where D + (D + ) denotes the upper (lower) right-hand Dini derivative.

Theorem 2 removes the possibility of periodic solutions that reside inside P syn , and therefore, we will continuously focus our attention on a fixed point. In this respect, Theorem 2 further implies that under Assumption M1, the existence of a fixed point in the interior of P syn is equivalent to the existence of an open positively invariant set in P syn . This positively invariant set becomes the domain of attraction due again to Theorem 2. An explicit sufficient condition for the existence of a fixed point in this context is given in the next subsection.

Existence of a fixed point p * in the interior of P syn

According to the illustration of the integrate-and-fire model ( 9) in Section 4, the only explicit knowledge that we have on the point after synchronous spiking of the entire network is that it is inside the set

i∈N [x mid i , x i ].
Therefore, the only explicit sufficient condition for the existence of a fixed point p * ∈ P syn , or equivalently, the existence of a positively invariant set P pos ⊆ P syn that we can conclude based on Properties 1-5 should be X ++ (p) ∈ P pos for all p ∈ i∈N [x mid i , x i ] with some P pos ⊆ P syn . Thus, the condition will be most relaxed for P pos = P syn .

In this regard, the most relaxed explicit sufficient condition for the existence of a fixed point p * ∈ P syn that we can conclude based on Properties 1-5 is

X ++ (p) ∈ p =col(x 1 ,...,x N )∈B {x j(p ) } × i =j(p ) [x i , x i ), ∀p ∈ i∈N [x mid i , x i ], (18)
which is based on the approximation of P syn given in (10). In particular, under this sufficient condition, there exists a unique fixed point p * in the interior of P syn . Now, (18) means that for each p ∈ i∈N [x mid i , x i ], if x ++ j (p) = x j , then there exists a point p = col(x 1 , . . . , x N ) ∈ B such that j(p ) = j and

x i ≤ x ++ i (p) < x i , ∀i = j.
In other words, τ j (x j ) = min i∈N τ i (x i ) and

0 = τ i (x i ) ≤ τ i (x i ) -τ j (x j ) < τ i (x i ), ∀i = j, ( 19
)
where p =: col(x 1 , . . . , x N ). Now, when x j becomes smaller as x mid j and x i become larger as x i for all i = j, we still have x ++ j (p) = x j , and this implies

∃p = col(x 1 , . . . , x N ) ∈ B s.t. j(p ) = j and τ i (x i ) -τ i (x i ) < τ j (x mid j ), ∀i = j, (20)
and τ j (x mid j ) ≤ τ i (x i ) for all i = j. Therefore, ( 18) is equivalent to (20) for all j ∈ N fast ⊆ N , where

N fast := {j ∈ N : τ j (x mid j ) ≤ min i∈N τ i (x i )}, (21)
because then for any p = col(x 1 , . . . , x N ) ∈ i∈N [x mid i , x i ] there exists j ∈ N fast satisfying x ++ j (p) = x j , and thus,

τ i (x i ) -τ i (x i ) ≤ τ i (x i ) -τ i (x i ) < τ j (x mid j ) ≤ τ j (x j ) ≤ τ i (x i ), ∀i = j.
A simple sufficient condition is (20) for all j ∈ N . In particular, under this sufficient condition, there exists a fixed point p * inside j∈N

{x j } × i =j [x i , x F i (j)) =: P pos ⊆ P syn (22)
with some x F i (j) > x i for i = j.

Note that this condition is robust to heterogeneity in the sense that one can preserve (20) under a slight perturbation on the values x i , x mid i , x i , and the function h i . Moreover, note that this robustness becomes stronger with respect to the monotonicity (Assumption M1), as when H i (•) is fixed, by decreasing c S→F or c F →S , the speed of travel right after the spike,

|h i (x)| at x ∈ [x mid i , x i ],
becomes faster relative to the speed of travel right before the spike,

|h i (x)| at x ∈ [x i , x i ) ⊂ [x i , x mid i ]. In- deed, if individual neurons are identical, in particu- lar, h(•) ≡ h i (•), x ≡ x i , x ≡ x i , x mid ≡ x mid i
, and x mid ≡ x mid i for all i ∈ N , then the sufficient condition (20) simply becomes

τ (x) -τ (x F ) < τ (x mid ),
where x F := min j∈N max p ∈B s.t. j(p )=j min i =j x i . This is equivalent to τ (x) -τ (x mid ) < τ (x F ), or further to

x mid x 1 h(x) dx < x x F 1 h(x) dx.
So, if the function h(•) is monotone in the mean:

max x∈[x,x F ] |h(x)| ≤ max x∈[x,x mid ] |h(x)| ≤ c min x∈[x mid ,x] |h(x)|,
then we only need the monotonicity to be strong enough:

c < x F -x x -x mid .
Now, in the next subsection, we emphasize the important role of the existence of an open neighborhood around the fixed point in leading to rapid global convergence.

Rapid global convergence

We first emphasize that under the sufficient condition given in the previous subsection, P syn also corresponds to a network behavior of persistent synchronous spiking. Now, this network behavior becomes globally stable if there is strong monotonicity: the speed of travel is sufficiently fast outside the neighborhood of the lower threshold. This is intuitive as the neurons that have not spiked yet 'wait' (by moving very slowly) for the neurons that have spiked, which significantly reduces the phase difference, eventually ending up synchronously firing.

Theorem 3 (Stability inside P sp ) Given the integrate-and-fire model (9) with Properties 1-5, in addition to Assumptions G1 and M1 assume that there exists a positively invariant set P pos of the form (22) with x F i (j) > x i , i = j (e.g., by assuming (20) for all j), and

• (Sufficient monotonicity) there is r i ∈ (0, 1) so that

τ i (x i ) -τ i ((1 -r i )x i + r i x F i ) < τ j (x mid j ), ∀i = j,
where

x F i := min j =i x F i (j), hence the fixed point p * = col(x * 1 , . . . , x * N ) satisfies x * i -x i < r i (x F i -x i )
for all i ∈ N , and there exists r η ≥ 0 and η * ∈ [0, min i∈N (1-

r i )[x F i -x i ]) so that (N -1)D P sp < (1 -c S→F rc F →S )d P syn -r η (d P sp -d P syn ),
where

D P sp := 2 max i∈N [τ i (x i )-τ i (x F i -η * )], d P syn := min i∈N [τ i (x F i -η * )-τ i ((1-r i )x i +r i x F i )], and
d P sp := max i∈N [τ i (x i ) -τ i (x i )]. 1
Then, for sufficiently small η: η

∈ [0, η * ] that satisfies 2 max i∈N max xi∈[x i ,xi] max x + i ∈[max{xi-η,x i },xi] h i (x i ) h i (x + i ) ≤ 1 + r η , (23)
any trajectory that starts from P sp enters the positively invariant set P syn of persistent synchronous spiking within a couple of cycles (rapid convergence).

PROOF. By ( 23), analogous to the argument in the proof of Theorem 1, we can still conclude

max i∈N δτ ++ i -min i∈N δτ ++ i ≤ (1 + r η ) max i∈N δτ i -min i∈N δτ i
at point p = col(x 1 , . . . , x N ) in P sp \ P syn , for any infinitesimal perturbation col(δx 1 , . . . , δx N ) that preserves N F (p) and an index j such that x j = x j .

This means that if there is a path that connects p and p inside P sp such that on the path, we preserve N F (p) and an index j such that x j = x j , then we can conclude

d(X ++ (X + (p)), X ++ (X + (p ))) ≤ (1 + r η )d(p, p ),
analogous to Theorem 2. Now, when this does not apply is when the path passes a point p where the index set of synchronously firing neurons changes. However, we still know that

lim p+→p + lim p-→p - d(X ++ (X + (p + )), X ++ (X + (p -)))
is bounded by a positive constant D P sp . In particular,

lim p+→p + lim p-→p - d(X ++ (X + (p + )), X ++ (X + (p -))) = lim p+→p + lim p-→p - max i∈N [τ i (x + i (p + )) -τ i (x + i (p -))] -lim p+→p + lim p-→p - min i∈N [τ i (x + i (p + )) -τ i (x + i (p -))] ≤ 2 max i∈N F (p+) N F (p-) [τ i (x i ) -τ i (x F i -η)] ≤ D P sp ,
where N + N -denotes the symmetric difference between N + and N -, i.e., (

N + \ N -) ∪ (N -\ N + ).
This finally implies that on the path that connects p * and p ∈ P sp \P syn , as in the proof of Theorem 2, for any η > 0, there exists p ∈ P syn such that d(p, p ) ≤ d(p * , p )d(p * , p) + η and d(p * , p) ≥ d P syn , which satisfies

d(p * , X ++ (X + (p ))) = d(X ++ (X + (p * )), X ++ (X + (p ))) ≤ d(X ++ (X + (p * )), X ++ (X + (p))) + d(X ++ (X + (p)), X ++ (X + (p ))) ≤ c S→F rc F →S d(p * , p) + (N -1)D P sp + (1 + r η )d(p, p ).
This is because, by Property 3, for a monotonically increasing (decreasing) path, the points where the index set of synchronously firing neurons changes are isolated, hence there are at most N -1 of them. By (Sufficient monotonicity), there exists a sufficiently small c ∈ (0, 1) such that

(N -1)D P sp < (1 + r η -c S→F rc F →S )d P syn -(c + r η )d P sp
and this implies with the former argument that

d(p * , X ++ (X + (p ))) ≤ -(1 + r η -c S→F rc F →S )d P syn + (N -1)D P sp + (1 + r η )d(p * , p ) + (1 + r η )η < (1 -c)d(p * , p ) + (1 + r η )η .
As this holds for any η > 0, we can conclude the proof. A significant difference from previous studies is our focus on the spiking map from P syn (local), contrary to the usual focus on the spiking map from P sp \ P syn (global). In particular, by departing from the assumptions of identical neurons and all-to-all coupling, we focused our attention on the behavior of persistent synchronous spiking itself, more precisely the existence of a hyperbolic fixed point. It is the robustness of this attractor in the singular limit that allows us to conclude from singular perturbation analysis away from the singular limit. As in previous studies [START_REF] Mauroy | On the dichotomic collective behaviors of large populations of pulse-coupled firing oscillators[END_REF], the monotonicity (in the mean) is found essential. In this subsection, we further remarked that global convergence property can be recovered by strong monotonicity.

The biophysical model illustrated in Section 3 is standard, and thus, the integrate-and-fire model that we proposed in Section 4 is well suited for the study of network behaviors in a neuronal network. For instance, it provides robustness to the phenomenon of absorptions, motivated by the biophysical model.

6 Main result: biophysical model

Discussion

(1) In the first part of Section 5, the assumption on the monotonicity in the mean (Assumption M1) has been shown to be sufficient for the local stability of the fixed point and the contraction inside P syn , which also gives uniqueness and the existence of a fixed point (given the positive invariance of P syn ). For the biophysical model of Section 3, h i (•) (H i (•)) in Assumption M1 is determined by the slow (fast) behavior on the lower (upper) branch. Thus, as first observed in (Somers & Kopell, 1993), the key property is the "rate of compression," which measures relative velocities of the slow variable just before and just after an ultrafastfast jump of the neuron, the parameters c S→F and c F →S in Assumption M1. In particular, Assumption M1 is the exact counterpart of the stability condition given in (Somers & Kopell, 1993), especially the weaker one mentioned in its Remark 2.5, while in (Somers & Kopell, 1993) systems with a twotime scale (slow=fast-ultrafast) are considered instead. In the proof of Theorem 1, we have explained how such a compression condition arises in an analytical treatment of simple linearization, even for heterogeneous networks. This property is shared by biophysical conductance-based models of neuronal systems (Somers & Kopell, 1993). (2) We emphasize that, for the biophysical model of Section 3, Assumption M1 is completely determined by the individual behavior of neurons: the shape of the individual nullcline, the time constant, and the robustness margin. In particular, a sufficient condition for Assumption M1 is that their slow trajectory has a "scalloped" shape (Somers & Kopell, 1993): that there exists a sufficiently small λ < 0 such that

h i (x) ≤ λ < 0, ∀x ∈ (x i (0), x i (M i )), ∀i ∈ N .
This means that the properties of the individual neurons determine the rate of convergence to synchrony, independent of the coupling strength and the network topology (however, they may affect the size or the positive invariance of P syn , see the next item).

(3) In Section 5.2, an explicit sufficient condition (20) for the existence of a fixed point has been found based on the approximation (10). According to Section 4.1, for the biophysical model illustrated in Section 3, especially by Property 2 , P syn can be exactly characterized with the collection B of points p B T indexed by a spanning tree T of G. Also, by Property 4 , x i can be relaxed as x i (g i d i ), and thus, the explicit sufficient condition (20) becomes ∃ spanning tree T of G having a root j s.t.

τ i (x i (g i d i )) -τ i (x i (g i d T i )) < τ j (x j (0)), ∀i = j, (24)
for all j ∈ N such that τ j (x j (0)) ≤ min i∈N τ i (x i (g i d i )). In particular, this is equivalent to whether for each j ∈ N the resulting set N T of the following algorithm is equal to the entire set N . (a) Let l = 1 and

N T 0 = {j}. (b) Let N T l = {i : τ i (x i (g i d i )) -τ j (x j (0)) < τ i (x i (g i d T i ))} \ N T l-1 , where d T i := k∈Ni∩N T l-1 α ik . (c) If N T l = ∅, then let l = l + 1 and repeat (b). Otherwise, let l * = l-1 and define N T := N T l * .
As noted in Section 5.2, this is easier to satisfy if the rate of compression is small.

Note that this can even be satisfied for arbitrarily weak coupling strength (when g i → 0), under the assumption that there exists T > 0 such that

T = τ i (x i (0)) -τ i (x i (0)) = τ i (x i (0)), ∀i ∈ N ,
meaning that all neurons should have an identical period for their decoupled limit cycle. This is also necessary. Under this additional condition, our sufficient condition (24) can be guaranteed for sufficiently small g coup (where g i = g coup g i ) if and only if ∃ spanning tree T of G having a root j s.t.

τ i (x i (0))x i (0)g i d i < τ i (x i (0))x i (0)g i d T i (25)
for all i = j and j ∈ N . Again this is easier to satisfy if the rate of compression is small. In fact, this is always satisfied for biophysical conductancebased models if the graph is sufficiently balanced, i.e., d i = min j =i min T ,root=j d T i =: d i , because, by monotonicity in the mean, we have

τ i (x i (0)) = - 1 h i (x i (0)) < - 1 h i (x i (0)) = τ i (x i (0))
and biophysical conductance-based models satisfy 3

x i (0) = E i -v i (0) -∂fi ∂x (x i (0), v i (0)) ≤ E i -v i (0) -∂fi ∂x (x i (0), v i (0)) = x i (0),
e.g., for conductance-based models as in Sections 2 and 7, we have (-∂f i /∂x)(x, v) = g i,K (0.7+v) > 0.

A similar argument shows that if the conditions on monotonicity and conductance-based models are extended for g ∈ (0, g coup ), then (24) is satisfied. An example of a sufficiently balanced graph with unitary weights, i.e., α ij = 1 for all (j, i) ∈ E, is the directed ring network, where N i = {i + 1} for all i = N and N N = {1} (d i = d i = 1). This conclusion supports the numerical observation at the end of Section 3.1 of (Somers & Kopell, 1993). Indeed, sparsity in connection reduces the difference between d i and d i , which further favors our sufficient condition (25). ( 4) We emphasize again that the strong contraction provided by the nodal excitability (given by Assumptions E1-E3) and the weak localized synaptic coupling (defined by Assumptions S1 and S2) is what gives rapid convergence by constructing

3 The equalities can be derived from the identities

(d/dm)f m i (x i (m), v i (m)) = (d/dm)f m i (xi(m), vi(m)) = 0 and (∂f m i /∂v)(x i (m), v i (m)) = (∂f m i /∂v)(xi(m), vi(m)) = 0.
an open set P syn of persistent synchronous spiking around the fixed point p * , and thus, maintains the convergence rate c S→F rc F →S even when the weak coupling strength g i is further reduced. These assumptions are robust to heterogeneity.

(5) Theorem 3 in Section 5.3, underlines the important role of the domain of attraction P syn having a volume in leading to rapid global convergence. In particular, the global convergence property is recovered by sufficient compression: the speed of slow travel is sufficiently fast outside the neighborhood of the left knee. This is consistent with the numerical observation in item 1 of Section 3.3 of (Somers & Kopell, 1993). ( 6) Now, the conclusion of Theorem 3 results in an almost global convergence towards the fixed point p * with the help of the first step of the convergence illustrated in the Introduction: rapid convergence (finite-time convergence for the = 0 limit) to P sp , which is obtained by preserving the biophysical nature of each individual neuron by weak coupling. See Appendix B for the complete argument. This leads to our main theorem illustrated in the next subsection.

Main theorem

Theorem 4 (Rapid and robust synchronization)

In addition to Assumptions E1-E3, S1-S2, G1, and M1, assume (24) for all j ∈ N , and

• (Sufficient compression) there is r i ∈ (0, 1) so that

τ i (x i (g i d i )) -τ i (x i (r i g i d i )) < τ j (x j (0)), ∀i = j,
hence the fixed point p * = col(x * 1 , . . . , x * N ) satisfies

x * i ∈ [x i (0), x i (r i g i d i
)) for all i ∈ N , and

(N -1)D P sp < (1 -c S→F rc F →S )d P syn ,
where

D P sp := 2 max i∈N [τ i (x i (g i d i )) -τ i (x i (g i d i ))] and d P syn := min i∈N [τ i (x i (g i d i )) -τ i (x i (r i g i d i ))].
Then, there exists * > 0 such that for each ∈ (0, * ), the network (3) with (4) has an almost semi-globally stable limit cycle, that corresponds to a persistent synchronous spiking network behavior, which has rapid convergence.

In particular,

(1) by preserving the biophysical nature of individual neurons by weak coupling (Assumption S1), the network rapidly converges to a neighborhood of S, i.e., each neuron rapidly converges to a neighbor of its limit cycle attractor. Then, after less than a cycle, it converges into a neighborhood of P sp . (2) After a couple of cycles, this results in a neighborhood of P syn , a positively invariant set that corresponds to a persistent synchronous spiking net-work behavior, as presynaptic neurons only affect the postsynaptic neuron when synchronously spiking (the third step of the convergence). (3) Finally, these network behaviors of persistent synchronous spiking converge to the unique limit cycle with a convergence rate that is independent of the coupling strength and the network topology (the second step of the convergence).

Theorem 4 follows from Theorems 2-3 and from the arguments in Sections 4.1, 4.2, 6.1, and Appendices A-B. This emphasizes that the network under study achieves rapid convergence to a persistent synchronous spiking network behavior that is robust to heterogeneity via weak synaptic coupling.

Remark 4 Theorem 4 proves the conjecture of Somers and Kopell in Section 3.4 of (Somers & Kopell, 1993).

Simulation

To illustrate Assumptions E1-E3, S1-S2, G1, M1, and especially the robustness of the existence of p * with respect to heterogeneity, we simulate a network consisting of the type of systems given as

τ (v i ) ẋi = -x i + x ∞ (v i ) vi = g i,L (-0.4 -v i ) + g i,Ca m ∞ (v i )(1 -v i ) + g i,K x i (-0.7 -v i ) + 0.4 + g coup (1 -v i ) j∈Ni S(v j )
where = 0.0004 and m ∞ (•) is given in Section 2, which resembles the standard Morris-Lecar excitable model in Section 2. Here, heterogeneous parameters are the maximal conductances g i,L , g i,Ca , and g i,K , as in the neuronal network. We take the parameters random with uniform distribution from the intervals [0.3, 0.75], [1.75, 2.25], and [2.75, 3.25], respectively. If x ∞ (•) is a sigmoidal function ranging in [0, 1], then we see that D i = [0, 1] × [-0.7, 1] is positively invariant for any m ≥ 0, as

f m i (x, v) = g i,L (-0.4 -v) + g i,Ca m ∞ (v)(1 -v) + g i,K x(-0.7 -v) + 0.4 + (1 -v)m < -1.4g i,L + 0.4 ≤ -0.02 < 0
for all x ≥ 0 and v > 1, while f m i (x, v) > 0.4 > 0 for all x ≥ 0 and v < -0.7.

The nullcline f

m i (x, v) = 0 has N -shape if [0.3g i,L + 1.7m + 0.4]/g i,Ca < 0.71 * 0.99 * m ∞ (0.01) -1.7m ∞ (0.01) = 1.4260, < 0.7 * 1 * m ∞ (0) -1.7m ∞ (0) = 1.4833.
This can be deduced by calculating the existence condition of two zero roots for the partial derivative of the nullcline f m i (x, v) = 0 with respect to v. In fact, this is the condition that guarantees (∂x m i /∂v)(0.01), (∂x m i /∂v)(0) > 0, where x m i (v) is given in Lemma 1. Thus, we take our robustness margin as Note that (∂x m i /∂v)(0.01), (∂x m i /∂v)(0) > 0 also tells us that the horizontal lines v = 0.01 and v = 0 are always located between the left knee and the right knee. Then, Assumptions E1 and E3 are satisfied if x ∞ (•) is 0 in the lower branch and 1 in the upper branch. This is because, we have

x i (M i ) < x Mi i (v := -0.1460) ≤ x 0 i (v := 0.0930) < x i (0),
where the second inequality holds if and only if

M i (1 -v )(0.7 + v ) ≤ -(0.3g i,L + 0.4)(v -v ) -g i,Ca m ∞ (v )(1 -v )(0.7 + v ) + g i,Ca m ∞ (v )(1 -v )(0.7 + v ). Moreover, 0 < x i (m) < x i (m) < 1 because f m i (0, v) ≥ min{0.3(-0.4 -v), 0.75(-0.4 -v)} + 1.75m ∞ (v)(1 -v) + 0.4 > 0, v ∈ [-0.7, 0], f m i (1, v) ≤ 0.3(-0.4 -v) + (2.25m ∞ (v) + 0.36)(1 -v) + 2.75(-0.7 -v) + 0.4 < 0, v ∈ [0.01, 1].
To satisfy the rest of the conditions in Assumptions E2, E3, and S2, we model τ (•), x ∞ (•), and S(•) as

x ∞ (v) = S(v) = 1, if v ≥ 0, 0, if v < 0, τ (v) =    q , if v ≥ E th , q + τ , if v ∈ [-0.258, E th ), q + τ , if v < -0.258,
with q = 1/2, E th = 0.01, and τ = 5 q . Note that the horizontal line v = -0.258 always has an intersection with the lower branch (we always have (∂x m i /∂v)(-0.258) < 0), hence the rate of compression can be approximated as

sup x∈[x i ,x m i (-0.258)) |h i (x)| inf x∈(x m i (-0.258),xi] |h i (x)| = τ τ = 5 τ .
We utilize a graph satisfying Assumption G1, where we selected 36 neighbors for each i as i + 1, . . . , i + 36 mod N , where the total number of neurons is N = 100. By selecting adjacency elements in {0, 1}, we get d i = 36 for all i ∈ N . Hence, we take our weak coupling strength as g coup = 0.01 to satisfy Assumption S1 with the robustness margin M i = 0.36. Then, we can see from the following simulation results that as the rate of compression c gets smaller: as τ > 0 gets larger, the conditions for the existence of p * and for the global convergence are further satisfied for the varying parameters g i,L , g i,Ca , and g i,K in a broad range.

Conclusion

We have illustrated how weak localized synaptic coupling combined with spiking oscillators can result in rapid spike synchronization in heterogeneous networks.

Traditionally, the analysis of synchronization in the vicinity of a limit cycle is achieved via a reduction to phase models (neuron-wise reductions). This reduction however fails in the singular limit of slow-fast attractors. We propose instead that direct linearization of a network-wise reduction of a state-space conductancebased model provides an elucidating connection to the classical linear theory of network synchronization by diffusive coupling. We showed the essential role of the "rate of compression" (monotonicity) for stability and robustness, which connects our analysis to the theory of synchronization in integrate-and-fire networks (Mirollo & Strogatz, 1990b;[START_REF] Mauroy | Clustering behaviors in networks of integrate-and-fire oscillators[END_REF][START_REF] Mauroy | On the dichotomic collective behaviors of large populations of pulse-coupled firing oscillators[END_REF]. 2020). Geometric singular perturbation theory beyond the standard form. Springer.

A A network-wise singular perturbation analysis

We first introduce a method to approximate singularly perturbed systems, in which the fast system is not uniformly stable, as in our network (3) with ( 4). The method is to first normalize the entire dynamics and then consider the = 0 limit. In particular, for our network (3) with (4), we get the parametrized normalized network as

ẋi = h i, (x i , vi ) n (x 1 , v1 , . . . , xN , vN ) vi = fi (x i , vi , {v j } Ni )/ n (x 1 , v1 , . . . , xN , vN ) , ∀i ∈ N (A.1)
where

h i, := -x i + x i,∞ (v i ) τ i, (v i ) , fi := f i (x i , vi ) + g i E i -vi j∈Ni α ij S ij (v j ), n := N i=1 h i, (x i , vi ) 2 + ( fi (x i , vi , {v j } Ni )/ ) 2 .
Now, by considering the = 0 limit, we obtain a discontinuous dynamical system given as ẋi

= h i (x i ) n S (x 1 , . . . ,
mi i (x i ) n F (x 1 , . . . , xN ) and vi = 0, ∀i ∈ N F , ẋi = 0 and vi = 0, ∀i / ∈ N F , n F (x 1 , . . . , xN ) := i∈N F H mi i (x i ) 2 , (A.4)
when col(x 1 , v1 , . . . , xN , vN ) ∈ F. Then, the Filippov solution of this network can be obtained by the corresponding differential inclusion [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF], which results in the sliding trajectory on S following the dynamics

ẋi = h i (x i ) nS (x 1 , v1 , . . . , xN , vN ) , vi = - (∂f i /∂x)(x i , vi ) (∂f i /∂v)(x i , vi ) h i (x i ) nS (x 1 , v1 , . . . , xN , vN ) =: R i (x i , vi ) h i (x i ) nS (x 1 , v1 , . . . , xN , vN ) , ∀i ∈ N , nS := N i=1 (1 + R i (x i , vi ) 2 ) h i (x i ) 2 .
(A.5)

Note that at jump points, i.e., points on P sp introduced in (7), if neuron i approaches the left knee, i.e, xi → x i (0), then we have R

i (x i , vi ) → -∞, hence R i (x i , vi ) h i (x i ) nS (x 1 , v1 , . . . , xN , vN ) = vi → 1
while ẋj , vj → 0 for j = i and ẋi → 0. This implies that we smoothly transit to the region UF, which follows (A.3). Similarly, the Filippov solution results in the sliding trajectory on F following the dynamics ẋi = H mi i (x i ) nF (x 1 , v1 , . . . , xN , vN ) , vi = -(∂f mi i /∂x)(x i , vi ) (∂f mi i /∂v)(x i , vi ) v1 , . . . , xN , vN ) =: R mi i (x i , vi ) v1 , . . . , xN , vN ) , i ∈ N F , ẋi = 0 and vi = 0, i / ∈ N F , nF := i∈N F

H mi i (x i ) nF (x 1 ,
H mi i (x i ) nF (x 1 ,
(1 + R mi i (x i , vi ) 2 )H mi i (x i ) 2 .

(A.6)

If neuron i ∈ N F approaches the right knee, i.e., xi → x i (m i ), then we have R mi i (x i , vi ) → -∞, hence R mi i (x i , vi ) v1 , . . . , xN , vN ) = vi → -1.

H mi i (x i ) nF (x 1 ,
This again implies that we smoothly transit to UF.

Notice that the obtained solution trajectory is a bit different from that of the network introduced in Sections 3.3, 4.1, and 4.2. In particular, (12) in Section 4.1 requires an infinite time of travel for a jump. However, what we only care about for this dynamics is the shape of the trajectory; the resulting point after a jump, and its effect on infinitesimal perturbation. In this respect, we could identify two networks, the one introduced in Sections 3.3, 4.1, and 4.2 and the one introduced in this subsection. This is because, the point before and after a jump is on the N -dimensional manifolds S and F, which is completely determined by x i 's. Note that x i 's and δx i 's do not change during this period of travel.

On the other hand, (6) (or ( 14)) is also different from (A.5) (or (A.6)), since (6) (or ( 14)) becomes infinitely faster than (A.5) (or (A.6)) as it approaches the knee (for some x i ). However, again we are only interested in its trajectory; the resulting point after slow (fast) travel, and its effect on infinitesimal perturbation.

In this respect, we could identify two networks, since we preserve the identity In particular, analysis performed in Sections 3.3-5.1 all applies to the network defined here. Now, therefore, from Theorems 1-3, we can conclude that the discontinuous dynamical system obtained in this subsection, also has a unique almost globally stable limit cycle, which is contained in a positively invariant set representing the behavior of persistent synchronous spiking, which has the convergence rate independent of the weak coupling strength and the network topology.

For sufficiently small , then we can find a slightly smaller (but extended in the direction of ultrafast convergence) positively invariant set containing the unique stable limit cycle, which has rapid convergence. Now, the semi-global convergence follows from the continuous dependence for any trajectory that has a finite length [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF]. The exception is made around the measure zero set ∪ N ⊂N X (N ) defined in Appendix B, which can be made arbitrarily small by reducing .

Therefore, we can prove the existence of the almost semiglobally rapidly stable limit cycle for the normalized network with a sufficiently small . Since, for fixed , we have finite velocity in the entire region, the shape and the rapid stability of the limit cycle are preserved.

B Fast-slow behavior in the rest of the region D \ (S ∪ U F ∪ F)

To complete the illustration on the network that we obtain by the time scale separation for the = 0 limit, we will detail the fast-slow behavior in the rest of the region D \ (S ∪ UF ∪ F). Recall that any trajectory that starts from the slow region S moves to the fast region F passing through the ultrafast region UF, and any trajectory that starts from the fast region F moves to the slow region S passing through the ultrafast region UF after possibly a finite number of iteration F → UF → F. Now, assume without loss of generality that

D i = [x i (0), x i (M i )]×[v sm i (0, x i (M i )), v lg i (M i , x i (0))], ∀i.
Then, the remaining region can be partitioned as

N ⊂N M(N ),
where M(N ) is the set of all points col(x 1 , v 1 , . . . , x N , v N ) such that

0 = f i (x i , v i ) + g i (E i -v i ) j∈Ni α ij S ij (v j ) =: f i (x i , v i ) + g i (E i -v i )d i (v 1 , . . . , v N ), i ∈ N ,
and v i = v md i (g i d i (v 1 , . . . , v N ), x i ) if and only if i / ∈ N , i.e., v i = v sm i (g i d i (v 1 , . . . , v N ), x i ) or v i = v lg i (g i d i (v 1 , . . . , v N ), x i ) for i ∈ N , where v md i (m, x) denotes the middle root of f m i (x, v) = 0 for x ∈ (x i (m), x i (m)) so that v sm i (m, x) < v md i (m, x) < v lg i (m, x).

A similar argument as earlier shows that neurons i in the lower branch or the upper branch, i.e., i ∈ N , does not move towards the middle branch, hence any trajectory that starts from M(N ) can only move to either S ∪ F or M(N ) with N ⊇ N , passing through the ultrafast region UF. In particular, note that any trajectory that starts from the set M(N ) ends up in either S or M S (N ) := {col(x 1 , v 1 , . . . , x N , v N ) ∈ M(N ) :

v i = v sm i (g i d i (v 1 , . . . , v N ), x i ), i ∈ N }
with some N ⊇ N , via the ultrafast-fast behavior.

Therefore, in the remaining of this section, we will show that M S (N ) is partitioned into a measure zero set X (N ) and the rest M S (N ) \ X (N ) where a trajectory that starts from it moves to either S or M S (N ) with some N N , via the ultrafast-fast behavior. In particular, let X (N ) := {col(x 1 , v 1 , . . . , x N , v N ) ∈ M S (N ) :

∃i / ∈ N s.t. x i = x i,∞ (v i )}, S M (N ) := {col(x 1 , v 1 , . . . , x N , v N ) ∈ M S (N ) :

v i ≤ E th i , ∀i / ∈ N },
M S,lef t (N ) := {col(x 1 , v 1 , . . . , x N , v N ) ∈ M S (N ) :

v i ≤ E th i or x i < x i,∞ (v i ), ∀i / ∈ N },
M S,right (N ) := M S (N ) \ (X (N ) ∪ M S,lef t (N )).

Then, X (N ) is a measure zero set and any trajectory that starts from M S,lef t (N ) \ S M (N ) moves either to S M (N ) via the fast behavior or to M S (N ) via the ultrafast-fast behavior with some N N . On the other hand, any trajectory that starts from M S,right (N ) moves to M S (N ) via the ultrafast-fast behavior with some N N . Finally, any trajectory that starts from S M (N ) moves either to M S,right (N ) via the slow behavior or to M S (N ) via the slowultrafast behavior with some N N .

Fig. 1 .

 1 Fig. 1. Rapid synchronization of two identical (a) and non-identical (b) excitable systems under weak excitatory synaptic coupling. (c) Poor synchronization of the same non-identical excitable systems under strong diffusive coupling.

Fig. 2 .

 2 Fig. 2. An illustration of Assumptions E1-E3 and S2 (see text for notations). Assumption E1 (Spiking phase portrait) The nullcline x = x i,∞ (v) has a monotonic shape and the nullcline f i (x, v) = 0 is N -shaped, as illustrated in Figure 2. The nullcline x = x i,∞ (v) is on the left (right) to

M

  i := 0.36 < 1.1003 = 1.4260 * 1.75 -0.3 * 0.75 -* g i,Ca -0.3 * g i,L -0.4 1.7 .

Fig. 3 .

 3 Fig.3. A simulation result of 100 neurons modeled in the manner explained in Section 7 is shown as a raster plot of spike times (i.e., times of the extreme peak for each cycle) of individual neurons, where the initial condition is also taken randomly with uniform distribution from Di. The top, middle, and bottom figures correspond to τ = 0.5, 5, and 30, respectively. For illustration, decoupled heterogeneous neurons are shown until t = 100 and the synaptic coupling becomes active at the red line.

  the argument in Sections 4.2 and 5.1 and the resulting point.

  The mechanism analyzed in the present paper is of general nature in excitable models. It would be of interest to extend the analysis to other types of synchronization, for instance, synchronization mechanisms in heterogeneous ensembles of bursting neurons.

	Systems Magazine, 34 (2), 50-74.
	Somers, D. & Kopell, N. (1993). Rapid synchronization
	through fast threshold modulation. Biological Cyber-
	netics, 68, 393-407.
	Somers, D. & Kopell, N. (1995). Waves and syn-
	chrony in networks of oscillators of relaxation and non-
	relaxation type. Physica D: Nonlinear Phenomena,
	89 (1-2), 169-183.
	Stan, G.-B. & Sepulchre, R. (2007). Analysis of inter-
	connected oscillators by dissipativity theory. IEEE
	Transactions on Automatic Control, 52 (2), 256-270.
	Strogatz, S. H. (2000). From Kuramoto to Crawford:
	exploring the onset of synchronization in populations
	of coupled oscillators. Physica D: Nonlinear Phenom-
	ena, 143 (1-4), 1-20.
	Wang, W. & Slotine, J.-J. E. (2005). On partial contrac-
	tion analysis for coupled nonlinear oscillators. Biolog-
	ical Cybernetics, 92, 38-53.
	Wechselberger, M. (

  , v1 , . . . , xN , vN ) ∈ S, ẋi = 0 and vi = fi (x i , vi , {v j } Ni ) nUF (x 1 , v1 , . . . , xN , vN ) , ∀i ∈ N ,

			when col(x 1 nUF :=	N	fi (x i , vi , {v j } Ni ) 2 ,
				i=1
				(A.3)
			when col(x 1 , v1 , . . . , xN , vN ) ∈ UF, and
			ẋi =	H
	xN )	and vi = 0, ∀i ∈ N ,
		N	(A.2)
	n S (x 1 , . . . , xN ) :=	h i (x i ) 2 ,
		i=1

Preprint submitted to Automatica

Note that DPsp /dPsyn can be made arbitrarily small if the monotonicity is further extended so that there is a sufficient difference in the speed in the region[(1-ri)x i +rix F i , x F i -η * ) and (x F i -η * , xi].

Note that the left-hand side of the inequality can be made arbitrarily close to 1 by decreasing η.