Mean-field approximation for networks with synchrony-driven adaptive coupling

A population of coupled neurons with synaptic plasticity may be modeled as an adaptive dynamical network consisting of a very large number of excitable nodes and/or oscillators. One approach to circumvent the complexity of the full network is to consider the mean-field dynamics. Here, we introduce a mean-field model for a population of $θ$-neurons that are adaptively coupled according to a phase-difference-dependent plasticity rule. In this context, the mean-field dynamics describe important properties of neural synchrony and conductance. We analyze the validity of the model and discover a rich tapestry of dynamics, including multiple routes to chaos, induced by the adaptive coupling used to represent synaptic plasticity.

Synaptic plasticity is the primary mechanism through which the brain learns, stores information, and adapts to new experiences. Extensive experimental evidence demonstrates that the strength of synaptic connections can be dynamically enhanced or diminished in response to experience and neural activity.^1–9 Understanding synaptic plasticity is essential for unraveling the neural basis of cognitive functions and developing effective therapeutic interventions for neurological disorders. Both invasive and non-invasive brain stimulation techniques have been shown to induce lasting modifications in neural circuits.¹⁰ Among these, transcranial magnetic stimulation (TMS) has demonstrated promising cognitive and behavioral benefits in post-stroke rehabilitation¹¹ and in the treatment of depression.¹² However, the precise mechanisms by which TMS elicits these effects and the optimization of its application remain areas of active investigation.

Mathematical models of synaptic plasticity date back to Donald Hebb, who in 1949 proposed that the repeated and persistent co-activation of two or more neurons strengthens their synaptic connections.¹³ Since then, a plethora of mathematical models have been developed to describe the mechanisms underlying the strengthening and weakening of synaptic connections (see review by Taherkhani et al.¹⁴). Of particular relevance for this study are adaptive network models, which serve as valuable proxies for understanding synaptic plasticity in the brain.¹⁵ In general, network adaptivity enriches the repertoire of dynamical behaviors exhibited by networks of oscillators, enabling phenomena such as multistability,¹⁶ frequency clustering,¹⁷ as well as the coexistence of incoherence and synchronization.¹⁸ Adaptive neuronal oscillator models typically employ phase-difference-dependent plasticity (PDDP), wherein coupling strengths are updated based on the similarity or dissimilarity of the oscillators’ phases. Duchet et al. recently demonstrated that in the network of Kuramoto oscillators, PDDP rules can be effectively equated with spike-timing-dependent plasticity rules.¹⁹ They also showed that a symmetric PDDP rule permits the use of the Ott–Antonsen Ansatz²⁰ to derive a low-dimensional mean-field description of the system. In this study, we extend the work of Duchet et al. to the more biologically realistic $θ$-neuron model with first order synaptic dynamics.

The $θ$-neuron model is known to admit an exact mean-field reduction.²¹ When synaptic dynamics are incorporated, this reduction yields a neural mass model driven by a dynamical variable representing population synchrony.^22–24 Neural mass models, also known as firing rate models, are commonly used to study coarse-grained neural activity. Traditionally, these models are phenomenological, offering limited connections to the underlying biological mechanisms. In contrast, this new form of the neural mass model, derived from a population of $θ$-neurons and dubbed the next generation neural mass model, provides an explicit link between the underlying spiking network and the measurable population-level activity. Similarly, Montbrió et al. showed that the quadratic integrate-and-fire (QIF) model also admits an exact low dimensional description and derived a conformal mapping that links their mean-field model to the next generation neural mass model.²⁵ 

As the $θ$-neuron and QIF models include accurate descriptions of the synaptic interactions between cortical cells, it is possible to include plasticity, and other neural mechanisms, in a biologically plausible manner. Laing added a spatial dependence to the neurons’ activity to derive an exact neural field model,²⁶ and later he extended this work to include and study the effect of gap junction coupling.²⁷ Taher et al. extended the model to include short-term synaptic facilitation and depression, exploring their role in working memory.²⁸ Similarly, Gast et al. incorporated short-term adaptation mechanisms to investigate bursting dynamics.^29,30

In Sec. II, we present the mean-field model for a population of $θ$-neurons incorporating a phase-difference-dependent plasticity (PDDP) rule to represent adaptive coupling strengths within the network. The accuracy of the resulting model is validated in Sec. III through comparisons with corresponding full-network dynamics. A bifurcation analysis of the mean-field model is conducted in Sec. IV, followed by a broader discussion of our findings in Sec. V.

The mean-field dynamics are given by (4)–(6) with $k= k ^$⁠. We note that (6) was derived previously by Duchet et al. but coupled to the Kuramoto model. Here, we combine it with the more biologically plausible $θ$-neuron network with first order synaptic coupling. The increased complexity of the $θ$-neuron model introduces a richer bifurcation structure and a wider range of dynamical behaviors compared to simpler oscillator models. These dynamics are explored in Sec. IV.

We study the $θ$-neuron network model (1) and (2), where the pairwise coupling strengths are updated according to (3) (i.e., network with pairwise or “local” updates). We compare these results with the same $θ$-neuron network where all coupling strengths are now updated using the mean evolution equation (6) (i.e., network with global updates), as well as the mean-field system (4)–(6) with $k= k ^$⁠. Figure 1 displays results from numerical simulations of a network of 1000 $θ$-neurons (see  Appendix B for computational details). Each column corresponds to a different value of $η 0$⁠, showcasing regimes with qualitatively distinct dynamical behaviors. The green, blue, and pink curves show the evolution of the dynamic variables for the network with pairwise updates, the network with global updates, and the mean-field model, respectively. The top row shows the Kuramoto order parameter (⁠ $z$⁠). The second, third, and fourth rows show the time courses for the mean coupling strength (⁠ $k ^$⁠), network synchrony (⁠ $R= |z |$⁠), and mean conductance (⁠ $s$⁠), respectively. The bottom row shows the coupling weight distributions of the network with pairwise updates.

For fixed point dynamics, the mean-field model is a good approximation for both the network with local updates and the network with global updates, especially when the fixed point lies close to the boundary $|z |=1$ [Fig. 1(a)]. When $|z |=1$⁠, the neurons are perfectly synchronized and the phase differences are zero. While for $|z |≈1$⁠, the heterogeneity in the phase differences is minimal and the distribution of coupling weights has a single narrow peak around the mean coupling strength of $∼1.9$ [Fig. 1(a)(v)]. A finite size analysis, conducted for $N=500$⁠, $1000$⁠, $2000$⁠, $5000$⁠, confirms that the distribution of coupling weights converges to a unimodal distribution, which peaks at $∼1.9$ as $N$ increases (see supplementary material, Fig. S1). This result highlights that the global coupling strength update rule closely approximates the local pairwise coupling strength update rule in this regime. Although the population-level dynamics settle to a fixed point, the neurons are not necessarily quiescent. As shown by Byrne et al.,²² the firing rate is positive for all $|z |<1$ and increases as the dynamics approach $z= e i π$⁠, where firing is maximal as all neurons cross threshold simultaneously (see Fig. 8 of Byrne et al.²²). In the regime where $|z |≈1$ and $ϕ≠π$ the majority of neurons are synchronized at rest, and as such, this fixed point corresponds to a population of predominantly quiescent neurons.

As the fixed point moves away from the $|z |=1$ boundary, the firing rate increases, and the population of neurons transitions to an active state. This shift away from $|z |=1$ also reduces the level of synchrony, introducing greater heterogeneity in the phases. Consequently, the dynamics of the networks with local and global updates begin to diverge [Fig. 1(b)]. Despite this divergence, networks with global updates still closely match the mean-field dynamics. The coupling weight distribution starts to become bimodal [Fig. 1(b)(v)], which suggests that the oscillators may be forming clusters. The finite size analysis confirms that, as $N$ increases, the distribution converges towards a bimodal distribution with peaks at $∼0.3$ and $∼1.9$ (see supplementary material Fig. S1). Although the mean coupling strength is higher in the network with local updates [Fig. 1(b)(ii)], the Kuramoto order parameter still fluctuates around the mean-field fixed point value. This is particularly apparent in the synchrony time course [Fig. 1(b)(iii)].

For limit cycle dynamics, the mean-field model is a good approximation for the network with global updates only [Fig. 1(c)]. Although the trajectories do not match the network with local updates, they lie in roughly the same region of the Kuramoto order parameter phase plane [Fig. 1(c)(i)]. The mean-field model slightly underestimates the mean coupling strength for the network with global updates and the oscillatory amplitude of the conductance, while slightly overestimating the oscillatory amplitude of synchrony. Despite these discrepancies, the dynamics remain qualitatively similar. As in the previous regime, the coupling weight distribution appears bimodal [Fig. 1(c)(v)] and a finite size analysis confirms that this bimodal distribution persists and stabilizes as $N$ increases (see supplementary material Fig. S1).

We note that although the dynamics for the network with global updates and the mean-field model are qualitatively similar across all three regimes, they are not exactly matched. The primary source of this mismatch is the variability in the background drives $η j$⁠, which are randomly sampled from a Lorentzian distribution. To illustrate this variability, we present the time course of the mean coupling strength across 15 different realizations of the network simulations (Fig. 2) at three network sizes, namely, $N=500$⁠, 1000, and $2000$⁠. Zoomed-in versions of these time courses, highlighting the long-term behavior of the system, can be found in the supplementary material (Fig. S2). As the network size increases, the variability across realizations decreases, and the dynamics of the network with global updates converges toward the mean-field dynamics (pink). We see a similar trend for both synchrony $R$ and mean conductance $s$ (see supplementary material, Figs. S3 and S4).

The variability can be reduced by decreasing the heterogeneity of the individual neurons, i.e., reducing $Δ$⁠. However, modifying the heterogeneity also changes the qualitative behavior of the system. A more homogeneous system (small $Δ$⁠) is more likely to synchronize and if operating in the active state (as is the case for $η 0=10$⁠), oscillations will emerge. As reported in Fig. 2, oscillatory states exhibit greater variability across realizations compared to non-oscillatory states. In contrast to the findings of Gast et al., who reported that a small adaptation time rate (analogous to the learning rate $ϵ$ in our model) improved the agreement between the network dynamics and the mean-field model dynamics for a network of quadratic integrate-and-fire neurons with spike-frequency adaptation,³⁰ we found that for our synchrony driven plasticity rule, changing $ϵ$ did not significantly affect the goodness of the fit (see supplementary material, Figs. S5, S6, and S7). These findings support our claim that, for a fixed network size, the goodness of the fit is determined by the type of dynamics (i.e., node, spiral, and limit cycle), rather than specific parameter values.

Interestingly, when comparing the dynamics of the mean-field system to the network with global updates, we observed transitions in the network simulations away from the mean-field dynamics (Fig. 3). This behavior indicates the presence of another stable attractor, which was confirmed through bifurcation analysis (see Sec. IV). The transition observed in Fig. 3 suggests that the finite size fluctuations in the network model can push the system out of the basin of attraction of the limit cycle attractor and into the basin of attraction of the stable fixed point where the network is more strongly coupled (light blue curve). Increasing the network size from 500 to 1000 reduces the finite size fluctuations, making it more likely for the network dynamics to remain within the basin of attraction of the limit cycle attractor (dark blue curve). For smaller networks (⁠ $N=500$⁠), transitions to the stable fixed point occurred in approximately $80%$ of realizations [Fig. 3(b)], while for larger networks (⁠ $N=1000$⁠), transitions only occurred for $∼30%$ of realizations [Fig. 3(c)]. This behavior can be attributed to the long tail of the Lorentzian distribution, which ensures a non-zero probability of extremely high or low values of $η j$⁠. In larger networks, the impact of these outlier neurons on the overall dynamics is minimal. However, in smaller networks, these neurons have more of an effect. To explore the potential for multistability and to better understand the complex dynamics of the system, we perform a bifurcation analysis in Sec. IV.

We perform a preliminary bifurcation analysis of the mean-field approximation of the network with PDDP (4)–(6) by means of numerical simulation and continuation. The goal of this analysis is to highlight the richness of the dynamics that emerge when the coupling is made adaptive. Specifically, we focus on the dynamics in terms of the mean background drive, $η 0$⁠, which is the bifurcation parameter commonly used in previous studies of the model without PDDP, as well as the parameter $α$⁠, which controls the strength of adaptiveness (i.e., how strongly the coupling strength depends on the state of the system). While this analysis is not exhaustive, it serves as a starting point, with further avenues for detailed exploration discussed in Sec. V. For the numerical analysis, we employ the fourth-order Runge–Kutta method to identify stable behaviors, including fixed points, periodic, quasi-periodic, and chaotic attractors, and generally exclude transient dynamics unless otherwise specified. We use the matlab-based continuation package MatCont³² to obtain unstable fixed point solutions and periodic orbits, as well as to calculate their stability properties and identify bifurcations.

To begin, we present a bifurcation diagram of the model without PDDP (i.e., $ϵ=0$ and fixed $k=1$⁠) to demonstrate the typical behavior of the model with fixed coupling. Figure 4(a) shows a thin curve representing fixed points, in terms of synaptic conductance $s$⁠, that are either stable (blue) with zero unstable eigendirections or saddle (red) with two unstable eigendirections. The thick blue curve shows the maximum $s$ value of stable periodic orbits emerging from the Hopf bifurcation, represented by the magenta circle. This bifurcation diagram is consistent with previous studies that have shown the system undergoes a supercritical Hopf bifurcation under appropriate parameter conditions.^21,23 This relatively simple behavior contrasts with the complex dynamics that arise when adaptive plasticity is incorporated into the model.

Figure 4(b) shows the bifurcation diagram for the same mean-field model, now with adaptive coupling (⁠ $α=30$ and $ϵ=0.5$⁠). The thin green curve represents fixed points with one unstable eigendirection. The thick green and red curves represent saddle periodic orbits with one and two unstable Floquet multipliers, respectively. Blue, black, and red circles indicate saddle-node, period-doubling, and torus bifurcations. Panel (c) shows the same bifurcation diagram in terms of mean coupling strength $k ^$⁠.

When the mean background drive is zero (⁠ $η 0=0$⁠), the only attractor is a fixed point representing a strongly coupled network. As $η 0$ increases, the coupling strength decreases until the branch of equilibria loses stability at a saddle-node bifurcation at $η 0≈10$⁠. A branch of saddle fixed points leads to stable fixed points, representing a weakly coupled network, via another saddle-node bifurcation. The lower branch of stable fixed points is short-lived since a supercritical Hopf bifurcation is encountered, producing a branch of stable periodic orbits. The periodic orbits become unstable twice as $η 0$ increases: they gain one unstable Floquet multiplier between a pair of period-doubling bifurcations and two unstable Floquet multipliers between a pair of torus bifurcations.

The numerous bifurcations evident in Figs. 4(b) and 4(c) spurred an investigation of the bifurcation set in the (⁠ $η 0$⁠, $α$⁠)-plane. We identified $α$ as a suitable second bifurcation parameter, since it relates to the effect of the adaptive coupling, as it represents an upper bound on the coupling strength. Figure 5(a) shows the saddle-node bifurcations of fixed points, observed in Figs. 4(b) and 4(c), that meet at a cusp bifurcation. Inside the two curves exist a region of bistability. For the parameter values considered here, the upper branch of fixed points remains stable, while the lower branch undergoes numerous bifurcations. Figure 5(b) shows the bifurcations that emerge from the lower branch in Figs. 4(b) and 4(c) continued in the (⁠ $η 0,α$⁠)-plane. In other words, the Hopf (magenta), period-doubling (black) and torus (red) bifurcation points in Figs. 4(b) and 4(c) now correspond to curves in Fig. 5(b). The dashed gray line indicates the value of $α$ used in Fig. 4(c). In addition, we have blue curves representing saddle-node bifurcations of periodic orbits, additional period-doubling curves in black, and shaded regions of the plane displaying maximal Lyapunov exponents (see  Appendix C for computation details). For clarity, we only show the positive maximal Lyapunov exponents.

Figure 5(b) provides valuable insights into the behavior of the system. First, we see that the pair of period-doubling bifurcations shown in Fig. 4(c) belong to the same bifurcation curve in Fig. 5(b). We show further period-doubling bifurcation curves to demonstrate the period-doubling cascade that gives rise to a region of chaotic behavior. We find that solutions within this region possess positive maximal Lyapunov exponents, providing evidence that these solutions are indeed chaotic. Note that the maximal Lyapunov exponents are calculated by incrementally increasing $α$ for fixed values of $η 0$ so it is plausible that some positive maximal Lyapunov exponents are overlooked due to multistability.

A more detailed bifurcation diagram in $η 0$ for fixed $α=30$⁠, which includes the period-doubling route to chaos, is shown in Fig. 6. The solutions with a positive maximal Lyapunov exponent are highlighted yellow, revealing the branch of chaotic attractors between the sets of period-doubling cascades.

Returning to Fig. 5(b), as $η 0$ is increased further, the periodic orbit born at the Hopf bifurcation loses stability at a torus bifurcation. Our simulations reveal that small invariant tori emerge beyond the torus bifurcation curve, where the periodic orbit is unstable. These tori persist in long simulation runs, suggesting that the torus bifurcations are supercritical (for the range of parameter values considered here). Aside from quasiperiodic solutions on the invariant tori, they will also have associated $p : q$ Arnold tongues, as labeled in Fig. 5(b). Generally, where the underlying torus is stable, the Arnold tongue will contain a pair of stable and saddle $p : q$ periodic orbits. Boundaries of the tongue are given by saddle-node bifurcations where these two periodic orbits meet and disappear. Here, we show three example Arnold tongues. In accordance with theory,^33,34 the $1 : 5$ Arnold tongue is rooted at the point on the torus bifurcation curve where its rotation number (or winding number) becomes $1/5$⁠, indicated by the red circle. This means that two curves of saddle-node bifurcations, which give the boundaries of the Arnold tongue, approach the torus bifurcation curve from the same side and meet (see the inset in Fig. 7 for a clearer depiction). In theory, there exist an infinite number of Arnold tongues, each rooted at an associated point along the torus bifurcation curve with a rational rotation number. The boundaries of the $1 : 3$ and $1 : 4$ Arnold tongues, on the other hand, do not appear to be rooted at their associated resonance points along the torus bifurcation curve. This suggests that these Arnold tongues possess a more complicated geometry, as discussed in Sec. V.

The maximal Lyapunov exponents in Fig. 5(b) also indicate regions of chaotic behavior due to dynamics on tori for larger values of $η 0$⁠. The emergence of chaotic regimes, resulting from the overlapping of Arnold tongues, has been observed in many systems throughout the literature; for example Refs. 35, 36, and 37. The chaotic attractor may appear and disappear via multiple routes. For example, in Fig. 7, we show three period-doubling curves that exist within the $1 : 5$ Arnold tongue, representing a cascade of period-doubling bifurcations that leads to chaotic behavior. In some cases, the disappearance of chaotic attractors appears to coincide with saddle-node bifurcations of periodic orbits. In such cases, our simulations reveal very long chaotic transients. This suggests that chaotic attractors may disappear at boundary crises, where a chaotic attractor meets the stable manifold of a saddle invariant set.³⁸ To provide further evidence of this, Fig. 8 displays the results of 3500 simulations for parameter values near where the chaotic attractor disappears (⁠ $η 0=8.94$ and $α=30$⁠; black cross in Fig. 7) using randomly chosen initial conditions. They are chosen from the former basin of attraction of the chaotic attractor. The black circles on the plot show the fraction of the 3500 simulations that remain on the ghost of the chaotic attractor as a function of time (see  Appendix D for computation details). The horizontal axis is scaled by the average chaotic transient time $τ ¯=1920$ and the vertical axis has a natural logarithmic scaling. We see that all trajectories eventually leave the ghost of the chaotic attractor. They all settle to the stable $1 : 5$ periodic orbit. The black circles are in good agreement with the gray curve, which indicates the theoretical (exponential) probability distribution of chaotic transient lengths in the vicinity of a boundary crisis.³⁸ 

In this work, we derive a low-dimensional description for a network of $θ$-neurons with first order synapses and adaptive coupling. We compare the dynamics of the reduced system to the $θ$-neuron network with both local pairwise coupling updates and with global population-wide coupling updates. Our results show that the reduced system accurately describes the evolution of the synchrony and the mean coupling strength for the global updates. For the local update rule, the reduced system only reproduces the network dynamics when the neurons do not cluster, and therefore the coupling weights form a unimodal distribution. The accuracy of the mean-field model is further influenced by the narrowness of this distribution, due to the $k= k ^$ assumption.

We also present a preliminary bifurcation analysis of the mean-field model to highlight the complex dynamics that arise once the coupling becomes adaptive. Figure 5(b) clearly shows that the more interesting behavior, compared to the fixed coupling case, only occurs when $α$ is sufficiently large. For smaller values of $α$⁠, only a Hopf bifurcation is encountered as $η 0$ is varied, such that the behavior is similar to the fixed coupling case. For larger values of $α$⁠, the periodic orbit created at the Hopf bifurcation undergoes further bifurcations, leading to two distinct regimes of chaotic behavior. One is due to a period-doubling cascade. The other is due to overlapping Arnold tongues created along the torus bifurcation curve.

Phase-difference-dependent plasticity (PDDP) rules are commonly used in adaptive network models as a proxy for spike-time-dependent plasticity (STDP).^16,31,39,40 Duchet et al. studied the validity of this approximation, constructing PDDP rules to fit both symmetric and asymmetric STDP rules, showing that the resulting dynamics and coupling weights closely match.¹⁹ In making the mean-field reduction, any notion of spike times is lost, and we are forced to rely on synchrony as a measure of how similar or dissimilar the activities are. For single-cluster states, where neurons exhibit similar activities, the reduced model accurately describes networks with STDP-like plasticity. If the plasticity induces clustering, population-level synchrony is no longer an accurate surrogate for spike-time difference or phase difference. However, the framework presented here is readily extendable to study multi-cluster states via the introduction of multiple populations and/or higher-order Kuramoto–Daido order parameters.

Although STDP has dominated the synaptic plasticity literature, particularly in the field of mathematical and computational neuroscience, it is not the only mechanism for activity dependent changes in synaptic coupling strengths.⁴¹ Synaptic scaling, a global plasticity mechanism, strengthens and weakens all synapses by the same factor. This homeostatic plasticity mechanism acts to stabilize and regulate the network activity, to ensure normal propagation of signals through the network.^42,43 Our global update rule, which takes the form of a synaptic scaling law, regulates the network activity based on the population synchrony. In its current form, the plasticity rule promotes synchrony, increasing the coupling strength when the neurons are synchronized, but one could easily promote asynchrony or partially synchronized states by adding a phase shift to the original PDDP rule. Previous studies have examined global update rules based on synchrony in Kuramoto oscillator networks,^18,44 but the coupling update in these studies is linear in synchrony, whereas our model has a quadratic dependence. The novelty of our approach lies in deriving this global update rule from a pairwise rule for adaptive oscillators.^15,16,31 It is not clear what form of pairwise rule would provide a linear dependence on synchrony. Comparing the network dynamics for the linear update rule with the dynamics for our quadratic update rule would be an interesting follow-up study.

We note that the PDDP rule considered here pays no recourse to the actual firing time and instead updates the coupling strengths based on the phase difference at all times. Duchet et al. showed that this simple version of the PDDP rule provided a reasonable approximation of the network and coupling weight dynamics for a network of Kuramoto oscillators with STDP.¹⁹ The authors also introduced an event-based PDDP rule whereby the coupling weights are only updated at the spiking phases. Since no distinctly different dynamics were observed with the event-based rule, we opted for the simpler version here, though the framework could easily be extended to incorporate the event-based rule. Another interesting avenue for future work is the plasticity rule proposed by Bergoin et al.⁴⁵ The authors proposed a $2π$-periodic exponential rule with asymmetric potentiation and depression phase intervals. Using the framework of Duchet et al.,¹⁹ one could compute the Fourier expansion of this rule and construct an approximate mean-field description for the network dynamics.

Our bifurcation analysis reveals similarities to the dynamics observed by So and Barreto⁴⁶ in a model of bimodal Kuramoto oscillators with periodic coupling in the thermodynamic limit. In contrast to the PDDP rule considered here, the time-dependent coupling strength considered in Ref. 46 was sinusoidal with an enforced timescale. In other words, the model was non-autonomous with periodic parametric forcing. For certain timescales of the coupling, the authors identified complicated dynamics, including period-doubling routes to chaos and attractor crises, which have also been found here.

There are several interesting avenues for future work on the analysis of the mean-field model considered here. First, our bifurcation analysis focuses on two parameters: namely, the intrinsic excitabilities, $η 0$⁠, and the maximum coupling strength, $α$⁠, which also relates coupling strength to the state of the system. This choice of parameters is made to directly compare the behavior that emerges when network plasticity is included with the behavior in the non-adaptive model, which has been studied in the past.^21,23 However, the emergence of complex dynamics in the adaptive system may also depend on other parameters, such as the timescale of plasticity, $ϵ$⁠, and $Δ$⁠, which encodes the heterogeneity of the intrinsic excitabilities of the neurons in the network. Previous work has highlighted the key role that $Δ$ plays in controlling system dynamics (e.g., Ref. 25). In fact, in bifurcation diagrams similar to Fig. 4, we find evidence of additional Hopf bifurcations due to Bogdanov–Takens bifurcations, when $Δ$ is varied. However, at this point, we leave a fully comprehensive bifurcation analysis of the model, which could also benefit from nondimensionalization, for a thorough treatment in the future, since this would be a significant body of work.

The preliminary bifurcation analysis above already hints at intriguing details to be studied in the future. It remains unclear how the $1 : 3$ and $1 : 4$ Arnold tongues shown in Fig. 5(b) emerge from the torus bifurcation curve. We detect curves of torus bifurcations and neutral saddles within these Arnold tongues that could belong to so-called Chenciner bubbles, suggesting that the Arnold tongues attach to the torus bifurcation curve via folds of tori that connect regimes of stable and saddle tori. To show this clearly would require a careful and intricate numerical analysis, as demonstrated by Keane and Krauskopf.⁴⁷ 

In Sec. IV, we provide evidence of a boundary crisis route to chaos (Fig. 8). To clearly visualize such crises in this system would require the calculation of higher-dimensional saddle manifolds using advanced numerical techniques; for example, the work of Hasan et al.⁴⁸ and Dankowicz et al.⁴⁹ A deeper appreciation of the role of saddle manifolds in this system could also be beneficial for understanding the limitations of the mean-field model. In particular, if an attractor of the mean-field model is near the boundary of its basin of attraction, then this mean-field solution may not correspond to a solution in the full-network system, since the variability around its mean state may force the system into a different basin of attraction. For example, in Fig. 3, we observe that the simulation of the mean-field model follows a periodic solution, while the full-network simulation appears to approximately follow the periodic solution for some time before jumping to a different attractor with a larger coupling strength. Figure 9 shows the same phase space projection as Fig. 3(a), containing the relevant invariant sets of the mean-field model. The thick blue curve represents the stable periodic orbit, the blue/green circle represents the stable/saddle fixed point. The Jacobian of the model evaluated at the saddle has the following eigenvalues: $(−3.016,−0.5240±2.224i,0.5915)$⁠. As such, the saddle has a one-dimensional unstable manifold and a three-dimensional stable manifold. In Fig. 9, the thin red curve represents the unstable manifold of the saddle fixed point, found by numerically integrating the model forward in time from the saddle. The thin blue curve is an approximation of the strong stable manifold, found by numerically integrating the model backward in time from the saddle. Since the boundary between the basins of attraction is given by the stable manifold of the saddle, Fig. 9 paints a convincing picture that the jump observed in Fig. 3 is due to the natural variability of the full-network model pushing the trajectory from the periodic orbit beyond the stable manifold of the saddle and, therefore, into the basin of attraction of the stable fixed point. However, the picture in Fig. 9 is somewhat misleading since the stable manifold of the saddle is three-dimensional, and it only shows the one-dimensional strong stable manifold. Also, it cannot be ruled out that the close vicinity of the periodic orbit to the strong stable manifold of the saddle in Fig. 9 is due to a fortunate projection from four dimensions onto two dimensions. Hence, providing a clear and precise picture of how trajectories interact with saddle manifolds would require their calculation in higher dimensions using more advanced techniques.

Building on the phase space analysis of mean-field models, Klinshov et al. showed that for a non-adaptive network of quadratic integrate-and-fire neurons, the finite size fluctuations can be included in the mean-field description as an additional shot noise term.^50,51 They observed noise-induced switching between attractors in the resulting stochastic model, similar to the transitions observed in Fig. 3. Their framework for calculating state lifetimes could be useful for understanding attractor transitions in our model.

Finally, it could be worthwhile to consider smaller values of $ϵ$⁠, representing a slower adaptive coupling rate. For sufficiently small $ϵ$⁠, a multiple time-scale approach using geometric singular perturbation theory^52,53 could provide a new perspective on the mean-field model, as demonstrated by Ciszak et al.⁴⁴ for Kuramoto models. It is likely that such an analysis could provide another perspective on the behavior of our mean-field model.

See the supplementary material for additional figures showing the distribution of coupling weights as a function of the network size (Fig. S1), a zoomed-in version of Fig. 2 (Fig. S2), zoomed-in time courses for the network synchrony (Fig. S3) and conductance (Fig. S4) for different network sizes, and the dynamics of the mean coupling strength (Figs. S5 and S7) and the Kuramoto order parameter (Fig. S6) for various learning rates ϵ and levels of heterogeneity Δ.

We thank Cris Hasan for helpful discussions about saddle manifolds in $R 4$⁠. We also thank Andrew Flynn for advice on computing Lyapunov exponents and Ruarai Goodfellow for assistance in spotting typos in the cumbersome Jacobian matrix of the mean-field model. We are grateful for the feedback from two anonymous reviewers that improved the manuscript. This publication has emanated from research conducted with the financial support of Science Foundation Ireland under Grant No. 18/CRT/6049. For the purpose of Open Access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission. R.L. acknowledges support from the EPSRC Grant Nos. EP/V013068/1, EP/V03474X/1, and EP/Y028872/1.

The authors have no conflicts to disclose.

N. Fennelly and A. Neff are joint first authors.

N. Fennelly: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). A. Neff: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). R. Lambiotte: Project administration (equal); Supervision (equal); Writing – review & editing (equal). A. Keane: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Á. Byrne: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.

We simulate two networks of $θ$-neuron oscillators (1) and (2): one with pairwise adaptive coupling (3) and the other with global adaptive coupling (6). To address the discontinuities introduced by the Dirac delta at each spike, we employ the Euler method with a fixed time step of $dt=0.0001$ for these simulations. Numerical stability tests confirm that this step size is sufficiently small (solutions are unchanged for smaller step size). Our primary simulations focus on populations of $N=1000$ neurons, while additional simulations with $N=500$⁠, $N=2000$⁠, and $N=5000$ are conducted to assess the validity of the mean-field reduction. Each neuron’s background drive $η j$ is drawn randomly from a Lorentzian distribution with mean $η 0$ and half-width at full maximum $Δ$⁠; simulations are made reproducible by setting random seeds for each run. As the mean-field dynamics are continuous, they are computed using the adaptive step-size solver solve_ivp from SciPy’s integrate package. We use the fourth-order Runge–Kutta method with relative tolerance $10 − 6$⁠.

In this section, we provide a brief overview of how Lyapunov exponents are computed for the purpose of this paper; for details, we refer to Refs. 54 and 55. Lyapunov exponents are a very useful tool from dynamical systems theory. They measure the exponential convergence or divergence of infinitesimal perturbations. More generally, one may consider the expansion or contraction rate of an $n$-dimensional volume and compute $n$ orthogonal Lyapunov exponents, referred to as a Lyapunov spectrum. In this paper, we only consider computing the maximal (largest) Lyapunov exponent. The maximal Lyapunov exponent, $λ m$⁠, may imply that an attractor is a chaotic attractor (⁠ $λ m>0$⁠), a stable limit cycle or at a bifurcation point (⁠ $λ m=0$⁠), or a stable fixed point (⁠ $λ m<0$⁠).

In Fig. 8, we demonstrate the exponential distribution of durations of transients near the disappearance of a chaotic attractor in order to provide evidence for a boundary crisis.³⁸ In this case, we consider the behavior of the mean-field model for parameter values where there is no chaotic attractor (i.e., $η 0=8.94$ and $α=30$⁠, shown by the black cross in Fig. 7), while a chaotic attractor does exist for a slightly larger value of $η 0$ (e.g., $η 0=8.96$⁠, indicated by the positive maximal Lyapunov exponents in Fig. 7). If we let $D∈ C× R 2$ denote the basin of attraction of the chaotic attractor when $η 0=8.96$⁠, then we expect that trajectories for $η 0=8.94$ with random initial conditions within $D$ will exhibit very long chaotic transients.³⁸ In other words, the trajectories remain on the “ghost” or “remnant” of the chaotic attractor for a very long time before finally converging to an attractor.

The duration of a transient could be approximated by pre-defining a region in phase space that encloses the ghost of a chaotic attractor and then determining how long it takes for the system to leave this region. Here, we take a slightly different approach. In this case, we know that the transient is over once the system approaches a stable $1 : 5$ periodic orbit. Therefore, we approximate the duration of the transient by determining how long it takes for the system to approach the stable $1 : 5$ periodic orbit.

This is done by numerically integrating the mean-field model with a random initial condition within $D$ and recording all the local maxima of the transients. We identify maxima based on the Euclidean norm of the state of the system $x=(z,s, k ^)$⁠. We say that the trajectory has settled on the stable $1 : 5$ periodic orbit once every fifth local maximum is (sufficiently close to being) identical.

To illustrate this procedure, Fig. 10 shows the local maxima taken from an example trajectory used to construct Fig. 8. Figure 10(a) shows each local maximum, denoted with a label $i$ (i.e., $x 1$ represents the first observed local maximum, $x 2$ represents the second, etc.). Panel (b) shows the difference between every fifth local maximum. The transient is considered over once the values shown in panel (b) converge to zero with a tolerance of 0.001.