Networks of piecewise linear neural mass models

Neural mass models are ubiquitous in large scale brain modelling. At the node level they are written in terms of a set of ODEs with a nonlinearity that is typically a sigmoidal shape. Using structural data from brain atlases they may be connected into a network to investigate the emergence of functional dynamic states, such as synchrony. With the simple restriction of the classic sigmoidal nonlinearity to a piecewise linear caricature we show that the famous Wilson-Cowan neural mass model can be analysed at both the node and network level. The construction of periodic orbits at the node level is achieved by patching together matrix exponential solutions, and stability is determined using Floquet theory. For networks with interactions described by circulant matrices, we show that the stability of the synchronous state can be determined in terms of a low-dimensional Floquet problem parameterised by the eigenvalues of the interaction matrix. This network Floquet problem is readily solved using linear algebra, to predict the onset of spatio-temporal network patterns arising from a synchronous instability. We consider the case of a discontinuous choice for the node nonlinearity, namely the replacement of the sigmoid by a Heaviside nonlinearity. This gives rise to a continuous-time switching network. At the node level this allows for the existence of unstable sliding periodic orbits, which we construct. The stability of a periodic orbit is now treated with a modification of Floquet theory to treat the evolution of small perturbations through switching manifolds via saltation matrices. At the network level the stability analysis of the synchronous state is considerably more challenging. Here we report on the use of ideas originally developed for the study of Glass networks to treat the stability of periodic network states in neural mass models with discontinuous interactions.


Introduction
The Wilson-Cowan model [47,48] is one of the most well-known neural mass models for modelling the activity of cortex, and for a historical perspective see [11]. Neural mass models generate brain rhythms using the notion of population firing rates, aiming to side-step the need for large-scale simulations of more realistic networks of spiking neurons. Although they are not derived from detailed conductance-based models, they can be motivated by a number of phenomenological arguments [9], and typically take the form of systems of non-linear ordinary differential equations (ODEs). The Wilson-Cowan neural mass model describes the dynamics of two interacting populations of neurons, one of which is excitatory and the other inhibitory. Interactions are mediated between the populations with the use of a non-linear sigmoidal firing rate function. In its most simple incarnation, it consists of two non-linear ODEs, and as such has been widely studied using techniques from phase-plane analysis and numerical bifurcation theory. At the network level, the model can either be posed on a graph or a continuous space, and since the 1970s there has been a large amount of attention devoted to the analysis of these models and their application in neuroscience [14]. Recent examples of their use include reconciling information from anatomical and functional data [49], understanding phase-amplitude coupling (whereby the amplitude of a higher frequency brain rhythm is modulated by the phase of lower frequency activity) [38], modelling epilepsy [36], and understanding the emergence of cortical resonant frequencies [32]. Indeed, there are many variants of the Wilson-Cowan neural mass model now in use for interpreting neuroimaging data [46], including those of Zetterberg et al. [50], Jansen and Rit [26], and Liley et al. [34]. Neural mass models are a key component of the Virtual Brain project that aims to deliver the first simulation of the human brain based on individual large-scale connectivity [41]. Such large-scale brain network models are especially relevant to understanding resting state networks [5], whereby different regions of the brain's sensorimotor system oscillate slowly and synchronously in the absence of any explicit task.
However, at heart it is well to note that from a mathematical modelling perspective, all neural mass models to date are essentially low-dimensional coupled ODEs with a sigmoidal firing rate non-linearity, exemplified by the Wilson-Cowan model. Using extensions of the techniques originally developed by Amari [2], the continuum or neural field [9] Wilson-Cowan model has been analysed when the choice of this firing rate non-linearity is a Heaviside function. This has been possible because of a smoothing of the firing rate with a spatial kernel representing anatomical connectivity. However, when posed on a graph, representing a network of interacting neural populations, no such smoothing arises. Surprisingly, there are hardly any mathematical results for such networks, as opposed to their continuum counterparts for which there are now a plethora ranging from the properties of localised states through to travelling waves, as reviewed in [8]. This discrepancy is really a reflection of the fact that there are many more techniques for studying smooth dynamical systems as opposed to non-smooth. However, the body of mathematical work in this area is rapidly growing, driven in part by its importance to engineering [7,15]. Given their relevance to large-scale brain dynamics, it is highly desirable to develop mathematical techniques for the analysis of Wilson-Cowan style neural mass models at the network level. Here, we advocate for the replacement of smooth sigmoidal non-linearities in neural mass models by more tractable functions, including piecewise linear (PWL) and piecewise constant functions. A PWL continuous choice has been used in several previous studies, including those of Hansel and Sompolinsky [24], and Kilpatrick and Bressloff [29], whilst the discontinuous Heaviside (piecewise constant) choice has proven especially popular since the seminal work of Amari [2]. In these instances, this has facilitated the construction of certain types of localised states in continuum neural field models. However, in a discrete neural network context, there is a major mathematical difference in the analysis of network states for the case of continuous versus discontinuous firing rates. As well as introducing a simple methodology to treat the construction of periodic orbits in idealised Wilson-Cowan networks, this is one of the major topics we wish to address in this paper.
First, in Section 2, we introduce the model for an isolated Wilson-Cowan node with a PWL firing rate. The description of dynamical states with reference to switching manifolds becomes very useful. We show how matrix exponentials can be used to patch together a periodic orbit, and that Floquet theory simplifies considerably to yield explicit formulas for determining solution stability. Next, in Section 3, we consider a network of PWL Wilson-Cowan nodes, with nodes arranged along a ring with distance-dependent interactions. This particular choice of coupling guarantees the existence of the synchronous state. We then develop a linear stability analysis of this state and show that this leads to a tractable variational problem of a very similar type to that for the single node, albeit now parameterised by the eigenvalues of the connectivity matrix. We use this to determine instabilities that can lead to the formation of spatio-temporal network patterns. Next, in Section 4, we consider the case that the firing rate is a Heaviside function, for which the techniques developed for studying PWL systems break down. Once again periodic orbits can be constructed using matrix exponentials, although standard Floquet theory must be now augmented to cope with the evolution of linearised perturbations through the switching manifolds. This is most readily achieved with the use of saltation matrices that have commonly been used for the study of non-smooth mechanical systems [33]. However, at the network level, the stability of the synchronous state is much harder to determine than for the continuous model. Here, we show that ideas from the study of Glass networks developed by Edwards [17] are particularly useful, and that stability is strongly influenced by the temporal order in which network components cross-switching manifolds, and that this in turn is determined by the choice of initial perturbation. Finally, in Section 6, we conclude with an overview of the new results about synchrony in networks of neural mass models, and discuss the natural extension of this work to treat non-synchronous states.

The Wilson-Cowan model and a piecewise linear reduction
For their activity-based neural mass model, Wilson and Cowan [47,48] distinguished between excitatory and inhibitory sub-populations. This seminal (space-clamped) model can be written succinctly in terms of the pair of coupled differential equations: Here, u = u(t) is a temporal coarse-grained variable describing the proportion of excitatory cells firing per unit time at the instant t. Similarly, the variable v represents the activity of an inhibitory population of cells. The constants w αβ , α, β ∈ {u, v}, describe the weight of all synapses from the αth population to cells of the βth population, and τ is a relative time-scale. The non-linear function F describes the expected proportion of neurons in population α receiving at least threshold excitation per unit time, and is often taken to have a sigmoidal form. Here, the terms I α represent external inputs (that could be time varying). For a historical perspective on the Wilson-Cowan model see [14], and for a more recent reflection by Cowan see [12]. To reduce the model to a mathematically tractable form, we consider the choice of a PWL firing rate function given by For appropriate choices of parameters the Wilson-Cowan model, with the firing rate given by (2.2), can support stable oscillations. An example is shown in Figure 1, where we also plot the four switching manifolds defined by the condition that arguments to the function F in (2.1) take on the values zero and . Away from the switching manifolds, the dynamics governing the evolution of trajectories is linear, and may be constructed using matrix exponentials. To simplify further analysis, it is first convenient to introduce new variables )/|W |, as well as the matrices With these choices, (2.1) transforms to d dt In the representation (2.4), we see that the four switching manifolds are simply defined by U = 0, U = , V = 0, and V = . The periodic orbit shown in Figure 1 (encircling an unstable fixed point) crosses each of these manifolds twice, so that the periodic trajectory is naturally decomposed into eight separate pieces. On each piece, we shall denote the time-of-flight for a trajectory to travel from one switching manifold to another by Δ i , i = 1, . . . , 8, so that the period of the orbit is given by Δ = 8 i=1 Δ i . As an explicit example of how to construct a trajectory between two switching manifolds, consider the region where 0 6 U 6 and V < 0. In this case, the solution of (2.4) is given by It is a simple matter to write down the trajectories in each of the remaining regions of phase space visited by a periodic orbit. We may then use these matrix exponential formulas to patch together solutions, setting the origin of time in each region such that initial data in one region comes from final data from a trajectory in a neighbouring region. We shall denote the periodic orbit by (U, V ) such that (U(t), V (t)) = (U(t + Δ), V (t + Δ)).
If we consider initial data with (U(0), V (0)) = (U 0 , 0) then the eight times-of-flight and the unknown U 0 are determined self-consistently by the nine equations V (Δ 1 ) = , U(Δ 2 ) = , The numerical solution of this non-linear algebraic system of equations can be used to construct periodic orbits such as the one shown in Figure 1. Note that the construction of periodic orbits that do not cross all of the switching manifolds can similarly be performed (requiring the simultaneous solution of fewer equations). To determine stability we can turn directly to Floquet theory for planar systems which tells us that the non-trivial Floquet exponent is given by where D(s) denotes the Jacobian of the system evaluated along the periodic orbit. In general this is a hard quantity to evaluate for systems where the periodic orbit is not available in closed form. However, for the PWL Wilson-Cowan model the Jacobian is piecewise constant and we have that where Thus a periodic orbit is stable if σ < 0. In Figure 2, we present a plot of σ as a function of τ to show that the periodic solution in Figure 1 is stable. Given the above method to construct and determine the stability of a periodic orbit, we next show how to extend this approach to treat synchronous solutions in networks of Wilson-Cowan oscillators. Parameters as in Figure 1. Periodic orbits emerge via a supercritical Hopf bifurcation as τ increases through τ Hopf = (w vv + )/(w uu − ) ∼ 0.3. We see that the branch of periodic orbits shown is stable, with stability decreasing to zero as the solution is lost with increasing τ. This loss of existence occurs because of a grazing bifurcation (coincident with a saddle-node bifurcation of periodic orbits) at τgraze ∼ 0.6 whereby part of the trajectory develops a point of inflection on the switching manifold v = (Iu + w uu u)/w vu (red solid line in Figure 1), such that beyond bifurcation the trajectory does not cross the switching manifold and instead is attracted to the stable fixed point at (u, v) = (0, 0).

A piecewise linear Wilson-Cowan network
The study of coupled oscillator networks in biology, physics, and engineering is now commonplace. Two particularly well-known tools for studying patterns of phase-locked states and their instabilities are the theory of weakly coupled oscillators [30], and the master stability function (MSF) [39]. The reduction of a coupled limit cycle network to a phase oscillator network has proven very useful for gaining insight into phenomena ranging from the synchronisation of flashing fireflies [19] to behaviours in social networks [4], and for a recent review see [16]. However, there is an obvious limitation to such an approach, namely the restriction to weak interaction (and near identical oscillators). The MSF approach (for identical oscillators) does not require any such restriction on coupling strength, and can be used to determine the stability of the synchronous state in terms of the eigenstructure of the network connectivity matrix. However, the numerical evolution of a system of dynamical equations, arising from a Floquet variational problem, must be performed. Importantly the MSF approach 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 [40,43]. Here, we favour the MSF approach and show how it simplifies considerably for a PWL choice of firing rate function. This allows us to improve upon previous mathematical studies of Wilson-Cowan networks, such as those by Campbell and Wang [6] (who treated networks with nearest neighbour coupling and established the condition for synchrony), Ueta and Chen [45] (who performed a numerical bifurcation analysis for small networks), and Ahmadizadeh et al. [1] (who used perturbation techniques and numerics to study synchrony in networks with diffusive coupling).
We now consider a network of Wilson-Cowan nodes given by These row-sum constraints are natural for networks arranged on a ring, and guarantee the existence of a synchronous orbit ( is given by the solution of (2.1).
It is now convenient to introduce a vector notation with Here, the symbol ⊗ denotes the usual tensor product for matrices, and 1 N is an N-dimensional vector with all entries equal to unity. This means that the switching manifolds can be succinctly described by Y i = 0 and Y i = , and the dynamics takes the form where J is given by (2.3) and I N is the N ×N identity matrix. If we denote the synchronous solution by ) and consider small perturbations such that Y = Y + δY , then these evolve according to where DF(Y ) is the Jacobian of F evaluated along the periodic orbit. Given the constraints on the matrices W αβ , with α, β ∈ {u, v} it is natural to take these to be circulant matrices with W αβ ij = W αβ |i−j| . In this case, the normalised eigenvectors of If we introduce the matrix of eigenvectors P = [e 0 e 1 . . . e N−1 ], then we have that (1), . . . , ν αβ (N − 1)), and Moreover, it is easy to establish that in the above notation (P ⊗ If we now consider perturbations of the form δZ = (P ⊗ I 2 ) −1 δY , then from (3.6), we find that the linearised dynamics is described by the system where D ∈ R 2×2 is the Jacobian of (F(U), F(V )), and is a piecewise constant matrix that is only non-zero if 0 < U(t) < or 0 < V (t) < . In the former case [DF] 11 = −1 with all other entries zero, and in the latter case [DF] 22 = −1 with all other entries zero. We see that (3.10) has a block structure where the dynamics in each of N 2 × 2 blocks is given by Thus, comparing to (2.4), we see that the variational equation for the network is identical to that for a single Wilson-Cowan unit with W replaced by Λ(p). We note that for p = 0 the variational problem is identical to that for an isolated node since Λ(0) = W (using ν αβ (0) = N−1 μ=0 W αβ μ = w αβ ). Thus, to determine the stability of the synchronous state, we only have to consider a set of N two-dimensional variational problems. Exploiting the fact that between switching manifolds the variational problem defined by (3.11) is time-independent, we may construct a solution in a piecewise fashion from matrix exponentials and write ξ(t) = exp[(A(p) + Λ(p)JD)t]ξ(0). We may then build up a perturbed trajectory over one period of oscillation in the form where (3.13) Thus if a periodic orbit of an isolated Wilson-Cowan node is stable, then the synchronous network solution will be stable provided all the eigenvalues of Γ (p), for p = 0, . . . , N − 1, lie in the unit disc (excluding the one that arises from time-translation invariance, with a value +1). For a fixed value of p, one of three bifurcations is possible, namely a tangent instability defined by det(Γ (p) − I 2 ) = 0, a period-doubling instability defined by det(Γ (p) + I 2 ) = 0, and a Neimark-Sacker bifurcation defined by det Γ (p) = 1. If there is a p = p c such that one of these instabilities occurs, then the excited network state will correspond to the eigenvector Re e pc .

Example: a ring network
By way of illustration of the above theory, let us consider a network of Wilson-Cowan nodes arranged on a ring with an odd number of nodes. Introducing a distance between nodes indexed by i and j as dist(i, j) = min(|i − j|, N − |i − j|), we can define a set of exponentially decaying connectivity matrices according to (3.14) Thus, we have a set of four circulant matrices parametrised by the four spatial scales σ αβ that respect the row-sum constraints N j=1 W αβ ij = w αβ . In Figure 3, we show a plot of the eigenvalues of Γ (p) for p = 0, . . . , N − 1 for two different parameter choices. In one case, all of the eigenvalues (excluding the one arising from time-translation invariance) lie within the unit disc, whilst in the other one leaves the unit disc along the negative real axis. This latter scenario predicts an instability of the synchronous state, and is consistent with direct numerical simulations. Moreover, by studying the spectrum under parameter variation, we can find the value of p = p c which goes unstable first. In Figure 4, we show time courses (obtained by direct numerical simulation) for the components u i (t) of the emergent network state just beyond the point of instability, as well as a plot of the real part of the spatial eigenvector e pc . We see that the spatial pattern of the network state is well-predicted by e pc , suggesting that the bifurcation is supercritical.

The Heaviside world
In a recent paper, Harris and Ermentrout [25] considered a single Wilson-Cowan population with a Heaviside non-linearity, where the firing rate in (2.1) takes the form The choice of a Heaviside firing rate has been very popular in mathematical neuroscience ever since the seminal work of Amari (for neural field models), as nicely exemplified by his recent article on the 'Heaviside World' [3]. A case in point is the work of Laing and Chow [31] for understanding binocular rivalry. They considered a neural mass network model with recurrent excitation, cross-inhibition, adaptation, and synaptic depression and showed that the use of a Heaviside non-linearity allowed the explicit calculation of the dominance durations of perceptions. A more recent use of the Heaviside firing rate has been by McCleney and Kilpatrick [35] for neural activity models with spike rate adaptation to understand the dynamics of up-down states. Using techniques from Filippov systems and differential inclusions, Harris and Ermentrout made a study of periodic orbits for a Heaviside firing rate using a boundary value problem approach. Here, we show that we can recover their results using the matrix exponential approach of Section 2. Moreover, we also extend their work on a single node by showing how to determine the stability of periodic orbits using a non-smooth version of Floquet theory. In the representation (2.4), with F = H, we see that the there are two switching manifolds defined by U = 0 and V = 0. If we introduce the indicator functions h 1 (U, V ) = U and h 2 (U, V ) = V , then we can define these manifolds (lines in this case) as then the U-nullclines are given by and the V -nullclines are given by An example set of nullclines is shown in Figure 5.
To discuss fixed points and their stability it is first necessary to complete the description of the dynamics on the switching manifolds. We do this using the Filippov convex method [20] and extend our discontinuous system into a convex differential inclusion. The Filippov extension of (2.4) is then d dt . (4.7) As illustrated in Figure 5, it is possible for two nullclines to intersect and create a fixed point (U ss , V ss ). In the example shown, this occurs for U < 0 and V < 0, so that (U ss , V ss ) = (I u , I v ). Linear stability analysis shows that this is a stable node (with eigenvalues of A, namely −1 and −1/τ). Moreover, this system also supports pseudo equilibria, where either a nullcline touches a switching manifold, or two switching manifolds intersect. A thorough exploration of the pseudo equilibria of (2.1) can be found in [25]. Here, we shall simply focus on the pseudo equilibrium at (U ss , V ss ) = (0, 0), and characterise its stability by considering trajectories around this point. In fact given the PWL nature of the dynamics, it is sensible to consider the construction of periodic orbits, and determine the stability of the pseudo equilibrium in terms of the stability of encircling small amplitude orbits.

Periodic orbits and their stability
A non-sliding periodic orbit around (0, 0) can be constructed in terms of the times-of-flight in each region D αβ . If we denote these four times by the symbols Δ αβ , then the period of the orbit is given by Δ = Δ ++ + Δ −+ + Δ −− + Δ +− . We may then use a matrix exponential solution to patch together solutions, setting the origin of time in each region such that initial data in one region comes from final data from a trajectory in a neighbouring region. We shall denote the periodic orbit by (U, V ) such that (U(t), V (t)) = (U(t+Δ), V (t+Δ)). To indicate which region we are considering, we shall simply add αβ subscripts to the formula in (4.8).
In this way, a periodic orbit that visits all four regions in turn can be parameterised by the five unknowns U ++ (0)  Here, DF(Y (t)) is the piecewise constant matrix described after (3.10). Consider, for example, the time-of-flight, t 1 ( ), between U = and U = 0. For small , we may estimate t 1 ( ) using the result that U(t) U(t 0 ) +U t=t0 (t − t 0 ), giving t 1 ( ) = − /U t=Δ++ . The corresponding change in state across this small time interval can be obtained by integrating (4.10) to give Thus, we obtain δY (T + ) = K 1 δY − , with the saltation matrix K 1 given by The other saltation matrices (describing the passage through -neighbourhoods of U = 0 and V = 0) are constructed in a similar fashion, and found to be It is straightforward to check that the saltation matrices (4.12)-(4.13) are equivalent to those defined by (4.9). Between switching events the perturbations evolve according to exp(A(t − T ))δY (T + ), for t > T , where δY (T + ) is the perturbation at the switching time. Thus, after one period of oscillation, we may put this all together to obtain (4.14) The periodic orbit will be stable if the eigenvalues of Γ 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 Γ = e σΔ × 1. Hence, det e AΔ+− det e AΔ−− det e AΔ−+ det e AΔ++ . (4.15) Using the fact that det e At = e Tr A t , we find . (4.16) A periodic orbit will be stable provided σ < 0. We shall say that the pseudo-equilibrium at (0, 0) is unstable (stable) if it is enclosed by a stable (unstable) periodic orbit of arbitrarily small amplitude. We shall say that there is a pseudo-Hopf bifurcation at (0, 0) when the pseudo-equilibrium changes stability, namely when σ = 0. A plot of σ = σ(τ) (not shown) for the parameters of Figure 2, shows very similar behaviour as for the steep PWL firing rate function. In essence, we may regard the second term on the right-hand side of (4.16) as a correction term to standard Floquet theory to cope with the non-smooth nature of the Heaviside firing rate.

An unstable periodic sliding orbit
The Wilson-Cowan node can also support an unstable periodic orbit that has a component which slides along the switching manifold U = 0 for V ∈ [V 1 , V 2 ], as depicted in Figure 5. The points V 1,2 are easily calculated by determining the points at which the U-nullclines touch the switching manifold, where U = 0, and are found to be V 1 = (A 11 I u + A 12 I V − w uu )/A 12 and V 2 = V 1 + w uu /A 12 . In reverse time initial data close to a sliding trajectory would be attracted to it. Thus we can think of constructing an unstable periodic sliding orbit, of the type shown in Figure 5, by breaking it into five pieces. All pieces of this orbit are constructed similarly to before (see above), except the component that slides. Using the Filippov method and equation (4.7), we find κ 2 = (A 11 I u − A 12 V + A 12 I v )/w uu , with the sliding dynamics prescribed by d dt In backward time, the periodic sliding orbit shown in Figure 5 would slide up along U = 0 until the point V = V 2 , where it would leave the switching manifold. We now turn our attention to networks built from Wilson-Cowan nodes with a Heaviside firing rate.

A network of Heaviside Wilson-Cowan nodes
As we have shown in Section 4, the replacement of a sigmoidal firing rate by a Heaviside function can lead to highly tractable models for which substantial analytical results can be obtained (with the use of matrix exponentials and saltation matrices). However, at the network level the mathematical differences between the treatment of smooth and nonsmooth firing rates are considerably amplified relative to those at the single node level. At the node level, it is well-known that regarding the Heaviside function as the steep limit of a sigmoidal function can lead to arbitrarily many different non-equivalent dynamical systems. This is simply due to the non-uniqueness of the singular limits by which smooth functions may tend towards discontinuities. For a recent perspective on this issue, see the work of Jeffrey [27]. Thus, there is no reason to assume that taking the limit → 0 for the PWL network considered in Section 3 will be relevant to a Wilson-Cowan network with a Heaviside non-linearity. Namely the approximation of a Heaviside function by a continuous function such that H(x) = lim →0 F(x), where F(x) is given by (2.2), may have little utility given that pointwise convergence need not imply distributional convergence.
We now return to the network introduced in Section 3, but replace the dynamics of each node with the Heaviside limit studied in the previous section. For the following analysis, it is convenient to rewrite (3.4) as The synchronous network state is given by (4.8) (remembering the row-sum constraint on the network connections). To study its linear stability, we consider values of the perturbed network state Y that are close to the synchronous network state at the unperturbed crossing times. Let T i denote the time that the synchronous state moves between one of the four quadrants (as illustrated in Figure 5). We then make the ansatz that the perturbed network state Y can be expressed with respect to the synchronous orbit at one of the switching times T i and write Y (t) = Y (T i ) + δY (t) with t in the neighbourhood of T i .
We first construct the saltation matrix through a switch, indexed by i = 1, . . . , 4. Suppose that the kth crossing occurs at a perturbed crossing time T i,k . The network states at two consecutive crossings are related via This equation is obtained by integrating (5.1) using the observation that F is constant between crossings. By linearising (5.2), we can relate the perturbations between crossing events as where e m is the mth canonical basis vector in R 2N . The saltation matrix for each of the four switches is then given by The ordering of matrix multiplications in (5.6) is determined by the iterative minimisation of the perturbations given by (5.4).
In the next step, we analyse how a perturbed network state is propagated between saltation events. Let T + i denote the time when the last node crosses between quadrants. Here, the superscript makes explicit that all nodes have crossed into the next quadrant. The next network event occurs when one of the nodes crosses into the subsequent quadrant. This happens at a time T − i+1 , where the superscript indicates that only one node has crossed. We will make the ansatz that from which we obtain after linearisation where we have used the fact that F(Y (T + i )) = F(Y (T − i+1 )), since F is constant between crossing events. Here, Y (t) denotes the differential of Y (t) with respect to t. As above, the component of δY (T − i+1 ) that corresponds to the node that switches first, say at position m, vanishes. Taking the mth component of (5.8) then yields an expression for the perturbation of the crossing time , (5.9) where the vector f i ∈ R 2N is given by e AΔi δY ( We again find the value of m by minimising (5.9) over all admissible values of m and refer to it as m i . This and . (5.11) Taken together, we obtain after one period The matrices Γ i act to propagate perturbations across a quadrant, and the L i propagate perturbations through a switch. At first sight, the definition of G i suggests that we have introduced a dependence of Γ i on δY (0) through the inclusion of δY (T + i ). This dependence can be avoided by noting that δT + i = δT − i + k δT i,k and the repeated use of (5.4), (5.5), and (5.9). The drawback of this approach is that the resultant operator does not lend itself to an interpretation of successive propagations and saltations, nor is it numerically advantageous. Moreover, this operator would only remove the explicit dependence of Ψ on δY (0). The minimisation steps that are necessary to determine the order in which nodes switch already leads to an implicit dependence of Ψ on δY (0). Changing δY (0) can lead to a different order of switching, and since matrix multiplication does not commute, Ψ can be different for different δY (0). This has profound implications for asserting linear stability. The usual argument that the eigenvalues of Ψ determine linear stability does not hold anymore. To see this, consider the propagation of δY (0) over multiple periods, i.e., (2) , . . . (5.13) so that 14) The eigenvalues of Ψ (i) and Ψ (j) can be different for i = j. For some value of i, Ψ (i) may have all eigenvalues in the unit disc, whilst for another value of i there may be some eigenvalues outside the unit disc. Over one period, perturbations can therefore grow or shrink. This entails that for a product of operators as in (5.14), δY (m) may be smaller than δY (0) , although some Ψ (i) might have some eigenvalues that lie outside the unit disc. Instead of looking at the eigenvalues of individual Ψ (i) , we could have studied the eigenvalues of the product of operators in (5.14). We would have come to the same conclusion since eigenvalues of the product operator move into and out of the unit disc as we increase m. Figures 6 and 7 illustrate the dependence of the spectra on random initial conditions δY (0). In both figures, the left panel shows the spectra for initial conditions when all eigenvalues of Ψ (0) lie within the unit disc. The middle panel displays spectra with some eigenvalues outside the unit disc, and the right panel is a blow-up of the middle panel around the unit disc. For Figure 6, we chose a value of σ such that the synchronous orbit of the PWL network, with a small values of = 0.001, is linearly stable. We observe that the eigenvalues of the Heaviside network cluster around those of the PWL network. While it appears that the majority of synchronous solutions are stable (for this parameter choice), some initial conditions lead to eigenvalues outside the unit disc. When zooming into the unit disc, we see some degree of clustering, although this is not as pronounced as for the stable solutions. For larger values of σ, the synchronous state of the PWL network becomes unstable (for small ). The left panel of Figure 7 shows that the eigenvalues of the Heaviside network that all fall into the unit disc exhibit only a weak association with the eigenvalues of the PWL network. In addition, it seems that more initial conditions lead to unstable synchronous solutions than stable ones. This mirrors the behaviour in Figure 6, where the majority of initial conditions gives rise to stable solutions. The blow-up in the right panel of Figure 7 illustrates that the eigenvalues of the Heaviside network form clusters around those of the PWL network. While the notion of linear stability in terms of eigenvalues of the propagator is lost for the Heaviside network, it appears that the clustering of these eigenvalues reflects the stability of the PWL system, at least for small values of (where the PWL firing rate becomes more switch like).

Conclusion
In this paper, we have shown that the combination of two popular approaches in dynamical systems, namely PWL modelling of low-dimensional oscillators and the MSF, can be combined to give insight into the behaviour of network states in neural mass network models. This is natural for this type of system since the sigmoidal non-linearity, ubiquitous throughout neuroscience modelling of large-scale brain dynamics, is wellcaricatured by a PWL reduction. We have focussed here on the bifurcation of the synchronous network state, and shown how this can be determined in terms of a set of lowdimensional Floquet problems, each of which can be solved using simple linear algebra. In essence, the PWL aspect of the model allows the variational problem for stability to be solved without recourse to the numerical solution of an ordinary differential equation. Closed form solutions are patched together, and although this may appear inelegant at first sight, it does lead to explicit formulas for Floquet exponents at the single node level, and is easily cast into algorithmic form for accurate numerical computations at the network level. This nicely highlights the benefits of PWL modelling. Importantly, the approach advocated here is not just limited to the construction and stability of the synchronous state. Pecora et al. [40] and Sorrentino et al. [43] have recently extended the MSF approach to treat more exotic states making extensive use of tools from computational group theory. Thus, the work presented here is readily extended to treat non-synchronous states, such as clusters, and for a further discussion see [37]. From a neuroscience perspective, it would also be important to treat delays, arising from the finite propagation speed of action potentials relaying signals between distinct brain regions [13]. In this case, we would hope to exploit the growing body of knowledge on PWL dynamics with time delay, as exemplified by [42].
From a mathematical perspective, we have also seen that there is an important difference between the analysis of a high gain continuous PWL sigmoid and that of a discontinuous switch-like Heaviside firing rate. Although this can be facilitated with the use of saltation matrices (to propagate perturbations through switching manifolds) there is no MSF style approach that reduces the study of synchrony to a set of sub-network Floquet problems. Moreover, in contrast to the linear stability analysis of continuous systems, there is now a new challenge of addressing the temporal order in which perturbations to network states pass through a switching manifold. To treat this, we have made use of ideas originally developed for Glass networks [17], though note that similar issues of ordering also arise in the analysis of pulse-coupled systems [23,28,44]. In essence, the analysis of a Wilson-Cowan network with a Heaviside firing rate must be performed carefully, and with non-standard tools, as its behaviour can differ from that of a similar network with a high gain PWL sigmoid.