Synchrony in networks of coupled non-smooth dynamical systems: Extending the master stability function

The master stability function is a powerful tool for determining synchrony in high-dimensional networks of coupled limit cycle oscillators. In part, this approach relies on the analysis of a low-dimensional variational equation around a periodic orbit. For smooth dynamical systems, this orbit is not generically available in closed form. However, many models in physics, engineering and biology admit to non-smooth piece-wise linear caricatures, for which it is possible to construct periodic orbits without recourse to numerical evolution of trajectories. A classic example is the McKean model of an excitable system that has been extensively studied in the mathematical neuroscience community. Understandably, the master stability function cannot be immediately applied to networks of such non-smooth elements. Here, we show how to extend the master stability function to non-smooth planar piece-wise linear systems, and in the process demonstrate that considerable insight into network dynamics can be obtained. In illustration, we highlight an inverse period-doubling route to synchrony, under variation in coupling strength, in globally linearly coupled networks for which the node dynamics is poised near a homoclinic bifurcation. Moreover, for a star graph, we establish a mechanism for achieving so-called remote synchronisation (where the hub oscillator does not synchronise with the rest of the network), even when all the oscillators are identical. We contrast this with node dynamics close to a non-smooth Andronov–Hopf bifurcation and also a saddle node bifurcation of limit cycles, for which no such bifurcation of synchrony occurs.


Introduction
Real world networks such as those found in the brain, heart, World Wide Web, geoeconomic structures and ecologies show wonderfully rich emergent behaviour. The observed complex patterns of network activity reflect both the connectivity and the nonlinear dynamics of the network components. The new science of networks [1] has proven especially fruitful in probing the role of connectivity, yet has very little to say about the role of dynamics, and this has mainly focused on the properties of synchronisation in complex networks [2,3]. This is perhaps not too surprising given the challenge of understanding even low-dimensional dynamical systems. However, there has been an appreciation for some time in the applied sciences, particularly in engineering and biology, of the benefits of studying caricature systems built from piece-wise linear (pwl) and possibly discontinuous dynamical systems. An example from mechanical engineering is the Lazer-McKenna suspension bridge model which has been used to explain the large oscillations seen in the Tacoma Narrows bridge collapse [4], using a pwl restoring force model for bridge stays. Indeed there is a very strong history of pwl modelling across all of engineering [5], and particularly in electrical engineering [6], that has now begun to pervade other disciplines, including the social sciences, finance and biology [7]. In a neuroscience context, the McKean model [8] is a classic example. This replaces the cubic nullcline of the FitzHugh-Nagumo model [9] for action potential generation with a pwl function that preserves the original shape, allowing explicit calculations that could not otherwise be performed in the smooth case. At heart, pwl modelling allows analytical insight to be obtained for a non-linear model by breaking down the phase space into zones where trajectories obey simple linear dynamical systems, and patching these together across the boundaries between the zones. The approach can also handle discontinuous dynamical systems, such as those that arise naturally when modelling impacting mechanical oscillators [10], integrate-and-fire models of spiking neurons [11] and caricatures of cardiac oscillators [12] with both state and time-dependent switching. Although a beautifully simplistic modelling perspective, the necessary loss of smoothness precludes the use of many results from the standard toolkit of smooth dynamical systems, and one must be careful to correctly determine conditions for existence, uniqueness and stability of solutions [13].
Given the considerable challenge of understanding the dynamics of real world networks, it is of interest to develop a new mathematical framework for understanding networks of pwl dynamical elements. The expectation being that the existing success of such an approach at the single node level, in disciplines ranging from engineering to biology, will reap similar benefits for understanding network states. As a first step in this direction, we will focus on networks of identical units with linear coupling. Specifically, we will show how to extend the powerful machinery developed by Pecora and Carroll [14] to analyse synchrony in an oscillatory network of coupled pwl non-smooth dynamical systems. Namely, we augment the master stability function (MSF) approach to treat networks of non-smooth dynamical elements.
In Section 2, we introduce the class of pwl linear planar dynamical models that we treat throughout the rest of this paper. This is sufficiently broad that we can discuss three distinct oscillator models that exemplify (i) a saddle node bifurcation of limit cycles, (ii) a non-smooth Andronov-Hopf bifurcation and (iii) a homoclinic bifurcation. The extension of standard Floquet theory is developed in Section 3, providing a method to determine the stability of periodic orbits in a non-smooth setting. A similar approach is used in Section 4 to construct the appropriate MSF for studying networks of non-smooth dynamical elements. As an application of this new form of MSF, we consider global and star linearly coupled networks in Section 5 and contrast the stability properties of the synchronous state for each of the oscillator models (i)-(iii). Finally, in Section 6, we discuss natural extensions of the work in this paper.

Piece-wise linear models
As well as providing useful caricatures of systems in engineering and biology, non-smooth pwl models can also be viewed as the uniform limit of smooth non-linearities. As such, the global dynamics of smooth models can sometimes be approximated by pwl models and vice versa. For a very readable review of pwl dynamical systems, we refer the reader to the tutorial by Ponce [15]. Here, we will treat only planar systems, though higher dimensional systems may be studied with a similar approach. Our main interest will be in the construction of periodic orbits.
Consider the system described by z = (v, w) ∈ R 2 withż = F(z) and Here, A L,R ∈ R 2×2 and c L,R ∈ R 2 . We define an indicator function h as so that the switch in the system occurs when h(z) = 0. We shall refer to the line v = a as the switching manifold, call F L (F R ) the left (right) half-system of (2.1) and call the region of the plane where the left (right) zone. If an equilibrium exists in the left (right) zone, then its stability is determined by the eigenvalues of A L (A R ). These are easily calculated as λ ± (A L ) (λ ± (A R )), where We may then classify each half-system of (2.1) in the trace-determinant plane, for which there are five types of linear system: saddle, attracting node, repelling node, attracting focus and repelling focus. Piece-wise linear systems of the type (2.1) may support bifurcations that do not exist in smooth systems [16]. For example, if an equilibrium crosses a switching manifold as a system parameter is continuously varied, we expect the eigenvalues to change discontinuously because of the discontinuity in the Jacobian. In this case, the equilibrium may disappear or its stability may change, and this is often termed a discontinuous bifurcation. Many bifurcations caused by discontinuities have been previously studied, and are often referred to as "C" bifurcations (after the Russian wordšvejnye for "sewing") [17]. Examples of C-bifurcations include discontinuous saddle-node bifurcations (that have an analogue in smooth systems), grazing bifurcations, where a periodic orbit touches a switching manifold, and sliding bifurcations, where a trajectory moves for some time along the switching manifold. A recent example of a sliding homoclinic orbit in a pwl neural mass model can be found in [18].
In general, it is hard to find closed form solutions for periodic orbits in non-linear dynamical systems. The exception to this rule being pwl systems, in which case solutions of local linear systems may be glued together to find global trajectories. In essence, we solve the system in each of its linear regimes and demand continuity of solutions to construct orbits of the full non-linear flow. To see how to do this, we denote a trajectory on S L (S R ) by z L (z R ) and integrate (2.1) to obtain z L (t, t 0 ) = z(A L , c L ; t, t 0 ) (z R (t, t 0 ) = z(A R , c R ; t, t 0 )), for t > t 0 , where and I 2 is the 2 × 2 identity matrix. If A has real eigenvalues λ ± (calculated from (2.4)), such that Aq ± = λ ± q ± with q ± ∈ R 2 , then we may diagonalise and write G(A; t) in the computationally useful form G(A; t) = P e Λt P −1 , where Λ = diag(λ + , λ − ), P = [q + , q − ], and q ± = [(λ ± − a 22 )/a 21 , 1] T , where a ij denote the components of A with i, j = 1, 2. If A has complex eigenvalues ρ ± iω, then the associated complex eigenvector is q such that Aq = (ρ + iω)q, q ∈ C 2 . In this case, G(A; t) = e ρt P R ωt P −1 , where with ω = ω/a 12 and ρ = (ρ − a 11 )/a 12 . Note that ρ and ω may be written using the invariance of Tr and det as ρ = (a 11 + a 22 )/2, ω 2 = a 11 a 22 − a 12 a 21 − ρ 2 > 0. Imagine now a closed orbit built from connecting two trajectories, starting from initial data z(0) = (a, w(0)), in each of the half-systems according to for some Δ > Δ R > 0. A periodic orbit may be found by demanding that z be Δ-periodic. Thus, we must simultaneously solve the three equations: for the three unknowns (Δ, Δ R , w(0)). This may be done using Matlab (which readily implements matrix exponentials using expm and can solve systems of algebraic equations with fsolve). The functions (v(t), w(t)) are constructed from (2.5) and (2.8) with z(Δ R ) = z R (Δ R , 0). We shall refer to Δ R and Δ L = Δ − Δ R as the time-of-flight on S R and S L , respectively. We shall now consider three different pwl models, each of which supports oscillatory behaviour, and utilise the above analysis to construct periodic orbits.

The McKean model
Like the FitzHugh-Nagumo [9] model that of McKean [8] is a model for action-potential generation -though its voltage nullcline has an infinitely steep middle branch. The equations of motion take the forṁ v = −γv + μH(v − a) − w,ẇ = bv, (2.10) where H is the Heaviside step function. Thus, the model may be written in the form The vector field of the model is not defined if v = a, and we have some freedom on how to extend the vector field to this surface. The widely accepted way to resolve this freedom is to consider a set valued extension (see e.g. [19]) and writė where Σ = {z | v = a} so that R 2 = S L Σ S R and co(S) denotes the smallest closed convex set containing S. In this particular case, The extension (or convexification) of a discontinuous system (2.10) into a convex differential inclusion (2.12) is known as the Filippov convex method [20]. The function Figure 1. For a > 0, the fixed point at (v, w) = (0, 0) has eigenvalues λ ± = (−γ ± γ 2 − 4b)/2. Hence, for b, γ > 0, the fixed point is stable. For γ 2 − 4b > 0, the fixed point is globally attracting. When γ 2 − 4b < 0, the fixed point is a focus, and periodic behaviour may be possible. As a trajectory spirals in to the fixed point, it may hit the section Σ. At this point, there are two possibilities: (a) the trajectory exits Σ and enters into either S L or S R ; (b) the trajectory remains in Σ and slides.
Using the approach described in Section 2, we may construct periodic orbits, such as those illustrated in Figure 2. Interestingly, stable periodic orbits can co-exist with the stable fixed point, separated by an unstable periodic orbit. In the case that this is comprised of the union of an unstable sliding mode on Σ together with a trajectory on S L , we may construct the orbit by considering evolution in backward time. Under the exchange t → −t, the signs ofv in Figure 1 are reversed and sliding motion along Σ is now stable. Given thatẇ = −ba (remembering that we are evolving backward in time) along Σ, then trajectories will slide down (withv = 0) until they hit the point where κ(w) = 0, at which point they will be subject to dynamics withv < 0 and will spiral out in S L until they return and hit Σ. If this is below the point where w = −γa + μ, then the trajectory will slide down again, resulting in a periodic orbit. The time-of-flight in S L can be calculated Hence, the unstable periodic orbit is the union of z L (t, 0) and the vertical line joining (a, w L (Δ L )) and (a, −γa). Given that the model can naturally support a co-existing stable and unstable periodic orbit, we should expect a saddle node bifurcation of limit cycles under parameter variation. We shall not pursue the construction of bifurcation diagrams here, though it is useful to mention that the annihilation of periodic orbits is very easily achieved if we allow the slopes of the voltage nullcline to vary on either side of Σ and use the ratio of these as a bifurcation parameter. For a further exposition of the model, we refer the reader to [21] and for an alternative approach to the construction of periodic orbits see [22]. It is worth emphasising here that the advantage of studying the McKean model, as opposed to the FitzHugh-Nagumo model, is that it can be solved explicitly. Moreover, the generalisation of the model to describe an axon with the inclusion of a spatial diffusion term also admits to a level of analysis for the shape, speed and stability of waves that would be otherwise be hard to achieve [23].

The absolute model
In contrast to the McKean model, the absolute model has a continuous v-nullcline and the equations of motion are written in the forṁ 14) with a > 0, 0 < g < 1 and w 0 + (a − v 0 )/g < 0. Thus, the model may be written in the form (2.1) with with stability determined by the eigenvalues of A R : Hence, the fixed point for v > a is an unstable focus. In contrast if the fixed point lies , then stability is determined by the eigenvalues of A L : In this case, the fixed point is a stable focus. The absolute model can thus support a non-smooth Andronov-Hopf bifurcation that arises when an equilibrium encounters a switching manifold. This occurs when the equilibrium moves from v < a to v > a and crosses the switching manifold in such a way that its eigenvalues "jump" across the imaginary axis. For a more thorough discussion of this bifurcation, we refer the reader to [24]. An example of a stable periodic orbit, that can emerge from a non-smooth Andronov-Hopf bifurcation, is shown in Figure 3, together with the nullcline structure for the model. It is an interesting exercise to check, that in contrast to a smooth Andronov-Hopf bifurcation, emergent periodic orbits have a period that is independent of g and whose amplitude scales linearly with g.

A model with a homoclinic loop
Consider a scenario with a saddle for z ∈ S L and a repelling focus for z ∈ S R . In this case, it is possible to construct a homoclinic orbit that is tangential to the stable and unstable eigen directions of the saddle in S L , and also manages to enclose the repelling focus in S R . An example of such an orbit is shown in Figure 4 for the pwl model defined by wherev andẇ are chosen to be continuous at v = a (by fixing a = 0). Here, we see that Tr A μ = τ μ and det A μ = δ μ for μ ∈ {L, R}. If we denote λ ± (A L ) = λ ± and λ ± (A R ) = α ± iβ, with λ − < 0 < λ + and β > 0, then a (non-degenerate) homoclinic exists if and only if Λ = 0 [25], where Moreover, the existence of a limit cycle or a homoclinic loop requires τ L τ R 6 0 [25]. These conditions are merely the result of demanding that straight line trajectories in the (v, w) plane leaving and terminating on the saddle in S L along its unstable and stable manifolds (which because of linearity are tangent to the eigendirections q + and q − , respectively) can form a continuous global trajectory after matching to curved trajectories in S R at the switching manifold, where v = a. It is also possible to construct homoclinic orbits when the fixed point in S R is a centre (α = 0), and see [26] for a further discussion. In this case, it is easy to see that Λ = 0 for the choice Tr A R = 0 = Tr A L and det A R = 1 = − det A L , although in this special case there are also an infinite number of periodic orbits residing within the homoclinic loop.

Floquet theory
One may naturally ask how the determination of stability of periodic orbits has been performed in the above. For the case the vector field is not differentiable or even continuous one must be careful not to invoke standard Floquet theory without further thought. For a smooth planar dynamical system, one of the Floquet exponents is equal to zero, corresponding to perturbations along the periodic orbit, and the other, σ smooth , is given by Tr DF(z(s))ds, (3.1) where z(t) denotes the Δ-periodic orbit, and DF is the Jacobian of the vector field. However, much like the construction of continuous periodic orbits can be readily achieved for pwl planar systems, so too does the Floquet theory simplify. In Appendix A, we show that the Floquet exponent for the pwl models considered here takes the form where T R = Δ R and T L = Δ are the switching times defined by h(z(T μ )) = 0. Here, T ± μ = lim x→0 + T μ ± x. We note that ifv is continuous across the switching manifold, as is the case for the absolute model, then σ reduces to σ smooth as expected. Periodic orbits are stable if σ < 0 and unstable if σ > 0. For a more general discussion of the stability of periodic orbits in planar Filippov system with exactly one switching line (not restricted to pwl models), see [27]. Interestingly, the result (3.2) was probably first discovered by Leonid Shilnikov during his Master's thesis sometime in the 1950s and subsequently discussed in papers with Neimark 1 [29,30].
For the McKean model, we have that Given that the green orbit in the left plot of Figure 2 has a nearly continuous gradient the formula above shows that it is stable, since σ −γ < 0. The stability of the other orbits is found by the numerical calculation of the full formula for σ. For the absolute model, we have that and for the homoclinic model (3.5)

Extending the master stability function
Synchrony is an important network state due to its relevance in both real world applications and the foundations of non-linear dynamics [31]. The MSF technique has proven a very powerful tool for understanding synchrony in coupled systems of identically coupled oscillators [14] in terms of the eigenstructure of the network connectivity matrix. However, as it assumes that the underlying oscillation (or even chaotic time-series) is generated by a smooth dynamical system, we cannot directly employ the MSF approach for the models considered here. The reason being that, as for the Floquet theory described above in Section 3, one must be careful when propagating perturbations through switching manifolds. However, we may easily incorporate the Floquet treatment of non-smooth systems within the MSF framework for oscillators. As such, it is first convenient to review the MSF approach for smooth systems.
To introduce the MSF formalism, it is convenient to consider N nodes (oscillators) and let z i = (v i , w i ) ∈ R 2 be the two-dimensional vector of dynamical variables of the ith node with isolated (uncoupled) dynamicsż i = F(z i ), with i = 1, . . . , N. The output for each node is described by a vector function H ∈ R 2 . For a given connectivity matrix with components w ij and a global coupling strength , the network dynamics of N coupled identical systems, to which the MSF formalism applies, is specified by Here, the matrix G ∈ R N×N with entries G ij has the graph-Laplacian structure G ij = −w ij + δ ij k w ik . The N − 1 constraints z 1 (t) = z 2 (t) = . . where ξ l is the lth (right) eigenmode associated with the eigenvalue λ l of G (and DF(s) and DH(s) are independent of the block label l). Since j G 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 G, the eigenvalues λ l may be complex, which brings us to consideration of the system For given s(t), the MSF is defined as the function which maps the complex number α to the greatest Floquet exponent of (4.3). The synchronous state of the system of coupled oscillators is stable if the MSF is negative at α l = λ l , where λ l ranges over the eigenvalues of the matrix G (excluding λ 1 = 0). For a further discussion about the use of the MSF formalism in the analysis of synchronisation of oscillators on complex networks, we refer the reader to [3,32], and for the use of this formalism in spiking neural networks of non-linear integrate-and-fire type see [11,33]. This approach has recently been extended to cover the case of cluster states by making extensive use of tools from computational group theory to determine admissible patterns of synchrony [34] in unweighted networks. We now restrict our attention to pwl node dynamics with linear coupling in the vcomponent so that H(z) = (v, 0), which is often referred to as diffusive coupling [35]. In this case, both Jacobians, DF(s) and DH(s), are piece-wise constant matrices. Thus, we may solve (4.3) using matrix exponentials, since time-ordering is easily handled, being careful to treat how perturbations evolve through the switching manifold. Using the techniques developed for constructing the Floquet exponent, as in Appendix A, we find that after one period ξ(Δ) = Γ (l)ξ(0), where and Here, the matrices G μ (l) describe the evolution of a system linearised around a periodic orbit between switching events, whilst the saltation matrices K μ describe how to map perturbed trajectories through a switching manifold. We note that for a continuous vector field the matrices K μ reduce to I 2 as expected. Thus, the synchronous solution will be stable if the periodic orbit of a single node is stable (in the absence of coupling) and the eigenvalues of Γ (l) for l = 2, . . . , N lie within the unit disc. It is worth noting that the MSF formalism can be applied to chaotic systems where, instead of Floquet exponents, one would calculate Liapunov exponents. In this case, one must also be careful in treating dynamics at switching events, and for a further discussion of the calculation of Liapunov exponents in discontinuous dynamical systems see [11,36].

Global and star graph connectivities
Although the MSF approach (and its extension to non-smooth systems) is valid for an arbitrary network, it is informative to consider an application to a globally coupled network with > 0. In this case, the eigenvalues of the graph-Laplacian are easily calculated since for the choice w ij = N −1 we have that G ij = −N −1 + δ ij , and the non-zero (N − 1 degenerate) eigenvalue is +1. In this case, we need only consider the eigenvalues of M 1 ( ) = K L G(A L ( ); Δ L )K R G(A R ( ); Δ R ), where G(A; t) is given in (2.6) and We note that the condition for an instability is independent of the system size, since neither A μ ( ) nor K μ depend on N.
The condition for eigenvalues to cross the unit circle along the real axis is det[M 1 ( ) ± I 2 ] = 0, and off of the real axis we have the condition det M 1 ( ) = 1. For the McKean and absolute models discussed earlier, and specifically the models with parameters as in Figures 2 and 3, we find that the synchronous state is stable for weak coupling and that this stability persists with increasing . However, for the homoclinic node model, the synchronous state is unstable for weak and can restabilise with increasing when det[M 1 ( ) + I 2 ] = 0, namely via an inverse period-doubling bifurcation. For the parameters of Figure 4, this occurs at = pd 2.36, and is found to be in excellent agreement with direct numerical network simulations. Below the point of instability where 0 < < pd simulations further show that the network typically exhibits quasi-periodic behaviour. This is expected since the associated network eigenvectors have components that sum to zero, describing perturbations that act to push the oscillations apart pairwise. To quantify this, we introduce the mean field variables (E(v), E(w)), where and trace the mean field signal for a typical large network simulation in Figure 5. Here, it can be seen that the mean field trajectory does not settle down, and moreover is repelled away from the synchronous solution (also plotted). For just less than pd , we see a period doubled orbit as expected. For large values of , it is also possible that a tangent bifurcation can occur, where det[M 1 ( ) − I 2 ] = 0, leading to the formation of a cluster state. For the parameters of Figure 4, this is found to occur at t 4.45, with an instability occurring as is increased through t . We note that for large N the system may support a splay state with a constant mean field. In this case, one may use the techniques in [37] to determine the period and stability of this network state, though we shall not pursue this here.
It is also informative to consider a hub and spoke model of N nodes, with a central oscillator connected to N − 1 other leaf nodes (which are not connected to each other), describing a star network. These arise very naturally in computer network topologies where one central computer acts as a conduit to transmit messages (providing a common connection point for all nodes through a hub). This star graph architecture is described by a connectivity matrix that we write in the form Figure 5. Simulation of a globally linearly coupled network of homoclinic node models with = 0.3 and N = 100. Other parameters as in Figure 4. The solid green line shows the unstable synchronous orbit, whilst the thin black line shows the mean field trajectory (E(v), E(w)).
If we make the convenient choice K = N − 1, then the eigenvalues of the graph-Laplacian are 2, 0 and 1 (with N − 2 degeneracy), then any instabilities will be independent of the system size. Thus, the spectrum of eigenvalues that determines stability is the same as that for the globally coupled system together with the eigenvalues of the matrix M 2 ( ) = K L G(A L (2 ); Δ L )K R G(A R (2 ); Δ R ). Thus, if synchrony is stable in a globally coupled network, then it will also be stable in a star network (since M 1 ( ) and M 2 ( ) are equivalent under a simple re-scaling of ). If we consider a network of homoclinic node models then the critical value of determined from M 2 ( ) for a period doubling instability is precisely pd /2. This would not be met before the inverse period doubling bifurcation at pd (defined by det[M 1 ( ) − I 2 ] = 0) when decreasing . Since the associated network eigenvector is (0, u 2 , . . . , u N ) with N i=2 u i = 0, the instability at = pd will excite a mode in which the leaf nodes desynchronise. However, a tangent bifurcation, where det[M 2 ( ) − I 2 ] = 0 is also possible. In this case, the associated network eigenvector is (−1, 1, . . . , 1) and we would expect the network to break into a cluster state in which the leaf nodes are synchronised with each other, but not with the hub. This is sometimes referred to as remote synchronisation [38]. For the parameters of Figure 4, a cluster bifurcation is found to occur at t 2.23, with an instability occurring as is increased through t . Thus, the synchronous state is always unstable, being unstable to both a tangent and period-doubling bifurcation for t < < pd and unstable only to a period doubling instability for < t and only unstable to a tangent bifurcation for > pd . An example of an emergent network state in each of these last two parameter regimes is shown for a star graph with six nodes in Figure 6. Given that, in this example, t and pd are so close we would expect to see the network state that is most close to synchrony in the parameter window ∈ ( t , pd ), and this is confirmed in simulations. It is interesting to note that although remote synchronisation, as illustrated in Figure 6 (right), has previously been seen in star networks of Stuart-Landau oscillators [39], this has required a frequency mismatch between the hub oscillator and the leaf ones. For the homoclinic node model considered here, no such frequency mismatch is required. Figure 6. Simulation of a linearly coupled star network of six homoclinic node models. Parameters as in Figure 4. The solid green line shows the unstable synchronous orbit. Blue lines are trajectories of the hub node, and other lines those of the leaf nodes (colour coded according to the star graph inset). Left: = 1, where the system is only unstable to a period doubling instability. Right: = 4, where the systems is only unstable to a tangent bifurcation. Here, we see an example of a cluster state exhibiting remote synchronisation.

Discussion
In this paper, we have extended the MSF approach for analysing synchrony to treat structured networks of identical planar non-smooth oscillators. To illustrate the application of the theory, we have considered global linearly coupled networks and used this to highlight the importance that local node dynamics can have on network states. In particular, we have shown that when the local dynamics is poised near a homoclinic bifurcation that synchrony can be unstable for weak to moderate values of the coupling strength for a globally coupled network.
As well as local node dynamics, it is now increasingly recognised that there is a major role to be played by motifs in the synchronisation of complex networks [40]. Motifs are specific patterns of network connections (either direct or indirect) that occur at a rate significantly higher than in randomised versions of the network [41]. In the framework presented here, the critical value of coupling strength to ensure the stabilisation of the synchronous state is determined in terms of the eigenvalues of a finite set of 2×2 matrices, each dependent on the eigenvalues of the graph-Laplacian for the network motif. We have utilised this approach for a star graph to establish that remote synchronisation can also arise when the choice of local node dynamics is that of an oscillator close to a homoclinic bifurcation.
The matrix formalism that we have employed is equally well suited to studying higher dimensional pwl node dynamics, and represents one natural extension of the work in this paper. We now reflect briefly on other extensions.
The synchronous state is a very special network solution, especially when compared to more exotic ones, such as phase-locked clusters, chimeras (mixing synchronous and asynchronous states), chaos and death (where the dynamics of some nodes is silenced), as recently reviewed in [42] from a neuroscience perspective. Although the tools of weakly coupled oscillator theory have shed some light on these states, as for example described in the book by Hoppensteadt and Izhikevich, [43], we currently lack a general theory of network dynamics valid for strong coupling. However, for the explicit forms of pwl models that we have discussed in this paper, we further expect to be able to develop an explicit theory valid for strong coupling. We have already realised this for synchronous states (in networks built from identical units), and we advocate for a push to develop this methodology further to treat more exotic network states, including those expected in the presence of noise [44]. Moreover, it seems very fruitful to classify emergent dynamics based upon the group of structural symmetries of the network [45], as has recently been done in [34] to determine cluster states.
Importantly, the framework that we have presented here, and illustrated for linear coupling, is particularly well suited to studying the type of event driven synaptic interactions that arise naturally in models of biological neural networks. Indeed the techniques for constructing the relevant saltation matrices for handling synaptic dynamics with pulsatile forcing as well as discontinuous reset in non-linear integrate-and-fire systems have recently been described by Ladenbauer et al. [33]. The combination of this with pwl caricatures of Izhikevich style integrate-and-fire models [46], such as described in [47], opens the door for a very thorough examination of network states in biologically realistic spiking networks.
The above are all topics of current interest and will be reported upon elsewhere.
Here, we have introduced the notation z(T ± μ ) = lim 0 z(T μ ± ), to make sure that derivatives are well defined. Using the fact that h(z(T μ )) = 0 = h( z( T μ )), we obtain ∇ z h(z(T μ )) · δz(T μ ) + z (T − μ )δT μ = 0. (A 4) Using the result that ∇ z h(z) = (∂ v , ∂ w )(v − a) = (1, 0), the above can be re-arranged to give the perturbation in the switching time in terms of the difference between the perturbed and unperturbed trajectories as (A 5) We now construct the deviation between the two trajectories at the perturbed switching time as If δT μ > 0, then the unperturbed trajectory will already have switched, in which case the two trajectories are on either side of the switching manifold. A similar argument holds for δT μ < 0. Thus, we may write Combining (A 5) and (A 7) gives We may write the above in matrix form as ) or equivalently as δz(T μ + δT μ ) = K(T μ )δz(T μ ) with From (A 1), between switching events at Δ R and Δ, the perturbations evolve according to G μ (t)δz 0 , where G μ (t) ≡ G(A μ ; t) = e A μ t and δz 0 is the perturbation at the switching time.
Thus, after one period of oscillation we may put this all together to obtain δz(Δ) = Mδz(0), The periodic orbit will be stable if the eigenvalues of M lie within the unit disc. Note that one of the Floquet multipliers is equal to one, corresponding to perturbations along the periodic orbit. Let us denote the other eigenvalue by e σΔ and use the result that det M = e σΔ × 1. Hence, Using the fact that det e At = e Tr A t , we find where T R = Δ R , T L = Δ and Δ L = Δ − Δ R .