The Dynamics of Networks of Identical Theta Neurons

We consider finite and infinite all-to-all coupled networks of identical theta neurons. Two types of synaptic interactions are investigated: instantaneous and delayed (via first-order synaptic processing). Extensive use is made of the Watanabe/Strogatz (WS) ansatz for reducing the dimension of networks of identical sinusoidally-coupled oscillators. As well as the degeneracy associated with the constants of motion of the WS ansatz, we also find continuous families of solutions for instantaneously coupled neurons, resulting from the reversibility of the reduced model and the form of the synaptic input. We also investigate a number of similar related models. We conclude that the dynamics of networks of all-to-all coupled identical neurons can be surprisingly complicated.


Introduction
Due to their analytical intractability, many investigations of large networks of model neurons involve extensive numerical simulation [1][2][3][4]. Without any detailed knowledge of the connectivity between neurons, one might assume the simplest form of for k = 1, 2, . . . , N, where κ is the strength of coupling (which could be positive or negative), and η is the input current to all neurons when uncoupled. The j th term in the sum (2) represents the pulse of current emitted by the j th neuron as it fires, i.e. θ j increases through π . The function (1 − cos θ) 2 is non-zero except when θ = 0, and thus this form of coupling may be regarded as non-physical, but the pulse can be localised more around θ = π by increasing the second power in (2) as will be discussed below. Note that each neuron is actually coupled to itself, but this term is only one out of N in the sum (2), so will be negligible for large N . We can write (1) as where ω = η + κI + 1 and H = i(η + κI − 1). Note that ω is real and H is imaginary. The WS ansatz [9,28] states that there is a transformation giving almost any solution of (1)-(2) (θ k ) in terms of N constants ({ψ k }, k = 1, 2, . . . , N) and three variables (ρ, Φ and Ψ ), where these variables satisfy the ODEs Thus, while it may seem possible that (for a given initial condition) a solution of (1)- (2) can explore the full N -dimensional phase space described by the N values of θ k , the solution is actually constrained to lie on a three-dimensional manifold with coordinates ρ, Φ and Ψ . The dynamics on this manifold will depend on the values of the constants {ψ k }.
There are N variables {θ k }, N constants {ψ k } and three variables ρ, Φ and Ψ , and thus one needs to specify three constraints so that there is a unique relationship between {θ i } and {ψ k , ρ, Φ, Ψ }, to determine initial conditions, for example. One way to do this is to set ρ(0) = Φ(0) = Ψ (0) = 0, so that {ψ k } = θ k (0). Then integrate (5)- (7) with ρ(0) = Φ(0) = Ψ (0) = 0 and use the solutions ρ(t), Φ(t) and Ψ (t) to reconstruct {θ k (t)} using (4). Here, the constraints are on ρ(0), Φ(0) and Ψ (0). Another set of constraints often taken is [28] N k=1 e iψ k = 0; Re N k=1 e 2iψ k = 0. (8) Given the N initial values θ k (0), one can usually uniquely determine the values of ρ(0), Φ(0) and Ψ (0) and {ψ k } such that both (4) and (8) hold [9]. "Usually" refers to solutions for which fewer than half of the neurons have exactly the same state. This rules out full synchrony, which is often a state of interest; however, we can understand the fully-synchronous state by simply considering the behaviour of one self-coupled neuron. We will use (8) unless specified otherwise. In order to use (5)- (7), we need to write I (and thus ω and H ) in terms of the new variables and the constants {ψ k }. Now [28] where overbar indicates complex conjugate, z = ρe iΦ , and We define and see that C 0 = 1, and from (8), C 1 = 0. Using a series expansion of [1 + ρe i(ψ k −Ψ ) ] −1 , we can write in the general case For the special case of ψ k = 2πk/N , i.e. evenly spaced ψ k , C n = 0 except when n is a multiple of N , when it equals 1. Then (14) is a geometric series and (note that this expression is given incorrectly in [17,21]).
In the general case, For evenly-spaced ψ k , Irrespective of ψ k , if ρ = 1 then γ = γ 2 = 1. If ρ < 1 and {ψ k } are uniformly distributed, then as N → ∞ we see that γ → 1 and γ 2 → 1. In these cases (γ = γ 2 = 1) I becomes independent of Ψ and (5) and (6) decouple from (7), i.e. the dynamics becomes two-dimensional. An alternative description of the dynamics of (1)- (2) when N = ∞ is given by writing the continuity equation governing the evolution of the probability density of θ s, p(θ, t) [26]. One can then use the Ott/Antonsen (OA) ansatz [29,30] to reduce the dynamics of this evolution equation to a single complex equation for the evolution of the order parameter z ≡ 2π 0 p(θ, t)e iθ dθ : where ω and H are as above and I = 3/2 − (z +z) + (z 2 +z 2 )/4. Substituting z = ρe iΦ into (20) and taking real and imaginary parts, we recover (5) and (6). Thus the OA ansatz corresponds to a special case of the WS ansatz: when N = ∞ and the constants {ψ k } are uniformly spread over [0, 2π] [28]. Note that while the OA ansatz gives dynamics on an invariant manifold in the space of all p(θ, t), if the neurons are identical, the manifold is not attracting, and thus the full dynamics is not given by (20) and must be described using the WS ansatz [17,28].

Infinite N , Equally-Spaced Constants
We now consider the case of N = ∞ with uniform density of ψ k . Thus we are interested in solutions of (20), written as where Note that (21)- (22) are invariant under (z, t) → (z, −t), i.e. simultaneous reflection in the real axis and reversal of time. This has a significant effect on the possible dynamics.

Excitatory Coupling
First consider κ = 1. From setting dρ/dt = 0 in (5) and recalling that H is imaginary, we find two sorts of fixed points of (21)-(22): (i) those with ρ = 1 and (ii) those with Φ = 0. From (6), those with ρ = 1 satisfy where These solutions are plotted in Fig. 1(a). Note that these solutions can be found directly from (1)- (2). For identical neurons, ρ = 1 corresponds to full locking, so all θ i are equal to Φ and a simple trigonometric identity gives (24) from (2). Solutions of type (ii) with Φ = 0 have z = ρ. From (6), they satisfy where I = 3/2 − 2ρ + ρ 2 /2. These solutions (restricted to −1 ≤ ρ ≤ 1 for physical reasons) are shown in Fig. 1(b). Two solutions are created in a saddle-node bifurcation as η is increased (at a negative value of η). One is a saddle whose eigenvalues sum to zero, and the other is a focus with purely imaginary eigenvalues; these properties result from the reversibility of the system [31]. These fixed points correspond to splay states in which all neurons follow the same trajectory but are equally displaced from one another in time such that average quantities such as z are constant [9,32,33]. Note that the solutions shown in Fig. 1(a) collide with the saddle solution in Fig. 1(b) at (η, z) = (0, 1). A selection of solutions and the fixed points are shown in Fig. 2(a) for η = −0.5. We see that for these parameter values, initial conditions either tend to the stable fixed point in the lower half plane (i.e. quiescence), or (if they are in the region enclosed by the homoclinic orbit to the saddle fixed point) follow one of a continuous family of periodic orbits. For η > 0, the only remaining fixed point is the focus, and the phase space (ρ ≤ 1) is filled with a continuous family of periodic orbits (see Fig. 2(b)), again, as a result of the system's reversibility.

Inhibitory Coupling
Now consider κ = −0.5. As with excitatory coupling, there are fixed points with ρ = 1, given by (23)- (24) and shown in Fig. 3(a). For η > 0, there is also a focus fixed point (a splay state), as shown in Fig. 3(b). This fixed point is surrounded by a continuum of periodic orbits, as in Fig. 2(b). (Note that for strong inhibition, i.e. κ large and negative, the situation shown in Fig. 3(a) can become more complex, with multiple stable fixed points. This is due to the finite width of the pulses in (2), and the region of multistability disappears as the pulses are made narrower.) In summary, the dynamics of (21)- (22) is non-generic due to their reversible nature, which can lead to the existence of continuous families of neutrally-stable periodic orbits.

Finite N , Equally-Spaced Constants
We now consider the case of finite N but with equally-spaced {ψ k }. Thus we consider (5)- (7) where I depends on Ψ via γ and γ 2 . First we point out the reversibility of this system under the transformation (ρ, Φ, Ψ, t) → (ρ, −Φ, −Ψ, −t). This transformation interchanges the products zγ andzγ , and z 2 γ 2 and its conjugate. This leaves I (and thus ω and H ) unchanged. Keeping in mind that H is imaginary, we see that this transformation leaves (5)-(7) unchanged, i.e. they are reversible.

Fixed Points
The fixed points of type (i) with ρ = 1 persist, independent of N , since these solutions have γ = γ 2 = 1. The values of Φ are given by solving (23)- (24). However, these fixed points have arbitrary values of Ψ since dΨ/dt = 0. Regarding the fixed points of type (ii) that exist for N = ∞ as analysed in Sect. 2.2.1, they had constant and generically non-zero dΨ/dt. Thus, for finite N , we expect these to appear as time-dependent orbits, with the amplitudes of fluctuations in ρ and Φ going to zero as N → ∞. To understand this, assume to a first approximation that ρ is constant. Then γ and γ 2 will have N periods of oscillation as Ψ goes through one period of oscillation. Thus I, ω and H will all have N oscillations in one period of Ψ and so will ρ and Φ. Thus the fixed points of type (ii) which exist for infinite N will appear for finite N as quasiperiodic orbits in which ρ and Φ undergo N oscillations for each full rotation in Ψ . An example for N = 4 is shown in Fig. 4, where Ψ decreases from π to −π .
The amplitude of this type of periodic orbit goes to zero exponentially as N → ∞; see Fig. 5. Even though N is physically an integer, the expressions for γ and γ 2 do not require it to be so. Thus we can, for example, continue the saddle-node bifurcation seen in Fig. 1(b) as a function of the continuous parameter N . The result is shown in Fig. 6, where we also show interpolated values at integer N . Interestingly, while the curve oscillates, the values at integer N are monotonic. The periodic orbit corresponding to the focus fixed point in Fig. 1 In terms of the original system (1)-(2), the periodic orbit shown in Fig. 4 corresponds to a periodic orbit with all Floquet multipliers having magnitude 1, i.e. completely neutrally stable. Two of the multipliers are a complex conjugate pair corresponding to the rotation seen in Fig. 2, while the remaining N − 2 are equal to 1.

Other Orbits
Consider the continuous family of periodic orbits which exist in the infinite-N case for κ = 1 and η sufficiently positive (see Fig. 2(b)). These persist as what seems to be a continuous family of quasiperiodic orbits. Some are shown in Fig. 7 where we plot the value of z when α ≡ Φ − Ψ increases through a multiple of 2π . Now consider the case κ = −0.5. The dynamics in this case seems much more complex. An example is shown in Fig. 8 for η = 0.6. For some initial conditions, the dynamics seems quasiperiodic, while for others the orbits appear to be chaotic (as indicated by a positive Lyapunov exponent, not shown). This mixture of quasiperiodic and chaotic behaviour has been previously observed in reversible systems [34] including a resistively-loaded series array of Josephson junctions also studied using the WS ansatz [22]. The overall trend for this system is that the dynamics becomes more regular as η is increased. We leave the investigation of this dynamics for a future publication.

Finite N , Non-Uniform ψ k
Now consider non-uniformly spaced constants ψ k . We follow [28] and distribute these values along two arcs each of length qπ , where 0 < q ≤ 1 is a parameter. Choosing N to be even, we have ψ k = (1 − q)π/2 − qπ/N + 2πqk/N for k = 1, 2, . . . , N/2 and ψ k = (3 − q)π/2 − qπ/N + 2πqk/N for k = N/2 + 1, 2, . . . , N. For q = 1, this is a uniform distribution, and as q → 0 the distribution tends to the two points ±π/2. The saddle-node bifurcation shown in Fig. 1(b) moves as a function q, as shown in Fig. 9. Varying q for κ < 0 also gives a variety of different dynamics, as shown in Fig. 10. Here we see a mixture of quasiperiodic and more complex behaviour, as seen in Fig. 8.

Summary
In summary, for instantaneous synapses of the form studied here, there may be continuous families of periodic orbits (for drive η > 0 and for some η < 0 if κ > 0) in the infinite-N uniformly distributed ψ k case. This is due to the reversibility of (21)- (22). For finite N , some of these orbits become quasiperiodic, with some initial conditions showing either quasiperiodic or more complex types of behaviour. This type of network can be thought of as having two sources of degeneracy: even if the constants {ψ k } are fixed, depending on parameters, there may be continuous families of neutrally stable periodic or quasiperiodic orbits. Choosing different initial conditions for the WS variables ρ, Φ and Ψ selects between these orbits. Secondly, even for fixed initial conditions of the WS variables, there are many continuous families of orbits that can be obtained by varying {ψ k }.
Note that if the system is bistable (for excitatory coupling and sufficiently small negative drive η), then even for fixed initial values of ρ, Ψ and Φ, changing either N (for evenly spaced constants ψ k ) or the distribution of ψ k (for fixed N ) can lead to very different outcomes, as these changes can move the system from one basin of attraction to another.

Synaptic Dynamics
We now consider including some synaptic processing, which amounts to delaying the synaptic input, still in the form of a current input. We thus replace (2) with where The network of neurons is still described by (5)- (7) but with the addition of where z, γ and γ 2 are as in Sect. 2.1.

Infinite N , Equally Spaced Constants
Repeating analysis as in Sect. 2.1, we now have Fixed points of this pair of equations are the same as those of (21)- (22), but the stability of some of them has changed. Fixed points of type (i) shown in Fig. 1(a) gain another negative eigenvalue. The focus fixed point shown in Fig. 1(b) becomes Fig. 11 Frequency of the stable periodic orbit of (29)- (30). τ = 1, κ = −0.5 stable, while the saddle fixed point in that figure remains a saddle, but with two stable directions. The continuum of periodic solutions shown in both panels of Fig. 2 is destroyed, and the system has either one stable fixed point, or two (in the region of bistability for η < 0 but not too negative). In terms of attractors, (29)- (30) show what one would expect from an excitatorily self-coupled network. We have stable quiescence for large negative drive, bistability between quiescence and an active splay state for small negative drive, and a stable splay state for positive drive.
For the case of inhibitory coupling, fixed points of type (i) shown in Fig. 3(a) gain another negative eigenvalue (as for excitatory coupling). The focus fixed point in Fig. 3(b) now has one stable direction and two unstable ones. The continuum of periodic orbits seen for η > 0 mentioned in Sect. 2.2.2 is now replaced by a single stable periodic orbit with ρ = 1. The frequency of this orbit increases from zero as η does, as shown in Fig. 11.

Finite N , Equally Spaced Constants
As was found in the case of instantaneous synapses, solutions with ρ = 1 are unaffected by this change. However, for κ = 1, the stable fixed point that exists for η > 0 and N = ∞ now has small amplitude oscillations, as discussed in Sect. 2.3.1. The size of these oscillations decays exponentially with N , as was seen in Fig. 5. The only difference with the case discussed in Sect. 2.3.1 is that this low-amplitude periodic orbit is stable, whereas the one discussed in Sect. 2.3.1 was neutrally stable.

Finite N , Unequally Spaced Constants
Here we distribute the constants ψ k as in Sect. 2.4 and investigate the effects of varying q on the amplitude of the oscillations in ρ of the stable periodic orbit that exists for κ = 1, η = 1, τ = 1. The results are shown in Fig. 12. For q = 1, we saw previously that the amplitude of oscillations decays exponentially to zero as N increases. However, for q < 1, the amplitude is always finite and has a limiting value as N → ∞. As q → 0, the amplitude increases and the value of N becomes less relevant. This can be seen by examining the coefficients C n (13). For q = 0, the first N/2 of ψ k equal π/2, while the last N/2 equal −π/2. Inserting this into (13), we find that C n = 0 for n odd, and C 2n = (−1) n , independent of N .

Summary
In summary, adding synaptic dynamics of the form examined in this section destroys the reversibility of (21) and thus the non-generic behaviour of solutions of

Other Similar Models
In this section we discuss several similar models and consider the case of identical oscillators.

Gap Junction Coupling
Here we consider a network of all-to-all gap junctionally coupled theta neurons. Laing [25] showed, using the equivalence of a theta neuron and a quadratic integrateand-fire neuron [35], that a network of N identical gap junctionally coupled theta neurons can be written as follows: where g is the strength of coupling, I is a constant input, and where 0 < 1. We can write (31) as where ω = 1 + I + gQ and H = g + i(I + gQ − 1). Thus the solutions of (31) can be described by the three ODEs (5)- (7) with the above definitions of ω and H .
Extending the results of Laing [25], we find that where z = ρe iΦ , "c.c." indicates the complex conjugate of the previous term, For N = ∞ and equally-spaced ψ k , one can show [28] that all γ m = 1 and thus (5)-(6) are sufficient to describe this system. Written in complex form, (5)- (6) are Note that unlike (21)- (22), equation (37) In a similar way to the synaptically coupled network, (37) has two fixed points with ρ = 1 for I < 0 (these correspond to all neuron states being equal). From (5)-(6), these satisfy which could also be derived directly from (31). These fixed points are shown in Fig. 13(b).
For large and negative I , the fixed points with ρ = 1 are a source and a sink, and there exists another unphysical saddle fixed point with ρ > 1 (Fig. 13(a)). As I is increased, the unphysical solution is involved in a transcritical bifurcation with the source fixed point having ρ = 1, and it enters the unit circle. As I is increased through zero, the synchronous fixed points are destroyed in a saddle-node on invariant circle (SNIC) bifurcation, leading to synchronous oscillations with ρ = 1, with a source fixed point inside the unit circle.
The only attractors for this system have ρ = 1, and for these it is clear that all γ m = 1 for finite N and any ψ k , and thus their dynamics will still be described by (37) in this case. Solutions with ρ < 1 (none of which are attracting) will, however, be described by (5)-(7), and their dynamics will depend on N and the distribution of ψ k .

Conductance Dynamics
Coombes and Byrne [36] considered a model in which synaptic input was in the form of a current, equal to the product of a conductance and the difference between the voltage of a quadratic integrate-and-fire neuron and a reversal potential. A particular case of their model can be written as follows: where η is a constant drive, V syn is the reversal potential (positive for excitatory coupling, negative for inhibitory one), and where k is the coupling strength (0 < k) and δ is the Dirac delta. Equation (39) can be written as follows: where ω = 1 + η + gV syn and H = g + i(η + gV syn − 1). Thus, as in Sect. 4.1, the solutions of (39) can be described by the three ODEs (5)-(7) with these definitions of ω and H . We have where z = ρe iΦ , κ = k/π , and γ m is given in (36). For infinite N and equally spaced ψ k , the system is described by where all γ m = 1. Synchronous fixed points have ρ = 1 and g = 0, and thus satisfy 0 = 1 + η + (η − 1) cos Φ, which we could obtain directly from (39). These fixed points only exist for η ≤ 0. For excitatory coupling (V syn = 2, κ = 1), the fixed points are qualitatively the same as in Fig. 1, although the pair with ρ = 0 do not have Φ = 0. This pair (analogous to those in panel (b) of Fig. 1) are a stable focus and a saddle, just as in Sect. 3.1, leading to a region of bistability. For inhibitory coupling (V syn = −2, κ = 1), the fixed points are qualitatively the same as in Fig. 3; although, again, the fixed point for η > 0 does not have Φ = 0, and it is a stable focus.
We note that this model does not have the reversibility of that in Sect. 2, and thus none of the non-generic behaviour seen there. The analysis of this model for finite N should be similar to that in Sect. 3, but we do not present the results here.

Winfree Model
The Winfree model of N pulse-coupled oscillators dates from 1967 [37][38][39] and is written for identical oscillators as where we choose Q(θ ) = sin β − sin (θ + β) and P (θ) = a n (1 + cos θ) n where a n is a constant such that 2π 0 P (θ) dθ = 2π . The function Q is the phase response curve of an oscillator, which can be measured experimentally or determined from a model neuron [40], and P (θ) is the pulsatile signal sent by a neuron whose state is θ . We can write (44) as where ω = Ω + εσ h and H = εe −iβ h and Thus the solutions of (44) can be described by the three ODEs (5)-(7) with the above definitions of ω and H . To be concrete, choose n = 2, in which case where z = ρe iΦ and γ and γ 2 are as given in (11)- (12). As for the network of theta neurons, in the case of N = ∞ and equally-spaced ψ k , γ = γ 2 = 1 and (7) decouples from (5)- (6), which are, for this system, where we have rescaled time so that Ω = 1, and h = 1 and (ii) up to two for which ρ = 1 (locked states), satisfying These fixed points and their stability are shown in Fig. 14 as a function of ε for β = 0.1. For small ε, the only attractor is a periodic orbit with ρ = 1 corresponding to synchronous oscillations. As ε is increased, there is a saddle-node on invariant circle (SNIC) bifurcation destroying the periodic orbit and leading to the creation of type (ii) fixed points, one of which is stable. As ε is increased further, the type (i) fixed point undergoes a transcritical bifurcation with the saddle fixed point on ρ = 1 and leaves the unit circle, thereby becoming unphysical. This sequence of bifurcations happens for all 0 ≤ β < π/2. Note the similarity between this scenario and that for the gap junction coupled neurons in Sect. 4.1: both have a SNIC bifurcation on the circle ρ = 1, with the saddle later becoming a source as another fixed point leaves the unit circle.
As was the case in Sect. 4.1, since the only attractors have ρ = 1, these will persist unchanged for finite N and any ψ k , as γ and γ 2 will both equal 1 in this case.

QIF Neurons
The correspondence between a theta neuron and a quadratic integrate-and-fire (QIF) neuron is well known [11,25,41]. Thus, some of the degenerate behaviour seen here in the networks of identical theta neurons should also appear in all-to-all coupled networks of identical QIF neurons [32]. Consider the network in [41]: for i = 1, 2, . . . , N, with the rule that when V i = ∞ (i.e. neuron i fires) V i is set to −∞. Here, I 0 is a constant, J is the strength of coupling between the neurons, and r(t) is the rate at which the network is firing: where t k j is the kth firing time of the j th neuron, and the sum over k is only over past firing times. Using the transformation V j = tan (θ j /2), (52) becomes which is the same form as (1). Thus, under the assumption of infinite N and equally spaced constants ψ k , the dynamics of (54) is governed by (see (21)). The authors [41] showed that defining equation (55) can be written and that the real part of w, divided by π , is the firing rate r, while the imaginary part of w is the mean of the voltages V i in the original network (52). Thus writing w = πr + i V , (57) can be written as two real equations: Note the reversibility: ( V , t) → (− V , −t). Equations (58)-(59) are the same as equations (12a) and (12b) in [41], once identical neurons are considered. Equations (58)-(59) will describe the dynamics of (52) in the limit N → ∞, and when the constants ψ k are equally spaced. We now briefly discuss initial conditions for (52) and their relationship to ψ k . Using (4) and V k = tan (θ k /2), we have Equally-spaced ψ k means ψ k = 2πk/N , k = 1, 2, . . . , N, and since tan is periodic with period π , we see that in the limit N → ∞, the value of Ψ becomes irrelevant in determining V k . (We expect this, as in this case-as we have seen a number of times-(5)-(6) decouple from (7).) For finite N , we will set Ψ = 0 and initialise V k as Varying Φ and ρ corresponds to moving through a two-parameter family of initial conditions for (52) which, in the limit N → ∞, corresponds to varying the initial conditions of (58)-(59). We implemented (52) just as in [41] but with N = 1000 and identical neurons. Setting I 0 = 0.2, J = 0.1 and using three different values of Φ and ρ to initialise V k , we obtained the results in Fig. 15(a). (Both the mean voltage and firing rate were smoothed by convolving with a Gaussian of standard deviation 0.07 time units.) We see three out of a continuous family of periodic orbits. They match very closely numerical simulations of (58)-(59) (not shown) even though (58)-(59) are only valid in the limit N → ∞. Conversely, if we choose V k (0) not using (61), for example, taking them randomly from a unit normal distribution, we obtain the results in Fig. 15(b). The solution is quasiperiodic and clearly not described by a planar system. To understand this, consider the transformation (4) with θ k (0) = 2 tan −1 (V k (0)). For large N , it is generally impossible to choose Φ(0), Ψ (0) and ρ(0) such that ψ k are uniformly distributed. Thus the system must be described by three coupled equations of the form (5)-(7) which, even in the limit N → ∞, will not reduce to a planar system due to the non-uniformity of the distribution of ψ k .
Note that while (52) is synaptically driven by the instantaneous rate r, the network (1) is driven by the average of the pulse-like functions in (2). However, replacing (1 − cos θ) 2 in (2) by a n (1 − cos θ) n , where a n is chosen so that the integral of this function over [0, 2π] is independent of n, and taking the limit n → ∞, one can obtain a drive equal (up to a scale factor) to the rate given in terms of delta functions in (53) [24,36].

Discussion and Conclusion
In conclusion, we studied finite and infinite networks of identical all-to-all coupled theta neurons, with both instantaneous and delayed synaptic interactions, and several related models (gap junction-coupled theta neurons, theta neurons with conductance dynamics, Winfree oscillators and quadratic integrate-and-fire neurons). For all models, the WS ansatz gives a reduced description of the model in terms of three ODEs and a set of constants. Changing these constants while keeping the initial conditions of the WS variables fixed gives different dynamics in the original models, and thus there is a continuum of such solutions. For instantaneous synapses, even keeping the constants fixed, there can be many solutions, due to either the reversibility of the dynamics of the WS variables, or the coexistence of many quasiperiodic orbits. It is only this case, of instantaneous synapses, which seems non-generic. Adding synaptic processing destroys this reversibility, and none of the other models (gap junction-coupled theta neurons, theta neurons with conductance dynamics or Winfree oscillators) have such a reversibility.
We have considered identical sinusoidally-coupled oscillators, for which the WS ansatz [9] gives the correct description in the form of the three ODEs (5)- (7), together with the constants {ψ k }. This is an idealisation, so we should discuss the case of non-identical oscillators. For an infinite number of non-identical oscillators, one can consider the continuity equation governing the evolution of the probability density of the phases [9,26,29]. Solutions of this decay onto the Ott/Antonsen (OA) manifold, on which their dynamics can be found using the OA ansatz [24-26, 29, 30, 37, 42, 43]. The OA equations are simpler than those resulting from the Watanabe/Strogatz ansatz: only a pair of real equations are needed to describe (1)-(2) when the value of η for each neuron is different, for example, rather than (5)-(7) and {ψ k } (although we must take N → ∞ for them to be valid). The OA equations typically have hyperbolic and isolated fixed points, unlike those obtained using the Watanabe/Strogatz ansatz, which can be perturbed to a nearby similar state by varying one of ψ k , for example. However, there remains a gap: is there a dimension reduction applicable to finite networks of non-identical sinusoidally-coupled oscillators? See [44] for ideas in this direction. Another approach is that of [28] who consider a number of subpopulations, each of which has identical oscillators although parameters are different between subpopulations, and the WS ansatz can be applied to each subpopulation.
There are several extensions that could be performed using the ideas presented here. The first is time-dependant forcing of the network [45]. As long as each neuron experiences the same force, the WS ansatz will apply. The second involves each neuron receiving common noise [46,47]. Again, the WS ansatz will apply. We could also consider introducing a discrete delay in the synaptic processing, i.e. replacing I (t) in (1) by I (t − τ ) for some delay τ > 0. Another extension involves investigating a pair of populations of neurons, one excitatory and one inhibitory [41]. Coupling these may result in a PING rhythm [48], and if neurons are identical within each population, the coupled network may show interesting dynamics.