Integrability of a globally coupled complex Riccati array: quadratic integrate-and-fire neurons, phase oscillators and all in between

We present an exact dimensionality reduction for dynamics of an arbitrary array of globally coupled complex-valued Riccati equations. It generalizes the Watanabe-Strogatz theory [Phys. Rev. Lett. 70, 2391 (1993)] for sinusoidally coupled phase oscillators and seamlessly includes quadratic integrate-and-fire neurons as the real-valued special case. This simple formulation reshapes our understanding of a broad class of coupled systems - including a particular class of phase-amplitude oscillators - which newly fall under the category of integrable systems. Precise and rigorous analysis of complex Riccati arrays is now within reach, paving a way to a deeper understanding of emergent behavior of collective dynamics in coupled systems.

The study of complex systems often involves describing them using a few relevant variables known as order parameters.Such dimensionality reduction techniques are highly useful but challenging to discover, and they may not exist for every system.In this regard, Watanabe and Strogatz [1] (WS) made significant progress by demonstrating how globally coupled phase oscillators can be effectively described using only three order parameters.This reduction arises from a specific subclass of complex Riccati equations [2,3] that govern the dynamics of such oscillatory arrays.The macroscopic variables are parameters of a Möbius transform that relates the dynamical state variables to constants determined by the initial states of the oscillator array [4].Remarkably, Ott and Antonsen [5] showed that for large ensembles of nonidentical oscillators, the dynamics even collapses to an effectively 2D manifold.These descriptions have proven immensely useful in studying the collective dynamics of coupled oscillatory systems [6] and found many applications, particularly in neuroscience [7,8].
In this Letter, we present a framework that generalizes these dimensionality reductions to a larger class of systems: an arbitrary array j = 1, ..., N of globally forced complex Riccati equations, ẋj = ax 2 j + bx j + c , a, b, c ∈ C , where all a, b, c ∈ C can be arbitrary complex functions of time and x j ∈ C can start with arbitrary complex values.The choice for a, b, c allows selecting from a range of coupled oscillator systems.Certain choices reproduce known models, such as the quadratic integrate-and-fire neurons (QIF) and phase oscillators; but the possibilities go far beyond that and cover a broad class of 2D systems.The flow of array variables x j in (1) is given by the following Möbius transformation: where the three global complex-valued variables Q, y, s evolve according to: and ξ j ∈ C are constants determined by the initial values of the oscillator array variables [9] (cf.Eqs. ( 4) in [10] and Eqs.(15)(16)(17) in [11]).Since these equations completely generate the flow of the original system (1), this implies that the dynamics of Eqs. ( 1) is effectively six dimensional.See Supplemental Material [12] for a proof of this result.

Choosing initial conditions.
Since we started with an array of N variables {x j } and now we describe the system with the three macroscopic variables Q, y, s and N constants {ξ j }, we have some freedom in choosing the relationship between these variables.An additional constraint of choice is set on the initial values, that defines the exact relationship between x j and Q, y, s, ξ j .We present two convenient options here.(I) Identity conversion: a straightforward way to determine initial conditions is requiring that variables x j initially coincide with constants ξ j (cf.Eq. (3.7) in [2]) (II) Möbius conversion: in general, the relation between x j and ξ j is a Möbius transform.We here use In different situations different initial condition choices are more appropriate -we will see examples of both options in the later examples.Other options for initial constraints are possible as well, e.g., j ξ j = 0 (cf.Eq. (4.12) in [2]), see Supplemental Material [12] for more details.Special case of real-valued arrays.We can consider the special case of real coefficients and real initial conditions: a, b, c, x j (0) ∈ R. Consequently, the flow of variables is real-valued, x j (t) ∈ R for all t ≥ 0. The dynamics is then three dimensional, which is easiest to see with constraint (4) which make Q, y, s initially real, and since they follow Eqs.(3) with real coefficients, they remain real with evolution.An example of this special case is a globally coupled array of identical QIF neurons [13,14].Individual voltages x j (t) obey ẋj = x 2 j + I , if x j > x thr then x j → x reset , (6) where the voltage threshold and reset values are: x thr = ∞ and x reset = −∞.The current I can have a constant component I 0 associated with intrinsic neuronal dynamics, but it can also represent external forcing (even noisy) or global coupling, e.g., the input current generated by N globally coupled QIF neurons with pulses P (u) is expressed as: Within our formalism, the voltage spikes naturally occur when the denominator 1+sξ j in (2) crosses zero.As a result, the voltage in that instance reaches +∞ upon which it is reset to −∞.What has to be considered, however, is that depending on the chosen initial constraint, Q might also diverge.Indeed, this is the case if one chooses initial conditions according to (4): in which case additional resetting of variables would be needed [15].However, diverging variables and additional resetting can be avoided by simply choosing the appropriate initial conditions (5): The constants ξ j then relate to x j via a Möbius transform: . Since x j (0) take real values, the constants ξ j are unitary |ξ j | = 1, and can thus be described by their argument ψ j : ξ j = e iψj .
For this specific case of a, b, c, x j ∈ R and initial constraint (5) the following simplification is true: which reduces the dynamical equations (3) to where ζ ∈ R is the argument of s = e iζ .Note that these dynamics are three dimensional since Q now takes on complex values (even though x j ∈ R).The transformation (2) is reduced to: The variable Q remains bounded for all times, no additional resetting is needed and the spikes occur when the denominator 1 + e i(ψj +ζ) crosses 0, i.e., when ζ = π − ψ j .
A numerical example of this special case is shown later in Fig. 1.One can arrive at the same dynamics by transforming the QIF neurons into theta neurons via the transformation x j = tan(θ j /2) and then employing the Watanabe-Strogatz theory for phase oscillators [1,2,4].For more on how our formalism relates to the Watanabe-Strogatz theory for phase oscillators see Supplemental Material [12], which includes Refs.[16,17].
Such a description of QIF neurons has already been considered in the continuum limit N → ∞, cf.Eqs.(42) in [18].If additionally the distribution of x j variables is Cauchy-Lorenztian, the dynamics further simplifies and variable Q completely determines the macroscopic state, its real and imaginary component correlating to the firing rate and mean voltage [19]. Examples.
We now present some specific systems where the dimensionality reduction can be applied and with numerical simulations validate its exactness.First (I), a known and relatable example of pulse-coupled realvalued QIF neurons.It is known that this system possesses low-dimensional dynamics since the QIF model can be transformed into a theta neuron to which the WS theory applies.Our formalism provides not only a new perspective of this fact, but also justifies the voltage resetting at infinity by viewing the model as a limiting case of the extended complex model.Next, we show two examples which generalize known models to complex values.Example (II), the complex generalization of the QIF model, where we simply allow the "voltage" variables to attain complex values.And example (III) concerns a complex-valued generalization of phase oscillators, specifically we choose overdamped Josephson junctions.An additional example is found in Supplemental Material [12] of how our description applies to infinite ensembles in the thermodynamic limit, and how particular integrals can simplify with set initial conditions.
For ensembles of pure phase oscillators the Kuramoto order parameter is defined as Z 1 = 1/N N j=1 e iφj where its amplitude |Z 1 | quantifies the order in the system.This is easily generalized to the full complex plane by simply invoking the mean [20]: which is neatly expressed with dynamical variables Q, y, s and constants ξ j , The sum 1/N j ξ j /(1 + sξ j ) can be seen as a function constant in time, that holds the information of initial conditions and only depends on the variable s.Example (II): complex QIF model.Let us consider a simple generalization of the QIF neurons to the complex plane.If the voltages x j start off the real axis, then they never diverge and there is no need for resetting conditions in (6).Let us consider such generalized QIF neurons, globally coupled via the first moment Z 1 (11), where I 0 is the intrinsic input current and x + 0 the positive-imaginary fixed point of individual neurons (which can be absorbed in the current I 0 → I 0 − ϵx + 0 ).In the form of the initial Riccati equation (1) this model translates to the parameters: a = 1, b = 0, c = I 0 + ϵ(Z 1 − x + 0 ).On the real line the behavior of an individual unit x j tends towards infinity in finite time and hence one needs a reset rule (6) as well as an implementation of a pulse during that event.However, if the dynamics occurs instead in the complex plane, the trajectory of this unit naturally oscillates around a fixed point.The uncoupled system is time-reversible, invariant under the involution: Re[x j ] → −Re[x j ], t → −t, which allows for a continuum of closed orbits symmetric with respect to the imaginary axis if I 0 > 0, see Supplemental Material for details [12].Indeed, we find two center fixed points of the single unit dynamics: x ± 0 = ±i √ I 0 .When we couple several units together, they remain oscillatory.In our example we use N = 8 oscillators with I 0 = 1, ϵ = −5 and initial conditions: {x j (0)} = {x + 0 + j 2 20 exp(i π 2N (j − 1))}, j = 1, ..., N .The dynamics settles into periodic motion with a nontrivial limit cycle, as shown in Fig. 2.
Complex generalization of coupled quadratic integrate-and-fire neurons (13).Unlike the real QIF model, the voltage resetting in this complex generalization is redundant since the dynamics everywhere (except the special case of real line with real coupling) loops back and stays finite; see also Example (I) in Fig. 1 where voltages diverge, but the resetting naturally results from the transformation (2).N = 8 units settle into periodic motion with a nontrivial limit cycle.Trajectories of oscillators on the limit cycle in the complex plane are depicted with colored lines.Initial conditions are marked with (color-coded) points.The fixed point x + 0 is emphasized with a black dot.
Example (III): complex generalization of phase oscillators.Now let us generalize phase oscillators to include a free amplitude.Consider the phase dynamics equation in complex exponential form: but allow that oscillators have amplitude r j different from 1: x j = r j e iφj .This choice results in a special family of phase-amplitude oscillators: Note that r j = 1 defines an invariant subspace that oscillators cannot cross: oscillators that start on the inside of the unit disk stay inside the disk forever.For this example we consider a complex generalization of coupled Josephson junctions (cf.Eq. (3.16) in [2]), where Z 1 is the generalized Kuramoto order parameter (11).We use N = 8 units with initial conditions: {x j (0)} = {−i sin( π N j) exp(i 2π N j)}, j = 1, ..., N .Just like in the case of pure phase oscillators, for this parameter choice we observe chaos, see Fig. 3, cf.Fig. 4(c) in [2] and Fig. 2 in [21].16) with complex initial conditions inside the unit disk (colored points) generalize pure phase oscillators bound on the unit circle, see Fig. 4(c) in [2] and Fig. 2 in [21].Trajectory of one oscillator in the complex plane is shown in orange.Small blue scatter points depict the value of Z1 on the Poincaré section where the oscillator on the unit circle passes phase π/2.

Discussion.
Our study presents a novel lowdimensional description that generalizes the well-known WS [1,2] theory to arbitrary arrays of complex Riccati equations, and so includes the real QIF model as a limiting case.This exact formalism enables the consideration of a whole new class of complex oscillatory models, including a special case of phase-amplitude oscillators.To showcase its correctness and applicability, we provide numerical simulations for several interesting examples.
The new formalism is effectively six dimensional, as compared to the three dimensional WS theory -this is not surprising as the generalization from real R to complex C variables implies twice the dimensionality.However, it should be noted that the generalization does not simply involve allowing the dynamical variables to take complex values; rather, the equations we obtain are fundamentally different from those described by WS theory.What is common with the WS theory is that the overarching motif is the Möbius transform between initial values of the dynamical variables x j and the constants ξ j .As was explored later [4], cross-ratios of dynamical variables C j = are invariant under the Möbius transform and thus constants of motion.Where in the WS context these cross-ratios are real (even though x j = e iφj ∈ C), here they can be complex: C j ∈ C. Just as was shown in [4], there are N − 3 independent ratios C j , j = 1, ..., N − 3. Since the initial problem contains N complex variables (1), and there are N − 3 complex constants of motion, this leaves three complex variables to describe dynamics (3), thus confirming that the description is six dimensional.
For real-phase oscillators in the thermodynamic limit the inclusion of heterogeneity [5,22] or noise [23] leads to terms of complex-valued effective frequencies.Our description (3) is clearly applicable to describe those scenarios as well, as we have shown in previous work [18,24].This remarkable equivalence between noise and heterogeneity, and complex valued frequencies can be studied further with our approach.
The generalization to complex numbers is substantial and provides room for qualitatively different dynamics.A complex extension of the Kuramoto model has recently been explored in [25]; but our description applies to a much broader class of coupled systems described by complex Riccati arrays defined in Eqs.(1).This unlocks a whole spectrum of 2D dynamical systems, including a special case of phase-amplitude oscillators (15) we showcased here.The examples considered here are simply complex generalizations of known models: Example (II) a generalization of QIF and Example (III) a generalization of phase oscillators.One can explore systems that go beyond just generalizing known models to complex initial conditions, and really consider complex units with complex coupling, thus tapping into the rich 2D dynamics of intrinsic units defined by Riccati equations (1).
Moreover, the new description provides a fresh perspective and insights on the relationship between phase oscillators and QIF neurons.In fact, the voltage resetting in QIF neurons arises naturally in our framework, providing additional motivation for its use in studying neuronal dynamics and development of new models.
Our findings have significant implications for both theoretical and practical research.The new description opens up many avenues for investigating the dynamics of complex systems.Several ideas for future work directly follow from our framework, and more are expected to arise from the research community.We briefly outline three ideas here.(A) It is well known that the general Riccati equation can be transformed into a linear second order ODE [26], which means that our formalism applies there as well.A more detailed study will be performed in a forthcoming work.(B) Most likely similar descriptions exist for higher dimensional systems as well; both for systems with larger spatial dimensionality [27], x j ∈ R n , n > 1, such as the generalizations to vectors [28] and matrices [11], as well as allowing for states x j that belong to higher number systems, such as quaternions or octonions.(C) For the special case of phase oscillators and QIF neurons, in the thermodynamic limit N → ∞ a particular initial state leads to an additional dynamical reduction as described by Ott and Antonsen [5].One would expect that similar reductions can be found for the more general case of a, b, c, x j ∈ C. It would be interesting to know which states allow for additional reductions and what are the reduced dynamics.(D) Throughout this work we strictly considered identical oscillators, i.e., at all times every oscillator felt the same global forcings exerted by a, b, c.Thus, the oscillators only differed by their states x j which are determined by the initial conditions.However, one may consider adding heterogeneity in the forces by assuming that a, b, c in some way differ between oscillators.Just like Ott and Antonsen incorporated Cauchy-Lorentzian inhomogeneities in frequencies [5], one can add them to our formalism in the thermodynamic limit as well.It is even likely that particular heterogeneity can be incorporated into finite arrays.

FIG. 1 .
FIG. 1. N = 8 quadratic integrate-and-fire neurons with global coupling interacting via Gaussian pulses.The dynamics is exactly described by the low dimensional Eqs.(9).Left panel: time series of the input current I(t), alongside the individual neurons' firing events (top).Right panel: trajectory of the macroscopic variable Q(t) determining voltages via (10).