Clusters in nonsmooth oscillator networks
Rachel Nicks Lucie Chambon Stephen Coombes









rachel.nicks@nottingham.ac.uk
Centre for Mathematical Medicine and Biology, School of Mathematical Sciences,

University of Nottingham,NG7 2RD, UK.




lucie.chambon@inria.fr

Centre de recherche INRIA Sophia-Antipolis Méditerranée, Borel building

2004, route des Lucioles – BP 93

06 902 Sophia Antipolis Cedex., France.




stephen.coombes@nottingham.ac.uk
Centre for Mathematical Medicine and Biology, School of Mathematical Sciences,

University of Nottingham,NG7 2RD, UK.





For coupled oscillator networks with Laplacian coupling the master stability function (MSF) has proven a particularly powerful tool for assessing the stability of the synchronous state.  Using tools from group theory this approach has recently been extended to treat more general cluster states.  However, the MSF and its generalisations require the determination of a set of Floquet multipliers from variational equations obtained by linearisation around a periodic orbit.  Since closed form solutions for periodic orbits are invariably hard to come by the framework is often explored using numerical techniques.  Here we show that further insight into network dynamics can be obtained by focusing on piecewise linear (PWL) oscillator models.  Not only do these allow for the explicit construction of periodic orbits, their variational analysis can also be explicitly performed.  The price for adopting such nonsmooth systems is that many of the notions from smooth dynamical systems, and in particular linear stability, need to be modified to take into account possible jumps in the components of Jacobians.  This is naturally accommodated with the use of saltation matrices.  By augmenting the variational approach for studying smooth dynamical systems with such matrices we show that, for a wide variety of networks that have been used as models of biological systems, cluster states can be explicitly investigated.  By way of illustration we analyse an integrate-and-fire network model with event-driven synaptic coupling as well as a diffusively coupled network built from planar PWL nodes, including a reduction of the popular Morris–Lecar neuron model.  We use these examples to emphasise that the stability of network cluster states can depend as much on the choice of single node dynamics as it does on the form of network structural connectivity.  Importantly the procedure that we present here, for understanding cluster synchronisation in networks, is valid for a wide variety of systems in biology, physics, and engineering that can be described by PWL oscillators.


05.45.Xt,05.45.-a,02.20.-a






Clusters in nonsmooth oscillator networks
    Stephen Coombes
    December 30, 2023
=========================================






§ INTRODUCTION


The study of synchrony in coupled oscillator networks has a very long history, dating all the way back to the work of Huygens on two interacting pendulum clocks.  Since his observations about the emergence of “an odd kind of sympathy" there have been countless other examples of synchrony discussed in the natural sciences and engineering ranging from the dynamics of populations of flashing fireflies to those of coupled Josephson junctions.  For an excellent review see Arenas et al. <cit.> and the recent focus issue on patterns of network synchronisation <cit.>.  However, perfect global synchronisation is just one of many states expected to emerge in structured oscillator networks.  Indeed instabilities of the synchronous state are generically expected to give rise to cluster states, in which sub-populations may synchronise, though not with each other.  This class of solutions has been relatively well explored for phase-oscillator networks, as in the work by Brown et al. <cit.> and Ashwin et al. <cit.>, though less so for networks of limit-cycle oscillators.  For this more general scenario Golubitsky and Stewart have made use of the Hopf bifurcation with symmetry to understand cluster states in ring networks, and provided a number of case studies of systems with nearest-neighbour coupling <cit.>.  Pogromsky and collaborators have also exploited symmetry under permutation of a given network of dynamical systems coupled through diffusion to determine the stability of cluster states, with the aid of an appropriate Lyapunov function <cit.>.  For the case of diffusively coupled oscillators, Belykh et al. have further provided a complete classification of cluster states, together with conditions that can be used to determine the coexistence of stable cluster states with an unstable synchronous state <cit.>.  More recently Pecora et al. <cit.>  and Sorrentino et al. <cit.>, have extended the master stability function (MSF) approach of Pecora and Carroll <cit.> making extensive use of tools from computational group theory.  This work paves the way for a systematic study of cluster states and their bifurcations for a very broad class of networks with state-dependent interactions, including arrays of electrochemical oscillators <cit.>, gene networks <cit.>, and the brain <cit.>.  Moreover, it will likely play a key role in addressing the growing interest in dynamics on networks <cit.>.

At heart the MSF approach exploits the fact that the Floquet theory for the stability of a synchronous network state decouples into a set of lower dimensional Floquet problems.  By studying just one of these lower dimensional Floquet problems, as a function of a single complex parameter, the linear stability of the synchronous network state can be determined.
Indeed the approach is sufficiently general that, for Laplacian coupling of identical nodes (so that the network synchronous state is guaranteed to exist as long as a single node oscillates), then the effect of different network choices is easily explored.  Namely for a given choice of node dynamics the MSF maps a complex number to the largest Floquet exponent of the reduced problem, and by choosing the complex number to run over all possible eigenvalues of the unnormalised graph Laplacian of the network connectivity then stability is guaranteed if the MSF is always negative.  The challenge of exploring different networks thus reduces to determining the eigenstructure of the graph Laplacian, with the caveat that the Floquet problem can be solved.  However, for a general nonlinear system (without special structure) it is unlikely that this problem can be solved analytically.  Thus it would be highly complementary to the MSF formalism, and its extensions to treat cluster states, if this barrier could be reduced or removed entirely.  It is this issue that we address in this paper.  We do so by restricting the choice of node dynamics to that of piecewise linear (PWL) systems.  Although at first sight this may appear overly restrictive, there has been an appreciation for some time in the mathematical community of the benefits of studying caricatures of complex systems built from PWL and possibly discontinuous dynamical systems <cit.>.
Indeed, there is a now growing perspective in the applied dynamical systems community that piecewise models are highly suitable for modern applications in science and can complement the smooth dynamical systems approach that has dominated to date <cit.>.  It has recently been shown how to extend the MSF formalism for synchrony to PWL models with state-dependent coupling <cit.>.  Here we go further and show how to treat discontinuous systems of integrate-and-fire (IF) type, and how to analyse the existence and stability of more general cluster states in a broad class of PWL network models with both state and time-dependent interactions.

In <ref> we give a brief recapitulation of the MSF for smooth dynamical systems as well as its extension to networks of nodes built from PWL oscillators.
In illustration of the utility of this combination of MSF and nonsmooth systems we show how to analyse synchrony in PWL IF networks with balanced synaptic coupling.  Next in <ref> we survey the techniques from computational group theory that allow the discovery of cluster states from the topology of the network.  For PWL systems we then  show, in  <ref>, how such orbits can be explicitly determined.  Similarly the variational equations that determine the stability of a cluster state are equally tractable.  By exploiting the PWL nature of the network dynamics we are able to determine the stability of cluster states without recourse to numerically evolving the variational equations.  Our approach also facilitates an explicit bifurcation analysis, which we illustrate for a simple five-node network with a PWL analogue of a Hopf normal form.
In <ref> we further explore a node dynamics that exhibits a PWL analogue of a homoclinic bifurcation, as well as a PWL IF model.   We use these examples to highlight that the stability of network cluster states can depend as much on the choice of single node dynamics as it does on the form of network structural connectivity.  Finally in <ref> we discuss natural extensions of the work in this paper.






§ THE MASTER STABILITY FUNCTION:  A RECAPITULATION


The MSF allows the determination of the linear stability of the synchronous state in a quite large class of smooth networks of identical nodes.
To describe the MSF formalism it is convenient to consider N nodes (oscillators) and let x⃗_i ∈^m be the m-dimensional vector of dynamical variables of the ith node with isolated (uncoupled) dynamics ẋ⃗̇_i = F⃗(x⃗_i), with i=1,…, N.  The output for each oscillator is described by a vector function H⃗∈^m.
For a given coupling matrix with components W_ij and a global coupling strength σ the network dynamics, to which the MSF formalism applies, is specified by

    tx⃗_i    = F⃗(x⃗_i) + σ∑_j=1^N W_ij [ H⃗(x⃗_j)-H⃗(x⃗_i)  ] ,
       ≡F⃗(x⃗_i) - σ∑_j=1^N 𝒢_ijH⃗(x⃗_j) .

Here the matrix 𝒢 with components 𝒢_ij has the unnormalised graph Laplacian structure 𝒢_ij = -W_ij + δ_ij∑_kW_ik.
The N-1 constraints x⃗_1(t)=x⃗_2(t)=… =x⃗_N(t) = s⃗(t) define the (invariant) synchronisation manifold, with s⃗(t) a solution of the uncoupled system, namely ṡ⃗̇=F⃗(s⃗).  To assess the stability of this state a linear stability analysis is performed by expanding a solution as x⃗_i(t) = s⃗(t) +δx⃗_i(t) to obtain the variational equation

    tδx⃗_i = D F⃗(s⃗) δx⃗_i - σ D H⃗(s⃗) ∑_j=1^N 𝒢_ijδx⃗_j .

Here D F⃗(s⃗) and D H⃗(s⃗) denote the Jacobian of F⃗(s⃗) and H⃗(s⃗) around the synchronous solution respectively.
The variational equation has a block form that can be simplified by projecting δx⃗ into the eigenspace spanned by the (right) eigenvectors of the matrix 𝒢.  If we organise these column eigenvectors into a matrix P then 𝒢P = P Λ, with Λ=diag(λ_1,…,λ_N), where λ_l are the corresponding eigenvalues for l=1,…,N.  If we collect the perturbations in a vector U⃗=(δx⃗_1, …, δx⃗_N) ∈^Nm, and introduce a new variable V⃗ according to the linear transformation V⃗ = (P ⊗ I_m)^-1U⃗, then we have that

    tV⃗ = (I_N ⊗ DF⃗(s⃗)) V⃗ - σ(Λ⊗ DH⃗(s⃗)) V⃗ .

Here the symbol ⊗ denotes the tensor (or Kronecker) product for matrices, and I_N is the N × N identity matrix.  Thus we have a set of N decoupled eigenmodes in the block form

    tξ⃗_l =  [ D F⃗(s⃗) - σλ_l D H⃗(s⃗)  ] ξ⃗_l,      l=1,…,N,

where ξ⃗_l is the lth (right) eigenmode associated with the eigenvalue λ_l of 𝒢 (and D F⃗(s⃗) and D H⃗(s⃗) are independent of the block label l).  Since ∑_j 𝒢_ij=0 there is always a zero eigenvalue, say λ_1=0, with corresponding eigenvector (1,1,…,1), describing a perturbation parallel to the synchronisation manifold.  The other N-1 transverse eigenmodes must damp out for synchrony to be stable.  For a general matrix 𝒢 the eigenvalues λ_l may be complex, which brings us to consideration of the system

    tξ⃗ =  [ D F⃗(s⃗) - β D H⃗(s⃗)  ] ξ⃗,     β = σλ_l ∈  ,

where ξ∈^m.
The MSF is defined as the function which maps the complex number β  to the largest Floquet exponent of (<ref>).
The synchronous state of the system of coupled oscillators is stable if the MSF is negative at  β=σλ_l  where  λ_l  ranges over the eigenvalues of the matrix 𝒢 (excluding λ_1=0).

We note that the Laplacian form of coupling in (<ref>) guarantees the existence of the synchronous state as long as a single uncoupled node has a periodic solution.  However, other forms of coupling, and in particular those described by adjacency matrices or their weighted counterparts are also common.  In this case (<ref>) would be replaced by the system

    tx⃗_i = F⃗(x⃗_i) + σ∑_j=1^N W_ijH⃗(x⃗_j) .

The synchronous solution (with an orbit inherited from the uncoupled node) is not a generic solution of (<ref>) unless there is a further constraint placed on the coupling matrix, namely that ∑_j=1^N W_ij = constant for all i.  If this constant value is 0, then we shall say that the system is balanced.  This scenario is ubiquitous in the modelling of many spiking neural networks <cit.>, and is relevant to understanding coding <cit.>, memory <cit.>, noise-robust neuronal selectivity <cit.>, and brain idling <cit.>.
Interestingly many spiking neural networks are modelled using integrate-and-fire (IF) neurons which are discontinuous dynamical systems owing to the hard reset of the voltage state variable upon reaching a spiking threshold <cit.>.  A naive application of the MSF formalism described above would lead to incorrect results since it assumes that the underlying dynamics is smooth.  Fortunately it is relatively straightforward to augment the MSF approach to handle firing discontinuities in IF networks using ideas from the study of impact oscillators.  This was originally done by Coombes to give meaning to Liapunov exponents for linear IF neurons <cit.>, and more recently by Ladenbauer et al. <cit.> for constructing the MSF in synaptically coupled networks of nonlinear (adaptive exponential) IF neurons. However, this latter work requires substantial numerical simulation as the periodic orbit for the synchronous solution is not available in closed form.  The perspective in this paper is to choose PWL caricatures of node dynamics to overcome this last barrier.  A case in point is the caricature of the adaptive exponential IF model developed in <cit.>,
which has a PWL subthreshold dynamics and an adaptive jump upon reset.  Since this is an exemplar PWL system we will describe it here in more detail to set the scene for the extension of the MSF formalism to cover nonsmooth systems in general and not just IF type models.  For a further perspective on the use of techniques from nonsmooth systems in neuroscience see <cit.>.




 §.§ MSF for PWL IF networks with synaptic coupling


There are now many variants on nonlinear integrate-and-fire neurons that are able to fit the spike trains of real neurons, such as those due to Gröbler et al. <cit.>, Izhikevich <cit.>, and Badel et al. <cit.>.  The planar adaptive exponential model has been particularly successful at fitting data from cortical fast spiking interneurons and regular spiking pyramidal neurons <cit.>.  Importantly it admits to a useful PWL reduction whereby its voltage nullcline is replaced by a PWL function <cit.>.  A similar approach is taken in <cit.> to develop an analytically tractable PWL IF model that can support both tonic (repetitive firing with a constant interspike interval) and bursting states (with regular alternations between short and long interspike intervals).  Here we subsume both models within a general description of the form ż⃗̇=F⃗(z⃗), z⃗∈^m with

    F⃗(z⃗) = 
    A_L z⃗ + c⃗_L,    z⃗∈ S_L 
    
    A_R z⃗ + c⃗_R,    z⃗∈ S_R
    ,

with A_L,R∈^m × m, c⃗_L,R∈^m, and ^m = S_L ∪Σ∪ S_R.
For simplicity we restrict attention to the situation that the phase-space of the model can be broken into two regions S_L,R separated by a boundary given implicitly by the zero of an indicator function h, such that Σ={z⃗ | h(z⃗)=0}.  It is a simple matter to break the phase-space up into more regions, to incorporate further switching manifolds, so we shall stick to describing the simplest situation (though will give examples in <ref> of systems with two switching manifolds).  Should the dynamics on the switching manifold not be continuous with that in regions S_L and S_R then it is common to adopt the Filippov convex method for dynamics on Σ <cit.>.

For a planar PWL IF model that caricatures the adaptive exponential IF model we set m=2, z⃗=(v,w), v≤ v_th, and

    A_L=[   a_L    -1; a_w/τ b_w/τ ],     
    A_R=[   a_R    -1; a_w/τ b_w/τ ],

with c⃗_L=(I,0)^𝖳=c⃗_R for some constant drive I.  The switching manifold is prescribed by the choice v=0.
Whenever the voltage variable v reaches the firing threshold v_th then the system is reset according to
z⃗→g⃗(z⃗) = (v_r, w +κ/τ).  This gives rise to another type of switching manifold that we shall refer to as the firing manifold. If we introduce the function h(z⃗;a)=v-a, then the indicator function for switching is given by h(z⃗;0), whilst that for firing by h(z⃗;v_th).
From now on we shall simply refer to the model described above as the PWL-IF model.  We note that the vector field defining the PWL-IF model is continuous at v=0 (so that there is no need to invoke the Filippov convex method).
A set of phase planes of the model (with nullclines and typical trajectories) for different parameter choices is shown in Fig. <ref>.  The PWL-IF model is able to support a tonic firing pattern, which itself is unstable to a period doubling bifurcation as seen in the top part of Fig. <ref>.  It can also support a burst firing pattern, which is unstable to spike-adding bifurcations as shown in the bottom part of Fig. <ref>.



Each of the periodic orbits shown in Fig. <ref> can be obtained without recourse to the numerical evolution of the underlying nonsmooth flow for the PWL-IF model.  Rather the shape of the orbits can be determined explicitly by exploiting the linearity of the model between switching and firing manifolds.  In any such region the model is given by a linear system of ordinary differential equations (ODEs) of the form ż=A z⃗+ c⃗ for appropriate choices of A and c⃗.  The explicit solution can be obtained in terms of initial data using matrix exponentials as

    z⃗(t) = ^Atz⃗(0) + A^-1[^At-I_m] c⃗ .



For example, to determine the tonic firing pattern shown in Fig. <ref>, which only ever visits the region of phase-space with z⃗∈ S_R we would simply use the formula for (<ref>) with A=A_R and c⃗=c⃗_R to determine the time-of-flight Δ according to v(Δ)=v_th with (v(0),w(0))=(v_r,w_0), with w_0 determined self-consistently according to w_0=w(Δ)+ κ/τ.
The resulting pair of nonlinear equations for (Δ,w_0) can in general be solved with a numerical root finding scheme (and choosing the solution with  the largest value of w_0).  Thus although we can eliminate the need to numerically solve an ODE, there is still some need for root-finding.  To determine the burst patterns shown in Fig. <ref>, with M spikes per burst, would require the simultaneous solution of M+2 equations (one for the value of w on the switching manifold, and the others to determine the times-of-flight between consecutive firing events, remembering that for this orbit the trajectory also visits the region S_L).  However, once the orbit is determined in this way the Floquet theory for stability simplifies considerably.  To see how this simplification arises it is enough to focus on the stability of simple tonic firing pattern.  In this case the Jacobian of the orbit is the constant matrix A_R, and is independent of the orbit shape.  Thus small perturbations to the periodic orbit δz⃗(t) can be simply constructed according to δz⃗(t)  = exp(A_R t) δz⃗(0).  The caveat being that this result only holds away from switching or firing events.  To propagate perturbations properly through switching and firing manifolds we make use of a saltation matrix <cit.>.  A derivation of the saltation matrix in the context of this paper is given in Appendix A.  Thus for our nonsmooth system the evolution of perturbations to the orbit over one period is given by δz⃗(Δ)  = Ψδz⃗(0), where Ψ=K(Δ) exp(A_R Δ) and K is the saltation matrix:

    K(t) = [         v̇(t^+)/v̇(t^-)                       0; (ẇ(t^+)-ẇ(t^-))/v̇(t^-)                       1 ] .

The orbit will be stable provided that the eigenvalues of Ψ lie within the unit disc.  Since one of these will be unity (reflecting time-translation invariance) we have that the orbit is stable provided Re r <0, where the Floquet exponent r can be calculated explicitly as

    r    = Tr A_R + 1/Δlnv̇(Δ^+)/v̇(Δ^-)
       =a_R+b_w/τ + 1/Δlna_R v_r-w(0)+I/a_R v_th-w(Δ)+I .

Here we have used the fact that Ψ = 1 ×λ, and defined the non-trivial eigenvalue λ of Ψ as λ = exp(r Δ).

In Fig. <ref> we show a plot of Re  r as a function of the parameter a_w.  This predicts the point of a period doubling instability, and is found to be in excellent agreement with direct numerical simulations.  The Floquet exponent for more complicated orbits (that burst and visit both S_L and S_R) can also be easily constructed using a revised structure for Ψ.  For example for a burst pattern with M spikes of the type shown in Fig. <ref> then we would have that

    Ψ   =  K(T_M+2)^A_L Δ_M+2K(T_M+1)^A_R Δ_M+1
       × K(T_M)^A_R Δ_M… K(T_2)^A_R Δ_2K(T_1)^A_R Δ_1 ,

where Δ_i indicate various times-of-flight, and the T_i the times at which firing or switching events occur.




We now consider a network of N synaptically coupled PWL-IF neurons with the introduction of an index i=1,…,N (labelling each node) and the replacement of the constant external drive by a time-dependent forcing such that
I → I + σ∑_j W_ij s_j(t).  Here the synaptic input from neuron j takes the standard event-driven form

    s_j(t) =  ∑_p ∈η(t-T_j^p),

where T_j^p denotes the pth firing time of neuron j and η describes the shape of the post-synaptic response.  Other forms of state-dependent synaptic coupling are also possible, such as the fast threshold modulation type used by Belykh and Hasler to study clusters in networks of bursting neurons <cit.>, though would require a further PWL reduction before they could be studied with the method presented here.
For a balanced network the existence of the synchronous network state is independent of any of the parameters describing the synaptic interaction.  This means that the form of synaptic coupling cannot induce any nonsmooth bifurcations, such as grazing, which occur if an orbit tangentially touches the firing threshold.
Here we shall focus on the common choice of a continuous α-function, so that η(t) = α^2 t ^-α tΘ(t), where Θ is a Heaviside step function.  In this case we may also write s_i(t) as the solution to an impulsively forced linear system:

    ( 1 + 1/αt ) s_i    = u_i, 
     ( 1 + 1/αt ) u_i    = ∑_p ∈δ(t-T_i^p).

Exploiting the linearity of the synaptic dynamics between firing events we may succinctly write the network model in the form ż⃗̇_i = F⃗(z⃗_i), where z⃗_i = (v_i,w_i,s_i,u_i), and F⃗ has the form of (<ref>) with

    A_L,R=[ a_L,R    -1     0     0; a_w/τ b_w/τ     0     0;     0     0    -α     α;     0     0     0    -α ],

and c⃗_L=(I,0,0,0)^𝖳=c⃗_R, with z⃗_i →g⃗(z⃗_i) = (v_r, w_i +κ/τ, s_i, u_i+ α) whenever h(z⃗_i;v_th)=0.  The vector function that specifies the interaction is given by H⃗(z⃗_i) = (s_i,0,0,0)^𝖳.
The MSF approach for a smooth network gives rise to a Floquet problem with a Jacobian D F⃗(s⃗) + σλ_l D H⃗(s⃗), where λ_l is an eigenvalue of the coupling matrix W.  The corresponding Jacobian for the PWL-IF network is A_L,R + σλ_l D H⃗, with the label L or R chosen according to whether the synchronous orbit is in S_L or S_R.  Moreover, D H⃗ is now a constant matrix with [D H⃗]_ij = 1 if i=1 and j=3, and is 0 otherwise.

The propagation of (linearised) trajectories through a switch or a firing event is achieved with the use of a saltation matrix K∈^m × m, and we write this as U⃗^+ = (I_N ⊗ K) U⃗^-, where U⃗^- and U⃗^+ are respectively the system states just before and just after the saltation.  In transformed coordinates (and see equation (<ref>)) it is simple to establish that V⃗^+ = (I_N ⊗ K) V⃗^-, so that saltation also acts blockwise with ξ_l^+ = K ξ_l^-.
Thus the synchronous solution is stable (following the same arguments as for a single node) if all of the eigenvalues of a set of matrices Ψ_l, l=1,…,N lie within the unit disc, excluding the one that arises from time-translation invariance (with a value +1).
For example for a synchronous tonic orbit of the type discussed above,

    Ψ_l = K(Δ) exp{ (A_R + σλ_l D H⃗) Δ} .

Here the saltation matrix is given by (and see Appendix A)

    K(t) = [           v̇(t^+)/v̇(t^-)                         0                         0                         0;   (ẇ(t^+)-ẇ(t^-))/v̇(t^-)                         1                         0                         0;   (ṡ(t^+)-ṡ(t^-))/v̇(t^-)                         0                         1                         0; (u̇(t^+)-u̇(t^-))/v̇(t^-)                         0                         0                         1;                           ] .

The periodic trajectory for z⃗(t) is subject to the constraints v(Δ)=v_th, w(0)=w(Δ)+κ/τ, s(0)=s(Δ), and u(0)=u(Δ)+α.
If we denote an eigenvalue of K(Δ) exp{ (A_R + β D H⃗) Δ}, β∈, by γ(β) then the MSF is the largest number in the set Re  ( lnγ(β) )/Δ, and the synchronous state is stable if the MSF is negative at all the points where β=σλ_l.
In Fig. <ref> we show a plot of the MSF for the synchronous tonic orbit.
This solution will be stable provided all of the eigenvalues of σ W lie within the shaded area  shown in Fig. <ref> (for a given α). For a positive semi-definite connectivity matrix then we see from Fig. <ref> that the synchronous solution is unstable for σ >0.  However, for σ<0 the same network can support a stable synchronous orbit for some sufficiently small value of |σ|.





We note here that for the choice of an exponential synapse as originally considered in <cit.>, with η(t) = α^-α tΘ(t), then synaptic interactions are discontinuous and the order in which perturbed components of the state vector cross the firing threshold becomes important <cit.>.  In this case the approach above, valid only for continuous interactions, must be modified (although this was not considered in <cit.>).  We shall treat this mathematically interesting case in another paper.







As a particular realisation of a network architecture that guarantees synchrony we choose a balanced ring network with N odd and W_ij = W(|i-j|), with distances calculated modulo (N-1)/2, and w(x) = (1-a|x|/d) ^-|x|/d.  Here the parameter a is chosen such that, for a given value of N and a scale d, ∑_j=1^N W_ij=0.  The circulant structure of this matrix means that it has normalised eigenvectors given by e⃗_l = (1,ω_l, ω_l^2, …, w_l^N-1)/√(N), where l=0,…, N-1, and ω_l=exp(2 π i l /N) are the Nth roots of unity.  The eigenvalues of the (symmetric) connectivity matrix are real and given by

    λ_l = ∑_j=0^N-1 W(|j|) ω_l^j.

We note that the balance condition enforces λ_0=0.  We also have that λ_N-l=λ_l for l=1,…, (N-1)/2, so that any excited pattern is given by a combination e⃗_m+e⃗_-m = 2 Re e⃗_m for some 1≤ m ≤ (N-1)/2.  Given the shape of the MSF function shown in Fig. <ref> the value of m is determined by λ_m= max_lλ_l.  In Fig. <ref> we compare simulations of a network against the predictions of the MSF.  When the network eigenvalues lie within the region where the MSF is negative, then small perturbations to synchronous initial data die away and the system settles to a synchronous periodic orbit as expected.  When one of the eigenvalues crosses the zero level set of the MSF (from negative to positive) we see two different types of instability emerge.  One is the emergence of a spatio-temporal pattern of spike-doublets, arising because an eigenvalue of Ψ_l leaves the unit disc at -1 (a period-doubling bifurcation), and another giving rise to periodic travelling wave (with asynchronous firing) because an eigenvalue of Ψ_l leaves the unit disc at +1 (tangent bifurcation).




§ CLUSTER STATES AND NETWORK STRUCTURE


We now move on to discuss the more general phenomenon of cluster synchronisation, where different groups of oscillators are exactly synchronised but there is no exact synchrony between the groups. We focus on networks with symmetry, where cluster states arise very naturally (although clustering can also occur in networks without symmetry where synchronous nodes have synchronous input patterns <cit.>). For networks of identical oscillators with dynamics described by (<ref>), a symmetry of the network is a permutation γ of the nodes which leaves the network equations unchanged. That is, if we denote x⃗= (x⃗_1, …, x⃗_N) ∈ℝ^Nm, M_γ the N× N permutation matrix for the permutation γ and 
    G⃗_i(x⃗) = F⃗(x⃗_i) - σ∑_j=1^N 𝒢_ijH⃗(x⃗_j),
 then the network dynamics (<ref>) are given by ẋ⃗̇ = G⃗(x⃗) for G⃗= (G⃗_1, …, G⃗_N)^𝖳 where the vector field G⃗ satisfies

    G⃗((M_γ⊗ I_m) x⃗) =(M_γ⊗ I_m) G⃗(x⃗).

 This results in the condition that γ is a symmetry of the network if M_γ𝒢 = 𝒢 M_γ (or equally M_γ commutes with the connectivity matrix). The network symmetries form a group Γ⊆ S_N. Many real-world networks have been shown to have a high degree of symmetry, arising from locally tree-like structures produced by natural growth of the network <cit.>. Determining symmetries of such large and complex networks is impossible by inspection. Indeed even networks with a small number of nodes can have a large symmetry group <cit.>. However, the computations required to determine the symmetry group of a given network are easily implemented using computational algebra routines <cit.>.



Recent work of Pecora et al. <cit.> and Sorrentino et al. <cit.> has demonstrated how the ideas behind the MSF can be combined with group theoretical techniques used in the study of symmetric dynamical systems to analyse the stability of cluster states within symmetric networks of dynamical units. The approach taken in <cit.> and <cit.> continues to focus on networks of coupled identical oscillators whose dynamics are given by (<ref>). This approach is equivalent to a restriction (to a particular form of admissible network equations (<ref>) in the special case of identical oscillators and one type of bidirectional coupling) of the more general theory of patterns of synchrony in coupled cell systems developed by Golubitsky and Stewart and collaborators and recently reviewed in <cit.>. Using terminology from the more general theory, here we review for the particular network dynamics (<ref>) how network structure can be used to determine a catalogue of possible cluster states. We also discuss techniques for determining the stability of any given cluster state, utilising methods from computational group theory. The exposition below is inspired by that of <cit.> and is presented as a recipe for how to extend the standard MSF approach to treat networks of PWL oscillators.

Many of the possible cluster states which a given network may support arise from the symmetry of the network. We first describe how these cluster states may be determined before highlighting the algorithm of Sorrentino et al. <cit.> which can be used to determine additional cluster states resulting from the particular choice of Laplacian coupling.





 §.§ Cluster states from network symmetries


Consider a network of N identical oscillators with one type of bidirectional coupling whose dynamics are given by (<ref>) and which has symmetry group Γ.  As a consequence of the dynamics satisfying the equivariance condition (<ref>), if x⃗(t) is a solution of (<ref>) (equilibrium or periodic state) then for any permutation γ∈Γ,

    tγx⃗ = γx⃗t = γG⃗(x⃗) = G⃗(γx⃗),
 so γx⃗(t) is also a solution. Here γ acts on x⃗ via the permutation matrix (M_γ× I_m). Thus solutions occur as group orbits, {γx⃗ :  γ∈Γ}. The isotropy subgroup of a solution is the subgroup Σ⊆Γ given by

    Σ = {γ∈Γ :  γx⃗(t) = x⃗(t),  ∀ t ∈ℝ}.
 This is the group of (spatial) symmetries of the solution x⃗(t) and the largest subgroup of Γ under which the solution is invariant. Solutions which lie on the same group orbit have conjugate isotropy subgroups and the same existence and stability properties. Given a subgroup Σ⊆Γ we can define its fixed-point subspace

    (Σ) = {x⃗ :  γx⃗ = x⃗, ∀γ∈Σ}.

Fixed-point subspaces are flow invariant: Let γ∈Σ and x⃗∈(Σ). Then

    γG⃗(x⃗) = G⃗(γx⃗) = G⃗(x⃗),
 so G⃗(x⃗) ∈(Σ).

Suppose that 
 Σ is any subgroup of Γ. The orbit under Σ of the node i is the set {γ(i)  :  γ∈Σ}. The orbits permute subsets of nodes among each other and in this way partition the nodes into clusters. Nodes which lie on the same orbit (in the same cluster) have synchronised dynamics, x⃗_γ(i)≡x⃗_i for any γ∈Σ - see <cit.>. Also, the synchronised state for each cluster is flow-invariant. Thus, to enumerate all possible synchronised cluster states of a given network which are due to network symmetries, we need to determine the network symmetry group Γ and all of its subgroups. There will be one type of cluster state for each isotropy subgroup for the action of Γ on the node space up to conjugacy.

Example: A five-node network

Consider the example five-node network studied by Sorrentino et al. <cit.> which has graph Laplacian matrix
 
    𝒢 = [r] 3    -1    0     -1     -1
     -1     3    -1     0     -1  
     0    -1    3     -1     -1
    -1     0    -1     3     -1  
     -1     -1     -1     -1     4.

  This network has symmetry group Γ≅𝐃_4 generated by the permutations
 
    ρ= (1234),   π=(24).
 The conjugacy classes of isotropy subgroups for this network and corresponding cluster states are given in Table <ref>. Nodes belonging to the same cluster are the same colour. These are all possible cluster states arising from network symmetries. Figure <ref> shows the lattice of isotropy subgroups (cluster states) where arrows indicate (up to conjugacy) inclusion of the symmetries of the cluster state at the tail of the arrow within the symmetry group of the cluster state at the head of the arrow. This gives an indication of likely symmetry-breaking bifurcations between these cluster states.  











Sorrentino et al. <cit.> provide an algorithm for computing the subgroups of the network symmetry group Γ which give rise to cluster states: Suppose that S={ p_1, …, p_k} is a set of permutations which generates Γ. Partition S into subsets S= S_1 ∪⋯∪ S_ν such that the set of vertices moved by permutations in S_i is disjoint from the set of vertices moved by permutations in S_j for i≠ j, i,j =1, …ν, and each S_i cannot itself be partitioned in this way. If we let ℋ_i denote the subgroup generated by S_i then we obtain the geometric decomposition of the group Γ into the direct product of subgroups 
    Γ = ℋ_1 ×⋯×ℋ_ν,
 <cit.>. Note that for each node which is not permuted by any p_i ∈ S, we include a factor of ℋ_i in the geometric decomposition where S_i contains only the identity. Thus for the five-node example above Γ = ℋ_1 ×ℋ_2 where S_1 = {ρ, π} and S_2 = { e}. One can use the geometric decomposition to determine all subgroups of Γ which may give rise to cluster states by considering all groups of the type Σ = 𝒦_1 ×⋯×𝒦_ν where 𝒦_i is a subgroup of ℋ_i. Efficient computational algebra routines <cit.> are available which will compute the geometric decomposition and subgroups of a given arbitrary permutation group, so the computations can be carried out even for large networks with tens of thousands of nodes <cit.>.



 §.§ Cluster states from Laplacian coupling


In the general theory of patterns of synchrony in coupled cell systems <cit.>, cluster states correspond to balanced colourings of the network nodes (where synchronous cells have equivalent inputs accounting for different types of cells and inputs and also self-coupling). Networks whose topology does not have any symmetries can support balanced colourings, while networks with symmetry can have additional balanced colourings alongside those which result from symmetry. When we consider networks with Laplacian coupling dynamics, as in (<ref>), we expect to see additional potential cluster states resulting from balanced colourings which arise from the self-coupling. These additional states will include global synchronisation which is guaranteed to exist in networks of identical oscillators with Laplacian coupling.

An algorithm for computing these additional potential cluster states is given in <cit.> which uses as building blocks the clusters which can be found from network symmetries. First choose an isotropy subgroup Σ = 𝒦_1 ×⋯×𝒦_ν⊆Γ. This gives a partition of the nodes into clusters (a potential cluster state from symmetry). Next, merge some of these clusters as possibilities for new cluster states before finally checking if the merged cluster state (which cannot be one of those which arises due to symmetry) is dynamically valid. By this we mean that if we set all of the x⃗_i to be equal for nodes in the merged clusters, the equations of motion are consistent. This checking for consistency in the network equations can be done by eye for small networks. However, a better way to complete this final step for large networks is to automate the checking process using computational algebra tools <cit.>. To do this we use  the following steps which are described in greater detail in <cit.>. If the clusters we wish to merge are synchronised then the dynamics are equivalent to those which would be obtained if nodes in the same merged cluster were not connected (since the feedback term, H⃗(x⃗_i),  will cancel coupling terms of nodes from the same merged cluster). This gives a dynamically equivalent graph Laplacian matrix 𝒢 which is the original graph Laplacian matrix 𝒢 with off-diagonal entries between nodes in the same merged cluster set to 0 and the diagonals set to the negative of the new row sums. If, when we compute the symmetry group Γ of the network with this graph Laplacian, one of its isotropy subgroups corresponds to our merged cluster state then the dynamics are flow-invariant and our merged cluster state is a dynamically valid synchronised state. Sorrentino et al. <cit.> call the cluster states found in this way Laplacian clusters. Note that there is no need to check mergings between subgroups of the same group ℋ_i.

Example: Laplacian clusters

For the same five-node network as discussed previously, we can see that if we take isotropy subgroup Σ = ℤ_2(ρ^2) which gives the clusters { 1,3}, { 2, 4} and {5} then one potential Laplacian merged cluster state may be { 1,3,5}, {2,4}. This merging gives new graph Laplacian matrix
 
    𝒢 = [r] 2    -1    0     -1     0
     -1     3    -1     0     -1  
     0    -1    2     -1     0
    -1     0    -1     3     -1  
     0     -1     0     -1     2.
 A network with this connectivity matrix has symmetry group Γ generated by S= { (13), (35), (24)}. Thus the isotropy subgroup Γ gives the clusters {1,3,5}, {2,4} and therefore this cluster state is dynamically valid. It is called L3 in Table <ref> (along with its conjugate cluster state).

Another potential Laplacian merged cluster state may be { 1,2,5}, {3,4} from merging clusters from the symmetry cluster state A2 (see Table <ref>). In this case the merging gives a new graph Laplacian matrix 
    𝒢 = [r] 1    0    0     -1     0
     0     1    -1     0     0  
     0    -1    2    0     -1
    -1     0    0     2     -1  
     0     0     -1     -1     2,
 which corresponds to a network with symmetry group Γ generated by S= { (12)(34)}. Here the only isotropy subgroups are Γ which gives the cluster state A2 and the trivial subgroup which gives the cluster state A5. Thus this merging of clusters does not give a dynamically valid state.

Table <ref> gives a list of all dynamically viable Laplacian cluster states for the five-node network.  


 

For Laplacian coupled networks with symmetry, Sorrentino et al. <cit.> argue that their algorithm for computing cluster states which do not arise directly from symmetry is much more computationally efficient than algorithms not based on symmetries which compute all balanced colourings for any network topology, in particular that of Kamei and Cock <cit.>. They are careful to point out that while the checking stage (computing geometric decompositions and subgroups) is computationally efficient, the number of possible cluster mergings which must be checked grows combinatorially with the number μ of cluster states from symmetry with an upper bound of B_μ the μ-th Bell number. Thus the algorithm of Sorrentino et al. <cit.> is much faster than that of Kamei and Cock <cit.> when the number of cluster states from symmetry is small compared to the size of the network, that is μ≪ N. Both Sorrentino et al. <cit.> and Kamei and Cock <cit.> note that finding all dynamically valid cluster states for networks with Laplacian coupling may be substantially more difficult than for symmetric networks without self-coupling where connectivity is described by an adjacency matrix and all cluster states arise due to network symmetries.




 §.§ Stability of cluster states


The methods outlined in sections <ref> and <ref> provide a catalogue of possible cluster states which may exist within a given network with graph Laplacian coupling (<ref>). In applications, which of these states may exist and be stable will depend on the particular choices we make for the local dynamics F⃗ and the coupling function H⃗ as well as the global coupling strength σ. The presence of symmetry within the system imposes constraints on the form of the Jacobian matrix which can be used to greatly simplify stability calculations. Here we briefly review well-established methods for stability calculations within symmetric systems which apply to the cluster states which arise from network symmetries <cit.>. We also summarise the results of Sorrentino et al. <cit.> which extend these techniques to Laplacian cluster states.

First consider a periodic cluster state arising from network symmetry that has isotropy Σ⊆Γ. The fixed-point subspace of this subgroup is Υ = (Σ) which is the synchrony subspace for the cluster state. The state consists of M clusters 𝒞_k, k=1, …, M, where M = dim( Fix (Σ))= (Υ). Letting s⃗_k(t) denote the synchronised state of nodes in cluster 𝒞_k, and using the notation of <ref>, the variational equation of (<ref>) about the cluster state is

    tU⃗ =     [ ∑_k=1^M E^(k)⊗ D F⃗(s⃗_k(t)) . 
       . - σ∑_k=1^M ( 𝒢 E^(k)) ⊗ DH⃗(s⃗_k(t))  ] U⃗ ,

where E^(k) is the diagonal N × N matrix such that E_ii^(k) = 1 if i ∈𝒞_k and E_ii^(k) = 0 otherwise. To determine the stability of the periodic cluster state we need to compute the Floquet exponents of (<ref>). This task can be greatly simplified due to the fact that the system of variational equations can be block-diagonalised using the symmetries present in the system.

The node space can be decomposed into a number of irreducible representations of the isotropy subgroup Σ. Some of these subspaces will be isomorphic to each other and we combine these to obtain the isotypic components of the node space. Each isotypic component is invariant under the variational equation (<ref>) so Floquet exponents can be found by considering the restriction of the equations to each isotypic component. Thus the decomposition puts the variational equations into block-diagonal form and we then compute Floquet exponents for each block to determine stability of the cluster state. The process of isotypic decomposition and its use in stability computations is derived in detail in <cit.>. Pecora et al. <cit.> give an explicit algorithm to determine the isotypic decomposition for a given cluster state from symmetry and compute the transformation matrix Q such that 𝒢^' = Q 𝒢 Q^-1 is block diagonal. Applying this transformation to the variational equation (<ref>) we obtain block diagonal system of equations

    tV⃗ =    [∑_k=1^M J^(k)⊗ D F⃗(s⃗_k(t)) . 
       .- σ∑_k=1^M ( 𝒢^' J^(k)) ⊗ DH⃗(s⃗_k(t)) ] V⃗ ,

where V⃗(t) = (Q ⊗ I_m) U⃗(t) and J^(k) = Q E^(k) Q^-1. In Table <ref> we give the block-diagonalised graph Laplacian matrix for each of the cluster states from symmetry. An important point to note is that the isotypic component of the trivial representation is (Σ) = Υ, the synchronisation manifold. The block corresponding to this component is the M × M block which appears in the top left in each of the examples in Table <ref>. This block corresponds to perturbations within the synchronisation manifold and will always have a Floquet exponent equal to 0.  The remaining blocks correspond to the isotypic components of other irreducible representations of Σ. When the node space representation contains l≥ 1 isomorphic copies of a particular irreducible representation this will result in a block of dimension l. For example, see the 2 × 2 block in the block-diagonalisation for cluster state A2 in Table <ref>. These blocks represent perturbations transverse to the synchronisation manifold and Floquet multipliers from these blocks will determine stability under synchrony-breaking perturbations. For the cluster state to be stable all Floquet exponents (except the one which is always 0) must have negative real part.

In the case of a periodic Laplacian cluster state, the synchronisation manifold is an invariant subspace, but it is not the fixed-point subspace of any subgroup of Γ. However, we can still find a block-diagonalisation of 𝒢 which has in the top left a block which corresponds to perturbations within the synchronisation manifold. We do this following the algorithm of Sorrentino et al. <cit.>. Suppose that we start with a cluster state from symmetry with isotropy Σ which has M clusters and whose variational equations are block diagonalised by the matrix Q. Now suppose that a Laplacian cluster state is obtained by merging together two of the clusters in this state. The dimension of the synchronisation manifold decreases by one, while the dimension of the transverse manifold increases by one. New coordinates on the synchronisation manifold are obtained by transforming the new synchronisation vector in the node space coordinates (this has a 1 in the position of every node in the new merged cluster and 0s elsewhere) into the coordinates of the block-diagonalisation of the cluster state with isotropy Σ. The orthogonal complement provides the new transverse direction. Normalising the resulting vectors and entering them as rows of an orthogonal matrix Q^' whose other rows have Q^'_ij=δ_ij, we find that the matrix χ=Q^' Q block-diagonalises 𝒢 to a matrix 𝒢^'' which has top left block of dimension (M-1) × (M-1). Thus the transformation matrix χ  will block diagonalise the variational equations for the Laplacian cluster state so we may more easily determine the m(M-1) Floquet exponents within the synchronisation manifold and the m(M+1) transverse Floquet exponents.

Example: Block diagonalisation

We revisit the example five-node network of Sorrentino et al. <cit.>. The block-diagonalisations for the cluster state A3 and the merging which gives cluster state L3 are covered in <cit.>. Here we consider the cluster states A4 and L4 for illustration of the algorithms outlined above. The cluster state A4 has (choosing one option from the conjugacy class as they have identical existence and stability properties) isotropy subgroup ℤ_2(π). This group has only two irreducible representations, the trivial representation and the representation where π acts as multiplication by -1. The node space representation contains 4 isomorphic copies of the trivial representation and one copy of the other representation. Following the algorithm given explicitly in <cit.> we obtain

    Q = 1/2[r] 2    0   0     0    0
     0     -√(2)   0     -√(2)    0  
     0    0    2     0     0
    0    0    0     0     2  
     0     -√(2)    0    √(2)    0,
 and 𝒢^' as in Table <ref>.

Now consider the Laplacian cluster state L4 which is found by merging {2,4} with {5} to give a two cluster state {2,4,5}, {1,3}. The new synchronisation vector in node coordinates will be v⃗_1 = (0,1,0,1,1)^T and new transverse vector is v⃗_⃗2⃗ = (0,1,0,1, -2)^T. Normalising Qv⃗_1 and Qv⃗_2 we see that we must put these vectors as rows 2 and 4 of the orthogonal matrix Q^' to obtain

    Q^' = 1/3[r] 3    0   0     0    0
     0     -√(6)   0    √(3)    0  
     0    0    3     0     0
    0    -√(3)   0     -√(6)    0  
     0     0      0    0      3,
 and 𝒢^''= Q^'𝒢^' (Q^')^-1 as in Table <ref>. 






§ ORBITS AND VARIATIONAL EQUATIONS FOR NONSMOOTH SYSTEMS



Using the algorithms in section <ref> one can, given a network structure, compute the catalogue of cluster states (from symmetry and Laplacian coupling) which may be observed for given local dynamics F⃗ and interaction function H⃗. It is also possible, to block-diagonalise the variational equations needed to compute the stability of a given cluster state. However in practice, to determine which of the cluster states are stable and to carry out any kind of bifurcation analysis, we observe from (<ref>) that we will require closed form solutions for the periodic cluster state. For many of the situations in which one would wish to apply the framework, this information simply is not available. Therefore, we turn our attention to networks of PWL oscillators where, as we will show here, it is relatively straightforward to construct the periodic orbits for the cluster state and additionally apply the required modifications to the Floquet theory to account for the lack of smoothness of the dynamics.

Here we demonstrate how the computations of periodic cluster states and their stability may be carried out for a 2-cluster state in a network with 2-dimensional local dynamics and one switching plane. The method outlined here can easily be extended to larger numbers of clusters and more complex local PWL dynamics (e.g. with more switching planes and jump discontinuities) and we consider such examples in section <ref>.

Consider the specific case of a two cluster state L3 in the five-node example network with graph Laplacian given by (<ref>). Let each of the nodes have two-dimensional local dynamics given by the absolute model, which is named because it has one nullcline described by the absolute value function (and see <cit.> for a further discussion), where

    A_L=[ -1 -1;  1 -g ],  A_R=[  1 -1;  1 -g ],  c⃗_L=[      0; gw̅-v̅ ]= c⃗_R,
 with 0<g<1, gw̅<v̅ and indicator function h(z⃗; 0) = v. This model supports a non-smooth Andronov–Hopf bifurcation which occurs when the equilibrium moves from v<0 to v>0 crossing the switching manifold. When the (unstable) equilibrium lies in S_R, a periodic orbit for the single node dynamics can be constructed by connecting two trajectories, one in each of the half planes. Starting from initial data z⃗(0) = (0, w(0))^𝖳, this requires that we simultaneously solve the three nonlinear equations v(T_R)=0, v(Δ)=0, w(Δ)=w(0) where Δ is the period of the orbit, T_R is the time of flight in S_R and thus Δ>T_R>0. We do this numerically using Matlab's fsolve routine.  The periodic orbit can be shown to be stable with Floquet exponent <cit.>

    r = -g - Δ-2T_R/Δ.
 The nullclines for the model together with the stable periodic orbit are shown in Fig. <ref>.










Now couple these nodes together into a network described by (<ref>) where 𝒢 is given by (<ref>) and the coupling function is given by H⃗(z⃗) = (v,0)^𝖳 (diffusive coupling). Numerical simulations of the network equations with the parameter values of Fig. <ref> indicate that a periodic L3 cluster state exists with coupling strength σ=-0.03. On this cluster state z⃗_1 = z⃗_3 =z⃗_5 = s⃗_1 and z⃗_2 = z⃗_4= s⃗_2. (Recall that the conjugate cluster state where z⃗_1 = z⃗_3 and z⃗_2 = z⃗_4= z⃗_5 has the same existence and stability properties.) On this invariant subspace the equations have the form ṡ⃗̇ = A_μ_1, μ_2s⃗ + c⃗_μ_1, μ_2 where s⃗ = (s⃗_1, s⃗_2)^𝖳,

    A_μ_1, μ_2 = [ A_μ_1 - 2 σ DH⃗         2 σ DH⃗;         3 σ DH⃗  A_μ_2 -3 σ DH⃗ ],  c⃗_μ_1, μ_2 = [ c⃗_μ_1; c⃗_μ_2 ],
 and 
    μ_i=   L,     v_i<0
     R,     v_i>0 .

This is simply a 4-dimensional PWL system with two switching planes, v_1=0 and v_2=0. The periodic orbit on the 4-dimensional synchronous manifold can again be constructed by connecting together trajectories which satisfy the linear ODEs in each region across the switching planes. Starting from initial data s⃗(0) = (0, w_1(0), v_2(0), w_2(0))^𝖳, we now have to solve a system of seven nonlinear algebraic equations for w_1(0), v_2(0), w_2(0) and the four switching times T_1,1, T_2,1, T_2,2 and T_1,2= Δ. See Fig. <ref>.





The initial data and switching times can now be used to compute explicitly the Floquet multipliers for the periodic orbit, making use of the block-diagonalisation of the variational equations (<ref>). Multipliers corresponding to perturbations within the synchronisation manifold can be computed without using the block-diagonalisation. We have

    tδs⃗ = A_μ_1, μ_2δs⃗ ,
 which can be solved using matrix exponentials, being careful to use saltation matrices to evolve perturbations through switching manifolds. After one period 
    δs⃗(Δ) = Ψ_sδs⃗(0) ,
 where from Fig. <ref>,

    Ψ_s =    K_12ℰ_L,R(Δ- T_1,1) K_11ℰ_R,R(T_1,1-T_2,2) 
       × K_22ℰ_R,L(T_2,2-T_2,1)K_21ℰ_R,R(T_2,1),

with saltation matrices

    K_ij = P_i ⊗ K_i(T_i,j),   
    P_1 = [ 1 0; 0 0 ],    P_2 = [ 0 0; 0 1 ],


    K_i(t) = [           v̇_i(t^+)/v̇_i(t^-)                             0; (ẇ_i(t^+)-ẇ_i(t^-))/v̇_i(t^-)                             1 ],    i=1,2,

and 
    ℰ_μ_1, μ_2(t) = ^A_μ_1, μ_2t,     μ_i ∈{ L, R}.

The Floquet multipliers are the eigenvalues of Ψ_s. One of these eigenvalues will always be 1 corresponding to perturbations along the periodic orbit.

In the directions transverse to the synchronisation manifold, the block-diagonalisation of 𝒢 results in the following three decoupled Floquet problems:

    V̇⃗̇_3    = (D F⃗(s⃗_1) -3σ D H⃗)V⃗_3,  
    V̇⃗̇_4    = (D F⃗(s⃗_2) -3σ D H⃗)V⃗_4,  
    V̇⃗̇_5    = (D F⃗(s⃗_1) -5σ D H⃗)V⃗_5,
 which again can be solved using matrix exponentials and saltation matrices: V⃗_i(Δ) = Ψ_t_iV⃗_i(0) where

    Ψ_t_3   = K_1(Δ) ℰ_L^3(Δ-T_1,1) K_1(T_1,1) ℰ_R^3(T_1,1)  ,  
    Ψ_t_4   =ℰ_R^3(Δ- T_2,2) K_2(T_2,2) ℰ_L^3(T_2,2-T_2,1)
       × K_2(T_2,1) ℰ_R^3(T_2,1) ,  
    Ψ_t_5   = K_1(Δ) ℰ_L^5(Δ-T_1,1) K_1(T_1,1) ℰ_R^5(T_1,1) ,

and 
    ℰ_μ^β(t) = ^(A_μ-βσ DH⃗)t.

The eigenvalues of the Ψ_t_j, j ∈3,4,5 give the Floquet multipliers for directions transverse to the synchronisation manifold.

Note that the change of basis from U⃗ coordinates to V⃗ coordinates has no effect on the action of the saltation matrices:
Recall V⃗ = (Q ⊗ I_m) U⃗. To progress U⃗ through a switch or discontinuity then U⃗^+ = K U⃗^-, where

    K = ∑_k=1^M E^(k)⊗ K_k.

Thus V⃗^+ =KV⃗ where

    K   = (T ⊗ I_m) K (T⊗ I_m)^-1  = ∑_k=1^M (T E^(k) T^-1⊗ K_k) 
        = ∑_k=1^M (J^(k)⊗ K_k).



Since the vector field for the absolute model is continuous, in this case all saltation matrices are the identity. Figure <ref> shows the Floquet multipliers for the choice of parameter values as in Fig. <ref>. As they all lie inside the unit disc (with the exception of the one which is forced to be 1) the periodic orbit, as shown in the left hand panel of Fig. <ref> is stable.






We can locate bifurcations of the periodic orbit by determining when Floquet multipliers leave the unit disc (noting that the ordering of times which the trajectories cross the switching planes may also change as parameters are varied). Treating all cluster types we can build up bifurcation diagrams. For the absolute model with the parameters for the local dynamics as given in Fig. <ref>, bifurcations which can be found upon varying coupling strength σ are shown in Fig. <ref>. All bifurcations from stable states are of tangent type where a Floquet multiplier passes through +1.










§ FURTHER APPLICATIONS AND EXAMPLES


Extending the techniques illustrated in section <ref>, here we discuss the cluster patterns and bifurcations which occur for different PWL local dynamics on the five node network (<ref>). We consider a diffusively coupled network of PWL Morris–Lecar oscillators (which has more than one switching plane) and the integrate-and-fire network model with event-driven synaptic coupling as discussed for the case of global synchrony in section <ref> (which has jump discontinuities within the vector field). We show that each model of local dynamics gives a very different variety of cluster state dynamics demonstrating that emergent network behaviour depends as much on node dynamics as network structure.




 §.§  Piecewise linear Morris–Lecar network


We now investigate dynamics of a network of planar PWL nodes with more than one switching plane. Consider the five-node network given by graph Laplacian coupling matrix (<ref>) where the local dynamics F⃗ are given by a PWL reduction of the Morris–Lecar model (PML) <cit.>, and the coupling function is again given by H⃗(z⃗) = (v,0)^𝖳 (diffusive coupling). Figure <ref> shows the phase plane for this model in the regime where there are type I oscillations (emerging with arbitrarily low frequency from a homoclinic bifurcation) and bistability between a periodic orbit and a stable fixed point.  The dynamics of the PWL model are defined by

    C v̇ =  f(v)-w +I,     ẇ =  g(v,w),

where

    f(v) = 
    -v,     v<a/2
    
    v-a,     a/2 ≤ v ≤ (1+a)/2
    
    1-v,     v > (1+a)/2
     ,

and

    g(v,w) = 
    ( v-γ_1 w + b^*γ_1 -b)/γ_1,     v<b
    
    ( v-γ_2 w + b^*γ_2 -b)/γ_2,     v≥ b
     ,

with a/2 ≤ b^* ≤ (1-a)/2 and a/2 ≤ b≤ (1+a)/2.





The periodic orbit (as shown in Fig. <ref>) consists of four pieces labelled by μ=1,…, 4. On each section z⃗_μ(t) = ℰ_μ(t) z⃗_μ(0) + 𝒦_μ(t) c⃗_μ
where ℰ_μ(t) = ^A_μ t, and 𝒦_μ(t) = A_μ^-1[ℰ_μ(t)-I_m] with
A_3=A_1, c⃗_3=c⃗_1 and

    A_1     = [   1/C  -1/C; 1/γ_2    -1 ],     A_2    = [  -1/C  -1/C; 1/γ_2    -1 ],
    
    A_4     = [   1/C  -1/C; 1/γ_1    -1 ] ,    c⃗_1    = [     (I-a)/C; b^* - b/γ_2 ], 
    c⃗_2    = [     (1+I)/C; b^* - b/γ_2 ],    c⃗_4    = [     (I-a)/C; b^* - b/γ_1 ].

The periodic orbit crosses two switching planes with indicator functions h_1(z⃗; 0) = v-b and h_2(z⃗; 0) = v-(1+a)/2. Therefore for initial data z⃗(0) = (b, w(0)), to find the periodic orbit we must solve the system of five nonlinear algebraic equations

    v(T_1)= (1+a)/2,       v(T_1+T_2) = (1+a)/2 ,
       v(T_1+T_2 + T_3)=b,       v(Δ) = b,       w(0) = w(Δ) ,
 for the switching times T_1, T_2, T_3, the period Δ and the initial value w(0). 
This periodic orbit can be shown to be stable with Floquet exponent <cit.>

    r = -1 + (Δ-2T_2)/(CΔ).

Note that other types of periodic solutions are possible, such as those that cross only v=b or v=(1+a)/2, or ones that do not cross either of the switching manifolds (as for a periodic orbit surrounding a linear center).


When nodes with PML local dynamics are diffusively coupled with the structure given by graph Laplacian (<ref>) a number of different stable cluster states are observed. Similarly to the case of node dynamics described by the absolute model as discussed in section <ref>, we can locate periodic orbits for cluster states and carry out linear stability calculations. We can then produce bifurcation diagrams such as that given in Fig. <ref>. Note that just as for the absolute model local dynamics, the PML model has a continuous vector field and therefore again all saltation matrices (<ref>) are the identity.
In practice, dealing with changes in the ordering of event crossings (defining a cluster) and to gain insight into what types of cluster state are possible (to provide good initial guesses for a numerical root-finding scheme) we also did simulations of a corresponding set of smooth ODEs, with a sufficiently steep switch to mimic the PML system.  






The bifurcation diagram in Fig. <ref> shows the presence of 6 of the possible 8 cluster states. All but the L4 branch possess a range of values of σ for which the periodic cluster state is stable. In the region 0.033 ≲σ≲ 0.071 where there are no stable branches of periodic cluster states, aperiodic behaviour such as that shown in the inset of Fig. <ref> dominates. There is also a stable multirhythm,
where nodes can oscillate at different frequencies (not shown in Fig. <ref>). This has x⃗_1(t) = x⃗_2(t+Δ/2), x⃗_3(t) = x⃗_4(t+Δ/2) and x⃗_5(t) = x⃗_5(t+Δ/2) so x⃗_5 has a period of oscillation which is half that of x⃗_i, i=1,…,4. The amplitudes of x⃗_1 and x⃗_5 are also very small resulting in a much smaller L_2-norm than the other solutions shown in Fig. <ref>.




 §.§  Piecewise linear Integrate-and-fire network


As our final example we consider a network of synaptically coupled PWL-IF neurons as discussed for the case of synchrony in section <ref>. We continue to use the same five-node structured network (<ref>). We also choose for the local dynamics parameter values as in the top left of Fig. <ref> that give tonic firing of each of the nodes when the coupling strength σ is 0. Setting the value of the synaptic rate parameter α=0.4 we observe a number of different stable cluster states depending on the value of the coupling strength σ. A selection of these cluster states is given in
Fig. <ref>. For these parameter values, the (tonic firing) synchronous state (L1) is stable for 0 ≲σ≲ 0.0334. This undergoes a period-doubling bifurcation at σ_pd≃ 0.0334 to an A3 cluster state of doublets such as that pictured in Fig. <ref> (a) for σ≃ 0.038. For 0.0395 ≲σ≲ 0.0485 this solution branch (as pictured for σ≃ 0.043 in Fig. <ref> (b)) is bistable with the periodic 4 cluster state A4 where x⃗_1 and x⃗_3 are not equal, but so close that it is hard to distinguish their time series and orbits (shown in blue and orange respectively in Fig. <ref> (c) for σ≃ 0.043). Further increasing σ beyond 0.0485 leads to a loss of periodicity of solutions (both aperiodic A4 and A5 solution branches are observed). At around σ=0.07 again we see a stable periodic A3 cluster state for which all orbits now pass through the switch at v_i=0. This is pictured in Fig. <ref> (d) for σ≃  0.08 where the solution is bistable with the two cluster L3 solution as in Fig. <ref> (e). The A3 cluster state remains stable for σ≳ 0.07 undergoing a series of spike adding bifurcations for large values of σ. Fig. <ref> (f) shows the solution branch when σ=8 where each cluster spikes three times per orbit.  We note that recent work by Chen et al. <cit.> has also considered cluster states in leaky IF networks with α-function synapses, although their results are restricted to weak-coupling in small all-to-all coupled networks.










§ DISCUSSION


In this paper we have shown how to extend the well known MSF formalism for networks of smooth dynamical systems to treat PWL models relevant to the modelling of coupled neural oscillators.  Moreover, we have also augmented the recent extension of the MSF approach that treats synchronous cluster states (exploiting spatial or pointwise network symmetries and balanced colourings).  We have illustrated the general theory using some relatively simple networks, though with various choices of node dynamics.  This has highlighted that the degree of synchronisation or type of clustering pattern observed in neuronal networks may have as much to do with the choice of connectivity as it has  with the choice of node dynamics.  This reinforces the notion that care must be taken when trying to recover structural connectivity from functional connectivity <cit.>.  Similarly it is known that the reconstruction of network dynamics from observables can be strongly influenced by symmetries <cit.>.  Although we have focused on example systems drawn from theoretical neuroscience, there are many PWL models in physics and engineering for which the framework we have presented can also be applied, such as electromechanical oscillator networks <cit.>.
Interestingly, recent work by Cho et al. <cit.> has made a connection between synchronised cluster states and chimeras <cit.>, where a sub-population of oscillators synchronises in an otherwise incoherent sea.  By exploiting the symmetries of the network they show how to construct cluster partitions, in which the synchronisation stability of the individual clusters in each set is decoupled from that in all the other sets.  Using this they show that a chimera can exist if the fully synchronous solution is unstable and there is a cluster state with a main stable group, and all other clusters unstable.  Thus it would be of interest to revisit this work, formulated in a smooth setting, using the PWL perspective presented here.

It is well to mention that we have only considered networks of identical oscillators, and it is well known that heterogeneity can disrupt synchrony <cit.>.  The MSF extension to nearly identical units is relatively straightforward <cit.>, though for PWL systems there is also the further possibility of the explicit construction of network states (by patching together local trajectories).  However, the MSF is not the only tool for analysing cluster states and it would be interesting to consider the augmentation of Lyapunov function methods, such as those developed by Belyk et al. <cit.>, to treat nonsnmooth systems with event driven interactions.  In a neural context it would also be important to treat delays (arising from the finite speed of action potential propagation), and there is a growing body of work in this area that could be revisited from a PWL perspective <cit.>.  Indeed delays are well known to lead to oscillatory instabilities, which at the network level may manifest as travelling waves.  This highlights an important mathematical extension of the work presented here, namely that periodic states have phase shift symmetries which must also be accounted for. When combined with the pointwise network symmetries these lead to spatio-temporal symmetries <cit.>.

A periodic solution x⃗(t) has a group K of spatial symmetries (the γ∈Γ such that γx⃗(t) = x⃗(t) for all t). It is these symmetries which we have discussed in this paper. However, periodic orbits also possess a group H of spatial symmetries under which the periodic orbit is invariant (that is the γ∈Γ such that {γx⃗(t)  :  t ∈ℝ} ={x⃗(t)  :  t ∈ℝ}). Thus there exists a phase shift ϕ∈𝕊^1 such that γx⃗(0) = x⃗(ϕ) and therefore γx⃗(t-ϕ) = x⃗(t) for all t. The pair (γ, ϕ) define a spatio-temporal symmetry of the periodic orbit. The group K is normal in H and H/K is a finite subgroup of 𝕊^1. If H=K then there are no phase-shift symmetries. The pair (H,K) determines the pattern of phase shifts (up to isomorphism). Patterns  of  phase  shifts  which may exist for a given network structure with symmetry group Γ can be determined using the H/K theorem <cit.>. This theorem gives conditions on subgroups K ⊆ H ⊆Γ such that the pair (H,K) could determine the spatio-temporal symmetries of a hyperbolic periodic solution of the Γ-equivariant network equations, allowing one to compute a catalogue of possible phase shifted states. These may include so called multirhythm states (such as that found here for the PWL Morris–Lecar local dynamics on the five-node network, see section <ref>) where some cells are forced to oscillate with a frequency that is a rational multiple of that of the other cells <cit.>. Other phase shifted cluster states may also arise due to additional constraints imposed by the choice of Laplacian or balanced coupling. There are also approaches exploiting spatio-temporal symmetries for computing stability of periodic states near Hopf bifurcation <cit.> which may also be extendable to phase-shifted cluster states from Laplacian or balanced coupling.  The extension of our work to treat spatio-temporal symmetries will be presented elsewhere.




The authors would like to thank Yi Ming Lai, Mason Porter, and Rüdiger Thul for helpful comments made on a first draft of this paper.
This work was supported by the Engineering and Physical Sciences Research Council [grant number EP/P007031/1].




§ SALTATION MATRICES





Saltation matrices are a useful way to augment the standard theory of smooth dynamical systems to treat the evolution of perturbations through regions in phase phase where there is a discontinuity in the vector field or the flow.  They have recently been used in a PWL network setting in <cit.>, and for nonlinear IF networks in <cit.>.  Here we combine these two perspectives to construct the saltation matrices for hybrid systems that may have switches in their dynamics as well as discontinuities arising from reset.  Given that these cannot occur at the same time we adopt a notation to track either of these types using an indicator function h(z⃗) such that an event (switching or reset) occurs when h(z⃗(T))=0 (with different choices of h for switching and reset).
The state of the system after an event is given by z⃗(T^+)=g⃗(z⃗(T)). Let us now consider an arbitrary trajectory of a dynamical system ż⃗̇ = F⃗(z⃗) and linearise around this trajectory.  Between switching or firing events small perturbations δz⃗ satisfy

    tδz⃗(t) = D F⃗(z⃗(t)) δz⃗(t),   δz⃗(0) = δz⃗_0 .

For the PWL models discussed in this paper then D F⃗(z⃗(t)) is a piecewise constant matrix, so that
δz⃗(t) = exp(A t) δz⃗_0 with A a constant matrix (such as A_L,R in <ref>).
A perturbed trajectory defined by z⃗(t) can be constructed as z⃗(t) = z⃗(t)+δz⃗(t). We can then define a corresponding perturbed event time T according to the condition h(z⃗(T))=0.  We shall denote the difference in the unperturbed and perturbed event times as δ T = T-T. For the case that δ T <0 we have that
z⃗(T^+)=z⃗(T^+-δ T) ≃z⃗(T^+) - ż⃗̇(T^+) δ T, and see Fig. <ref>.
Hence,

    z⃗(T^+) ≃g⃗(z⃗(T^-)) - ż⃗̇(T^+)δT 
       ≃g⃗(z⃗(T^-)+ż⃗̇(T^-)δT) - ż⃗̇(T^+)δT 
       ≃g⃗(z⃗(T^-)+δz⃗(T^-)+ż⃗̇(T^-)δT) - ż⃗̇(T^+)δT .

By performing a Taylor expansion of g⃗ we obtain:

    z⃗(T^+)    ≃z⃗(T^+) + D g⃗ (z⃗(T^-))
    [δz⃗(T^-)+ż⃗̇(T^-)δT] 
       - ż⃗̇(T^+)δT ,

where we have used the result that g⃗(z⃗(T^-))=z⃗(T^+).  To calculate δ T we Taylor expand the indicator function for the perturbed trajectory as

    h(z⃗(T^-)) = h(z⃗(T^-+δ T) + δz⃗ (T^-+δ T) ) 
       ≃ h(z⃗(T^-)+ ż⃗̇ (T^-)δ T + δz⃗ (T^-) ) 
       ≃ h(z⃗(T^-)) + ∇_z⃗ h(z⃗(T^-)) · [ ż⃗̇ (T^-)δ T + δz⃗ (T^-)].

Using the fact that h(z⃗(T^-)) = 0=h(z⃗(T^-)) we have that

    δ T = -  ∇_z⃗ h(z⃗(T^-)) ·δz⃗ (T^-)/∇_z⃗ h(z⃗(T^-)) ·ż⃗̇ (T^-) .

Now making use of the approximations ż⃗̇(T^-) ≃ż⃗̇(T^-) and ż⃗̇(T^+) ≃ż⃗̇(T^+) (and see Fig. <ref>), we may use (<ref>) and (<ref>) to construct δz⃗(T^+) in the form

    δz⃗(T^+) ≃ K(T) δz⃗(T^-),

where K(T) is the saltation matrix:

    K(T) = D g⃗(z⃗(T^-)) 
        +    [ż⃗̇(T^+) - D g⃗(z⃗(T^-)) ż⃗̇(T^-)  ]  [ ∇_z⃗ h(z⃗(T^-))  ]^𝖳/∇_z⃗ h(z⃗(T^-)) ·ż⃗̇ (T^-) .

The saltation matrix allows us to relate perturbations just before an event to those just after.  A similar calculation with δ T >0 gives precisely the same formula for K(T).  We note that for a switching event (where reset does not occur) then g⃗(z⃗)=z⃗ and D g⃗(z⃗) is the identity matrix.
Moreover, if in this case the vector field is continuous then ż⃗̇(T^+)=ż⃗̇(T^-) and K(T) reduces to the identity matrix as expected.

For the single PWL-IF node z⃗∈^2 with ∇_z⃗ h(z⃗;v_th)=(1,0)^𝖳, and

    D g⃗(z⃗) = [ 0 0; 0 1 ].

For the PWL-IF network  z⃗_i ∈^4 with ∇_z⃗ h(z⃗_i;v_th)=(1,0,0,0)^𝖳, and

    D g⃗(z⃗_i) = [ 0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1 ].





57
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Arenas et al.(2008)Arenas,
  Díaz-Guilera, Kurths, Moreno, and Zhou]Arenas2008
  author author A Arenas, author A Díaz-Guilera, author J Kurths, author Y Moreno,   and author C Zhou, title title Synchronization in complex networks,  @noop  journal journal Physics Reports  volume 469, pages 93–153 (year 2008)NoStop
[Abrams et al.(2016)Abrams,
  Pecora, and Motter]Abrams2016
  author author Daniel M.  Abrams, author Louis M.  Pecora,  and author Adilson E.  Motter, title title
  Introduction to focus issue: Patterns of network synchronization,  @noop  journal journal Chaos: An
  Interdisciplinary Journal of Nonlinear Science volume
  26, pages 094601 (year 2016)NoStop
[Brown et al.(2003)Brown,
  Holmes, and Moehlis]Brown03
  author author E Brown, author P Holmes,  and  author J Moehlis, title Perspectives and Problems in Nonlinear Science: A
  Celebratory Volume in Honor of Larry Sirovich,  (publisher
  Springer, address New York, year 2003) Chap.  chapter Globally coupled oscillator networks, pp. pages 183–215NoStop
[Ashwin et al.(2007)Ashwin,
  Orosz, Wordsworth, and Townley]Ashwin07
  author author P Ashwin, author G Orosz,
  author J Wordsworth,  and  author S Townley, title title Dynamics on networks of cluster states
  for globally coupled phase oscillators, @noop  journal journal SIAM Journal on Applied Dynamical Systems  volume 6, pages 728–758 (year 2007)NoStop
[Golubitsky and Stewart(1986)]Golubitsky1986
  author author M Golubitsky and author I N Stewart, title Multiparameter bifurcation
  theory,  (publisher American Mathematical Society,  address Providence RI, year 1986) Chap. chapter Hopf bifurcation with dihedral group symmetry: coupled nonlinear
  oscillators, pp. pages 131–173NoStop
[Pogromsky et al.(2002)Pogromsky, Santoboni, and Nijmeijer]Pogromsky2002
  author author A Pogromsky, author G Santoboni,
   and author H Nijmeijer,  title title Partial synchronization:
  from symmetry towards stability, @noop  journal
  journal Physica D volume 172,  pages 65–87 (year 2002)NoStop
[Pogromsky(2008)]Pogromsky2008
  author author A Y Pogromsky, title title A partial
  synchronization theorem, @noop  journal journal Chaos volume 18, pages
  037107 (year 2008)NoStop
[Belykh et al.(2000)Belykh,
  Belykh, and Hasler]Belykh2000
  author author V N Belykh, author I V Belykh,   and author M Hasler, title title Hierarchy and stability of partially
  synchronous oscillations of diffusively coupled dynamical systems,  @noop  journal journal Physical Review
  E volume 62, pages 6332–6345
  (year 2000)NoStop
[Pecora et al.(2014)Pecora,
  Sorrentino, Hagerstrom, Murphy, and Roy]Pecora2014
  author author L M Pecora, author F Sorrentino,
  author A M Hagerstrom, author T E Murphy,  and author R Roy, title
  title Cluster synchronization and isolated
  desynchronization in complex networks with symmetries, @noop 
  journal journal Nature Communications volume 5 (year 2014)NoStop
[Sorrentino et al.(2016)Sorrentino, Pecora, Hagerstrom,
  Murphy, and Roy]Sorrentino2016
  author author F Sorrentino, author L M Pecora, author A M Hagerstrom, author T E Murphy,  and author R Roy,  title title Complete characterization
  of the stability of cluster synchronization in complex dynamical
  networks, @noop  journal journal
  Science Advances volume 2, pages
  e1501737 (year 2016)NoStop
[Pecora and Carroll(1998)]Pecora1998
  author author L M Pecora and author T L Carroll, title title Master
  stability functions for synchronized coupled systems, @noop 
  journal journal Physical Review Letters  volume 80, pages 2109–2112 (year 1998)NoStop
[Kiss et al.(2007)Kiss,
  Rusin, Kori, and Hudson]Kiss2007
  author author I Z Kiss, author C G Rusin,
  author H Kori,  and author J L Hudson, title
  title Engineering complex dynamical structures:
  Sequential patterns and desynchronization, @noop  journal journal Science volume
  316, pages 1886–1889 (year 2007)NoStop
[Zhang et al.(2009)Zhang,
  Yuan, and Zhou]Zhang2009
  author author J Zhang, author Z Yuan,  and  author T Zhou, title title Synchronization and clustering of
  synthetic genetic networks: A role for cis-regulatory modules, @noop
   journal journal Physical Review E volume 79, pages 041903 (year
  2009)NoStop
[Zemonová et al.(2006)Zemonová, Zhou, and Kurths]Zemonova2006
  author author L Zemonová, author C Zhou,   and author J Kurths, title title Structural and functional clusters of
  complex brain networks, @noop  journal journal Physica D volume 224, pages
  202–212 (year 2006)NoStop
[Porter and Gleeson(2016)]Porter2016
  author author M A Porter and author J P Gleeson, @noop  title Dynamical Systems on
  Networks: A Tutorial., series Frontiers in Applied Dynamical
  Systems: Reviews and Tutorials, Vol. volume 4 (publisher Springer International Publishing, address
  Switzerland, year 2016)NoStop
[di Bernardo et al.(2008)di Bernardo, Budd, Champneys, and  Kowalczyk]diBernardo2008
  author author M di Bernardo, author C Budd,
  author A R Champneys,  and  author P Kowalczyk, @noop
   title Piecewise-smooth Dynamical Systems: Theory and
  Applications, Applied Mathematical Sciences (publisher
  Springer, year 2008)NoStop
[Glendinning(2015)]Glendinning2015
  author author P Glendinning, title The Princeton Companion to
  Applied Mathematics,  (publisher Princeton University
  Press, year 2015) Chap. chapter II.18 Hybrid
  Systems, pp. pages 103–104NoStop
[Coombes and Thul(2016)]Coombes2016
  author author S Coombes and author R Thul,  title title Synchrony in networks of
  coupled non-smooth dynamical systems: Extending the master stability
  function, @noop  journal journal
  European Journal of Applied Mathematics volume
  27, pages 904–922 (year 2016)NoStop
[van Vreeswijk and Sompolinsky(1996)]VanVreeswijk1996
  author author C van
  Vreeswijk and author H Sompolinsky, title title Chaos in
  Neuronal Networks with Balanced Excitatory and Inhibitory Activity,  @noop  journal journal Science volume 274, pages 1724–1726 (year
  1996)NoStop
[Denève and Machens(2016)]Deneve2016
  author author S Denève and author C K Machens, title title Efficient codes
  and balanced networks, @noop  journal journal Nature Neuroscience volume 19, pages 375–382 (year 2016)NoStop
[Roudi and Latham(2007)]Roudi2007
  author author Y Roudi and author P E Latham, title title A balanced
  memory network, @noop  journal journal
  PLOS Computational Biology volume 3, pages 1–22 (year 2007)NoStop
[Rubin et al.(2017)Rubin,
  Abbott, and Sompolinsky]Rubin2017
  author author R Rubin, author L F Abbott,   and author H Sompolinsky,  title title Balanced excitation and
  inhibition are required for high-capacity, noise-robust neuronal
  selectivity, @noop  journal journal
  Proceedings of the National Academy of Sciences volume 114, pages E9366–E9375 (year
  2017)NoStop
[Doiron and Litwin-Kumar(2014)]Doiron2014
  author author B Doiron and author A Litwin-Kumar, title title Balanced
  neural architecture and the idling brain, @noop  journal journal Frontiers in Computational Neuroscience  volume 8, pages 56 (year
  2014)NoStop
[Coombes et al.(2012)Coombes, Thul, and Wedgwood]Coombes2011a
  author author S Coombes, author R Thul,  and  author K C A Wedgwood,  title title Nonsmooth dynamics in
  spiking neuron models, @noop  journal journal Physica D volume 241, pages
  2042–2057 (year 2012)NoStop
[Coombes(1999)]Coombes99a
  author author S Coombes, title title Liapunov
  exponents and mode-locked solutions for integrate-and-fire dynamical
  systems, @noop  journal journal
  Physics Letters A volume 255, pages
  49–57 (year 1999)NoStop
[Ladenbauer et al.(2013)Ladenbauer, Lehnert, Rankoohi,
  Dahms, Schöll, and Obermayer.]Ladenbauer2013
  author author J Ladenbauer, author J Lehnert,
  author H Rankoohi, author T Dahms, author
  E Schöll,  and author
  K Obermayer., title title Adaptation controls synchrony and cluster states of
  coupled threshold-model neurons,, @noop  journal
  journal Physical Review E volume 88,  pages 042713 (year 2013)NoStop
[Naud et al.(2008)Naud,
  Marcille, Clopath, and Gerstner]Naud2008
  author author R Naud, author N Marcille,
  author C Clopath,  and author W Gerstner, title
  title Firing patterns in the adaptive exponential
  integrate-and-fire model, @noop  journal journal Biological Cybernetics volume 99,  pages 335–347 (year 2008)NoStop
[Gröbler et al.(1998)Gröbler, Barna, and Érdi]Grobler98
  author author T Gröbler, author G Barna,   and author P Érdi, title title Statistical model of the hippocampal
  CA3 region I. The single-cell module: bursting model of the pyramidal
  cell, @noop  journal journal
  Biological Cybernetics volume 79, pages 301–308 (year 1998)NoStop
[Izhikevich(2003)]Izhikevich03
  author author E M Izhikevich, title title Simple model
  of spiking neurons, @noop  journal journal IEEE Transactions On Neural Networks volume
  14, pages 1569–1572 (year 2003)NoStop
[Badel et al.(2008)Badel,
  Lefort, Berger, Petersen,
  Gerstner, and Richardson]Badel08
  author author L Badel, author S Lefort,
  author T K Berger, author C C H Petersen, author
  W Gerstner,  and author
  M J E Richardson, title
  title Extracting nonlinear integrate-and-fire models
  from experimental data using dynamic I-V curves, @noop 
  journal journal Biological Cybernetics volume 99, pages 361–370 (year
  2008)NoStop
[Coombes and Zachariou(2009)]Coombes09a
  author author S Coombes and author M Zachariou, title Gap junctions and emergent
  rhythms, in @noop  booktitle Coherent
  Behavior in Neuronal Networks, series Springer Series in
  Computational Neuroscience, Vol. volume 3, editor
  edited by editor K Josić,
  editor J Rubin, editor
  M Matias,  and editor
  R Romo (publisher Springer,  year 2009) pp. pages 77–94NoStop
[Filippov(1988)]Filippov1998
  author author A F Filippov, @noop  title Differential Equations
  with Discontinuous Righthand Sides (publisher Kluwer
  Academic Publishers, address Norwell, year
  1988)NoStop
[Aizerman and Gantmakher(1958)]Aizerman1958
  author author M A Aizerman and author F R Gantmakher, title title On the
  stability of periodic motions, @noop  journal
  journal Journal of Applied Mathematics and Mechanics volume 22, pages 1065–1078 (translated from
  Russian) (year 1958)NoStop
[Belykh and Hasler(2011)]Belykh2011
  author author I Belykh and author M Hasler,  title title Mesoscale and clusters of
  synchrony in networks of bursting neurons, @noop  journal journal Chaos: An Interdisciplinary Journal of
  Nonlinear Science volume 21, pages
  016106 (year 2011)NoStop
[Coombes(2000)]Coombes2000
  author author S Coombes, title title Chaos in
  integrate-and-fire dynamical systems, in @noop  booktitle Proceedings of Stochastic and Chaotic Dynamics in the Lakes,  Vol. volume 502 (organization American Institute
  of Physics Conference Proceedings, year 2000) pp. pages 88–93NoStop
[Golubitsky and Stewart(2016)]Golubitsky2016
  author author M Golubitsky and author I N Stewart, title title Rigid patterns
  of synchrony for equilibria and periodic cycles in network dynamics,  @noop  journal journal Chaos: An
  Interdisciplinary Journal of Nonlinear Science volume
  26, pages 094803–14 (year 2016)NoStop
[MacArthur et al.(2008)MacArthur, Sanchez-Garcia, and Anderson]MacArthur2008
  author author B D MacArthur, author R J Sanchez-Garcia,  and author J W  Anderson, title title
  Symmetry in complex networks, @noop  journal
  journal Discrete and Applied Mathematics volume 156, pages 3525 (year
  2008)NoStop
[GAP()]GAP4
  GAP,  title
  GAP – Groups, Algorithms, and Programming, Version 4.8.7, organization The GAP Group (year 2017)NoStop
[Developers(2017)]sagemath
  author author The Sage  Developers, @noop  title
  SageMath, the Sage Mathematics Software System (Version 7.6)
  (year 2017), note http://www.sagemath.orgNoStop
[Kamei and Cock(2013)]Kamei2013
  author author H Kamei and author P J A Cock, title title Computation of
  balanced equivalence relations and their lattice for a coupled cell
  network, @noop  journal journal SIAM
  Journal on Applied Dynamical Systems volume 12,  pages 352–382 (year 2013)NoStop
[Golubitsky et al.(1988)Golubitsky, Stewart, and Schaeffer]Golubitsky1988
  author author M Golubitsky, author I Stewart,
   and author D G Schaeffer,  @noop  title Singularities and Groups in
  Bifurcation Theory, Volume II (publisher Springer Verlag,  year 1988)NoStop
[Coombes(2008)]Coombes2008
  author author S Coombes, title title Neuronal
  Networks with Gap Junctions: A Study of Piecewise Linear Planar Neuron
  Models, @noop  journal journal SIAM
  Journal on Applied Dynamical Systems volume 7,  pages 1101–1129 (year 2008)NoStop
[Chen et al.(2017)Chen,
  Engelbrecht, and Mirollo]Chen2017
  author author Bolun Chen, author Jan R. Engelbrecht,  and author Renato  Mirollo, title title
  Cluster synchronization in networks of identical oscillators with
  -function pulse coupling, @noop  journal journal Physical Review E volume 95, pages 022207 (year
  2017)NoStop
[Stroud et al.(2015)Stroud,
  Barahona, and Pereira]Stroud2015
  author author J Stroud, author M Barahona,   and author T Pereira, title Dynamics of cluster synchronisation in modular networks:
  Implications for structural and functional networks, in @noop 
  booktitle Applications of Chaos and Nonlinear Dynamics in
  Science and Engineering, Vol. volume 4, editor
  edited by editor S Banerjee  and editor L Rondoni (publisher Springer International Publishing, year 2015)  pp. pages 107–130NoStop
[Alderisio et al.(2017)Alderisio, Fiore, and di Bernardo]Alderisio2017
  author author F Alderisio, author G Fiore,   and author M di Bernardo,  title title Reconstructing the structure
  of directed and weighted networks of nonlinear oscillators, @noop 
  journal journal Physical Review E volume 95, pages 042302 (year
  2017)NoStop
[Whalen et al.(2015)Whalen,
  Brennan, Sauer, and Schiff]Whalen2015
  author author A J Whalen, author S N Brennan,
  author T D Sauer,  and  author S J Schiff, title title Observability and controllability of
  nonlinear networks: The role of symmetry, @noop  journal journal Physical Review X volume 5, pages 011005 (year
  2015)NoStop
[Shiroky and Gendelman(2016)]Shiroky2016
  author author I B Shiroky and author O V Gendelman, title title Discrete
  breathers in an array of self-excited oscillators: exact solutions and
  stability. @noop  journal journal
  Chaos volume 26, pages 103112
  (year 2016)NoStop
[Cho et al.(2017)Cho,
  Nishikawa, and Motter]Cho2017
  author author Y S Cho, author T Nishikawa,   and author A E Motter,  title title Stable chimeras and
  independently synchronizable clusters, @noop  journal journal Physical Review Letters volume to appear (year 2017)NoStop
[Kuramoto and Battogtokh(2002)]kuramoto-battogtokh-2002
  author author Y Kuramoto and author D Battogtokh, title title Coexistence
  of coherence and incoherence in nonlocally coupled phase oscillators,  @noop  journal journal Nonlinear
  Phenomena in Complex Systems volume 5, pages 380–385 (year 2002)NoStop
[Abrams and Strogatz(2004)]Abrams2004
  author author D Abrams and author S Strogatz, title title Chimera states
  for coupled oscillators, @noop  journal journal Physical Review Letters volume 93,  pages 1–4 (year 2004)NoStop
[Denker et al.(2002)Denker,
  Timme, Diesmann, Wolf, and  Geisel]Denker2002
  author author M Denker, author M Timme,
  author M Diesmann, author F Wolf,  and author
  T Geisel, title title Breaking synchrony by heterogeneity in complex
  networks, @noop  journal journal
  Physical Review Letters volume 94, pages 074103 (year 2002)NoStop
[Sun et al.(2009)Sun,
  Bollt, and Nishikawa]Sun2009
  author author J Sun, author E M Bollt,  and  author T Nishikawa, title title Master stability functions for coupled
  nearly identical dynamical systems, @noop  journal
  journal EPL (Europhysics Letters) volume
  85, pages 60011 (year 2009)NoStop
[Belykh et al.(2003)Belykh,
  Belykh, Nevidin, and Hasler]Belykh2003
  author author I Belykh, author V Belykh,
  author K Nevidin,  and author M Hasler, title
  title Persistent clusters in lattices of coupled
  nonidentical chaotic systems, @noop  journal
  journal Chaos: An Interdisciplinary Journal of Nonlinear
  Science volume 13, pages 165–178
  (year 2003)NoStop
[Flunkert et al.(2010)Flunkert, Yanchuk, Dahms, and  Schöll]Flunkert2010
  author author V Flunkert, author S Yanchuk,
  author T Dahms,  and author E Schöll, title
  title Synchronizing distant nodes: A universal
  classification of networks, @noop  journal journal Physical Review Letters volume 105,  pages 254101 (year 2010)NoStop
[Dahms et al.(2012)Dahms,
  Lehnert, and Schöll]Dahms2012
  author author T Dahms, author J Lehnert,   and author E Schöll,  title title Cluster and group
  synchronization in delay-coupled networks, @noop  journal journal Physical Review E volume 86, pages 016202 (year
  2012)NoStop
[Huddy and Sun(2016)]Huddy2016
  author author S R Huddy and author J Sun,  title title Master stability islands for
  amplitude death in networks of delay-coupled oscillators, @noop 
  journal journal Physical Review E volume 93, pages 052209 (year
  2016)NoStop
[Buono and Golubitsky(2001)]Buono2001
  author author P L Buono and author M Golubitsky, title title Models of
  central pattern generators for quadruped locomotion: I. Primary gaits,  @noop  journal journal Journal of
  Mathematical Biology volume 42, pages
  291–326 (year 2001)NoStop




