Insights into oscillator network dynamics using a phase-isostable framework

Networks of coupled nonlinear oscillators can display a wide range of emergent behaviours under variation of the strength of the coupling. Network equations for pairs of coupled oscillators where the dynamics of each node is described by the evolution of its phase and slowest decaying isostable coordinate have previously been shown to capture bifurcations and dynamics of the network which cannot be explained through standard phase reduction. An alternative framework using isostable coordinates to obtain higher-order phase reductions has also demonstrated a similar descriptive ability for two oscillators. In this work we consider the phase-isostable network equations for an arbitrary but finite number of identical coupled oscillators, obtaining conditions required for stability of phase-locked states including synchrony. For the mean-field complex Ginzburg-Landau equation where the solutions of the full system are known, we compare the accuracy of the phase-isostable network equations and higher-order phase reductions in capturing bifurcations of phase-locked states. We find the former to be the more accurate and therefore employ this to investigate the dynamics of globally linearly coupled networks of Morris-Lecar neuron models (both two and many nodes). We observe qualitative correspondence between results from numerical simulations of the full system and the phase-isostable description demonstrating that in both small and large networks the phase-isostable framework is able to capture dynamics that the first-order phase description cannot.

Networks of coupled nonlinear oscillators can display a wide range of emergent behaviours under variation of the strength of the coupling.Network equations for pairs of coupled oscillators where the dynamics of each node is described by the evolution of its phase and slowest decaying isostable coordinate have previously been shown to capture bifurcations and dynamics of the network which cannot be explained through standard phase reduction.An alternative framework using isostable coordinates to obtain higher-order phase reductions has also demonstrated a similar descriptive ability for two oscillators.In this work we consider the phase-isostable network equations for an arbitrary but finite number of identical coupled oscillators, obtaining conditions required for stability of phase-locked states including synchrony.For the mean-field complex Ginzburg-Landau equation where the solutions of the full system are known, we compare the accuracy of the phase-isostable network equations and higher-order phase reductions in capturing bifurcations of phase-locked states.We find the former to be the more accurate and therefore employ this to investigate the dynamics of globally linearly coupled networks of Morris-Lecar neuron models (both two and many nodes).We observe qualitative correspondence between results from numerical simulations of the full system and the phase-isostable description demonstrating that in both small and large networks the phase-isostable framework is able to capture dynamics that the first-order phase description cannot.
The utility of the classical technique of phase reduction for describing the dynamics of networks of coupled oscillators is limited by the assumption that the dynamics for each node remain on the stable limit cycle of the uncoupled system.We here investigate reduced equations for networks of arbitrary finite size where the dynamics of each node is described in terms of its phase and the slowest decaying isostable coordinate, allowing for representation of trajectories away from (but near) the limit cycle.Specifically, we consider conditions for the existence and stability of phase-locked states generalising existing results for phase-reduced equations.We show that phase-isostable network equations provide the most accurate qualitative description of dynamics of the mean-field Ginzburg-Landau equation when compared with alternative frameworks such as higher-order phase reduction 1,2 .We further demonstrate the power of the general framework by considering networks of neural oscillators.We observe phenomena including the emergence of quasiperiodic behaviour that cannot be captured using first-order phase reduction.The results are shown to be in good qualitative agreement with the dynamics of the original network through numerical simulations and bifurcation analysis.

I. INTRODUCTION
Oscillations observed in biological, physical and chemical systems are often due to the presence of attracting limit cycles within the high dimensional dynamics 3,4 .The classical technique of first-order phase reduction [5][6][7][8][9][10] provides a rigorous way of describing the dynamics of weakly perturbed oscillators in terms of a single phase variable using the notion of isochrons that extend the phase variable for a limitcycle attractor to its basin of attraction 11,12 .The power of this approach has been demonstrated through its ability to reveal complex dynamics of weakly forced oscillations and emergent behaviours in weakly coupled oscillator networks in a variety of relevant systems 8,[13][14][15][16] .However, phase reduction assumes that the dynamics remain close to the unperturbed limit cycle and therefore requires that interactions are weak and convergence to the limit cycle is fast.The ability of phase reduced equations to accurately portray the dynamics of the full system diminishes and eventually breaks down with increasing interaction strength.
In recent years various strategies have been proposed to overcome the limitations of first-order phase reduction.The most widely employed approach is to include additional variables in directions transverse to the limit cycle which allows for a phase-amplitude description of transient trajectories away from the limit cycle.While some authors have used moving orthonormal coordinate systems to define the amplitudes 17,18 , there is a growing body of work that utilises the notion of isostable coordinates [19][20][21] .These coordinates represent level sets of the slowest decaying eigenmodes of the Koopman operator 20,22 and in the absence of perturbations, have exponential decay to zero within the basin of attraction of the limit cycle at rates given by the Floquet exponents of the linearisation of the flow about the limit cycle.First defined for systems with a stable fixed point 20 , isostable coordinates are now a well established concept for the analysis of perturbed trajectories near periodic orbits [23][24][25][26][27] and also for the control of oscillators (see 28 for a recent review).
A dimension reduction can be achieved for highdimensional systems by assuming that isostable coordinates which rapidly decay (associated with negative Floquet exponents of large magnitude) are always zero 27,29 .This results in a phase-amplitude reduced description of single oscillator dynamics in response to perturbations which can be taken to arbitrary orders of accuracy in the isostable coordinates [29][30][31] .In turn, through high order asymptotic expansion of interactions in terms of the phase and isostable coordinates 1 , it is possible to build network equations describing the evolution of phase and isostable coordinates of each node in a coupled oscillator network.To date, analysis of these network equations has been relatively limited and mainly restricted to networks with just two oscillators 1,[32][33][34] .
Two different but related frameworks for the use of the phase-isostable network differential equations for the study of oscillator network dynamics have emerged.The most frequently employed approach has been to use the isostable dynamics to derive a phase equation to higher order in the coupling strength.The isostable dynamics are assumed to be slaved to the phases so that at each order in the coupling strength the isostables are described by a linear ordinary differential equation 1,32,35,36 .This gives a higher-order description of the phase dynamics when the solution is substituted into the equation for the time evolution of the phase.The higher-order phase reduced equations for networks of more than two nodes contain non-pairwise phase interactions (i.e., terms involving the phases of three or more oscillators, noted in 32,35,36 , but overlooked in 1 ), despite the interactions between the nonlinear oscillators being pairwise.This is a typical result for higher-order phase reduced network equations 37 .Similar non-pairwise phase interactions were identified in 2 where an alternative strategy was used to compute phase reduced equations to second and third order in the coupling strength for the mean-field complex Ginzburg-Landau equation relying on the explicit expression for the isochrons (which can be obtained in this case).
The alternative framework retains both phase and one isostable coordinate for each node in the coupled oscillator network.One work to have adopted this approach is that of Ermentrout and Wilson 33 .Similarly to 1,32 , they consider as an example a pair of synaptically coupled thalamic neurons.All three works 1,32,33 show that their approaches can reveal bifurcations leading to stable phase-locked states at higher coupling strength as observed in full system simulations.Coombes et al., 34 also use network equations for the dynamics of phases and isostable coordinates to consider two linearly coupled planar piecewise-linear caricatures of the Morris-Lecar neuron model.Accounting for the nonsmoothness of the dynamics using results from 30 , a bifurcation diagram for the dynamics in the phase-isostable framework under variation of coupling strength reveals restabilisation of synchrony, the existence of stable phase-locked states other than synchrony and antisynchrony and also stable quasiperiodic states, none of which can be captured using first-order phase reduction of the dynamics (which predicts only stable antisynchrony and unstable synchrony for weak coupling).All of the dynamics revealed by the phase-isostable network analysis are shown to be in qualitative agreement with numerical simulations of the full system, however comparison with results using the more accurate master stability function 38 shows that the phase-isostable approach underesti-mates the value of the coupling strength at which synchrony restabilises 34 .
Previous work therefore indicates that for networks of two oscillators both phase-amplitude network equations and higher-order phase reduction based on phase-amplitude network interactions have the capability to reveal dynamics which cannot be explained using standard phase reduction.
It remains an open problem to fully explore the capabilities of both approaches and draw comparisons between the two in terms of the accuracy with which they can capture the existence and stability of network states and to extend analysis of the approaches to larger networks.
In this paper we go some way to address these challenges.Following discussion of the necessary background on isostable coordinates in section II, we derive the phaseisostable network differential equations and also the corresponding higher-order phase-reduced equations for a network of N identical oscillators in sections II and IV respectively.Using first-order averaging we obtain a system of phaseisostable network differential equations linear in the isostable coordinates and with six pairwise phase difference interaction functions.These equations are used to determine conditions for the existence and stability of various phase-locked states in section III.We are then able in section IV to compare results regarding the location of stability boundaries for full synchrony and the asynchronous splay state arising from each reduction method for the mean-field complex Ginzburg-Landau equation where it is also possible to derive the exact location of the stability boundaries for the full system.We observe that retaining the isostable coordinate through using the phase-isostable network equations gives the greatest accuracy in approximating the qualitative bifurcation structure.
We therefore proceed to use the phase-isostable network equations to approximate the dynamics of both small and large networks of neural oscillators in section V.For a two node network of linearly coupled Morris-Lecar neurons 39 we compare the bifurcation diagrams for the full model with that for the phase-isostable network equations.We observe that the phase-isostable description qualitatively agrees with the full model in capturing bifurcations of branches of synchronous and antisynchronous solutions in addition to the existence of stable phase-locked states.This marks a vast improvement on the descriptive power of the first-order phase reduction which cannot capture any of these phenomena.We also investigate a larger network of many globally linearly coupled Morris-Lecar neurons where we find that as for the full model considered in 40 there is an interval in the coupling strength where synchrony and the uniform incoherent (splay) state are both unstable.In this interval we observe stable cluster states in both the full dynamics and the phase-amplitude reduced equations.Finally in section VI we discuss the limitations of the framework as well as possible extensions of our work to further refine the accuracy of the description of the network dynamics.We also highlight other phenomena which may be revealed in driven phase oscillator dynamics through using phase-amplitude coordinates.

II. PHASE AND ISOSTABLE COORDINATES AND REDUCTIONS
We consider the dynamics of N identical pairwise coupled oscillators where x i ∈ R n .The coupling function G : R n × R n → R n describes the pairwise interactions which take the same form for each pair of oscillators.The overall strength of interactions is given by ε while w i j modulates the weight of connectivity from node j to node i, effectively describing the topology of the oscillator network.The uncoupled dynamics (for ε = 0) are described by F : R n → R n , where we assume that has a T-periodic hyperbolic limit cycle, γ and we denote points on γ by x γ (t).Let φ : R × R n → R n denote the flow induced by (2) so that φ (t, x 0 ) represents the solution of (2) with initial condition x(0) = x 0 .Then the monodromy matrix M(x γ (t * )) = ∂ φ (T, x)/∂ x| x=x γ (t * ) (time-T flow linearised about a point x γ (t * ) on the orbit) has eigenvalues λ m , m = 0, . . ., n − 1 called the Floquet (or characteristic) multipliers and we assume that 1 = λ 0 > |λ 1 | ≥ . . .≥ |λ n−1 | so that the limit cycle γ is stable.For simplicity we further assume that the Floquet multipliers λ m are positive, real and simple (and see 30 for discussion on how to proceed when this is not the case).The non-zero Floquet exponents are then κ m = log(λ m )/T , m = 1, . . ., n − 1 which are all real and negative.
A. First-order phase reduction and weakly coupled phase oscillators The periodic orbit γ of (2) can be parameterised by a phase θ ∈ [0, 2π) with an arbitrary point x γ 0 ∈ γ having θ (x γ 0 ) = 0 and θ (φ (t, x γ 0 )) = ωt where ω = 2π/T .The notion of phase can be extended to points x * in the basin of attraction of the limit cycle B(γ) by defining the asymptotic phase as the unique ( The level sets of θ (x * ) are called isochrons and contain all points in B(γ) with the same asymptotic phase 3,4,11,12 .The phase of each uncoupled node trajectory then evolves as dθ /dt = ω both on and off the limit cycle.Note that we can then use x γ (t) and x γ (θ ) to denote points on cycle where time is parameterised as t = θ /ω.In terms of phase variables, (1) becomes where is the phase response curve which quantifies the effects of a perturbation on the phase of the oscillator.If the coupling is weak, then the dynamics stay in a neighbourhood of the limit cycle and evaluating on the limit cycle we have to first order in ε, where θ is the gradient of the phase variable θ evaluated on the limit cycle and is known as the infinitesimal phase response curve (iPRC).While Z may be found directly 4,6,10 , it most commonly computed using the adjoint method 6,16,41 whereby Z is the T -periodic solution of the adjoint equation satisfying the normalisation condition Z(0) Here T denotes the transpose and J := DF(x γ (θ )) is the Jacobian of the vector field F evaluated on cycle.
A further simplification of ( 5) can be made by transforming the system to a rotating frame through θ i = φ i + ωt where, assuming weak coupling, φ i will slowly drift.First order averaging 16,41,42 then gives the phase dynamics back in the original variables as the simpler to analyse phase difference equations where is known as the phase interaction function.A solution of the averaged equations ( 7) is ε-close to a solution of the unaveraged equations (5) for times of O(ε −1 ).The existence and stability of phase-locked states in the first-order phase regime can be studied using (7) and specified in terms of the properties of H 1 and the choice of network structure given by W = (w i j ) [43][44][45][46][47] .In section II B we will show how the theory of isostables can be used to derive averaged phase-isostable network equations with six interaction functions H k , k = 1, . . ., 6 and in section III we investigate the existence and stability of phase locked states for larger magnitude interactions in terms of properties of the H k , and network structure.Equation (5) gives what we will refer to as the first-order phase reduction of the network dynamics.For further details on the phase reduction of oscillator networks see the recent reviews 8,48 .Phase-reduced equations have provided the standard framework for understanding dynamics of weakly coupled oscillator networks for the last four decades and can be extremely instructive.For instance, the Kuramoto model 3 has simply H 1 (χ) = sin(χ) yet is able to capture the basic mechanisms underlying synchronisation in many biological 49 , chemical 50 and physical 51 oscillator networks.However, the phase-reduction ( 5) is only valid provided the state of each oscillator remains close to its underlying limit cycle.This requires that interactions are weak (|ε| ≪ 1) and therefore (5)  cannot be used to describe network dynamics with interactions of larger magnitude which could potentially lead to oscillators spending significant time away from the limit cycle.In order to capture transient dynamics away from the limit cycle we introduce amplitude variables in directions transverse to the limit cycle.

B. Isostables and phase-isostable reduction of network equations
In this paper we use the concept of isostables to define coordinates off limit cycle 10,[19][20][21]26,27,30 . Isostabls for the uncoupled system (2) identify initial conditions with the same relaxation rate to the limit cycle and therefore approach the limit cycle together 26 .For each Floquet exponent κ m , m = 1, . . . ,  − 1, a set of isostables representing an amplitude degree of freedom can be defined and therefore for any point x ∈ B(γ) we can associate isostable coordinates ψ m , m = 1, . . ., n − 1.The isostables can be defined as level sets of certain eigenfunctions of the Koopman operator 20,25,26,52 and a constructive definition may be given for the slowest decaying isostable ψ 1 27,31 : Denote by w T the left eigenvector of M(x γ 0 ) associated with the eigenvalue where t k θ =0 is the time of the kth transversal of the θ = 0 isochron.The isostable coordinates are amplitudes which satisfy ψm = κ m ψ m in the absence of coupling or perturbations.
An alternative perspective is taken by the authors of 19,23,24 who note that for each point x ∈ B(γ) there is a coordinate transform given by an analytic map K such that x = K(θ , ψ 1 , . . .ψ n−1 ), and determine the map K using a parameterization method.
For oscillator i in the coupled network (1) we can now describe the phase and isostable dynamics by where the gradients are the phase and isostable response curves respectively that quantify the effects of a perturbation on the phase and amplitude coordinates of the oscillator.We now make the simplifying assumption that κ 2 , . . ., κ n−1 are all sufficiently close to zero that perturbations in the directions of isostables (ψ i ) 2 , . . ., (ψ i ) n−1 may be ignored.That is we consider a single isostable coordinate (ψ i ) 1 := ψ i for each oscillator with corresponding Floquet exponent κ 1 := κ.This assumption is similar to that made in previous work 1,27,32 .We have therefore (for node dynamics of dimension n > 2) made a reduction in order from Nn equations to 2N equations for the network dynamics, retaining only the most slowly decaying isostable dynamics for each oscillator.
Following 1,31 we can take asymptotic expansions of solutions away from the limit cycle and also the phase and isostable response curves about the limit cycle to arbitrary order in the O(ε) isostable coordinate.That is we can write where g (1) (θ ) is the Floquet eigenfunction (right eigenvector of M(x γ (θ )) associated with eigenvalue e κT ) and g (k) (θ ), k > 1 are higher order analogues.The gradient of the phase and amplitude coordinates evaluated on the limit cycle at phase θ are the iPRC Z (0) (θ ) and the infinitesimal isostable response curve (iIRC) I (0) (θ ) respectively.The terms Z (k) (θ ) and I (k) (θ ), k > 0, are higher order correction terms to Z (0) (θ ) and I (0) (θ ).Note that this generalises the linear order expansions developed in 27,29,30 .All of the T -periodic vector functions g (k) , Z (k) , I (k) can be computed using appropriately normalised adjoint equations 31 : where I n is the n × n identity matrix, e i are the unit basis vectors and are the vectors defined in (A7) and (A12) respectively.See 31 for the derivation of these equations for all N isostable coordinates and Appendix A for the derivation of the compact notation (12) when a single isostable coordinate is considered.The appropriate normalisations for (12) are also given in Appendix A. In general the solutions of the hierarchy of equations ( 12) must be computed numerically 1,10 .
In order to obtain the network equations (10) in terms of phase and isostable coordinates, we also require the appropriate expansion for the coupling function G : R n × R n → R n which can be expressed to arbitrary order following 1 .Here we will require only the expansion up to and including second order derivatives of G and therefore we may use more convenient notation than that in 1 T , the Taylor expansion of G can be expressed to second order as where DG is the Jacobian of G with respect to X and H q is the Hessian matrix of second order derivatives of G q where and all derivatives are evaluated at X γ .Since where is the Jacobian of G(x i , x j ) with respect to its kth argument and we can rewrite (13) as Using the expansions of ∆x i and ∆x j from (11a) and collecting terms up to quadratic order in ψ we have where We now have all of the expansions required to express (1) in terms of phase and isostable coordinates.Including terms up to O(ε 2 ) the phase-amplitude network equations are where the six interaction functions h k (θ j , θ j ) are given by The expansion of the network equations may be computed to higher order, but note that to O(ε m ) there are m(m + 1) interaction functions.Here we explore the capabilities of the O(ε 2 ) expansion ( 18) to describe network behaviours and discuss the limitations imposed by the truncation at this order in section VI.The network equations (and therefore analysis) may be simplified to a phase difference system by first-order averaging 33,42,53 yielding where H k are the 2π-periodic functions In the weak coupling limit where the isostable coordinates may be assumed to be zero, we recover the first-order phase reduced equation (7).The equations (20) generalise to networks of arbitrary size and structure those given in 33 and contain higher order terms than those derived in 32 .Analysis of existence and stability of phase-locked states in ( 20) is investigated next in Section III.Note that since (20) are averaged equations, the solutions of ( 18) and ( 20) are ε-close for times of O(ε −1 ).

III. PHASE-LOCKED STATES IN PHASE-AMPLITUDE NETWORK EQUATIONS
Having defined a system of phase-amplitude network equations (20), we now investigate conditions for the existence and stability of certain phase-locked states within the averaged equations, generalising well known results for phase-reduced networks [43][44][45]47,54 .
A 1 : 1 phase-locked state in the network of N identical oscillators (20) is defined by θ i = φ i + Ωt where the φ i are constant phase lags and Ω is the collective frequency of the oscillators.We denote such a state by Φ = (φ 1 , . . ., φ N ).In phase-reduced systems all node orbits are assumed to coincide with the stable limit cycle, γ of the node dynamics.By including information about the dynamics off-cycle through the isostable coordinates ψ i we allow each node to have a different trajectory in the phase-isostable phase space.That is, for a solution (θ 1 , . . ., θ N , ψ 1 , . . ., ψ N ) of ( 20), the projection ) can be different for each node.Substituting into the O(ε 2 ) truncated averaged equations (20) we have for i = 1, . . ., N, For a given fixed set of relative phases Φ and connectivity W = (w i j ), equations (22a) are N linear equations from which we determine that ψ i (t) are constants depending on Ω.That is, any phase-locked solution of (20) has node orbits which coincide with an isostable.Note that this result is particular to the truncation of the network equations (20) at O(ε 2 ) which includes only linear terms in the isostable coordinates and would not be expected to hold for higher order truncation.Denoting the constant isostable coordinates ψ i (t) = Ψ i we find from ( 22) that a phase locked state with relative phases Φ exists with collective frequency Ω only if where 1 N is the column vector of 1s, p = (p 1 , . . ., p N ) T , q = (q 1 , . . ., q N ) T with and P and Q are matrices with entries In this case the constant isostable values for the nodes are given by Ψ = (Ψ 1 , . . ., Ψ N ) T = Q −1 q.We now consider conditions for the stability of a general phase-locked state Φ = (φ 1 , . . ., φ N ) in the O(ε 2 ) truncated phase-isostable network equations (20).We linearise about the phase-locked solution by setting where ∆θ i (t) and ∆ψ i (t) are perturbations along and transverse to the periodic orbit Γ respectively.The linearisation of (20) gives where i j (Φ)∆θ j + H (2) and hence the Jacobian is the 2N × 2N block matrix where the N × N matrices H (q) (Φ) are given in Appendix B. The phase-locked state Φ is stable if the eigenvalues of J all have negative real part with the exception of the eigenvalue which is forced to be zero corresponding to perturbations along the periodic orbit.In general the H (q) are not in graph-Laplacian form as they are for the phase-only case 47,54 and therefore the eigenvalues of J cannot be related directly to the eigenvalues of W = (w i j ).Nonetheless we may still derive conditions required for the existence and stability of particular phase-locked solutions.

A. Phase-locked states with a common orbit
We first consider phase-locked states where all oscillators have a common orbit with isostable coordinate Ψ i = Ψ for all i = 1, . . ., N. Such states include full synchrony (φ i = 0 for all i), the splay (or uniform incoherent) state (φ i = 2πi/N) and balanced clusters where the N = Mq nodes lie in M clusters, each containing q nodes and phase difference between the clusters is 2π/M.In Section III B we consider unbalanced two-cluster states where each cluster is a different size and has a different orbit.

Full synchrony
Synchrony is the most widely studied example of a 1 : 1 phase-locked state.We extend established results regarding the stability of synchrony in weakly coupled (phase-reduced) networks by using the phase-isostable network equations (20).Here all nodes share a common phase θ 1 (t) = . . .= θ N (t) with θi = Ω and orbit with isostable coordinate Ψ.Without loss of generality, let φ i = 0. Then from ( 22) for all i = 1, . . ., N. Denoting the row sum ∑ N j=1 w i j = c i , we find that The values of Ω and Ψ must be identical for all i to guarantee the existence of the synchronous state.There are two ways this may be achieved.
a.If the row sums c i are independent of i (i.e.c i ≡ c for all i).This is the case for global coupling where w i j = 1/N which we will consider further later.
We consider the conditions required for linear stability in both cases.a. Connectivity matrix with constant row sums Taking c i = c for all i, (B1a)-(B1d) give the Jacobian for the synchronous state as where L is the graph Laplacian matrix with L i j = −w i j + δ i j ∑ N k=1 w ik , I N is the N × N identity matrix and ⊗ is the Kronecker (tensor) product.The eigenvalues of (33) depend on the particular choice of W with constant row sums.In the case of global coupling where w i j = 1/N so that c = 1 we see that J has non-zero eigenvalues κ + ε(H 5 (0) + H 6 (0)) and µ ± each of multiplicity N − 1 where µ ± are the eigenvalues of Therefore synchrony is stable when κ + ε(H 5 (0) + H 6 (0)) < 0, Trace(M (Ψ)) < 0 and det(M (Ψ)) > 0. Reducing back to the phase-only description by taking H 2 , . . ., H 6 ≡ 0 we recover the result for phase oscillators that synchrony is stable if −εH ′ 1 (0) < 0 43,47,54 .b. Interaction functions with H 1 (0) = 0 and H 4 (0) = 0 If H 1 (0) = 0 and H 4 (0) = 0 then Ψ = 0 and therefore the synchronous orbit coincides with the stable uncoupled node orbit and Ω = ω.The Jacobian is given by For the case of diffusive coupling which is linear in however there are no general results concerning the eigenvalues of J .In the case where coupling is global so that L i j = −1/N + δ i j the non-zero eigenvalues of J are κ and µ ± each of multiplicity N − 1 where µ ± are the eigenvalues of M (0) in (34) so that synchrony is stable for global diffusive coupling when κ + εH 5 (0

The splay state
Another important example of a 1 : 1 phase-locked solution is the splay state, sometimes referred to as a rotating wave or the uniform incoherent state (UIS).Here the phases are uniformly distributed around the unit circle and without loss of generality we take θ i+1 − θ i = 2π/N and Φ = (φ 1 , . . ., φ N ) where φ i = 2πi/N.The equations for the collective frequency Ω and collective orbit Ψ, given by ( 22) must be satisfied for all values of i for the existence of the splay state.In the case of global coupling, w i j = 1/N, we find that the splay state exists with where Considering now the stability of the state we see that the Jacobian takes the form (29) where for and where are the eigenvalues of the circulant matrix and hence Having found the eigenvalues of each of the blocks H (p) it remains to determine the eigenvalues µ of the Jacobian (29).Since each block H (p) is circulant, they can each be diagonalised as N−1 ), and therefore they commute.By 55 (Theorem 1), the eigenvalues µ then satisfy The characteristic equation for the 2N eigenvalues of J is therefore Hence the eigenvalues are µ ± q for q = 0, . . ., N − 1 where µ ± q are the eigenvalues of The splay state is stable if Re(µ ± q ) < 0 for all q = 0, . . ., N − 1 with the exception of the zero eigenvalue corresponding to perturbations around the limit cycle.
a.The splay state in the large N limit In the thermodynamic limit N → ∞ for global coupling, w i j = 1/N, the network averages are effectively Riemann sums and so may be replaced by time averages: where H k has a Fourier series representation Consequently Ω and Ψ from (37) become (1) q =ε (( The eigenvalues of the Jacobian J for the splay state in the large N limit are given by the eigenvalues of (45) with λ (p) q , p = 1, 2, 3, 4 now given by (48).The splay state is stable if Re(µ ± q ) < 0 for all non-zero q.It is straightforward to show that reducing back to phase-only by letting Ψ = 0 and taking |ε| ≪ |κ| we recover the result that for globally coupled phase oscillators in the thermodynamic limit, the splay state is stable if εqIm((H 1 ) −q ) < 0 for all non-zero q 43,45,47,54 .

Balanced cluster states in globally coupled networks
Consider now the more general states sharing a common orbit in a globally coupled network, w i j = 1/N, where the N = Mm nodes form M clusters, each with m nodes.Symmetry imposes that the phase difference between any of the M clusters is an integer multiple of 2π/M and without loss of generality we order the nodes such that where each distinct integer value occurs m times.From (23)  we find that the collective frequency and isostable coordinate of the collective orbit are where Considering stability of such a state, the Jacobian is given by (29), where each block H (p) (Φ) has block circulant structure: where 1 m denotes the m × m matrix of ones and where 47,56,57 By an application of the matrix determinant lemma 58 the eigenvalues λ (p) and λ Furthermore, each H (p) can be diagonalised as where with λ (p) repeated N − M times.As for the splay state, the H (p) therefore commute and so that the characteristic equation of J is q λ (3) q .
Hence the eigenvalues of J are the eigenvalues µ ± of each of multiplicity N − M (corresponding to intracluster perturbations) and the eigenvalues µ ± q of Λ q as in (45) for q = 0, . . .M − 1 where λ (p) q are as in (42) replacing N with M (the intercluster eigenvalues, which are the eigenvalues of the splay state in an M node network).Note that Λ 0 has a zero eigenvalue (the purely rotational eigenvalue).It is clear that the splay state is a special case where M = N, m = 1.Synchrony is the special case where M = 1, m = N and we note that the stability conditions for synchrony in globally coupled networks are recovered since here the non-zero eigenvalue of Λ 0 is κ + ε(H 5 (0) + H 6 (0)) and Λ = M(Ψ) as in (34).

B. Two-cluster states in globally coupled networks
We now consider network states where the N nodes in a globally coupled network (w i j = 1/N) lie in two clusters denoted C A and C B containing N A ≤ N/2 and N B = N −N A nodes respectively.We label nodes such that i ∈ C A for i = 1, . . ., N A and i ∈ C B for i = N A +1, . . ., N. Without loss of generality we may assume that phase of nodes in C A is θ A = Ωt and nodes in C B have phase θ B = Ωt + χ where Ω is the collective frequency of the solution and χ is the phase difference between the clusters.The clusters will generally have different orbits (coinciding with differing isostables of the uncoupled nodes) such that for i ∈ C A , ψ i (t) = Ψ A and for i ∈ C B , ψ i (t) = Ψ B where Ψ A and Ψ B are constants.
Using (22), equations for the collective frequency, phase difference and cluster orbits can be determined.The four equations defining Ω, χ, Ψ A and Ψ B are Equations ( 56b) and (56d) can be solved for Ψ A and Ψ B in terms of functions of χ.Substituting the result into (56a) and (56c) we see that χ is a root of a nonlinear periodic function in terms of functions of χ.The roots can be determined, typically using numerical root finding schemes.For given values of N A , N B and ε there can be many possible solutions (Ω, χ, Ψ A , Ψ B ). Synchrony, (Ω, 0, Ψ, Ψ), is a solution for all values of N A , N B and ε where Ψ and Ω are given by ( 31) and ( 32) respectively with c i = 1 for all i.However we are interested in states with χ ̸ = 0. Linearising around a a two-cluster state we find that the Jacobian J is of the form (29) where Φ = (0, . . ., 0, χ, . . .χ).
Here each of the N × N block matrices H (p) (Φ), p = 1, . . ., 4 themselves have a block structure where, letting 1 N×M denote the N × M matrix of ones, The coefficients a p , b p , c p , d p , s A p and s B p are given in Appendix B where they are arranged into 2 × 2 matrices A , B, C , D, S A and S B .The block structure of J facilitates the computation of its eigenvalues.As shown in Appendix B the eigenvalues corresponding to intracluster perturbations are found to be the eigenvalues, µ ± A , of S A each with multiplicity N A − 1 and the eigenvalues µ ± B of S B each with multiplicity N B − 1.For stability with respect to intracluster perturbations it is required that Trace(S A ) < 0, det(S A ) > 0, Trace(S B ) < 0 and det(S B ) > 0. The remaining 4 eigenvalues correspond to intercluster perturbations and are the eigenvalues of the 4 × 4 block matrix Row and column operations show that one of the eigenvalues of J M is the expected zero eigenvalue due to purely rotational eigenmodes.The others, corresponding to changes in the phase difference between clusters, and changes in the two isostable values Ψ A , Ψ B , are required to have negative real part for stability of the cluster state.They satisfy a cubic characteristic polynomial µ 3 + q 1 µ 2 + q 2 µ 3 + q 3 = 0.By Routh's stability criterion, all roots of this cubic polynomial lie in the left half-plane when q 1 > 0, q 2 > 0, q 3 > 0 and q 1 q 2 > q 3 .Reducing to the phase-only case, we recover known results for phase-reduced networks 44,47,59 .

IV. A COMPARISON OF HIGHER ORDER PHASE AND PHASE-ISOSTABLE REDUCTIONS
In the previous section we have used the phase-isostable network equations (20) to determine with greater accuracy conditions for the existence and stability of phase-locked states.Other authors have used an alternative approach based on isostable coordinates to describe network phenomena which cannot be observed using first-order phase reduction.Rather than work with the phase and isostable equations, it is possible to use isostable coordinates to obtain higher order phase reduction.The process to obtain this reduction, given by Park and Wilson in 1 , involves expanding the isostable coordinate as ψ(t) = ε p (1) . A hierarchy of linear differential equations can be found and solved for the p (k) (t) to any order.Substituting back into the phase equation gives the higher order phase reduced equation including terms up to O(ε 3 ) and after approximation using first-order averaging as where H 1 , . . ., H 8 are defined in Appendix C. Wilson and Ermentrout 32 previously explicitly derived these equations up to order ε 2 , while Park and Wilson 1 indicate how the equations may be derived to any order, giving explicit equations to order ε 3 .However, the form of the higher order phase reduction given in 1 differs from (60).In line with other work on higher order phase equations 2,37 we observe in (60) that for networks of more than two nodes the higher order terms have non-pairwise phase interactions (i.e., terms involving the phases of three or more oscillators), despite the interactions between the nonlinear oscillators in (1) being pairwise.Park and Wilson 1 however arrive at a higher order phase reduced equation with pairwise interaction terms which we find to be erroneous.We give details of our calculations leading to (60)  in Appendix C.More recently Park and Wilson 35 have updated their calculations allowing for non-identical oscillators.We therefore see that there are two related but different frameworks, both using isostable coordinates, which extend standard first-order phase reduction; the phase-isostable network equations (20) and the higher order phase reduction (60).An obvious question is which approach can most accurately capture the dynamics of the full system (1)?In this section we answer this question for the archetypal example of the meanfield complex Ginzburg-Landau equation where the linear stability boundaries for synchrony and the splay state are known and may be compared with approximations of the boundaries found using the phase-isostable network equations (20) and higher order phase equations (60).

A. The mean-field complex Ginzburg-Landau equation
The normal form of the Hopf bifurcation (or Stuart-Landau oscillator) is a ubiquitous example of an oscillatory system.When globally diffusively coupled the result is the mean-field complex Ginzburg-Landau equation (MF-CGLE).In real coordinates the governing equations are of the form (1) where and c 1 and c 2 are two real-valued parameters.When ε = 0 each uncoupled node has a stable limit cycle x γ (t) = (cos(c 2 t), − sin(c 2 t)) T with period T = 2π/c 2 and Floquet exponent κ = −2.We define the phase θ j on the limit cycle as θ j = c 2 t + φ j so that θ j ∈ [0, 2π) and ω = c 2 .In the full system (61), the synchronous state, where x j = x γ (t) for j = 1, . . ., N, exists for all parameter values.The splay state is given by x j = (r cos(ϕ j ), r sin(ϕ j )) T where r = √ 1 − ε and ϕ j = (−c 2 + ε(c 2 − c 1 ))t + 2π j/N and exists when ε < 1.The linear stability analysis of these solutions gives closed formulae for the marginal stability of synchrony (denoted ε s ) and the splay state (denoted ε 0 when N ≥ 3 and ε a when N = 2) for fixed values of c 1 and c 2 as 60,61 Here we compute and compare the approximations to the marginal stability curves for these states using both the phaseisostable network equations ( 20) and the higher order phase reduction (60) in order to compare their accuracy to the true curves.For completeness the results are also compared with those obtained by León and Pazó 2 who obtained an alternative higher order phase reduction based on the isochrons which are known in closed form for this system 10,23 .We note that both methods based on isostable reduction will generalise to any system, however the approach of 2 relies on having a closed formula for the isochrons.We choose to compare the phaseisostable network equations to first-order in ψ with the second and third order phase reductions since these reductions require comparable computational effort in terms of calculating the required terms in the expansions (11).Using (12), normalising as appropriate and denoting r(t) = (cos(c 2 t), − sin(c 2 t)) T and φ(t) = (sin(c 2 t), cos(c 2 t)) T we obtain g (1) (t) = A(r(t) + c 2 φ(t)), (t) = −3r(t) + c 2 φ(t), where A = ∥g (1) (0)∥ = (1 + c 2 2 ) −1/2 .Calculating the expansions for the coupling function ( 16), we find that J 2 = −J 1 and since the coupling is linear, all Hessian terms H q vanish.

The phase-isostable network equation transformation
Considering first the transformation to phase-isostable network equations (20), we calculate that the averaged coupling functions as in ( 21) are From (34) with Ψ = 0, we find that synchrony is stable in the globally coupled network when ε > −1 and 2ε(1 ) > 0. The linear stability boundary for synchrony using the phaseisostable network equations, denoted ε s,PI is therefore identical to the closed formula for the full dynamics (62a) ε s,PI = ε s .
For the splay state (that exists only when ε < 1), we find that the collective orbit satisfies Ψ = ε/(2A(ε − 1)) and Ω = c 2 − ε(c 2 − c 1 ) from (37).It can be seen that this agrees precisely with the orbit in the full system (61) as it can be computed that for the MF-CGLE, isostables are circles of radius r(ψ) = k/(ψ + k) where I (0) = −2k r(t) 10,23 .Since we have chosen to normalise such that k = −(2A) −1 we find that the orbit for the splay state in the phase-isostable transformed system is r(Ψ) = √ 1 − ε.Using the results of section III A 2 we can compute the conditions for linear stability of the splay state.When N = 2 (so the splay state is the antisynchronous solution) the eigenvalues of the Jacobian are 0, −2(1 − ε) and the eigenvalues of Λ 1 given by (45).Since the matrix Λ 1 has real entries in this case, the eigenvalues have negative real part when Tr(Λ 1 ) < 0 and Det(Λ 1 ) > 0. This gives stability boundaries ε = 1/2, ε = 0 and ε = ε a,PI where ε a,PI satisfies When N ≥ 3 and finite and also in the large N limit the non trivial marginal linear stability boundary for the splay state, denoted ε = ε 0,PI , can be shown to be given by (see Appendix D), (65)

The higher-order phase reduction
The synchronous solution of ( 60) with θ i = θ for i = 1, . . ., N in the case of global coupling w i j = 1/N is guaranteed to exist.The linearisation about the synchronous state has the Jacobian where L is the graph Laplacian for global coupling given by L i j = −1/N + δ i j .This has an (N − 1 degenerate) eigenvalue of +1 and therefore synchrony is stable in the higher-order phase reduction if ξ > 0. Calculating each of the functions H 1 , . . ., H 8 for the MF-CGLE as described in Appendix C we observe that Therefore, keeping only the terms up to ε 2 we obtain the linear stability boundary for the second-order phase reduction as which matches identically with the result from the second order phase reduction of León and Pazó 2 whose reduction is based on the isochrons for the MF-CGLE rather than the isostable coordinates used here.The third-order phase reduction gives the values for marginal linear stability for fixed values of c 1 and c 2 as ε s,3 satisfying whereas the third-order reduction of León and Pazó 2 gives the value as ε * s,3 satisfying Considering next the stability of the splay state for N ≥ 3, we find that the Jacobian has the form where and ϒ = ∑ N k=1 A k .The circulant structure allows us to determine that the eigenvalues of J are λ p = ν p − ϒ where ν p = ∑ N j=1 A j exp(2πi j p/N) for p = 0, . . ., N − 1.Then Since λ 0 = 0 the splay state is stable if Re(λ p ) < 0 for all p ̸ = 0.When N ≥ 3, the second-order phase reduction gives the boundary for linear stability of the splay state as which again agrees with the second-order phase result of 2 .
To third-order, we find the approximation to the boundary is given by which differs from the result of 2 , who find the approximation When N = 2 we find that for the antisynchronous solution the stability boundaries for the second and third order phase reductions are respectively

Comparison of approaches
In Figure 1 we illustrate for comparison the stability boundary for each of the approaches (exact, phase-isostable network equations, phase interactions of first, second and third order and the approximations of 2 ) of each phase-locked state (synchrony, antisynchrony when N = 2, splay state for N ≥ 3) when ε is positive.We choose to fix values of c 2 (the intrinsic parameter for the node dynamics determining the angular velocity) and plot the marginal stability curves in the (c 1 , ε) plane, following 2,61 .Figure 1 shows the boundaries for c 2 = 0.5 (left), c 2 = 1.1 (centre) and c 2 = 3 (right).Note that we avoid making the choice c 2 = 1 since in this case (74) and (75) both have a factor of 1 + c 1 giving c 1 = −1 as a stability boundary for both higher-order phase approximations.
For synchrony we have already observed exact agreement of the stability boundary for the MF-CGLE and the phaseisostable network approximation.The top row of Figure 1 shows that the third order phase approximation (ε s,3 in dashed purple), provides a slightly better approximation than the third order result of León and Pazó (ε * s,3 , in dotted green) 2 , in a neighbourhood of (c 1 , ε) = (−1/c 2 , 0).
For the antisynchronous solution when N = 2 (middle row of Figure 1) we see that phase-isostable approximation (ε a,PI of (64), shown in dash-dotted light blue) most closely matches the curve for the MF-CGLE over a range of values of c 1 > −1/c 2 and ε > 0.
For the splay state where N ≥ 3 (bottom row of Figure 1), while ε * 0,3 (dotted green) may provide the closest approximation to the curve for the MF-CGLE in a neighbourhood of (c 1 , ε) = (−1/c 2 , 0) we observe that this curve (and those for higher-order phase shown by dotted curves) blow-up for values of c 1 away from −1/c 2 , while the curve for the phaseisostable approximation ε 0,PI , shown in dash-dotted light blue, is the only one which has a similar shape to the curve for the MF-CGLE, providing a bifurcation locus for all values of c 1 < −1/c 2 .Therefore the phase-isostable approximation provides the most accurate description of the qualitative behaviour (i.e., the existence of a bifurcation) over a range of values of c 1 and moderate values of ε although it does not accurately predict the value of ε at which the bifurcation will occur, underestimating in all cases shown in Figure 1.
We further note that in the MF-CGLE, below the critical value of c 2 = √ 3 at which the boundaries ε s and ε 0 become tangent at ε = 0, there are regions of bistability between synchrony and the splay state 2 , as shown in Figure 2A for Both the full model and the phase-isostable approximation have a parameter region where neither synchrony nor the splay state are stable when c 2 = 1.1.For c 1 = −2 we compare the stable behaviour in this region for the full model with that predicted by the phase-isostable approximation for a network of N = 3 nodes.The bifurcation diagrams are as in Figure 3.We observe that the phase-isostable approximation captures the bifurcations of synchrony and the splay state (UIS) at moderate values of ε which cannot be seen with a first-order phase reduction.Furthermore, the phase-isostable approximation correctly shows the loss of stability of synchrony to a periodic two-cluster state as ε decreases through 0.48.It also correctly predicts that as ε increases the splay state losses stability at a Hopf bifurcation to a quasiperiodic nonuniform incoherent state (NUIS) where the three nodes all have different trajectories.Again we observe that the phase-isostable approximation captures the qualitative dynamics, despite discrepancies in the precise values of ε for which some of the bifurcations occur and the precise branching structure.For instance in the full model the Hopf bifurcation of the splay state occurs at ε = 0.393 which in the full model this occurs at ε = 0.444.The phase-isostable approximation shows the stable quasiperiodic NUIS bifurcating supercritically directly from the splay state, while in reality the Hopf bifurcation is subcritical and the stable quasiperiodic NUIS bifurcates from a periodic two-cluster state.The phase-isostable approximation also predicts increasingly complex branching of unstable solutions for increasing ε > 0.6 is not found in the full model.However, for values of ε this large the phase-isostable approximation is not expected to be valid.
Park and Wilson 35 show analysis of a similar three node network of diffusively coupled complex Ginzburg-Landau models using the higher-order phase reduction.They also find the loss of stability of the splay state at a Hopf bifurcation to a quasiperiodic NUIS with the third-order phase-reduction though they do this through numerical simulation and do not provide any bifurcation diagrams.Our analysis highlights that the interaction functions can be calculated analytically and our computations of stability conditions for phase-locked states in section III allows for explicit location of bifurcations.
We conclude that the phase-isostable network transforma-FIG.2. Stability regions for synchrony (blue) and splay state (red) for (A) the MF-CGLE, (B) the phase-isostable approximation, (C) the third-order phase approximation when c 2 = 1.1.The phaseisostable network equations (B) can qualitatively reproduce the regions of bistability observed for the full system (A), but higher-order phase approximations cannot (C).
tion better captures the qualitative bifurcation structure of the MF-CGLE when compared with second and third order phase reductions.Since it is the normal form for globally diffusively coupled Hopf normal form oscillators, it is natural to suppose that this framework is the most appropriate to use to study linearly coupled networks of oscillatory nodes near Hopf bifurcation.
It is noted by Park and Wilson 1 that the higher-order phase approximation to the analytical stability boundaries for synchrony and the antisynchronous solution for a network of N = 2 nodes does improve with increasing order.At 10th order good quantitative and qualitative agreement of the stability boundaries is achieved for moderate values of ε, however this does come at huge computational expense in terms of determining all required terms in the PRC and IRC expansions.Using the phase-isostable network equation approximation to first order may be an acceptable compromise, to be able to indicate the qualitative network behaviour using a simpler lower order approximation.

V. NETWORKS OF MORRIS-LECAR NEURONS
The comparison of higher-order phase reductions and the phase-isostable network equations for the MF-CGLE in section IV has shown that the phase-isostable network equation approach is able to predict the qualitative network behaviour for moderate values of interaction strength ε, which the higher-order phase reductions of similar computation expense cannot.We now use the phase-isostable framework to reveal network behaviours in small and large networks of planar Morris-Lecar neurons 39 which cannot be described using traditional first-order phase reduction.The Morris-Lecar model is planar model of neuronal excitability in which oscillations can occur through Hopf, SNIC or homoclinic bifurcations.The equations describing the node dynamics and our choices of parameter values are given in Appendix E. Our parameter choices place the node dynamics near the homoclinic bifurcation.The stable periodic orbit has κ = −0.4094and is indicated in Figure 4 along with isochrons and selected isostables normalised so that negative values of the isostable coordinate correspond to points inside the limit cycle.
The Floquet eigenfunction g (1) , and the iPRC and iIRC, Z (0) and I (0) and their first order correction terms Z (1) and I (1) are numerically determined as appropriately normalised solutions of ( 12) and the first (voltage) component of the response vectors are depicted together with the limit cycle in Figure 5A.
We consider networks diffusively coupled through the voltage variable v so that x i = (v i , w i ) T and G(x i , x j ) = (v j − v i , 0) T and w i j = 1/N in (1). Figure 5B shows the resulting six interaction functions H 1 , . . ., H 6 given by (21).In section V 1 we compare the dynamics of the phase-isostable approximation with that of the full model for two coupled nodes and consider a larger network of 200 nodes in section V 2.

Two coupled nodes
For a network of two nodes, bifurcation diagrams using the interaction strength ε as the bifurcation parameter are computed for both the full model (Figure 6A) and the phaseisostable network equations (20) (Figure 6B).For the full model the bifurcation diagram, Figure 6A, shows the maximum of the voltage variable for one of the two nodes on the vertical axis for relatively small values of the coupling strength ε.In the corresponding diagram for the phaseisostable approximation (Figure 6B) the vertical axis shows χ = θ 2 − θ 1 , the phase difference between the two nodes, so that periodic orbits of the full model correspond to equilibrium points of the phase-isostable reduced system.Since we choose to work with phase-isostable network equations ( 20) that are linear in the isostable coordinates, we find unique branches of synchronous and antisynchronous solutions with the isostable coordinate given by ( 31) and (37) respectively in Figure 6B.These branches also exist for all values of ε, except where the value of the isostable coordinate asymptotes to positive or negative infinity.There may be more than one phase-locked state with χ ̸ = 0, π for a given value of ε depending on the interaction functions H 1 , . . ., H 6 .The existence, stability of solution branches and locations of bifurcations in Figure 6B all agree with the explicit calculations in Section III.Taking Figure 6A as ground truth, we here comment on the extent to which the phase-isostable framework approximates the dynamics observed in the full system as the strength ε of the linear coupling increases.61) with the angular difference between two nodes on the vertical axis.(B) shows the approximation using the phaseisostable network equation framework (20) with the phase difference θ 2 − θ 1 on the vertical axis.Solid (dotted) lines indicate stable (unstable) solutions.Shown are two of the three symmetric branches of each type.We do not include the branches with θ 2 = θ 1 since these would obscure the synchronous branch.Parameter values are c 1 = −2 and c 2 = 1.1.The bifurcation diagrams are produced using the implementation of AUTO in XPPAUT 62 and bifurcation locations in (B) are in agreement with stability computations is section III.FIG. 4. Phase portrait for the Morris-Lecar model (E1) indicating the stable limit cycle (in white).Parameter values are as given in Appendix E. These parameters place the dynamics near the homoclinic bifurcation so that solutions on the limit cycle spend a significant portion of the period near the saddle point.Also shown are the boundary of the basin of attraction of the limit cycle, a selection of isostables and colours indicate isochrons connecting points with the same asymptotic phase.These are computed numerically using the constructive definition ( 9) and the method of Fourier averages 63 respectively.
Considering first the synchronous state (in light red) we see that in the full model this solution is unstable in the weak coupling regime for positive ε, but restabilises at ε = 0.499 where it meets the branch of phase-locked solutions (shown in purple).Synchrony is also stable for a small interval of negative ε, before undergoing a period doubling bifurcation at ε = −0.015.The phase-isostable approximation, as shown in Figure 6B, correctly predicts the stability type of synchrony in the neighbourhood of ε = 0 (as does a first-order phase approximation).The phase-isostable approximation is additionally able to capture bifurcations of synchrony for both positive FIG. 5. (A) The first (voltage) components of response vectors Z (0) , I (0) , Z (1) and I (1) together with the voltage component v(θ ) of the limit cycle, and (B) the interaction functions H 1 , . . ., H 6 given by ( 21) for the Morris-Lecar model (E1) with parameter values as in Appendix E and linear coupling.and negative ε, showing synchrony gaining stability where it meets a phase-locked state at ε = 0.0934 and losing stability at a Hopf bifurcation at ε = −0.407.We note that while this reproduces the qualitative behaviour of the synchronous FIG. 6. Bifurcation diagrams for two linearly coupled Morris-Lecar neurons with interaction strength ε as the bifurcation parameter.(A) shows the full model equations (E1) with the maximum of the voltage variable for one of the two nodes on the vertical axis.(B) shows the approximations using the phase-isostable network equation framework (20) with the phase difference χ = θ 2 − θ 1 on the vertical axis.Stable (unstable) periodic solutions are shown by solid (dashed) lines with synchrony in light red, antisynchrony in blue, phase locked states with χ ̸ = 0, π in purple and other periodic states in green ((A) only).In (B), stable (unstable) quasiperiodic solutions are indicated in yellow (dark red).Black circles/stars/diamonds indicate period doubling/torus/Hopf bifurcations respectively (noting that periodic solutions are given by equilibrium points of the phase-isostable approximation when working with the phase difference χ.) Black squares show limit points where the modulus of an isostable coordinate asymptotes to ±∞.The phase-isostable approximation captures many of the qualitative features of the full model dynamics as discussed in the text.Parameter values for the node dyanamics for both diagrams are as given in Appendix E. The bifurcation diagrams are produced using the implementation of AUTO in XPPAUT 62 .
branch, the phase-isostable approximation underestimates the precise value of ε at which the bifurcations occur, just as for the MF-CGLE in section IV A.
In the full system antisynchrony (shown in blue in Figure 6) is unstable for all negative values of ε, but is stable for 0 < ε < 0.0831 where the orbit lies inside, but close to the synchronous orbit γ.For 0.00491 < ε < 0.0831 the stable large amplitude antisynchronous solution co-exists with an unstable smaller amplitude antisynchronous solution.The qualitative stability of the branch of antisynchronous solutions of the full model with orbit closest to the uncoupled node orbit is also captured by the phase-isostable reduction, including the fact that it exists for ε < ε ∞ where the phase isostable framework estimates that ε ∞ = 0.0961 (compared to 0.0831 in the full system).This point is labelled as a limit point (black square) in Figure 6.As ε → ε ∞ from below the stable solution branch has Ψ → −∞ corresponding to a shrinking amplitude orbit in the (v, w) coordinates approaching the unstable fixed point inside the uncoupled node limit cycle.We also find a branch of unstable antisynchronous solutions for ε > ε ∞ with Ψ positive.As ε → ε ∞ from above this branch has Ψ → ∞ (corresponding to a growing amplitude orbit approaching the outer edge of the basin of attraction of the uncoupled node limit cycle, see Figure 4).There is no bifurcation at ε ∞ ; the change in stability here is coincidental.We have not found a corresponding antisynchronous state with orbit outside of γ in the full model since it may always be unstable.The phase-isostable reduction at the given order does not find the unstable small amplitude antisynchronous oscillations since, as previously noted, it gives a unique value of Ψ for each ε.It may be that if the equations (20) were taken to higher order in ε, then this other branch could be revealed.We note that the phase-isostable framework predicts a small region (0.0934 < ε < ε ∞ ) where both synchrony and antisynchrony are stable.This does not occur in the full system.
In the full coupled system there are also periodic solutions for which the nodes do not share a common orbit and therefore for each of these there is a symmetric solution under x 1 ↔ x 2 .For clarity, Figure 6A shows only the branch for which v 1 has the largest maximum value.A phase-locked solution with x 1 performing larger and x 2 smaller amplitude oscillations inside the synchronous orbit exists for −0.0218 < ε < 0.575 (shown in purple in Figure 6).Its stability varies along the branch as it undergoes a series of saddle node and torus bifurcations.Many of these occur within the region near ε = 0, but most obvious in Figure 6A are the torus bifurcations at ε = 0.0472, ε = 0.148 and ε = 0.538 marked by black stars.There is also an isola of periodic solutions for 0.375 < ε < 0.403 an example of which is depicted in Figure 7B,D.This solution is not shown in Figure 6A to avoid confusion with the synchronous branch since they share similar values of max(v 1 ).For values of ε ∈ (0.148, 0.375) ∪ (0.403, 0.499) numerical simulations indicate that the stable solutions are quasiperiodic with oscillations fluctuating around the phase-locked state, as indicated in Figure 7A,C.In the phase-isostable approximation, branches of stable phase-locked states with χ ̸ = 0, π bifurcate from synchrony at ε = 0.0934 losing stability at a Hopf bifurcation at ε = 0.0484.This gives rise to stable solutions where χ, Ψ 1 , Ψ 2 are all periodic corresponding to behaviour observed in numerical simulations of the full model (Figure 7), however, here these solutions exist over a much smaller interval of ε values.We also note that the stable phase-locked state has isostable coordinates of opposite sign for the two nodes corresponding to one node orbiting outside of γ and the other inside.As the limit point at ε = 0.0339 is approached, one of the isostable coordinates asymptotes to positive infinity while the other asymptotes to negative infinity.Beyond the limit point there is another (unstable) branch of phase-locked states with the signs of the isostable coordinates reversed.This branch of solutions persists for larger values of ε, but is unstable and always has one node with a large isostable coordinate.
The branch of periodic solutions shown in green in Fig- ure 6A is created at a homoclinic bifurcation at ε = 0.0381 and meets the unstable fixed point at ε = −0.0618.This solution has x 1 oscillating near the synchronous orbit γ, while x 2 performs very small amplitude oscillations near the fixed point outside of the basin of attraction of γ.Since phase-isostable coordinates are valid only within the basin of attraction of γ, the phase-isostable network equations are not able to capture this oscillatory solution.
We conclude that the phase-isostable approximation is able to capture the qualitative stability of synchrony and antisynchrony away from ε = 0. Furthermore, it is able to reveal phase-locked states and quasiperiodic solutions, capturing far more of the full dynamics than the first-order phase reduction of the model.Since approximation using phase-isostable coordinates can describe qualitative dynamics which first-order phase reduction cannot for a two node network, we now use the approximation to investigate the dynamics of larger networks of Morris-Lecar models in section V 2. We note that, while the phase-isostable framework can be instructive concerning network behaviours for relatively weak coupling, since the expansions rely on the assumption that ψ = O(ε) predictions for larger values of ε must be interpreted with caution.

Networks of many Morris-Lecar neurons
We next consider a network of 200 diffusively coupled Morris-Lecar neurons.We again use the parameters values in Appendix E so that the response and interaction functions remain as in Figure 5.Using ( 36) and ( 45) we calculate the stability of the synchronous and splay states respectively over a range of small positive and negative coupling strength ε.The results are shown in Figure 8.The maximal real part of the eigenvalues for synchrony is positive (and therefore the state is unstable) when 0 < ε < 0.0934 and the splay state is unstable for all values of ε ̸ = 0.The discontinuity in the maximum real part of the eigenvalues for the splay state at ε = 0.0664 corresponds to the discontinuity in the isostable coordinate of the splay state when the denominator in (37) becomes zero.The splay state orbit lies inside the node orbit (Ψ < 0) when ε < 0.0664 and for ε > 0.0664 the splay state has Ψ > 0.
We take ε = 0.065 which lies in the range where both synchrony and the splay state in the phase-isostable approximated framework are unstable and numerically simulate the phase-isostable network equations to investigate the stable behaviour predicted in this region.We initialise the system with the 200 nodes in a tightly packed group with phase coordinates θ i ∈ [0.283725, 0.283735] and isostable coordinates ψ i ∈ [2.9794, 2.9798].Figure 9 shows the evolution in time of the isostable coordinates.We see the group of nodes initially move towards the isostable ψ = 0 and remains near the synchronous network state for some time.The group then rapidly desynchronises before settling into a stable 2-cluster state.The analysis in section III B confirm the existence of the observed cluster state with N A = 28 and N B = 172 and also that it has phase difference χ = 2.1407 between the clusters which have orbits coinciding with the Ψ A = −0.1694and Ψ B = −0.1868isostables, also in agreement with the simulation.The stability analysis further confirms that this state is stable for ε = 0.065.
For the same value of the coupling strength ε = 0.065 we also simulate the full network equations from an initial state with (v i , w i ) in a group close to (v, w) ≈ (-0.1,0.07) which has phase and isostable coordinates (θ , ψ) = (0.28373, 2.9796).FIG. 9.The isostable coordinates of the 200 nodes in the phaseisostable approximation of the diffusively coupled network of Morris-Lecar neurons evolving with time.The initially near synchronous group moves towards the synchronous isostable where ψ i = 0 for all i before desynchronising at t ≈ 80 and settling into a stable two cluster state from t ≈ 130.In this cluster state, there are 28 nodes in the cluster with ψ = −0.1694and 172 nodes in the cluster with ψ = −0.1868.Interaction functions are as in Figure 5 and ε = 0.065.
Snapshots from the simulation are shown in Figure 10, which also indicates the isostable ψ = 0 which coincides with the synchronous orbit and the isostable corresponding to the splay state for these parameter values (ψ = −19.3) in the (v, w) space.Initially the nodes are close to synchronous and are first drawn towards and orbit near to ψ = 0.However, since this state is unstable the cluster soon breaks up and the nodes spread out and oscillate nearer to the splay state isostable.This state is also unstable in the full network and therefore the system oscillates between the splay and synchronous states for a time before settling at a stable 3-cluster state.The simulations of the phase-isostable approximation and the full model for equivalent initial conditions show the same progression from near synchrony, subsequent desynchronisation eventually settling on a cluster state.
For stronger diffusive coupling (ε = 0.4) with the node parameters still as in Appendix E, Han et al. 40 observed that in the 200 node network initialised at a near synchronous state, the nodes would at first behave as we have observed at lower coupling strength, moving to the synchronous orbit before desynchronising.However at this larger coupling strength it was observed that the desynchronised nodes spiral in towards the unstable node inside the limit cycle, before again moving out towards the synchronous orbit.This behaviour repeats resulting in mean field voltage traces showing large amplitude fluctuations.Unfortunately this type of behaviour cannot be captured using the phase-isostable network equations at the order of (20) as the regime where ε = 0.4 appears to be beyond the scope of its predictive power.Nonetheless, we have been able to demonstrate that at smaller coupling strengths the phase-isostable framework is able to capture the formation of stable cluster states for parameter values where both the synchronous and splay states are unstable in agreement with simulations of the full network equations.

A. Remarks
In section III we observed that in the phase-isostable network equations truncated at linear order in the isostable coordinates (or equivalently O(ε 2 )) (20), the orbits for individual nodes all coincide with isostables since they have constant isostable coordinate in the most slowly decaying direction (which is the only direction we retain here).When the node dynamics are two-dimensional the phase-isostable coordinates are a transformation from the original coordinates of the node model.That is, for each x ∈ B(γ) there is a unique corresponding (θ , ψ).Therefore, since trajectories of the node dynamics cannot intersect, isostables also have no self intersections.However, for periodic orbits of the full network dynamics the projection of the dynamics of each node onto the node phase space can have shared orbits which intersect.For example, Figure 6A indicates a period doubling bifurcation of synchrony at ε = −0.015.The period doubled branch (not shown in Figure 6A) at ε = −0.0129 is a periodic two-cluster phase-locked state which is bistable with synchrony.The shared node orbit is shown in Figure 11A and we observe that this has a self-intersection.The linear order truncation of the phase-isostable network equations used here  6 A. Since this periodic orbit projected into the two-dimensional node phase space has a self intersection, this orbit cannot be described by the linear phase-isostable network equations whose phase-locked solutions all coincide with isostables which do not have self intersections.(B) shows the corresponding time series for v 1 , v 2 over a single period.cannot capture such solutions, but it is possible that a higher order expansion would be capable of recovering such orbits.Higher order truncation in the network equations ( 20) would also likely improve the quantitative accuracy of the predictions of the phase-isostable dynamics with the dynamics of the full system.

VI. CONCLUSIONS
There are two emerging frameworks providing improved understanding of the behaviour of coupled oscillator networks through the use of isostable coordinates to capture dynamics off-cycle.While many authors are choosing to use isostable dynamics to obtain higher-order phase reductions 1,32,35,36 we instead focus on the approach which retains the isostable dynamics in the slowest decaying direction.We have obtained existence and stability criteria for phase-locked states in networks of identical nodes and used these to demonstrate that for the MF-CGLE the phase-isostable network equations have far superior capabilities at capturing the qualitative shape of the stability boundaries of the synchronous and splay states in networks of arbitrary size when compared with phasereduction at up to cubic order in the interaction strength.For a network of size N = 3 the phase-isostable framework can also correctly identify the stable solution types beyond the loss of stability of synchrony and the splay state.We then use the framework to study networks of 2 and 200 diffusively coupled Morris-Lecar neurons.For the smaller network, comparison of bifurcation diagrams for both the full and phase-isostable approximations of the dynamics show that the phase-isostable framework can capture qualitative changes in stability of synchrony and the antisynchrnous state as well as the existence of other phase-locked states and quasiperiodic solutions which cannot be described by first-order phase reduction.However, it is not able to accurately predict the values of the coupling strength at which bifurcations occur, tending to underestimate these values.For a network of 200 nodes, for values of the coupling strength where both synchrony and the splay state are unstable in the full system and the phase-isostable approximation, the phase isostable network equations predict stable cluster states which are also seen in numerical simulations of the full equations in this regime.The phase-isostable network approximation (or reduction in the case of node dynamics of dimension three or higher) appears to be a useful tool to indicate qualitative network behaviour beyond that which can be revealed using first-order phase reduction whilst also keeping the computational complexity low since only six interaction functions need to be computed, requiring only four response functions.
For expansions of comparable order and computational effort, we have observed that retaining the notion of distance from cycle through the isostable dynamics outperforms using the isostable dynamics to refine phase-reduced dynamics in locating bifurcations of phase-locked solutions.Furthermore it allows for the analysis of off cycle dynamics through the tracking of the isostable coordinate.In section IV A we found that the phase-isostable network transformation best captures the qualitative bifurcation structure of the MF-CGLE.Since it is the normal form for globally diffusively coupled Hopf normal form oscillators, it is natural to suppose that this framework is the most appropriate to use to study linearly coupled networks of oscillatory nodes near Hopf bifurcation.We do however note that recent work of Bick et al. 64 generalises the techniques of 2 to consider Stuart-Landau like systems with phase-dependent amplitude emulating the deformed limit cycles expected away from Hopf bifurcation.They show that a second-order phase reduction is able to accurately predict the stability properties of the synchronized and splay orbits when all terms of up to second order in the coupling strength are included.It would be interesting to compare with corresponding results using a phase-isostable approach.
Although we have chosen to use the orbit of an uncoupled oscillatory node as the reference for setting up the phaseamplitude coordinate system, other choices are possible.This has been recognised by Wilson as a way to to handle perturbations that take one far from cycle with the development of an adaptive reduction strategy to construct a more useful reference limit cycle (so that distances to this reference cycle remain small) 65 .Alternatively, one could more simply take a phase-isostable reduction about a stable periodic network state (e.g.synchrony or another phase-locked state) 29 .This approach has been used to investigate control strategies for desynchronising populations of neural oscillators under periodic stimulation 28,66 .
In the present work we have investigated finite size networks with linear (state-dependent) coupling and twodimensional node dynamics.There are obvious extensions to this work, incorporating other forms of coupling such as event driven interactions, along the lines described for synaptic interactions e.g. in 67 , and also to investigate networks with higher dimensional node dynamics where retaining a single isostable coordinate dynamics in the most slowly decaying direction represents a dimension reduction for the system.The extension to treat dense graphs could naturally build on the work in 68,69 for Kuramoto networks using the notion of a graphon.Moreover, the analysis of continuum models with non-local interactions using phase-amplitude coordinates seems feasible by generalising the approach in 32 developed for reaction-diffusion equations.Finally, let us mention the challenge of dealing with delays.At the node level these can induce oscillations, and at the network level can strongly influence patterns of phase-locked states and their bifurcations.A method for constructing the infinitesimal phase response for delay induced oscillations has already been developed in 70,71 , and it would be interesting to use this functional setting to develop the corresponding amplitude response, as well as to incorporate delayed interactions within a phase-amplitude network setting.All of the above are topics of ongoing work and will be reported upon elsewhere.(1) i h 5 (θ i , θ j ) + p (1) j h 6 (θ i , θ j ) ,

FIG. 3 .
FIG.3.Bifurcation diagrams for the MF-CGLE when N = 3 with interaction strength ε as the bifurcation parameter.(A) shows the result for the full model equations(61) with the angular difference between two nodes on the vertical axis.(B) shows the approximation using the phaseisostable network equation framework(20) with the phase difference θ 2 − θ 1 on the vertical axis.Solid (dotted) lines indicate stable (unstable) solutions.Shown are two of the three symmetric branches of each type.We do not include the branches with θ 2 = θ 1 since these would obscure the synchronous branch.Parameter values are c 1 = −2 and c 2 = 1.1.The bifurcation diagrams are produced using the implementation of AUTO in XPPAUT62 and bifurcation locations in (B) are in agreement with stability computations is section III.

FIG. 7 .
FIG. 7. Direct numerical simulations of two linearly coupled Morris Lecar neurons for values of ε where all solutions indicated in Figure 6A are unstable.(A,C) ε = 0.25, shows quasiperiodic behaviour for both nodes, as is typical for ε ∈ (0.148, 0.375) ∪ (0.403, 0.499) (B,D) ε = 0.39, shows periodic behaviour for both nodes, as is typical for ε ∈ (0.375, 0.403).Panels (A,B) show network activity in the (v, w) plane (synchronous orbit γ is indicated in black), while the corresponding time series for v 1 , v 2 are shown in (C,D).In D a single period is shown.Node parameters are as in Appendix E.

FIG. 8 .
FIG. 8.The maximal real part of eigenvalues of phase-locked states in a 200 node, globally-coupled, phase-isostable Morris-Lecar network plotted against the coupling strength ε.The synchronous state is shown in black and the splay state is shown in red.

FIG. 10 .
FIG. 10.A series of snapshots of a simulation of the full network of 200 diffusively coupled Morris-Lecar neurons showing positions of each node (black dots), the uncoupled limit cycle at ψ = 0 (blue) and the splay state orbit at ψ = −19.3(red).Parameter values are as in Appendix E with ε = 0.065.

FIG. 11 .
FIG. 11. (A)A stable periodic phase-locked state in the network of two coupled Morris-Lecar neurons occurring at ε = −0.0129 on the solution branch bifurcating from the periodic doubling bifurcation indicated at ε = −0.015 in Figure6A.Since this periodic orbit projected into the two-dimensional node phase space has a self intersection, this orbit cannot be described by the linear phase-isostable network equations whose phase-locked solutions all coincide with isostables which do not have self intersections.(B) shows the corresponding time series for v 1 , v 2 over a single period.