Networks of coupled nonlinear oscillators can display a wide range of emergent behaviors under the 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 the 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 generalizing 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 behavior 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–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 limit-cycle 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 behaviors in weakly coupled oscillator networks in a variety of relevant systems.^{8,13–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 utilizes the notion of isostable coordinates.^{19–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 linearization 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–27} and also for the control of oscillators (see Ref. 28 for a recent review).

A dimension reduction can be achieved for high-dimensional 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–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–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 Refs. 32, 35, and 36 but overlooked in Ref. 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 Ref. 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} Similar to Refs. 1 and 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 non-smoothness of the dynamics using results from Ref. 30, a bifurcation diagram for the dynamics in the phase-isostable framework under variation of coupling strength reveals restabilization 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 underestimates the value of the coupling strength at which synchrony restabilizes.^{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 Sec. II, we derive the phase-isostable network differential equations and also the corresponding higher-order phase-reduced equations for a network of $N$ identical oscillators in Secs. II and IV, respectively. Using first-order averaging, we obtain a system of phase-isostable 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 Sec. III. We are then able in Sec. 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 Sec. 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 Ref. 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 Sec. 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

### A. First-order phase reduction and weakly coupled phase oscillators

^{3,4,11,12}The phase of each uncoupled node trajectory then evolves as $ d\theta / dt=\omega $ both on and off the limit cycle. Note that we can then use $ x \gamma (t)$ and $ x \gamma (\theta )$ to denote points on cycle where time is parameterized as $t=\theta /\omega $. In terms of phase variables, (1) becomes

^{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

^{16,41,42}then gives the phase dynamics back in the original variables as the simpler to analyze phase difference equations,

^{43–47}In Sec. 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,\u2026,6$, and, in Sec. 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(\chi )=sin\u2061(\chi )$ yet is able to capture the basic mechanisms underlying synchronization 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 ( $ |\epsilon |\u226a1$) 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

^{10,19–21,26,27,30}Isostables 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 $ \kappa m$, $m=1,\u2026,n\u22121$, a set of isostables representing an amplitude degree of freedom can be defined, and, therefore, for any point $ x\u2208 B(\gamma )$, we can associate isostable coordinates $ \psi m$, $m=1,\u2026,n\u22121$. 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 $ \psi 1$:

^{27,31}Denote by $ w T$ the left eigenvector of $M( x 0 \gamma )$ associated with the eigenvalue $ \lambda 1=exp\u2061( \kappa 1T)$. The corresponding right eigenvector is $v$ with $\Vert v\Vert = v \u22c5 v=1$ and $ w Tv=1$. Then,

An alternative perspective is taken by the authors of Refs. 19, 23, and 24 who note that for each point $ x\u2208 B(\gamma )$, there is a coordinate transform given by an analytic map $K$ such that $ x=K(\theta , \psi 1,\u2026 \psi n \u2212 1)$, and determine the map $K$ using a parameterization method.

^{1,27,32}and its validity is dependent on a separation of scale of the magnitude of $ \lambda 1$ from $ \lambda m$, $m=2,\u2026,n\u22121$. 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.

^{31}

^{,}

^{1,10}

^{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 Ref. 1. For $ X= [ x i T , x j T ] T\u2208 R 2 n$, we write $ G( X)= [ G 1 ( X ) \u2026 G n ( X ) ] T$, where $ G q( X)\u2208 R$. Defining $ X \gamma = [ x \gamma ( \theta i ) T , x \gamma ( \theta j ) T ] T$ and $\Delta X= [ \Delta x i T , \Delta x j T ] T$, the Taylor expansion of $ G$ can be expressed to second order as

^{33,42,53}yielding

## III. PHASE-LOCKED STATES IN PHASE-AMPLITUDE NETWORK EQUATIONS

Having defined a system of phase-amplitude network Eq. (20), we now investigate conditions for the existence and stability of certain phase-locked states within the averaged equations, generalizing well known results for phase-reduced networks.^{43–45,47,54}

^{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 $ \Psi i=\Psi $ for all $i=1,\u2026,N$. Such states include full synchrony ( $ \varphi i=0$ for all $i$), the splay (or uniform incoherent) state ( $ \varphi i=2\pi 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\pi /M$. In Sec. III B, we consider unbalanced two-cluster states, where each cluster is a different size and has a different orbit.

#### 1. Full synchrony

If the row sums $ c i$ are independent of $i$ (i.e., $ c i\u2261c$ for all $i$). This is the case for global coupling where $ w i j=1/N$ which we will consider further later.

If $ H 1(0)=0$ and $ H 4(0)=0$.

##### 1. Connectivity matrix with constant row sums

Therefore, synchrony is stable when $\kappa +\epsilon ( H 5(0)+ H 6(0))<0$, $ Trace( M(\Psi ))<0$ and $det( M(\Psi ))>0$. Reducing back to the phase-only description by taking $ H 2,\u2026, H 6\u22610$, we recover the result for phase oscillators that synchrony is stable if $\u2212\epsilon H 1 \u2032(0)<0$.^{43,47,54}

##### 2. Interaction functions with *H*_{1}(0) = 0 and *H*_{4}(0) = 0

#### 2. The splay state

##### 1. The splay state in the large *N* limit

^{43,45,47,54}

#### 3. Balanced cluster states in globally coupled networks

^{47,56,57}

^{58}the eigenvalues $ \lambda ( p )$ and $ \lambda q ( p )$ of $ H ( p )$ satisfy

### 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\u2264N/2$ and $ N B=N\u2212 N A$ nodes, respectively. We label nodes such that $i\u2208 C A$ for $i=1,\u2026, N A$ and $i\u2208 C B$ for $i= N A+1,\u2026,N$. Without loss of generality, we may assume that phase of nodes in $ C A$ is $ \theta A=\Omega t$ and nodes in $ C B$ have phase $ \theta B=\Omega t+\chi $, where $\Omega $ is the collective frequency of the solution and $\chi $ 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\u2208 C A$, $ \psi i(t)= \Psi A$ and for $i\u2208 C B$, $ \psi i(t)= \Psi B$, where $ \Psi A$ and $ \Psi B$ are constants.

^{44,47,59}

## IV. A COMPARISON OF HIGHER-ORDER PHASE AND PHASE-ISOSTABLE REDUCTIONS

^{32}previously explicitly derived these equations up to order $ \epsilon 2$, while Park and Wilson

^{1}indicate how the equations may be derived to any order, giving explicit equations to order $ \epsilon 3$. However, the form of the higher-order phase reduction given in Ref. 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 equation (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 mean-field 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 equation (20) and higher-order phase equation (60).

### A. The mean-field complex Ginzburg–Landau equation

^{60,61}

^{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 generalize to any system; however, the approach of Ref. 2 relies on having a closed formula for the isochrons. We choose to compare the phase-isostable network equations to first-order in $\psi $ 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).

#### 1. The phase-isostable network equation transformation

For the splay state (that exists only when $\epsilon <1$), we find that the collective orbit satisfies $ \Psi =\epsilon /(2A(\epsilon \u22121))$ and $\Omega = c 2\u2212\epsilon ( c 2\u2212 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(\psi )= k / ( \psi + k )$ where $ I ( 0 )=\u22122k r ^(t)$.^{10,23} Since we have chosen to normalize such that $k=\u2212 ( 2 A ) \u2212 1$, we find that the orbit for the splay state in the phase-isostable transformed system is $r( \Psi )= 1 \u2212 \epsilon $.

#### 2. The higher-order phase reduction

^{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 $ \epsilon s , 3$ satisfying

^{2}gives the value as $ \epsilon s , 3 \u2217$ satisfying

#### 3. Comparison of approaches

In Fig. 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 Ref. 2) of each phase-locked state (synchrony, antisynchrony when $N=2$, splay state for $N\u22653$) when $\epsilon $ 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,\epsilon )$ plane, following.^{2,61} Figure 1 shows the boundaries for $ c 2=0.5$ (left), $ c 2=1.1$ (center), 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=\u22121$ 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 phase-isostable network approximation. The top row of Fig. 1 shows that the third-order phase approximation ( $ \epsilon s , 3$ in dashed purple) provides a slightly better approximation than the third-order result of León and Pazó ( $ \epsilon s , 3 \u2217$, in dotted green),^{2} in a neighborhood of $( c 1,\epsilon )=(\u22121/ c 2,0)$.

For the antisynchronous solution when $N=2$ (middle row of Fig. 1), we see that phase-isostable approximation ( $ \epsilon a , P I$ of (64), shown in dashed-dotted light blue) most closely matches the curve for the MF-CGLE over a range of values of $ c 1>\u22121/ c 2$ and $\epsilon >0$.

For the splay state where $N\u22653$ (bottom row of Fig. 1), while $ \epsilon 0 , 3 \u2217$ (dotted green) may provide the closest approximation to the curve for the MF-CGLE in a neighborhood of $( c 1,\epsilon )=(\u22121/ 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 $\u22121/ c 2$, while the curve for the phase-isostable approximation $ \epsilon 0 , P I$, shown in dashed-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<\u22121/ c 2$. Therefore the phase-isostable approximation provides the most accurate description of the qualitative behavior (i.e., the existence of a bifurcation) over a range of values of $ c 1$ and moderate values of $\epsilon $ although it does not accurately predict the value of $\epsilon $ at which the bifurcation will occur, underestimating in all cases shown in Fig. 1.

We further note that in the MF-CGLE, below the critical value of $ c 2= 3$ at which the boundaries $ \epsilon s$ and $ \epsilon 0$ become tangent at $\epsilon =0$, there are regions of bistability between synchrony and the splay state,^{2} as shown in Fig. 2(a) for $ c 2=1.1$. The phase-isostable approximation can reproduce the bistability regions qualitatively [see Fig. 2(b)], while the phase-based approximations and those of Ref. 2 cannot. See Fig. 2(c) for the case of the third-order phase approximation when $ c 2=1.1$, which has very different regions of bistability.

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=\u22122$, we compare the stable behavior 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 Fig. 3. We observe that the phase-isostable approximation captures the bifurcations of synchrony and the splay state (UIS) at moderate values of $\epsilon $, 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 $\epsilon $ decreases through $0.48$. It also correctly predicts that as $\epsilon $ 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 $\epsilon $ 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 $\epsilon =0.393$, which in the full model this occurs at $\epsilon =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 $\epsilon >0.6$ is not found in the full model. However, for values of $\epsilon $, this large 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 Sec. III allow for the explicit location of bifurcations.

We conclude that the phase-isostable network transformation 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 tenth order, good quantitative and qualitative agreement of the stability boundaries is achieved for moderate values of $\epsilon $; 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 behavior 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 Sec. IV has shown that the phase-isostable network equation approach is able to predict the qualitative network behavior for moderate values of interaction strength $\epsilon $, which the higher-order phase reductions of similar computation expense cannot. We now use the phase-isostable framework to reveal network behaviors 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 the 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 $\kappa =\u22120.4094$ and is indicated in Fig. 4 along with isochrons and selected isostables normalized 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 normalized solutions of (12) and the first (voltage) component of the response vectors are depicted together with the limit cycle in Fig. 5(a).

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 \u2212 v i , 0 ) T$ and $ w i j=1/N$ in (1). Figure 5(b) shows the resulting six interaction functions $ H 1,\u2026, H 6$ given by (21). In Sec. V A, 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 Sec. V B.

### A. Two coupled nodes

For a network of two nodes, bifurcation diagrams using the interaction strength $\epsilon $ as the bifurcation parameter are computed for both the full model [Fig. 6(a)] and the phase-isostable network equation (20) [Fig. 6(b)]. For the full model the bifurcation diagram, Fig. 6(a) 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 $\epsilon $. In the corresponding diagram for the phase-isostable approximation [Fig. 6(b)], the vertical axis shows $\chi = \theta 2\u2212 \theta 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 Eq. (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 Fig. 6(b). These branches also exist for all values of $\epsilon $, except where the value of the isostable coordinate asymptotes to positive or negative infinity. There may be more than one phase-locked state with $\chi \u22600,\pi $ for a given value of $\epsilon $ depending on the interaction functions $ H 1,\u2026, H 6$. The existence, stability of solution branches, and locations of bifurcations in Fig. 6(b) all agree with the explicit calculations in Sec. III. Taking Fig. 6(a) 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 $\epsilon $ of the linear coupling increases.

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 $\epsilon $, but restabilizes at $\epsilon =0.499$, where it meets the branch of phase-locked solutions (shown in purple). Synchrony is also stable for a small interval of negative $\epsilon $, before undergoing a period doubling bifurcation at $\epsilon =\u22120.015$. The phase-isostable approximation, as shown in Fig. 6(b), correctly predicts the stability type of synchrony in the neighborhood of $\epsilon =0$ (as does a first-order phase approximation). The phase-isostable approximation is additionally able to capture bifurcations of synchrony for both positive and negative $\epsilon $, showing synchrony gaining stability where it meets a phase-locked state at $\epsilon =0.0934$ and losing stability at a Hopf bifurcation at $\epsilon =\u22120.407$. We note that while this reproduces the qualitative behavior of the synchronous branch, the phase-isostable approximation underestimates the precise value of $\epsilon $ at which the bifurcations occur, just as for the MF-CGLE in Sec. IV A.

In the full system, antisynchrony (shown in blue in Fig. 6) is unstable for all negative values of $\epsilon $, but is stable for $0<\epsilon <0.0831$, where the orbit lies inside but close to the synchronous orbit $\gamma $. For $0.00491<\epsilon <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 $\epsilon < \epsilon \u221e$ where the phase-isostable framework estimates that $ \epsilon \u221e=0.0961$ (compared to $0.0831$ in the full system). This point is labeled as a limit point (black square) in Fig. 6. As $\epsilon \u2192 \epsilon \u221e$ from below the stable solution branch has $\Psi \u2192\u2212\u221e$ 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 $\epsilon > \epsilon \u221e$ with $\Psi $ positive. As $\epsilon \u2192 \epsilon \u221e$ from above, this branch has $\Psi \u2192\u221e$ (corresponding to a growing amplitude orbit approaching the outer edge of the basin of attraction of the uncoupled node limit cycle; see Fig. 4). There is no bifurcation at $ \epsilon \u221e$; the change in stability here is coincidental. We have not found a corresponding antisynchronous state with orbit outside of $\gamma $ 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 $ \Psi $ for each $\epsilon $. It may be that if equation (20) were taken to higher order in $\epsilon $, then this other branch could be revealed. We note that the phase-isostable framework predicts a small region ( $0.0934<\epsilon < \epsilon \u221e$), 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\u2194 x 2$. For clarity, Fig. 6(a) 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 $\u22120.0218<\epsilon <0.575$ (shown in purple in Fig. 6). Its stability varies along the branch as it undergoes a series of the saddle node and torus bifurcations. Many of these occur within the region near $\epsilon =0$, but most obvious in Fig. 6(a) are the torus bifurcations at $\epsilon =0.0472$, $\epsilon =0.148$, and $\epsilon =0.538$ marked by black stars. There is also an isola of periodic solutions for $0.375<\epsilon <0.403$, an example of which is depicted in Figs. 7(b) and 7(d). This solution is not shown in Fig. 6(a) to avoid confusion with the synchronous branch since they share similar values of $max( v 1)$. For values of $\epsilon \u2208(0.148,0.375)\u222a(0.403,0.499)$, numerical simulations indicate that the stable solutions are quasiperiodic with oscillations fluctuating around the phase-locked state, as indicated in Figs. 7(a) and 7(c). In the phase-isostable approximation, branches of stable phase-locked states with $\chi \u22600,\pi $ bifurcate from synchrony at $\epsilon =0.0934$ losing stability at a Hopf bifurcation at $\epsilon =0.0484$. This gives rise to stable solutions where $\chi , \Psi 1, \Psi 2$ are all periodic corresponding to behavior observed in numerical simulations of the full model (Fig. 7); however, here these solutions exist over a much smaller interval of $\epsilon $ values. We also note that the stable phase-locked state has isostable coordinates of opposite signs for the two nodes corresponding to one node orbiting outside of $\gamma $ and the other inside. As the limit point at $\epsilon =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 $\epsilon $ but is unstable and always has one node with a large isostable coordinate.

The branch of periodic solutions shown in green in Fig. 6(a) is created at a homoclinic bifurcation at $\epsilon =0.0381$ and meets the unstable fixed point at $\epsilon =\u22120.0618$. This solution has $ x 1$ oscillating near the synchronous orbit $\gamma $, while $ x 2$ performs very small amplitude oscillations near the fixed point outside of the basin of attraction of $\gamma $. Since phase-isostable coordinates are valid only within the basin of attraction of $\gamma $, 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 $\epsilon =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 Sec. V B. We note that while the phase-isostable framework can be instructive concerning network behaviors for relatively weak coupling, since the expansions rely on the assumption that $\psi =O(\epsilon )$, predictions for larger values of $\epsilon $ or solution branches where the value of the isostable coordinate becomes large must be interpreted with caution. Additionally, the range of values of coupling strength $\epsilon $ where we expect the phase-isostable network equations to be predictive of the dynamics of the full network will be limited by the fact that first-order averaging has been used in the derivation of (20). An approximate upper bound on predictive values of the coupling strength $\epsilon $ can be shown to be $1/T$ where $T$ is the period of the uncoupled oscillations.^{64} However, for the Morris–Lecar node dynamics $1/T=0.122$ and we have observed from the direct comparison of full and approximated dynamics above that the phase-isostable network equations are able to predict qualitative network behavior beyond this estimated upper bound. If necessary, higher-order averaging techniques^{64} could be employed to improve the validity region of the averaged equations.

### B. Networks of many Morris–Lecar neurons

We next consider a network of 200 diffusively coupled Morris–Lecar neurons. We again use the parameter values in Appendix E so that the response and interaction functions remain as in Fig. 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 $\epsilon $. The results are shown in Fig. 8. The maximal real part of the eigenvalues for synchrony is positive (and, therefore, the state is unstable) when $0<\epsilon <0.0934$ and the splay state is unstable for all values of $\epsilon \u22600$. The discontinuity in the maximum real part of the eigenvalues for the splay state at $\epsilon =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 ( $\Psi <0$) when $\epsilon <0.0664$ and for $\epsilon >0.0664$, the splay state has $\Psi >0$.

We take $\epsilon =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 behavior predicted in this region. We initialize the system with the 200 nodes in a tightly packed group with phase coordinates $ \theta i\u2208[0.283725,0.283735]$ and isostable coordinates $ \psi i\u2208[2.9794,2.9798]$. Figure 9 shows the evolution in time of the isostable coordinates. We see the group of nodes initially move toward the isostable $\psi =0$ and remain near the synchronous network state for some time. The group then rapidly desynchronizes before settling into a stable 2-cluster state. The analysis in Sec. III B confirms the existence of the observed cluster state with $ N A=28$ and $ N B=172$ and also that it has phase difference $\chi =2.1407$ between the clusters that have orbits coinciding with the $ \Psi A=\u22120.1694$ and $ \Psi B=\u22120.1868$ isostables, also in agreement with the simulation. The stability analysis further confirms that this state is stable for $\epsilon =0.065$.

For the same value of the coupling strength $\epsilon =0.065$, we also simulate the full network equations from an initial state with $( v i, w i)$ in a group close to $( v \xaf, w \xaf)\u2248$ (-0.1,0.07) which has phase and isostable coordinates $(\theta ,\psi )=(0.28373,2.9796)$. Snapshots from the simulation are shown in Fig. 10, which also indicates the isostable $\psi =0$ that coincides with the synchronous orbit and the isostable corresponding to the splay state for these parameter values ( $\psi =\u221219.3$) in the $(v,w)$ space. Initially the nodes are close to synchronous and are first drawn toward and orbit near to $\psi =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 desynchronization eventually settling on a cluster state.

For stronger diffusive coupling ( $\epsilon =0.4$) with the node parameters still as in Appendix E, Han *et al.*^{40} observed that in the 200 node network initialized 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 desynchronizing. However, at this larger coupling strength, it was observed that the desynchronized nodes spiral in toward the unstable node inside the limit cycle, before again moving out toward the synchronous orbit. This behavior repeats resulting in mean-field voltage traces showing large amplitude fluctuations. Unfortunately, this type of behavior cannot be captured using the phase-isostable network equations at the order of (20) as the regime where $\epsilon =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 synchronous and splay states are unstable in agreement with simulations of the full network equations.

#### 1. Remarks

In Sec. III, we observed that in the phase-isostable network equations truncated at linear order in the isostable coordinates (or equivalently $O( \epsilon 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\u2208 B(\gamma )$, there is a unique corresponding $(\theta ,\psi )$. 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, Fig. 6(a) indicates a period doubling bifurcation of synchrony at $\epsilon =\u22120.015$. The period doubled branch [not shown in Fig. 6(a)] at $\epsilon =\u22120.0129$ is a periodic two-cluster phase-locked state, which is bistable with synchrony. The shared node orbit is shown in Fig. 11(a) and we observe that this has a self-intersection. The linear order truncation of the phase-isostable network equations used here 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 Eq. (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 an improved understanding of the behavior 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 that 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 phase reduction 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 the 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, the comparison of bifurcation diagrams for both full and phase-isostable approximations of the dynamics shows that the phase-isostable framework can capture qualitative changes in stability of synchrony and the antisynchronous 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 behavior beyond that which can be revealed using first-order phase reduction while 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 the 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 Sec. 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 the recent work of Bick *et al.*^{65} generalizes 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.

We note that the restriction to a single slowly decaying isostable direction described in this work is most relevant in cases where there is a significant separation in the magnitude of the real, negative Floquet exponents $ \kappa 1$ and $ \kappa 2$, (i.e., $ \kappa 2\u226a \kappa 1<0$). In this case, decay to limit cycle in directions corresponding to $ \kappa m$, $m=2,\u2026,n\u22121$ may be assumed to occur instantaneously so that $ ( \psi i ) m\u22610$ is a suitable approximation. In cases where node dynamics have distinct Floquet exponents $ \kappa 1$, $ \kappa 2$ of comparable magnitude, isostable dynamics in both corresponding isostable directions are likely to have comparable influence on the node dynamics and should both be included in the phase-isostable network equation approximation, requiring an extension of the work presented here.

Although we have chosen to use the orbit of an uncoupled oscillatory node as the reference for setting up the phase-amplitude coordinate system, other choices are possible. This has been recognized by Wilson as a way to handle perturbations that take one far from the 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).^{66} 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 desynchronizing populations of neural oscillators under periodic stimulation.^{28,67}

In the present work, we have investigated finite size networks with linear (state-dependent) coupling and two-dimensional 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 Ref. 68, 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 Ref. 69 and 70 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 generalizing the approach in Ref. 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 Refs. 71 and 72 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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**R. Nicks:** Conceptualization (equal); Formal analysis (lead); Investigation (equal); Methodology (equal); Project administration (equal); Software (supporting); Supervision (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **R. Allen:** Formal analysis (supporting); Investigation (equal); Software (lead); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal). **S. Coombes:** Conceptualization (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

### APPENDIX A: HIGHER-ORDER TERMS IN PRC AND IRC EXPANSIONS

When the $n$-dimensional system (2) has been reduced to the dynamics of the phase $\theta $ and a single isostable coordinate $\psi $ as in (10a) and (10b), the expansion to any order of the solution about the periodic orbit, the PRC and the IRC are given by (11a)–(11c). Following Ref. 31 and simplifying for the case of a single isostable coordinate, we here derive the adjoint equations (12a)–(12c), which must be satisfied by the vector functions $ g ( k )$, $ Z ( k )$, and $ I ( k )$.

#### 1. Asymptotic expansion of eigenfunctions

^{31,73}The Taylor expansion of (2) about $ x \gamma (\theta )$ then gives

^{31}We can obtain an explicit expression for these differential equations by noting that

#### 2. Asymptotic expansion of the PRC and IRC

^{31}we now consider the expansions of $ Z(\theta ,\psi )$ and $ I(\theta ,\psi )$ of (11b) and (11c), respectively. Using standard arguments, Wilson

^{31}shows that $ Z$ satisfies the adjoint equation

##### a. Normalization of terms in PRC and IRC expansions

### APPENDIX B: STABILITY OF PHASE-LOCKED STATES

#### 1. General phase-locked states

#### 2. Two-cluster states

^{58}says that for invertible $n\xd7n$ matrix $P$ and invertible $m\xd7m$ matrix $Q$, if $U$ and $V$ are $n\xd7m$ matrices, then

### APPENDIX C: HIGHER ORDER PHASE EQUATION DERIVATION

^{1}to derive the higher-order phase equation (C6) explicitly, including terms up to cubic order in the interaction strength $\epsilon $. Park and Wilson

^{1}take the isostable coordinate as an $O(\epsilon )$ term which may be expressed as $\psi (t)=\epsilon p ( 1 )(t)+ \epsilon 2 p ( 2 )(t)+\cdots $ where $ p ( k )(t)$ are $O(1)$. Using this additional expansion in (10b) and matching terms at different powers of $\epsilon $ results in a hierarchy of linear first-order differential equations for the $ p ( k )(t)$.

^{1}The equations at $O(\epsilon )$ and $O( \epsilon 2)$ are

### APPENDIX D: SPLAY STATE LINEAR STABILITY BOUNDARY FOR $N\u22653$

### APPENDIX E: THE MORRIS–LECAR MODEL

^{39}with two ionic currents: an outward non-inactivating potassium current and an inward, non-inactivating calcium current. Assuming that the calcium dynamics operate on a much faster time scale than the potassium, the dynamics of an isolated node are given by

^{40}$\varphi =1.15$, $ g C a=1$, $ g K=2$, $ g L=0.5$, $ E C a=1$, $ E K=\u22120.7$, $ E L=\u22120.5$, $ V 1=\u22120.01$, $ V 2=0.15$, $ V 3=0.1$, $ V 4=0.145$, and $ C m=1$. For these parameter values, a limit cycle of the system arises at $ I b=0.0730$ through a homoclinic connection. We choose $ I b=0.075$ to place the system near to the homoclinic bifurcation. For these parameter values, we find that the periodic orbit has Floquet exponent $\kappa =\u22120.4094$ and period $T=8.1654$.

## REFERENCES

*Chemical Oscillations, Waves, and Turbulence*, Springer Series in Synergetics (Springer-Verlag, Berlin, 1984).

*The Geometry of Biological Time*, 2nd ed., Springer Study Edition (Springer-Verlag, New York, 2001).

*Mathematical Foundations of Neuroscience*, Interdisciplinary Applied Mathematics Vol. 35 (Springer, New York, 2010).

*Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting*

*Weakly Connected Neural Networks*, Applied Mathematical Sciences (Springer-Verlag, New York, 1997).

*Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields*, Applied Mathematical Sciences Vol. 42 (Springer-Verlag, 1990).

*Neurodynamics: An Applied Mathematics Perspective*, Texts in Applied Mathematics Vol. 75 (Springer, 2023).

*Averaging Methods in Nonlinear Dynamical Systems*

*Perspectives and Problems in Nonlinear Science: A Celebratory Volume in Honor of Lawrence Sirovich*, edited by E. Kaplan, J. E. Marsden, and K. R. Sreenivasan (Springer New York, New York, NY, 2003), pp. 183–215.

*Simulating, Analyzing, and Animating Dynamical Systems—A Guide to XPPAUT for Researchers and Students*