Università degli studi di Firenze Dottorato in dinamica non lineare e sistemi

complessi XXV ciclo

Collective Dynamics in Complex Neural Networks

Simona Olmi

Tutore

dott. Antonio Politi

Cotutore

dott. Alessandro Torcini

Coordinatore

prof. Stefano Ruffo

Settore disciplinare FIS/03

2

Contents

Introduction 1

1 Excitatory globally coupled neural networks 6 1.1 The Network Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2 Collective Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.2.1 The Splay State . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.2 Partial Synchronization . . . . . . . . . . . . . . . . . . . . . . . 10

1.3 Network Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.1 Linear Stability Analysis of the Splay State . . . . . . . . . . . 11 1.3.2 Bifurcation Characterization . . . . . . . . . . . . . . . . . . . . 13

1.4 Examined Topologies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2 Stability of the splay state in fully coupled networks 16 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2 Model and Main Results . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3 Event Driven Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4 Splay State Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.4.1 Zeroth Order Approximation . . . . . . . . . . . . . . . . . . . . 24 2.4.2 First Order Approximation . . . . . . . . . . . . . . . . . . . . . 25 2.4.3 Second Order Approximation . . . . . . . . . . . . . . . . . . . 25 2.4.4 Third Order Approximation . . . . . . . . . . . . . . . . . . . . 26

2.5 Linear Stability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.5.1 Continuum Limit . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.5.2 Zeroth Order Approximation . . . . . . . . . . . . . . . . . . . . 28 2.5.3 First Order Approximation . . . . . . . . . . . . . . . . . . . . . 29 2.5.4 Second Order Approximation . . . . . . . . . . . . . . . . . . . 30

2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3 Collective oscillations in massively connected neural networks 32 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.2.1 Event-driven Map . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2.2 Linear Stability Analysis . . . . . . . . . . . . . . . . . . . . . . 35 3.3 Collective Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Lyapunov Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3.4.1 Globally Coupled Network . . . . . . . . . . . . . . . . . . . . . 40 3.4.2 Diluted Network . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4 Coherent periodic activity in Erdös-Renyi neural networks 46 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.2 Model and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4.2.1 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2.2 The Connectivity Matrix . . . . . . . . . . . . . . . . . . . . . . 49 4.2.3 Characterization of Macroscopic Attractors . . . . . . . . . . . . 50 4.2.4 Lyapunov Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.3 Phase Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.3.1 Erdös-Renyi networks with constant probability . . . . . . . . . 51 4.3.2 Erdös-Renyi networks with γ-dependent probability . . . . . . . 54

4.4 Network Homogeneity . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.5 Chaotic vs Regular Dynamics . . . . . . . . . . . . . . . . . . . . . . . 56 4.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5 Collective dynamics and extensive chaos in sparse networks 60 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.2 Coupled Maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.3 Stuart-Landau Oscillators . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.4 Leaky Integrate and Fire Neurons . . . . . . . . . . . . . . . . . . . . . 64 5.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

6 Collective chaos and chimera states in a pulse-coupled network of networks 67 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.2 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.3 Phase Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 6.4 Collective Chaos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

A Numerical algorithms 79 A.1 Simulation of Single Networks . . . . . . . . . . . . . . . . . . . . . . . 79 A.2 Simulation of Network of Networks . . . . . . . . . . . . . . . . . . . . 80

B Linear Stability Analysis of the Splay State 82 B.1 Fixed-point expansion (general case) . . . . . . . . . . . . . . . . . . . 82 B.2 Fixed-point expansion (LIF model) . . . . . . . . . . . . . . . . . . . . 84

B.2.1 Zeroth Order . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 B.2.2 From first to third order . . . . . . . . . . . . . . . . . . . . . . 86

B.2.3 Fourth Order . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 B.3 Expansion in tangent space around the fixed point . . . . . . . . . . . . 87

B.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 B.3.2 N → ∞ limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

B.4 Comparison with previous results . . . . . . . . . . . . . . . . . . . . . 91

Introduction

In several scientific fields, e.g. in physics, chemistry and biology, it is common to distinguish between a ”microscopic” and a ”macroscopic” level. Intermediate or ”mesoscopic” levels are less commonly considered, but are becoming increasingly in focus [1, 2, 3, 4]. For example, in lasers dynamics or in certain chemical reactions, spatio-temporal patterns emerge at scales much larger than the constituent atoms or molecules. The formation of such patterns, which is highly relevant to mesoscopic physics, chemistry and biology has been intensely studied [5, 6]. In the last years the investigation of mesoscopic patterns is becoming a relevant topic also in the field of neural systems [7].

Mesoscopic brain dynamics is intermediate between microscopic and macroscopic neurodynamics. The microscopic level could refer to calcium waves, ion channels or single neurons. The macroscopic scale, on the other hand, could correspond to the largest scales possible to measure with regard to brain activity. This would be the dynamics related to systems, such as cortico-thalamic, or cortico-cortical interactions, usually measured with positron emission tomography (PET), functional magnetic resonance imaging (fMRI), or other brain imaging techniques, capturing the dynamics associated with blood flows and metabolism. Also scalp electroencephalography (EEG) and magnetoencephalography (MEG) may capture this kind of dynamics [8].

Mesoscopic brain dynamics refers typically to the dynamics of neural populations, networks or columns within cortical areas. It is characterized by its high complexity, often involving oscillations of different frequencies and amplitudes, interrupted by (chaotic or pseudo-chaotic) irregular behaviours [9]. The mesoscopic brain dynamics is influenced by the activity over other scales. For example, it is often mixed with noise, generated at a microscopic level by spontaneous activity of neurons and ion channels. It is also affected by macroscopic activity, such as slow rhythms generated by cortico-thalamic circuits or neuromodulatory influx from different brain regions.

A fundamental question in neuroscience concerns the functional significance of mesoscopic brain dynamics, including the observed phase transitions between various oscillatory states and chaotic or noisy states. Most cognitive and mental functions presumably involve macroscopic brain areas or even networks of interconnected cortical areas, but the mesoscopic dynamics of such an area may still reflect some aspect or part of the mental activity. For example, different conscious states, such as sleep or alertness may result in similar type of dynamics, such as alpha or gamma waves across several areas, and detectable at mesoscopic spatial scales. Similarly, some aspects of a Figure 1: Modelization of the olfactory cortex with a simplified network model. The attractor dynamics during recognition of an unknown input pattern. After an initial chaotic-like transient, the system converges to a near-limit-cycle attractor (40 Hz), corresponding to the recognized odor. The activity of three arbitrary network nodes are plotted against each other in arbitrary units. Figure taken from [15].

macroscopic phenomenon such as face recognition or eating behaviour may be reflected into the microscopic activity of single neurons (a prominent example is represented by mirror neurons [10]).

In order to elucidate the significance of mesoscopic brain dynamics, computational approaches are used to supplement experimental methods. For example, there is strong computational, as well as experimental evidence for a population (relational) coding in cortical networks, where mesoscopic brain dynamics apparently play a functional role [11, 12]. Such a coding principle implies that information is contained not only in the activation level of individual neurons but also in the relations between the activities of distributed neurons.

Numerical investigations with cortical neural network models support the view that complex dynamics makes neural information processing more efficient, providing a fast and accurate response to external stimuli or in associative memory tasks [13, 14]. For example, with an initial chaotic-like state, sensitive to the input signal, the system can rapidly converge to a limit cycle attractor memory state, as shown both from modelization (see Fig. 1 [15]) and from data analysis of experimental results [16].

Numerical simulations show that mesoscopic network dynamics can be driven into, or out of, different oscillatory states by small changes in the densities of the ion-channels or by altering connection strengths in a network model [17]. It is demonstrated that the blocking of specific ion channels, as a possible effect of some anesthetics, can modify the global brain activity from high-frequency states (awake) to low-frequency ones (anesthetized), as apparent in the recorded and simulated EEG (See Fig. 2)

It can further be demonstrated that microscopic noise can induce global syn-Figure 2: Simulated effect of anesthetics. The dynamics of a network, where K channels are increasingly blocked by anesthetics. The two upper time series show the activity of single excitatory and inhibitory neurons, respectively, while the lower time series is the network mean activity [17]

chronous oscillations in cortical networks and shift the system dynamics from one dynamical state to another [13, 18, 19]. This in turn can change the efficiency of the information processing of the system, e.g. system performance can be maximized at an optimal noise level, analogous to the case of stochastic [20, 21] and coherence resonances [22], furthermore spontaneous activity can facilitate learning and associative memory [23].

All in all, mesoscopic neurodynamics can be seen as resulting from the dynamic balance between opposing processes at several scales, from the influx and efflux of ions along channels in the membrane, from the concurrent activity of synaptic inhibition and excitation, from the competition of population and neural dynamics etc. Such interplay between opposing processes often results in (transient or continuous) oscillatory or more complex behaviour. Indeed, brain activity is constantly changing, due to neural information processing, intrinsic fluctuations, neuromodulation, sensory input, and internal state shifts. An essential feature of mesoscopic brain dynamics is the emergence of coherent patterns of activity, appearing at the collective level of a large number of neurons. Waves of activity move across the surface of sensory cortices, with oscillations at various frequency bands. This kind of activity is often associated with mental processes and states, which can be characterized by the specific pattern involved.

A combination of mathematical analysis and computational modeling can serve as an essential complement to clinical and experimental methods in furthering our understanding of neural and mental processes. In particular, when concerning the inter-relation between structure, dynamics and function of the brain and its cognitive functions, this method may be the best way to make progress. The study of transitions among different collective states in mesoscopic brain dynamics could be one of the most fruitful approaches in this respect [24, 25, 26]. In fact, relating different spatial and temporal scales in the nervous system, and linking them to mental processes may be seen as one of the greatest challenges to modern neuroscience.

In the present dissertation we use numerical and analytical techniques, usually employed in dynamical systems analysis and statistical mechanics, in order to characterize the emergence of collective solutions in neural networks of different topologies. In the following we refer to these collective states as occurring at macroscopic scales and not at a mesoscopic level as it would seem more appropriate from the previous discussion. This because with our modelization we are not attempting to reach the largest brain scales, which refers to brain areas, but we limit our analysis to a single or two interacting neural populations. Therefore we can identify only two levels at which the dynamics takes place: the microscopic single neuron evolution and the macroscopic dynamics of the whole considered network. In particular, we consider pulse-coupled networks, where neurons are described by simplified models, which however retains the essential features of neural dynamics. Neurons are usually schematized in this dissertation in terms of the so-called leaky integrate-and fire model [27], apart in Chapter 2 where more general phase-like models are analyzed.

In more details, the first Chapter is devoted to the collective regimes spontaneously emerging in globally coupled excitatory networks. Two distinct collective states can be identified: the splay state and the partial synchronization. In the splay state single neurons display periodic motions, while the macroscopic field associated to the network activity is constant. Furthermore the motion in this case is perfectly asynchronous. On the other hand, partial synchronization is characterized by a periodic macroscopic field and a quasi-periodic motion of each membrane potential.

In Chapter 2 we analytically investigate the stability of the splay state for generic pulse-coupled phase-like neurons in a fully coupled network, for large system sizes. Single membrane potentials are continuously driven by a velocity field from reset to threshold. We develop a perturbative technique which allows determining the Floquet exponents for generic velocity fields and implement the method for a given pulse shape. We find that in the case of discontinuous velocity fields, the Floquet spectrum scales as 1/N^{2} with the number N of neurons and the stability is determined by the sign of the jump at the discontinuity. Altogether, the form of the spectrum depends on the pulse shape, but it is independent of the velocity field.

In the remaining Chapters we analyze how the macroscopic dynamics of the network is influenced by the topology. In particular, we consider random networks of different topologies, ranging from massively connected to sparse networks. The last Chapter is devoted to the study of a network of networks, composed of two interacting families of neurons.

More specifically, in Chapter 3 we investigate the onset of collective oscillations in a random diluted network of pulse-coupled leaky integrate-and-fire neurons in the presence of quenched and annealed disorder. We find that the disorder induces a weak form of chaos vanishing for diverging system sizes. In addition to this, in the thermodynamic limit, the random network dynamics reduces to that of a fully homogeneous system with a suitably scaled synaptic coupling strength. In Chapter 4, we generalize the previous study to directed Erdös-Renyi with average connectivity scaling sub-linearly with the number of neurons in the network. In particular, we try to understand if the onset of collective oscillations can be related to a minimal average connectivity. We find that a finite minimal connectivity is sufficient to observe coherent activity in systems of any size irrespective of the kind of considered network: massively connected or sparse. However, this value depends on the nature of the disorder quenched or annealed. For annealed disorder, the critical value required to observe the onset of macroscopic behaviors is noticeably smaller than for quenched networks. The microscopic and macroscopic dynamics of sparse networks is investigated in details in Chapter 5. By simulating chaotic maps, Stuart-Landau oscillators, and leaky integrate-and-fire neurons, we show that a finite connectivity (of the order of a few tens) is able to sustain a nontrivial collective dynamics even in the thermodynamic limit. In addition to this, the microscopic evolution is extensive; this property is highly nontrivial since the sparse network dynamics is non additive and the random graph cannot be approximated with the juxtaposition of almost independent sub-structures.

In the final Chapter we study the dynamics of two symmetrically coupled populations of leaky integrate-and-fire neurons. Upon varying the coupling strength, we find symmetry-breaking transitions that lead to the onset of various chimera states as well as to a new regime, where the two populations are characterized by a different degree of synchronization. Symmetric collective states of increasing dynamical complexity are also observed. The computation of the finite-amplitude Lyapunov exponent allows us to establish the degree of chaoticity of the (collective) dynamics in a finite region of the phase plane. The further numerical study of the standard Lyapunov spectrum reveals the presence of several positive exponents, indicating that the microscopic dynamics is high dimensional.

Finally, we briefly summarize the contents of this dissertation, while the algorithms implemented for the reported numerical simulations are described in Appendix A and the details of the calculations outlined in Chapter 2 are reported in Appendix B.

Chapter 1

Excitatory globally coupled neural networks

1.1 The Network Model

To analyze in detail the neural network dynamics we should take in account: (i) the topology of the connections, (ii) the synaptic plasticity, (iii) the single neuron dynamics and (iv) the presence of intrinsic noise. Such a complex mathematical model, able to treat all the main characteristics of a real network, would be extremely difficult or impossible to simulate. On the other hand, if we reduce the complex network dynamics back to few ingredients, it is possible to analyze numerically the system behaviors by using a simplified model capable to reproduce the main features of a real neural system. In the following we use a simplified single neuron model, the leaky integrate-and-fire model (LIF), in order to treat analytically and numerically the network evolution. In particular, we restrict the analysis to excitatory pulse-coupled networks. Previous numerical and analytical results demonstrated that globally coupled excitatory networks of LIF neurons show stable, collective, non trivial behaviors which emerge spontaneously [28, 29]. In a network of N neurons, the membrane potential u_{i}(t) of the i − th oscillator evolves according to the following ordinary differential equation

u̇_{i}(t) = a − u_{i}(t) + g E(t)  i = 1, ⋯, N  (1.1)

where all variables and parameters are expressed in adimensional rescaled units. According to the above equation, the membrane potential u_{i} relaxes towards the value a + gE(t), but as soon as it reaches the threshold value u_{i} = 1, it is reset to u_{i} = 0 and a spike is simultaneously sent to all neurons (see below for the spike definition). This resetting procedure is an approximate way to describe the discharge mechanism operating in real neurons. The parameter a > 1 is the suprathreshold input DC current and g > 0 gauges the synaptic coupling strength of the excitatory interaction with the neural field E(t). This field represents the synaptic current injected in each neuron and is given by the superposition of all the pulses emitted by the network in the past. By following Ref. [25], we assume that the shape of a pulse emitted at time t = 0 is given by s(t) = \frac{α^{2}t}{N}e^{−αt}, where 1/α is the pulse–width. For this choice of the pulse shape it is easy to show that the field evolution is ruled by the following second order differential equation

Ë(t) + 2αĖ(t) + α^{2}E(t) = \frac{α^{2}}{N} ∑_{n|t_{n}<t} δ(t − t_{n}) .  (1.2)

In other words, s(t) represents the stereotyped shape of a single pulse emitted at time t_{n} by a neuron reaching the threshold value. The solution E(t) for a generic time t_{n} < t < t_{n+1} between two spike emissions is the linear combination of such pulse shapes and represents a macroscopic variable reproducing the network activity.

Then the dynamics is described by N first order differential equations for the membrane potentials (1.1) and by the second order differential equation (1.2), thus totalizing N + 2 degrees of freedom. It is possible to simplify the analysis by transforming the differential equations into a discrete-time event-driven map; this can be achieved by integrating between two successive pulse emissions Eqs. (1.1,1.2) [30]. The map for the field E and the membrane potentials reads as

E(n + 1) = E(n)e^{−ατ(n)} + P (n)τ(n)e^{−ατ(n)}  (1.3)
P (n + 1) = P (n)e^{−ατ(n)} + \frac{α^{2}}{N}  (1.4)
u_{i}(n + 1) = u_{i}(n)e^{−τ(n)} + a(1 − e^{−τ(n)}) + gH(n) ,  (1.5)

where τ(n) = t_{n+1} − t_{n} is the interspike time interval and, for the sake of simplicity, we have introduced the new variable P := αE + Ė. In addition to this

H(n) = \frac{e^{−τ(n)} − e^{−ατ(n)}}{α − 1} (E(n) + \frac{P (n)}{α − 1}) − \frac{τ(n)e^{−ατ(n)}}{(α − 1)} P (n) .  (1.6)

The interspike time interval τ(n) is determined by the equation for the largest membrane potential (identified by the label q(n)) reaching the threshold value u_{q} = 1 at time t_{n+1}, as follows

τ(n) = ln [\frac{a − u_{q}(n)}{a + gH_{q}(n) − 1}] ,  (1.7)

Since the map connects network states always characterized by a neuron crossing the threshold, at least one membrane potential value is constant and we loose a degree of freedom, thus reducing the degrees of freedom of the system to N + 1. The discrete time equations ( 1.3, 1.4, 1.5) are an exact and efficient method to simulate the network dynamics. In order words the set of Eqs. (1.3,1.4,1.5,1.6, 1.7) defines a discrete-time map that is fully equivalent to the original set of ordinary differential equations. Altogether, it turns out that a network of N identical neurons is described by N + 1 equations: Two of them account for the dynamics of E(n), while the remaining N − 1 equations describe the evolution of the neurons (u_{q} is no longer a variable being, by definition, equal to 0).

In order to quantify the network degree of synchronization we use the so-called Kuramoto order parameter r defined [31] as

r(t) = \frac{1}{N} ∑_{j=1}^{N} e^{iθ_{j}(t)} = R(t)e^{iψ(t)},  θ_{j}(t) = 2π \frac{(t − t_{j,n})}{t_{q,n+1} − t_{q,n}}  j = 1, …, N (1.8)

where θ_{j} is the phase of the j-th neuron, t_{j,n} is the last time (the n-th spike) at which a spike is emitted by the j-th neuron and q identifies the index of the last spiking neuron. In the following we use the modulus of the Kuramoto order parameter R(t) ∈ [0, 1] to characterize the level of synchronization in the network. There are two limiting cases: on one hand all the oscillators phases are equispatially distributed on the unitary circle and R = 0 (asynchronous state); on the other hand all phases are equal and R = 1 (synchronized state).

1.2 Collective Solutions

In excitatory pulse-coupled LIF networks two distinct collective states can be identified: the splay state and the partial synchronization. In this section we briefly describe the peculiar characteristics of each state, starting from the microscopic point of view corresponding to the membrane potential dynamics and finishing with the macroscopic point of view by associated to the behavior of the field E.

1.2.1 The Splay State

The splay state is a collective mode emerging in non linear fully coupled oscillator networks. In this case all the oscillators evolve in the same way but their phases are

Figure 1.1: Raster plot of a network with N = 200 neurons and a = 1.3, g = 0.4 for (a) α = 3 and (b) α = 9.

equispatially distributed in the interval [0, 2π]. Each oscillator can be described by a wave-function as follows

x_{k}(t) = X(t + \frac{kT}{N}) = A cos (ωt + \frac{2πk}{N} )  k = 1, ..., N  (1.9)

where T is the collective oscillation period and ω = 2π/T. As it was previously shown in [32], in fully coupled neural networks neurons reach the threshold in an ordered manner and this order never changes in time. So, in order to characterize the network dynamics, it is convenient to order the neurons according to their potential values and then to plot the index of the neuron that reaches the threshold as a function of the emission time of the pulses (see Fig. 1.1(a)). This spike time diagram (“raster plot”) shows that the interspike interval between two consecutive neurons reaching the threshold is constant, independently of the neurons, in a splay state regime.

At a macroscopic level, the field E(t) remains constant in time, thus indicating a constant average network activity. In addition to this, in the splay state regime, the interspike interval between two successive spike emissions in the entire network is constant and equal to T/N. The j-th neuron phase is θ_{j}(t) = 2π[1 − (j − 1)/N] while the order parameter R is exactly zero, as it can be verified by inserting the previous expression for the phase in Eq. (1.8). In other words the membrane potentials are subjected to a complete asynchronous collective motion, where neurons reach the threshold, one after the other, at constant time intervals T/N. In addition to this it has been demonstrated [33] that, for an excitatory neural network, there exist a critical line α_{c}(a, g) in the parameter space (α, g) which separates the stable splay state region from the unstable region (see Fig. 1.2). In the instability region the splay state solution is replaced by another collective solution which is described in detail in the following Section.

Figure 1.2: Critical curve α_{c} in the parameter space (g, α) for a = 1.3 in the N → ∞ limit [33].

Figure 1.3: Return map for τ_{n} (black curve) with parameter values α = 3, g = 0.4 e a = 1.3; Return map for τ_{n} (red curve) with parameter values α = 9, g = 0.4 e a = 1.3 .

1.2.2 Partial Synchronization

Above the critical line α_{c} the stable collective state, termed “partial synchronization”, is characterized by a periodic field E and a quasi-periodic motion of each membrane potential. By analyzing the Fig. 1.1(b) is it possible to see that a group of neuron reaches the threshold almost simultaneously; however the neurons participating to this almost synchronized group change in time. If we plot the return map of the single neuron interspike interval^{1} (see Fig. 1.3) we find, in this regime, a closed curve, as expected for a quasi-periodic motion; for a periodic motion, as in the splay state regime, we would have expected a fixed point since the membrane potentials are periodic (as shown in Fig. 1.3). The recombination of the quasi-periodic microscopic motion into a macroscopic periodic motion is absolutely not trivial and it is still object of studying. Finally, the Kuramoto parameter R is different from zero and oscillates periodically in time with the same period as the macroscopic field E, while the oscillation period of the membrane potentials is always smaller than the macroscopic field one.

1.3 Network Dynamics

Since in networks of identical fully coupled neurons the order of the potentials u_{i} is preserved, it is convenient to redefine the membrane potentials in a comoving frame, i.e. x_{j}(n) ≡ u_{j−n}(n). In this frame, the label of the closest-to-threshold neuron is constant and can be chosen equal to 1, without prejudice of generality. The updating equation writes,

1The interspike interval is the time between subsequent action potentials (also known as spikes) of a neuron.

x_{j−1}(n + 1) = x_{j}(n)e^{−τ(n)} + 1 − x_{1}(n)e^{−τ(n)}  j = 1, …, N − 1 ,  (1.10)

with the boundary condition x_{N} = 0 and where,

τ(n) = ln [\frac{x_{1}(n) − a}{1 − gH(n) − a}] .  (1.11)

In this framework, the periodic splay state reduces to a fixed point that satisfies the following conditions,

τ(n) ≡ \frac{T}{N} ,  (1.12)
E(n) ≡ Ẽ , P (n) ≡ P̃ ,  (1.13)
x̃_{j−1} = x̃_{j}e^{−T/N} + 1 − x̃_{1}e^{−T/N} ,  (1.14)

where T is still the time elapsed between two consecutive spike emissions of the same neuron. A simple calculation yields,

P̃ = \frac{α^{2}}{N} (1 − e^{−αT/N})^{−1} , Ẽ = T P̃ (e^{αT/N} − 1)^{−1} .

The solution of Eq. (1.14) involves a geometric series that, together with the boundary condition x̃_{N} = 0 , leads to a transcendental equation for the period T . Not surprisingly, the result turns out to be independent of the pulse–width α. For simplicity, we report only the leading terms for N ≫ 1 ,

x̃_{N−j} = \frac{aT + g}{T} (1 − e^{−j T/N}) ,  (1.15)
T = ln [\frac{aT + g}{(a − 1)T + g} ] .  (1.16)

By assuming that the uncoupled neurons are in the repetitive firing regime, i.e. by taking a > 1, the period T is well defined in the excitatory case (g > 0) only for g < 1 (T → 0, when g approaches 1), while in the inhibitory case (g < 0), a meaningful solution exists for any coupling strength (T → ∞ for g → −∞).

1.3.1 Linear Stability Analysis of the Splay State

In order to perform the stability analysis of the splay state, it is necessary to linearize Eqs. (1.3,1.4,1.10) around the fixed point (1.12,1.13,1.14),

δE(n + 1) = e^{−αT/N}δE(n) + T e^{−αT/N}δP (n)  (1.17)
− (αẼ − NP̃e^{−αT/N}) δτ(n) ,
δP (n + 1) = e^{−αT/N}δP (n) − αP̃e^{−αT/N}δτ(n) ,  (1.18)
δx_{j−1}(n + 1) = e^{−T/N}[δx_{j}(n) − δx_{1}(n)] + e^{−T/N}(x̃_{1} − x̃_{j})δτ(n) , (1.19)

where we have introduced δτ(n) whose expression can be obtained by linearizing Eqs. (1.11,1.6)

δτ(n) = τ_{x}δx_{1}(n) + τ_{E}δE(n) + τ_{P} δP (n) ;

where τ_{x} := ∂τ/∂x_{1} and analogous definitions are adopted for τ_{E} and τ_{P}.

The boundary condition x_{N} ≡ 0 imposed in the comoving frame yields δx_{N} = 0 . In practice, the stability problem is solved by computing the spectrum of Floquet multipliers {μ_{k}}, k = 1, …, N + 1 associated with the eigenvalue problem of the set of linear equations (1.17,1.18,1.19). In general the explicit computation has to be performed numerically.

For the sake of clarity, it is first convenient to discuss the trivial case of vanishing interaction, i.e. g = 0 . After simple calculations the eigenvalues turn out to be μ_{k} = exp(iϕ_{k}), where ϕ_{k} = \frac{2πk}{N}, k = 1, …, N − 1, and μ_{N} = μ_{N+1} = exp(−αT/N) . The last two exponents concern the dynamics of the coupling field E(t), whose decay is ruled by the time scale α^{−1} .

As soon as the coupling is on, small amplitude fluctuations ∼ O(g/N) affect the neuron dynamics and the spectrum of Floquet multipliers takes the general form

μ_{k} = e^{iϕ_{k}}e^{T (λ_{k}+iω_{k})/N} , ϕ_{k} = \frac{2πk}{N} , k = 1, …, N − 1 ,  (1.20)
μ_{N} = e^{T (λ_{N} +iω_{N})/N} , μ_{N+1} = e^{T (λ_{N+1}+iω_{N+1})/N} .

where, λ_{k} and ω_{k} are the real and imaginary parts of the Floquet exponents. As an example, in Fig. 1.4(a) we show the spectrum of the Floquet multipliers of the splay state for excitatory coupling (g > 0) and finite values of N . The multipliers with k = 1, …, N − 1 are very close to the unit circle, while the two isolated multipliers μ_{N} and μ_{N+1} are located very close to the real axis inside the unit circle. Already for g/N ≈ O(10^{−2}) the multipliers of the coupled case can be viewed as small “perturbations” of the uncoupled one. As another example we report the Floquet exponents for the parameter values g = 0.4, a = 1.3 (see Fig. 1.4(b)-(c)). For the splay state there is a stability region in the space (g, α) for α < α_{c}. Increasing α, for α > α_{c} the maximum Floquet exponent becomes positive, thus indicating that the splay state is no more stable. In particular splay state looses stability via a Hopf bifurcation.

The variable ϕ_{k} plays the same role of the wave-number in the linear stability analysis of spatially extended systems, therefore we can say that λ_{k} characterizes the stability of the k–th mode. In the following analysis it is convenient to distinguish between modes characterized by ϕ_{k} ≈ 0, mod(2π) + O(1/N), and the other modes. Actually they identify two spectral components, that require a different mathematical treatment. The first component corresponds to the condition ||μ_{k} − 1|| ∼ N^{−1} and is referred to as long wavelengths (LWs), while the second one corresponds to ||μ_{k} − 1|| ∼ O(1) and is referred to as, short wavelengths (SWs).

The stability analysis of the splay state in networks of generic phase-like neurons will be discussed in details in Chapter 2. In particular in that Chapter we develop a perturbative technique which allows determining the Floquet exponents for a generic velocity field [34].

Figure 1.4: (a) Unitary circle and exact Floquet multipliers spectra for the complete maps for N = 20 (red circles) and N = 10 (blue crosses) for excitatory coupling. The parameters are a = 3.0, g = 0.4, and α = 30.0. (b)-(c) Floquet exponents as a function of their angular index ϕ_{k} respectively for α = 3 and α = 9. In both cases N = 200.

Figure 1.5: (a) Minima and maxima of E(t) are shown increasing α and fixing the parameter values g = 0.4, a = 1.3. (b) Amplitude δ of E(t) oscillations as a function of the distance from the bifurcation point α_{c} (magenta circles). The black line represents the function δ = κ\sqrt{α − α_{c}} in the case κ = 0.81 e α_{c} = 8.32.

1.3.2 Bifurcation Characterization

To identify the bifurcation one can observe the behavior of the field by varying the inverse of the pulse width α and by maintaining fixed g and a. At the critical point α_{c} the field E starts to oscillate, while for α < α_{c} it is constant ( see Fig. 1.5(a)). In Fig. 1.5(b) it is reported the amplitude δ of the oscillations of the field E(t) as a function of the distance |α − α_{c}| for α > α_{c}. The behavior of δ is reproducible by using a square Figure 1.6: (a) Floquet multipliers as a function of α (α > α_{c}). (b) Synchronization indicator R as a function of α. (c) Limit cycles in the phase plane (E, P) for different α values. In all cases g = 0.4, a = 1.3.

root fit κ\sqrt{α − α_{c}}. The square root scaling with α − α_{c} is a peculiar characteristic of a Hopf bifurcation . As an additional analysis of the bifurcation we consider the Floquet multipliers calculated around the critical point α_{c}. By increasing α, a pair of complex conjugated multipliers leaves the unit circle, thus indicating that the splay state is no more stable and once more indicating that the bifurcation is of the Hopf type (see Fig. 1.6(a)). In other words, approaching the critical values α_{c}, the system looses his asynchrony and gains a partial synchronization. The level of synchronization increases by increasing α and tends to be a complete synchronization in the α → ∞ limit. As shown in previous papers [28, 29], the Hopf bifurcation gives rise to collective periodic behaviour, as it can be appreciated by measuring the Kuramoto order parameter, (reported in Fig. 1.6(b)). Above the bifurcation the collective solution is represented by a periodic attractor of increasing amplitude. This behavior is absolutely non trivial since the single neurons oscillate in a quasi-periodic manner. In Fig 1.6(c) is shown the macroscopic attractor in the phase plane E, P: the curves get bigger by increasing α.

1.4 Examined Topologies

Starting from the previous results, obtained in excitatory globally coupled networks, we will analyze, in the following chapters, how the presence of disorder in the connections can influence the emergence of coherent periodic activity. In order to discuss the dependence of the dynamical behavior of the network on its topology here we specify the definition of two frequently used extreme cases [35]:

• Massively connected network. In this case the average number of synaptic inputs per neuron k_{ij} varies proportionally to the size of the system. In order to avoid that the total synaptic input to each neuron will diverge in the thermodinamic limit, the strength of the synapse connecting neuron j to neuron i should be rescaled by N as follows

G_{ij}^{syn} = \frac{g_{ij}}{N}

where the conductances g_{ij} are independent of N.

• Sparse network. In this case the average number of synaptic inputs per neurons, k_{i,j}, remains O(1) compared to the size of the network and the coupling should in principle not be rescaled in the thermodinamic limit. However we want that the dynamics of the network would have a definite value also in the limit of large k_{ij}. Therefore, also in this case, the conductances between cells are rescaled by the number of afferent synapses

G_{ij}^{syn} = \frac{g_{ij}}{k_{ij}} .

Stability of the splay state in fully coupled networks

2.1 Introduction

The first objective of a (neural) network theory is the identification of the asymptotic regimes. The previous research activity led to the discovery of fully- and partiallysynchronized states, clusters and splay or asynchronous states in pulse-coupled networks [28, 25, 35, 36]. It has also been made clear that ingredients such as disorder (diversity of the neurons and structure of the connections) are very important in determining the asymptotic behaviour, as well as the possible presence of delayed interactions and plasticity [37, 38]. However, even if one restricts the analysis to identical, globally-coupled oscillators, there are very few theoretical results: they mostly concern fully synchronized regime or specific types of neurons (e.g. the leaky integrate-and-fire model) [35, 39, 40]. In particular the splay state regime has been mainly investigated in leaky-integrate-and-fire (LIF) neurons [28, 25, 33, 41, 42], but some studies have been also devoted to the θ-neurons [43] and to more realistic neuronal models [44]. All in all, splay states are important since they provide the simplest instance of asynchronous behaviour and can be thereby used as a testing ground for the stability of a more general class of dynamical regimes.

In this chapter, we develop a perturbative analysis for the stability of splay states in ensembles of N globally pulse-coupled identical neurons. Our model neurons are characterized by a membrane potential u that is continuously driven by the velocity field F (u), from the resetting value u = 0 towards the threshold u = 1 (see the next section for a more precise definition). As threshold and resetting value can be identified with one another and thereby u interpreted as a phase, it will be customary to refer to the case F (1) ≠ F (0) as to that of a discontinuous velocity field. Additionally, we assume that the post-synaptic potential (PSP) has a stereotyped shape, the so-called α-pulse, that is characterized by an identical rise and decay time 1/α [25]. The linear stability analysis reveals that the eigenvectors are characterized by different spatial frequencies (when moving from the neuron with the smallest to that one with the largest membrane potential). It is therefore convenient to use the frequency ϕ (scaled to the average phase separation 1/N) to parametrize the Floquet spectrum. As already discussed in [33] and mentioned in Chapter 1, there exist two components, namely short (SWs) and long (LWs) wavelengths. SWs vary on “microscopic” scales, i.e. correspond to finite ϕ values: they are typically marginally stable in the thermodynamic limit (N → ∞). LWs vary on scales of order O(1), i.e. correspond to vanishing frequencies: they have been studied in the context of a mean-field theory, i.e. by analyzing a suitable functional equation for the probability distribution of the membrane potential u [25, 45]. By developing an approach that is valid for arbitrary coupling strength and is perturbative in the inverse system-size 1/N, here we prove that the Floquet spectrum scales as 1/N^{2} and is proportional to F (1) − F (0). We are also able to determine the spectral shape and find it to be independent of the structure of the velocity field. The transition from SWs to LWs is signalled by the occurrence of a singularity in the spectrum for the frequency ϕ → 0. In the crossover region ϕ = k/N (where k is large but small compared to N) we show that the exponents remain finite and coincide with those determined in the weak coupling limit by Abbott and van Vreeswijk [25] with their mean field approach. This result is non trivial, since it is not a priori obvious that the “macroscopic” description discussed in [25] is fully contained in the “microscopic” description derived here, as they refer to two different levels of description.

More specifically, we first build a suitable event-driven map and expand it in powers of 1/N (a posteriori, we have verified that it is necessary to reach the fourth order). Afterwards, an expression of the splay state is determined: this task corresponds to finding a fixed point of the event-driven map in a suitably moving reference frame - analogously to what previously done in Section 1.3 for the LIF and in other specific contexts [33, 46, 47]. In practice this task is carried out by first taking the continuum limit for the various orders and then obtaining suitable differential equations. The solution of such equations allows proving that all finite-size corrections for both the period T and the membrane potential vanish up to the third order. Next, the stability analysis is carried out to determine the leading term of the Floquet spectrum. This task involves the introduction of a suitable Ansatz to decompose each eigenvector into the linear superposition of a slow and a rapidly oscillating component. The following continuum limit shows that the two components satisfy an ordinary and a differential equation, respectively.

Altogether, the proof of our main result requires determining all terms up to the third order in the 1/N expansion of the splay state solution, while some third order terms are not necessary for the tangent space analysis. Going beyond discontinuous fields would require extending our analysis to account for higher order terms and this might not even be sufficient to characterize analytic velocity fields. In fact, previous numerical simulations [46] suggest that the Floquet exponents scale with higher powers of 1/N that depend on which derivatives of F (u) are eventually discontinuous. Moreover, it is worth recalling that in the case of a strictly sinusoidal field, a theorem proved by Watanabe and Strogatz [48] implies that N − 3 Floquet exponents (N − 2 for a splay solution) vanish exactly for any value of N.

From the analysis of the SW spectra, one can conclude that the splay state is stable in excitatory (inhibitory) networks whenever F (0) > F (1) (F (0) < F (1)). These conditions can be extended also to the crossover region, where our results coincide with those obtained in [25] (in the limit of small coupling). Our analytical studies cannot, however, predict the behaviour of the LW component that may be responsible for the emergence of new collective solutions in excitatory networks [28, 29]. The overall scenario is partially reminiscent of the stability of synchronous and clustered regimes that is determined by the sign of the first derivative dF/du of the velocity-field averaged on the interval [0, 1] (this latter problem has been investigated in excitatory pulse-coupled integrate-and-fire oscillators subject to δ-pulses [36, 49]).

2.2 Model and Main Results

We consider a network of N identical neurons (rotators) coupled via a mean-field term. The dynamics of the i-th neuron writes as

u̇_{i}(t) = F (u_{i}) + gE(t) ≡ F_{i}(t)  i = 1, …, N ,  (2.1)

where u_{i}(t) represents the membrane potential, E(t) is the forcing field, and g is the coupling constant. When the membrane potential reaches the threshold value u_{i}(t) = u_{th} = 1, a spike is sent to all neurons (see below for the relationship between the single spikes and the global forcing field E(t)) and it is reset to u_{i}(t) = u_{R} = 0. The resetting procedure is an approximate way to describe the discharge mechanism operating in real neurons. The function F represents a velocity field for the isolated neuron and it is assumed to be everywhere positive (thus ensuring that the neurons repetitively fire, since they are supra-threshold), while F_{i} is the velocity field seen by the neuron i in the presence of a coupling with other neurons. While we consider both excitatory (g > 0) and inhibitory networks (g < 0), it is easy to show that F remains always positive to ensure the existence of splay states. For the simple choice

F (u) = a − u ,  (2.2)

the model reduces to the well known case of LIF neurons. The evolution of a membrane potential for a LIF suprathreshold neuron (a > 1) is reported in Fig. 2.1(a).

The field E is the linear superposition of the pulses emitted in the past when the membrane potential of each single neuron has reached the threshold value. By following Ref. [25], we assume that the shape of a pulse emitted at time t = 0 is given by s(t) = \frac{α^{2}t}{N}e^{−αt}, where 1/α is the pulse–width. This is equivalent to saying that the total field evolves according to the Eq. (1.2) here reported for simplicity

Ë(t) + 2αĖ(t) + α^{2}E(t) = \frac{α^{2}}{N} ∑_{n|t_{n}<t} δ(t − t_{n}) .  (2.3)

where the sum in the r.h.s. represents the source term due to the spikes emitted at times t_{n} < t.

Figure 2.1: (Color online) (a) Temporal course of the membrane potential for a suprathreshold LIF neuron in the absence of synaptic stimuli. The pulses s are exactly collocated in correspondence of the firing times. u_{R} represents the reset value and u_{th} the threshold value. (b) “Raster plot” for a splay state in a network of N = 100 neurons: neuron indices are reported as a fuction of firing times. The dots correspond to the times when the neurons emit an action potential. Indices are ordered according to the (first) time the neurons have reached the threshold.

It is convenient to transform the continuous-time model into a discrete-time mapping. We do so by integrating the equations of motion from time t_{n} to time t_{n+1} (where t_{n} is the time immediately after the n-th pulse has been emitted). The resulting map for the field variables are given by Eqs. (1.3,1.4), here rewritten for clarity

E_{n+1} = [E_{n} + τ_{n}P_{n}]e^{−ατ_{n}} ,  (2.4)
P_{n+1} = P_{n}e^{−ατ_{n}} + \frac{α^{2}}{N} ,

where τ_{n} = t_{n+1} − t_{n} is the interspike time interval and, for the sake of simplicity, we have introduced the new variable P := αE + Ė.

In this paper we focus on a specific solution of the network dynamics, namely on splay states, which are asynchronous states, where all neurons fire periodically with period T and two successive spike emissions occur at regular intervals τ_{n} ≡ T/N. A typical “raster plot” for this state is reported in Fig. 2.1(b) (see also Section 1.2.1 for more details).

In the large-N limit it is natural to consider 1/N as a smallness parameter and thereby to expand the evolution equations in powers of 1/N. In order to perform this expansion, the unique condition to require is the differentiability of the velocity field F (u) in the definition interval ]0, 1[. The only exception is represented by the boundaries of the interval, where a discontinuity is allowed.

The first result of this paper is that under the assumption that the velocity field F (u) is differentiable at least four times, the dependence of the period T onto the size N is of order o(1/N^{3}). In the specific case of LIF neurons, we show in Appendix B that the leading correction δT to the infinite size result is indeed of order O(1/N^{4}) and, more precisely,

δT = \frac{K(α) − 6}{720} \frac{[a(1 − e^{−T}) − 1]}{ge^{−T} + a (T + 1 − e^{−T}) − 1} \frac{T^{5}}{N^{4}} ,  (2.5)

where K(α) encodes the information on the pulse dynamics (see Eq. (B.12)). We did not dare to estimate the quartic contribution for generic velocity fields, not only because the algebra would be utterly complicated, but also since our main motivation is to determine the leading contributions in the stability analysis, and it turns out that it is sufficient to determine the splay state up to the third order.

The study of the stability requires determining the Floquet spectrum, i.e. the complex eigenvalues of a given periodic orbit of period T. With reference to a system of size N and by following [33], the Floquet multipliers can be written as in Eq. (1.20), where we have properly replaced τ_{n} = T/N. Notice that since the total number of exponents is N + 1 (the zero exponent has been removed by taking the Poincaré section), we are going to miss one of them. Furthermore, as shown in [33] the N-th and N +1-th exponents are associated to the field evolution and they will be not considered in our analysis.

In the following we prove that the leading term of the spectrum is,

λ_{k} = \frac{gα^{2}}{12} \frac{F (1) − F (0)}{F(1)F(0)} (\frac{6}{1 − cos ϕ_{k}} − 1) \frac{1}{N^{2}} k = 1, …, N − 1 ,  (2.6)

i.e., for discontinuous velocity fields, the real part of the spectrum scales as 1/N^{2}, while the imaginary part is of even higher order.

For continuous fields, it has been numerically observed that the scaling of the spectrum is at least O(1/N^{4}) [46]. In other words the shape of the spectrum is universal, apart from a multiplicative factor that vanishes if and only if F (1) = F (0), i.e. for true phase rotators where u = 0 coincides with u = 1.

In the limit ϕ_{k} → 0, Eq. (2.6) is singular and λ_{k} seems to diverge. However, one should recall that ϕ_{k} = 2πk/N, i.e. it is a discrete variable in a finite system. By expanding for k ≪ N, one finds that

λ_{k} = \frac{gα^{2}}{4π^{2}k^{2}} \frac{F (1) − F (0)}{F(1)F(0)} .  (2.7)

This expression, which holds in the limit of (1 ≪ k ≪ N), can be compared with the results obtained in [25], in the small-coupling limit (|g| << 1), for sufficiently large k (for more details see Appendix B.4). The equivalence of the methods has two implications: (i) the crossover component is “universal”, as it is valid also for large coupling constants; (ii) the “macroscopic” stability is fully contained in our “microscopic” analysis. In fact, numerical studies reveal a perfect correspondence also for the LW component that is not amenable to an analytic treatment [46, 50].

From Eqs. (2.6) and (2.7), it follows that the stability of the splay state can be inferred, for arbitrary coupling strength, from the sign of F (1) − F (0): in excitatory (inhibitory) networks, the state is stable whenever F (0) > F (1) (F (0) < F (1)). The same result was previously reported in the weak coupling limit in [25] ^{1}. It is, however, necessary to point out that such condition(s) do not account for instabilities that can arise in the LW component. This is, e.g., the case of the onset of partial synchronization via a supercritical Hopf bifurcation [28, 29].

2.3 Event Driven Map

By following Ref. [46, 26], it is convenient to pass from a continuous to a discrete time evolution rule, by introducing the event-driven map which connects the network configuration at subsequent spike emissions occurring at time t_{n} and t_{n+1}. The membranepotential value u_{i}^{−}(t_{n+1}) just before the emission of the (n+1)-st spike can be obtained by formally integrating Eq. (2.1)

u_{i}^{−}(t_{n+1}) − u_{n,i}(t_{n}) = ∫_{t_{n}}^{t_{n+1}} dtF (u_{i}(t)) + g ∫_{t_{n}}^{t_{n+1}} dt[E_{n} + (t − t_{n})P_{n}]e^{−α(t−t_{n})}
≡ A_{1}(u_{i}) + A_{2}(E_{n}, P_{n}) ,  (2.8)

where the minus superscript means that the map construction has not yet been completed. This task is accomplished by ordering the membrane potentials from the largest (j = 1) to the smallest (j = N) value and by passing to a comoving frame that advances with the firing neuron, i.e. by shifting the neuron index by one unit u_{n+1,j−1} = u_{j}^{−}(t_{n+1}), where the first subscript indicates that the variable is determined at time t_{n+1}. As a result, the evolution is described by an event-driven map that is composed of Eqs. (1.3,1.4), together with the following recursive relation

u_{n+1,j−1} = u_{n,j} + A_{1}(u_{j}, τ_{n}) + A_{2}(E_{n}, P_{n}, τ_{n}) ,  (2.9)

where τ_{n} = t_{n+1} − t_{n}. The time at which the (n + 1)-st spike is emitted can be determined implicitly from Eq. (2.9) by setting j = 1, since, by definition of the model, u_{n,0} ≡ 1. In this reference frame the splay state corresponds to a fixed point of the map.

Now we perform a perturbative expansion of both terms A_{1} and A_{2} (see Eq. (2.8)) in powers of δt = t − t_{n}. This is justified by the smallness of δt which is O(1/N). More precisely, the first integral appearing on the rhs of Eq. (2.8) is solved perturbatively by introducing a polynomial expansion of u_{i}(t) around t = t_{n},

u_{j}(t) = u_{n,j} + u̇_{n,j}δt + \frac{1}{2}ü_{n,j}δt^{2} + \frac{1}{6}u _{n,j}^{···}δt^{3} + O (δt^{4}) .  (2.10)

1Please notice that while Eq. (5.4) in [25] is consistent with our results and applies to both inhibitory and excitatory coupling, inequality (5.5) holds only for positive coupling.

Explicit expressions for the time derivates of u_{j} can be obtained from Eq. (2.1) and its time derivatives,

ü_{n,j} = F′ (u_{n,j}) u̇_{n,j} + gĖ_{n} ,
u _{n,j}^{···} = F′′ (u_{n,j}) u̇_{n,j}^{2} + F′ (u_{n,j}) ü_{n,j} + gË_{n} ,

where one can further eliminate Ë_{n} with the help of Eq. (1.2).

By inserting the expansion (2.10) into the expression of A_{1}, expanding the function F (u) around u_{n,j}, and performing the trivial integrations, one obtains

A_{1} = F_{n,j}τ_{n} + F′_{n,j}F_{n,j} \frac{τ_{n}^{2}}{2} + \{[F′′_{n,j}F_{n,j} + F_{n,j}^{′2}]F_{n,j} + gĖ_{n}F′_{n,j}\} \frac{τ_{n}^{3}}{6} + \{F′′′_{n,j}F_{n,j}^{3} (2.11)

+4F′_{n,j}F′′_{n,j}F_{n,j}^{2} + F_{n,j}^{′3}F_{n,j} + g [(3F′′_{n,j}F_{n,j} + F_{n,j}^{′2} − αF′_{n,j})Ė_{n} − αF′_{n,j}P_{n}]\} \frac{τ_{n}^{4}}{24} + O(τ_{n}^{5}) ,

where we have introduced the short-hand notation F_{n,j} for F (u_{n,j}) (and analogously for F).

In the case of A_{2}, it is possible to obtain an exact expression for the integral, which can be then expanded in powers of τ_{n}, as follows

A_{2} = \frac{g}{α}E_{n}(1 − e^{−ατ_{n}}) − \frac{g}{α}P_{n}τ_{n}e^{−ατ_{n}} + \frac{g}{α^{2}} P_{n}(1 − e^{−ατ_{n}})  (2.12)
= gE_{n}τ_{n} + gĖ_{n} \frac{τ_{n}^{2}}{2} − gα(Ė_{n} + P_{n})\frac{τ_{n}^{3}}{6} + gα^{2}(Ė_{n} + 2P_{n})\frac{τ_{n}^{4}}{24} + O(τ_{n}^{5}) .

By finally assembling Eqs. (2.9,2.11,2.12), we obtain a perturbative expression for the evolution rule of the membrane potential,

u_{n+1,j−1} = u_{n,j} + F_{n,j}τ_{n} + [gĖ_{n} + F′_{n,j}F_{n,j}] \frac{τ_{n}^{2}}{2} + \{F′_{n,j}[F′_{n,j}F_{n,j} + gĖ_{n}] (2.13)

+F′′_{n,j}F_{n,j}^{2} − gα[P_{n} + Ė_{n}]\} \frac{τ_{n}^{3}}{6} + \{−gα(Ė_{n} + P_{n})F′_{n,j} + 4F′_{n,j}F′′_{n,j}F_{n,j}^{2}

+F_{n,j}^{′3}F_{n,j} + g [(3F′′_{n,j}F_{n,j} + F_{n,j}^{′2})Ė_{n} + α^{2}(Ė_{n} + 2P_{n})] + F′′′_{n,j}F_{n,j}^{3} \} \frac{τ_{n}^{4}}{24} + O(τ_{n}^{5}) .

2.4 Splay State Solution

The splay state is a fixed point of the event-driven mapping with a constant interspike interval τ = T/N. Since the fixed point solutions do not depend on the index n they are denoted as,

E_{n} ≡ Ẽ , P_{n} ≡ P̃ , u_{n,j} ≡ ũ_{j} .  (2.14)

Substituting Eq. (2.14) into Eqs. (??) one obtains

P̃ = \frac{α^{2}}{N} \frac{1}{(1 − e^{−αT/N})} , Ẽ = \frac{T}{N} \frac{P̃}{(e^{αT/N} − 1)} .  (2.15)

By introducing Eq. (2.15) in Eq. (2.9) and eliminating the n dependence, we obtain a recursive equation for the variable ũ_{j}. The fixed point solution corresponds to the “trajectory” that, starting from ũ_{N} = 0, ends in ũ_{0} = 1 and can be found by tuning the “parameter” T/N. The existence of one or more solutions is related to the dependence on T. Simple calculations reveal that A_{2} = g/N, i.e. it is independent of T; moreover since A_{1} is the integral from time 0 to time T/N of a positive defined function (F > 0), it is a monotonically increasing function of T which vanishes in zero. Accordingly for any function F, the minimal value of u_{0} is g (obtained for T = 0), while the maximal value is unbounded from above. Therefore there exists one and only one solution, provided that g < 1

The variables Ẽ and P̃ depend explicitly on N, and it is natural to expect that the period T itself varies with N. By replacing the expansion of T in powers of the smallness parameter 1/N

T = ∑_{h=0}^{4} \frac{T^{(h)}}{N^{h}} + O (\frac{1}{N^{5}} ) .  (2.16)

in the Eq. (2.15), we obtain perturbative expressions for P̃, Ẽ, and ˙ Ẽ, namely Eqs. (B.1, B.2, B.3) reported in Appendix B.1. As for the membrane potential, it is necessary to introduce the formal expansion

ũ_{j} = ∑_{h=0}^{4} \frac{ũ_{j}^{(h)}}{N^{h}} + O (\frac{1}{N^{5}} ) .  (2.17)

Substituting the expansions of T, ũ_{j}, Ẽ, and P̃ into Eq. (2.13) (after having dropped the n dependence and expanded the “F ” functions) one obtains an equation for ũ_{j}^{(h)} and T^{(h)}, namely

∑_{h=0}^{4} \frac{ũ_{j−1}^{(h)} − ũ_{j}^{(h)}}{N^{h}} = ∑_{h=1}^{4} \frac{Q^{(h)}}{N^{h}} + O (\frac{1}{N^{5}} ) ,  (2.18)

where the Q variables are defined in Appendix B.1. Notice the dependence on T^{(h)} variables is hidden in Q terms.

In the large-N limit, one can introduce the continuous spatial coordinate x = j/N. In practice, this is tantamount to write,

U^{(h)}(x = j/N) = lim_{N→∞} ũ_{j}^{(h)} ,  h = 0, ⋯, 4 .  (2.19)

It is important to stress that the event-driven neuronal evolution in the comoving frame implies that U(0) = 1, i.e. the first neuron will fire at the next step, and U(1) = 0, i.e. the membrane potential of the last neuron has been just reset to zero. This implies that U^{(0)}(0) = 1 and U^{(0)}(1) = 0, while U^{(h)}(0) = U^{(h)}(1) = 0 for any h > 0.

Furthermore, by expanding U^{(h)}(x) around x = j/N, one obtains

U^{(h)}(x − 1/N) = U^{(h)}(x) + ∑_{m=1}^{4} \frac{1}{m!} (\frac{−1}{N})^{m} \frac{d^{m}}{dx^{m}} U^{(h)}(x) + O (\frac{1}{N^{5}} ) ,  (2.20)

By inserting this expansion into Eq. (2.18), we obtain an equation that can be effectively split into terms of different order that will be analyzed separately. Notice that by retaining terms of order h, it is possible to determine the original variables at order h − 1.

2.4.1 Zeroth Order Approximation

By assembling the first order terms, we obtain the evolution equation for the zeroth order membrane potential, namely

\frac{dU^{(0)}}{dx} = −g − T^{(0)}F (U^{(0)}) .  (2.21)

This equation is equal to the evolution equation of the membrane potential for a constant field E, with x playing the role of a (inverse) time. Please notice that, up to first order, Ẽ = 1/T^{(0)} (see Eq. (B.2)). An implicit and formal solution of Eq. (2.21) is,

1 − x = ∫_{0}^{U^{(0)}} \frac{dv}{g + T^{(0)}F (v)} ,  (2.22)

where we have imposed the condition U^{(0)}(1) = 0. However, there is a second condition to impose, namely U^{(0)}(0) = 1. This second condition transforms itself in the equation defining the interspike time interval T^{(0)}, when N → ∞ (i.e. in the thermodynamic limit)

1 = ∫_{0}^{1} \frac{dU^{(0)}}{g + T^{(0)}F (U^{(0)})} .  (2.23)

This result is, so far, quite standard and could have been easily obtained by just assuming a constant field E in equation (2.1). If we introduce the formal relation F′[U^{(0)}(x)] = \frac{dF (U^{(0)})}{dU^{(0)}} in Eq. (2.21) we obtain

\frac{dF (U^{(0)})}{g + T^{(0)}F (U^{(0)})} = −F′[U^{(0)}(x)]dx ,

which can be easily integrated

∫_{F (U^{(0)}(0))}^{F (U^{(0)}(1))} \frac{dF (U^{(0)})}{g + T^{(0)}F (U^{(0)})} = − ∫_{0}^{1} F′[U^{(0)}(x)]dx ,

giving the following relation (already derived in [29], by following a different approach)

\frac{e^{−T^{(0)}H(0)}}{F(U(0))} = \frac{e^{−T^{(0)}H(1)}}{F(U(1))} ,  (2.24)

where, for later convenience, we have introduced

H(x) = ∫_{0}^{x} F′[U^{(0)}(y)]dy ,  (2.25)

and where, for the sake of simplicity, the prime denotes derivative with respect to the variable U^{(0)} and the dependence of F and F′ on U^{(0)} has been dropped.

2.4.2 First Order Approximation

By collecting the terms of order 1/N^{2}, one obtains

\frac{dU^{(1)}}{dx} = −T^{(0)}F′U^{(1)} + \frac{1}{2} \frac{d^{2}U^{(0)}}{dx^{2}} − F T^{(1)} − \frac{1}{2}(T^{(0)})^{2}F′F − \frac{g}{2}T^{(0)}F′ .  (2.26)

An explicit expression for the second derivative of U^{(0)}(x) appearing in Eq. (2.26) can be computed by deriving Eq. (2.21) with respect to x. This allows rewriting Eq. (2.26) in a simplified form, namely

\frac{dU^{(1)}(x)}{dx} = −U^{(1)}T^{(0)}F′ − T^{(1)}F .  (2.27)

By imposing U^{(1)}(1) = 0, one obtains the general solution of Eq. (2.27),

U^{(1)}(x) = ∫_{x}^{1} du T^{(1)}F [U^{(0)}(u)] exp [T^{(0)}(H(x) − H(u))] ,

where H(x) is defined by Eq. (2.25). The further condition to be satisfied, U^{(1)}(0) = 0, implies T^{(1)} = 0 and thereby we have U^{(1)}(x) ≡ 0, i.e. first-order corrections vanish both for the period and the membrane potential.

2.4.3 Second Order Approximation

The second order corrections can be estimated by assembling terms of order 1/N^{3} and by imposing the previously determined conditions T^{(1)} = 0 and U^{(1)}(x) = 0,

\frac{dU^{(2)}}{dx} = −T^{(0)}F′U^{(2)} − \frac{1}{6} \frac{d^{3}U^{(0)}}{dx^{3}} − F T^{(2)} − \frac{g^{2}}{6} T^{(0)}F′′

−\frac{g}{6}(T^{(0)})^{2} [2F F′′ + F^{′2}] − \frac{(T^{(0)})^{3}}{6} [F′′F^{2} + F^{′2}F ] .

Once evaluated d^{3}U^{(0)}/dx^{3} from Eq. (2.21), the above ODE reduces to

\frac{dU^{(2)}}{dx} = −U^{(2)}T^{(0)}F′ − T^{(2)}F ;

which has the same structure as Eq. (2.27). Since one has also to impose the same boundary conditions as for the first order, namely U^{(2)}(0) = U^{(2)}(1) = 0, we can conclude that T^{(2)} = 0 and, consequently, U^{(2)}(x) ≡ 0. Therefore, second order corrections are absent too.

2.4.4 Third Order Approximation

By assembling terms of order 1/N^{4}, once imposed that first and second order corrections vanish, one obtains

\frac{dU^{(3)}}{dx} = −T^{(0)}F′U^{(3)} + \frac{1}{24} \frac{d^{4}U^{(0)}}{dx^{4}} − F T^{(3)} − \frac{g^{3}}{24}T^{(0)}F′′′ − \frac{g^{2}}{6} (T^{(0)})^{2}F′F′′

−\frac{g^{2}}{8} (T^{(0)})^{2}F F′′′ − \frac{g}{24}(T^{(0)})^{3} [F^{′3} + 8F F′F′′ + 3gF^{2}F′′′]

−\frac{(T^{(0)})^{4}}{24} F F′ [F^{′2} + 4F F′′ + F^{3}F′′′] .  (2.28)

By replacing d^{4}U^{(0)}/dx^{4} with its expression derived from Eq. (2.21), equation (2.28) takes the same form as in the two previous examined cases, namely

\frac{dU^{(3)}}{dx} = −U^{(3)}T^{(0)}F′ − T^{(3)}F .

Therefore, we can safely conclude that third order terms vanish too.

The LIF model can be solved exactly for any value of N, starting from the asymptotic value (N → ∞). As shown in Appendix B, it turns out that the leading corrections are of fourth order for both the period T and the membrane potential.

2.5 Linear Stability Analysis

The fixed-point analysis has revealed that the finite-size corrections to the stationary solutions are of order o(1/N^{3}). Since such deviations do not affect the leading terms of the linear stability analysis (as it can be verified a posteriori) they will be simply neglected. Therefore, for the sake of simplicity, from now on, T^{(0)} and ũ_{j}^{(0)} will be simply referred to as T and ũ_{j}.

The evolution rule in the tangent space is obtained by differentiating Eqs. (1.3,1.4) and then by expanding in powers of τ (this is equivalent to expanding in powers of 1/N, as the dependence of τ on N would only generate higher order terms),

δP_{n+1} = Z^{(1)}δP_{n} + P̃Z^{(2)}δτ_{n} ,  (2.29)
δE_{n+1} = Z^{(1)}δE_{n} + Z^{(3)}δP_{n} + [P̃Z^{(4)} − αẼZ^{(1)}] δτ_{n} ,  (2.30)

where the Z variables are polynomials of τ defined in Appendix B.3.1.

By further differentiating Eq. (2.13) around the fixed point solution, one obtains

δu_{n+1,j−1} = W^{(1)}δu_{n,j} + W^{(2)}δE_{n} + W^{(3)}δĖ_{n} + W^{(4)}δτ_{n} .  (2.31)

Finally, δτ_{n} can be determined by differentiating Eq. (2.13) for j = 1

δτ_{n} = −\frac{δE_{n}}{\overline{F}_{1}} W^{(5)} − \frac{δP_{n}}{\overline{F}_{1}} W^{(6)} − \frac{δu_{n,1}}{\overline{F}_{1}} W^{(7)} ,  (2.32)

where the auxiliary W variables are defined in Appendix B.3.1.

As usual, the eigenvalue problem can be solved by introducing the Ansatz,

δu_{n,j} = μ_{k}^{n}δu_{j}  δP_{n} = μ_{k}^{n}δP  δE_{n} = μ_{k}^{n}δE  δτ_{n} = μ_{k}^{n}δτ ,  (2.33)

where μ_{k} labels the eigenvalues, which must also be expanded as,

μ_{k} = e^{iϕ_{k}}e^{(λ_{k}+iω_{k})T/N} = e^{iϕ_{k}} (1 + ∑_{h=1}^{3} \frac{Γ^{(h)}}{N^{h}} + O (\frac{1}{N^{4}} )) .  (2.34)

where Γ^{(h)} is, in principle, a complex number and, for the sake of simplicity, we have dropped its dependence on k. Finally, as already shown at zeroth order, the eigenvalues correspond to a pure rotation (specified by ϕ_{k}) with no expansion or contraction, i.e. Γ^{(0)} = 0.

By inserting the above Ansätze in the map expression (2.29,2.30,2.31,2.32), one obtains, after eliminating δP, δE and δτ, a closed equation for δu_{j},

e^{iϕ_{k}} (1 + \frac{Γ^{(1)}}{N} + \frac{Γ^{(2)}}{N^{2}} + \frac{Γ^{(3)}}{N^{3}} ) δu_{j−1} = \{1 + \overline{F} ′_{j} \frac{T}{N} + [\overline{F} ′′_{j}\overline{F}_{j} + \overline{F} _{j}^{′2}] \frac{T^{2}}{2N^{2}}  (2.35)

+ [\overline{F} ′′′_{j}\overline{F}_{j}^{2} + 4\overline{F} ′_{j}\overline{F} ′′_{j}\overline{F}_{j} + \overline{F} _{j}^{′3}] \frac{T^{3}}{6N^{3}} \} δu_{j} − \{\overline{F}_{j} + \overline{F} ′_{j}\overline{F}_{j} \frac{T}{N} + [\frac{\overline{F} ′′_{j}\overline{F}_{j}^{2}}{2} + \frac{\overline{F} _{j}^{′2}}{2} \overline{F}_{j}

+\frac{g}{T} α^{2} \frac{e^{2iϕ_{k}} + 10e^{iϕ_{k}} + 1}{12(e^{iϕ_{k}} − 1)^{2}} (\begin{matrix}\overline{F} & j \\ \overline{F} & 1 \end{matrix} − 1)] \frac{T^{2}}{N^{2}} + [\frac{\overline{F} _{j}^{′3}}{6} \overline{F}_{j} + \frac{2}{3}\overline{F} ′_{j}\overline{F} ′′_{j}\overline{F}_{j}^{2} + \frac{\overline{F} ′′′_{j}}{6} \overline{F}_{j}^{3}

+ \frac{gα^{2}}{3T^{2}} Γ^{(1)} \frac{2e^{2iϕ_{k}} − 3e^{iϕ_{k}}}{(e^{iϕ_{k}} − 1)^{3}} (\begin{matrix}\overline{F} & j \\ \overline{F} & 1 \end{matrix} − 1) + \frac{g}{T} α^{2} \frac{e^{2iϕ_{k}} + 10e^{iϕ_{k}} + 1}{12(e^{iϕ_{k}} − 1)^{2}} \begin{matrix}\overline{F} & j \\ \overline{F} & 1 \end{matrix} (\overline{F} ′_{j} − \overline{F} ′_{1})

+\frac{gα^{2}}{T} \frac{5e^{iϕ_{k}} + 1}{12(e^{iϕ_{k}} − 1)^{2}} (\overline{F} ′_{1} \begin{matrix}\overline{F} & j \\ \overline{F} & 1 \end{matrix} − \overline{F} ′_{j}) + \frac{gα^{3}}{T} \frac{e^{iϕ_{k}}(e^{iϕ_{k}} + 1)}{(e^{iϕ_{k}} − 1)^{3}} (1 − \begin{matrix}\overline{F} & j \\ \overline{F} & 1 \end{matrix} )] \frac{T^{3}}{N^{3}} \} \frac{δu_{1}}{\overline{F}_{1}} ,

that is the object of our investigation. The overline means that the function is evaluated in ũ_{j}^{(0)}, corresponding to the infinite N limit.

2.5.1 Continuum Limit

Similarly to the splay-state estimation, it is convenient to take the continuum limit. However, at variance with the previous case, now one should take in to account also the presence of fast scales associated to the “spatial” dependence of ϕ_{k}.

Therefore, the correct Ansatz is slightly more complicated and we have to separate slowly and rapidly oscillating terms,

δu_{j} = π_{j} + ϑ_{j}e^{iϕ_{k}j},  (2.36)

where the complex exponential term accounts for the fast oscillations of the eigenvectors, while

π_{j} = ∑_{h=0}^{3} \frac{π_{j}^{(h)}}{N^{h}} + O (\frac{1}{N^{4}} ) , ϑ_{j} = ∑_{h=0}^{3} \frac{ϑ_{j}^{(h)}}{N^{h}} + O (\frac{1}{N^{4}} ) .

are slowly varying variables. On the one hand, the existence of the slow component π_{j} follows from the analogy with the real-space evolution. On the other hand, the presence of the rapidly oscillating terms e^{iϕ_{k}j}, first noticed in Ref. [25] in the uncoupled limit, suggests the presence of a second slow field, namely ϑ_{j}. Anyway, the correctness and uniqueness of Ansatz (2.36) is ensured a posteriori, by the consistency of the equations that are obtained for the various perturbation orders.

Nest, we can finally introduce the continuous variable x = j/N, as previously done in real space (see Eq. (2.19)),

Π_{j}^{(h)}(x = \frac{j}{N}) = π_{j}^{(h)} ,  Θ_{j}^{(h)}(x = \frac{j}{N}) = ϑ_{j}^{(h)} ,  (2.37)

where h = 0, ⋯, 3. This allows expanding δu_{j−1} = π_{j−1} + ϑ_{j−1}e^{iϕ_{k}(j−1)} around x = j/N, similarly to what done in Eq. (2.20). At variance with the computation of the fixed point, now there are also terms like U(1/N) and δU(1/N), whose computation requires a similar expansion but around x = 0. By incorporating all the expansion terms within Eq. (2.35), we have finally an equation, where terms of different orders are naturally separated from one another. The calculations are summarized in Appendix C and the final equation is (B.18). By separately treating the different orders, we obtain differential and ordinary equations for the Θ and Π variables. It turns out that it is necessary to consider in parallel different orders in the fast and slow terms to obtain Θ and Π to the same order. As a consequence, we will see that it is sufficient to expand δU(1/N) up to order O(1/N^{3}).

2.5.2 Zeroth Order Approximation

By assembling terms of order O(1/N) in Eq. (B.18), multiplied by the fast oscillating factor e^{iϕ_{k}j}, we obtain a first-order linear differential equation for Θ^{(0)}, namely

\frac{dΘ^{(0)}}{dx} = −Θ^{(0)}(T F′(U(x)) − Γ^{(1)}) ,  (2.38)

where Γ^{(1)} is the first order correction to the Floquet exponent which should be determined. It is important to remind that the prime denotes derivative with respect to the variable U^{(0)}, which has been simply redefined U, as previously mentioned. The solution is

Θ^{(0)}(x) = K^{(0)} exp [Γ^{(1)}x − T H(x)] ,  (2.39)

where we made use of the definition (2.25) and K^{(0)} is a suitable integration constant.

By assembling now the slow terms of zeroth order and reminding the definition of F(U(x)), we find the following algebraic equation

Π^{(0)}(x)(e^{iϕ_{k}} − 1) = − [e^{iϕ_{k}}Θ^{(0)}(0) + Π^{(0)}(0)] \frac{F(U(x))}{F(U(0))} .

With the help of Eq. (2.39), we obtain

Π^{(0)}(0) = −Θ^{(0)}(0) = −K^{(0)}e^{−T H(0)} ,
Π^{(0)}(x) = −K^{(0)}e^{−T H(0)} \frac{F(U(x))}{F(U(0))} .

We can now impose the boundary condition δU^{(0)}(x = 1) = Θ^{(0)}(1)+Π^{(0)}(1) = 0. This implies that

\frac{e^{−T H(1)+Γ^{(1)}}}{F(U(1))} = \frac{e^{−T H(0)}}{F(U(0))} .

From Eq. (2.24), we can conclude that Γ^{(1)} = 2πim, for 0 ≤ m < N (values outsided this range give identical solutions). Since the vector δu_{j} was assumed to be characterized by the phase jϕ_{j} (see Eq. (2.36)), the phase factor Γ^{(1)} = 2πim, would imply (through Eq. (2.39)) that ϕ_{k} has to be shifted by some amount, contrary to the initial assumption. Accordingly, the only solution consistent with the original Ansatz is m = 0, i.e. Γ^{(1)} = 0 and, from Eq. (2.39),

Θ^{(0)}(x) = K^{(0)}e^{−T H(x)},

i.e. the eigenvectors are independent of the phase ϕ_{k} and equal to one another. In other words the degeneracy has not been lifted.

2.5.3 First Order Approximation

By assembling the fast terms of order 1/N^{2} and by setting Γ^{(1)} = 0, we find that Θ^{(1)} satisfies the following first order differential equation,

\frac{dΘ^{(1)}}{dx} = Γ^{(2)}Θ^{(0)} − Θ^{(1)}T F′(U(x)) ,

whose solution is

Θ^{(1)}(x) = (Γ^{(2)}K^{(0)}x + K^{(1)}) e^{−T H(x)} ,  (2.40)

where K^{(1)} is an integration constant associated with the solution of the previous equation.

By collecting the slow terms of order 1/N in Eq. (B.18), one obtains the algebric equation

Π^{(1)}(x)(e^{iϕ_{k}} − 1) = − [e^{iϕ_{k}}Θ^{(1)}(0) + Π^{(1)}(0)] \frac{F(U(x))}{F(U(0))} ,

whose solution is,

Π^{(1)}(0) = −Θ^{(1)}(0) = −K^{(1)}e^{−T H(0)} ,
Π^{(1)}(x) = −K^{(1)}e^{−T H(0)} \frac{F(U(x))}{F(U(0))} .

By imposing the boundary condition δU^{(1)}(x = 1) = Θ^{(1)}(1)+Π^{(1)}(1) = 0, it is possible to evaluate Γ^{(2)},

Θ^{(1)}(1) + Π^{(1)}(1) = (K^{(0)}Γ^{(2)} + K^{(1)})e^{−T H(1)} − K^{(1)}e^{−T H(0)} \frac{F(U(1))}{F(U(0))} = 0 .

Again from Eq. (2.24) and using the same argument as in the previous section, we find that Γ^{(2)} = 0 and thereby (from Eq. (2.40))

Θ^{(1)}(x) = K^{(1)}e^{−T H(x)} .

Altogether, we can conclude that the second order correction to the Floquet exponent vanishes as well, and one cannot lift the degeneracy among the eigenvectors.

2.5.4 Second Order Approximation

By assembling fast terms of order 1/N^{3} appearing in Eq. (B.18) and by setting Γ^{(1)} = Γ^{(2)} = 0, the following first order differential equation for Θ^{(2)} can be derived

\frac{dΘ^{(2)}}{dx} = Γ^{(3)}Θ^{(0)} + Θ^{(2)}T F′(U(x)) ,

whose solution is

Θ^{(2)}(x) = (Γ^{(3)}K^{(0)}x + K^{(2)}) e^{−T H(x)} ,

where K^{(2)} is an integration constant associated with the solution of the previous differential equation.

Furthermore, by collecting the slow terms of order 1/N^{2}, we obtain the algebric equation,

Π^{(2)}(x)(e^{iϕ_{k}} − 1) = gα^{2}T Θ^{(0)}(0) \frac{e^{2iϕ_{k}} + 10e^{iϕ_{k}} + 1}{12(e^{iϕ_{k}} − 1)} \frac{F (U(0)) − F (U(x))}{[F(U(0))]^{2}}

− [e^{iϕ_{k}}Θ^{(2)}(0) + Π^{(2)}(0)] \frac{F(U(x))}{F(U(0))} .

By imposing that the above equation is satisfied for x = 0, it reduces to

Π^{(2)}(0) = −Θ^{(2)}(0) = −K^{(2)}e^{−T H(0)} ,
Π^{(2)}(x) = gα^{2}T Θ^{(0)}(0) \frac{e^{2iϕ_{k}} + 10e^{iϕ_{k}} + 1}{12(e^{iϕ_{k}} − 1)^{2}} \frac{F (U(0)) − F (U(x))}{[F(U(0))]^{2}} − Θ^{(2)}(0)\frac{F(U(x))}{F(U(0))} .

Finally, by imposing the boundary condition δU^{(2)}(x = 1) = Θ^{(2)}(1) + Π^{(2)}(1) = 0, it is possible to determine Γ^{(3)},

Γ^{(3)} = \frac{gα^{2}}{12} T \frac{F (U(0)) − F (U(1))}{F(U(0))F(U(1))} (\frac{6}{1 − cos ϕ_{k}} − 1) .  (2.41)

Accordingly, Γ^{(3)} is real and depends on the difference between F (U(x = 1)) ≡ F (0) and F (U(x = 0)) ≡ F (1), confirming the numerical findings in [46]. Therefore, the imaginary terms ω_{k} are smaller than 1/N^{2}.

In the specific example of a leaky integrate-and-fire neuron the expression for Γ^{(3)} reduces to

Γ^{(3)} = \frac{gα^{2}}{12} T (2 − e^{T} − e^{−T} ) (\frac{6}{1 − cos ϕ_{k}} − 1) ,  (2.42)

since, by using the equations that characterize LIF neurons, the following relation holds

\frac{F (U(1)) − F (U(0))}{F(U(1))F(U(0))} = \frac{1}{(a + \frac{g}{T})^{2}e^{−T}} = \frac{1}{(\frac{1}{1−e^{−T}})^{2}e^{−T}} = (e^{T} + e^{−T} − 2) .

All in all, Eq. (2.41) generalizes the expression found for the LIF model Eq. (2.42) [46]^{2}.

2.6 Discussion

We have derived analytically the short-wavelength component of the Floquet spectrum of the splay solution in a fully coupled network composed of generic suprathreshold pulse-coupled phase-like neurons in the large-N limit. Our analysis has revealed that, for discontinuous velocity fields, the spectrum scales as 1/N^{2} and the stability is controlled by the sign of the difference between the velocity at reset and at threshold. The shape of the spectrum is otherwise independent of the velocity field. It would be interesting to investigate the role of the pulse shape as well. Numerical studies suggest that the scaling behaviour is truly universal, while the shape of the Floquet spectrum depends on the pulse shape [50]. It would be interesting to discover whether and which pulse shapes may give rise to SW instabilities.

Networks of LIF neurons coupled via δ-like pulses are characterized by a finite (in)stability of the whole spectrum [33]. The difference with the case of α-pulse is so strong that it cannot be reconciled even by taking the limit α → ∞ (zero pulsewidth), indicating that the limits N → ∞ and zero pulse-width do not commute [33]. This reveals that even the development of a general stability theory of the simple splay states requires some further progress.

Finally, notice that although our analytical approach is able to cover the entire SW component and the crossover region, does not cover the truly long wavelengths which require going beyond a perturbative approach.

The results reported in this Chapter have been published in Ref. [34].

2In comparing with [46] one should pay attention to the different normalization used here to define μ_{k} in Eq. (2.34)

Collective oscillations in massively connected neural networks

3.1 Introduction

One of the most general and relevant dynamical phenomenon observed in the mammalian brain is the rhythmic coherent behavior involving different neuronal populations [9]. Collective oscillations are commonly associated with the inhibitory role of interneurons [51]; however, coherent activity patterns have been observed also in “in vivo” measurements of the developing rodent neocortex and hyppocampus for a short period after birth [52], despite the fact that at this early stage the nature of the involved synapses is essentially excitatory [53].

Independently, theoretical studies of fully coupled excitatory networks of leaky integrate-and-fire (LIF) neurons have revealed the onset of macroscopic collective periodic oscillations [28, 29] (CPOs). In this model the observed collective oscillations are a manifestation of the partial synchronization among the neuron dynamics, which arises from the destabilization of the splay state regime, as described in the Section 1.2.2.

Since real neural circuits are not expected to be fully connected [54], it is important to investigate the role of dilution for the occurrence and the stability of CPO. Therefore we investigate an excitatory network of LIF neurons with 20% of missing links in two different setups: quenched disorder, where the network topology is fixed; annealed disorder, where the active connections are randomly and independently chosen at each pulse emission. The physiological motivation of the latter topology relies on the fact that the synaptic transmission of signals is a stochastic, or unreliable [55], process. As a first step, we rewrite the dynamical equations as a suitable event driven map, by extending the approach developed in [33]. We do so by introducing a pair of variables for each neuron, to account for the evolution of the local electric field. This step is particularly important for the computation of the Lyapunov exponents, as it allows expressing the evolution equations into a “canonical” form and thereby allows the implementation of standard dynamical-system tools.

We find that the regime of CPOs is robust against the presence of dilution, both in the quenched and annealed setup. However, at variance with the homogeneous fully-coupled case, the dynamics of finite disordered networks turns out to be chaotic, although the degree of chaoticity decreases with the number N of neurons. In fact, the maximum Lyapunov exponent goes to zero as 1/N^{β}. The exponent β is smaller in the quenched setup, indicating that finite-size effects are stronger. In the thermodynamic limit, the random network dynamics reduces to that of a fully homogeneous system with a suitably rescaled coupling strength. In the homogeneous case, we are able to determine the full Lyapunov spectrum for sufficiently large numbers of neurons. As a result, we find that the first band of the spectrum scales as that of the splay state (as shown in the previous Chapter) [25, 46], namely as 1/N^{2}.

3.2 The Model

We study a diluted network of N neurons, whose individual dynamics is modeled as a LIF oscillator. Following the Eq. (1.1) and introducing the dilution, the membrane potential u_{i}(t) of the i − th oscillator evolves according to the differential equation

u̇_{i}(t) = a − u_{i}(t) + g E_{i}(t)  i = 1, ⋯, N  (3.1)

where a > 1 is the suprathreshold input current and g > 0 gauges the coupling strength of the excitatory interaction with the neural field E_{i}(t). At variance with the fullycoupled network, where all neurons depend on the same “mean field” E(t) (see Eq. 1.1), here we consider a general setup, where neurons have different connectivities inside the network. As a result, it is necessary and sufficient to introduce an explicit dependence of the neural field on the index i. The discharge mechanism operating in real neurons is modeled by assuming that when the membrane potential reaches the threshold value u_{i} = 1, it is reset to the value u_{i} = 0, while a pulse s(t) is transmitted to and instantaneously received by the connected neurons. The field E_{i}(t) is represented as the linear superposition of the pulses s(t) received by neuron i at all times t_{n} < t: the integer index n orders the sequence of the pulses emitted in the network. Each pulse s(t) is weighted according to the strength of the connection C_{j,i} between the emitting (j(n)) and the receiving (i) neuron. The matrix entries are chosen randomly with a constant probability: C_{j,i} = 1 with a probability p and C_{j,i} = 0 with a probability (1−p). This amounts to consider random networks with the same average connectivity as in the undirected Erdös-Renyi networks. In general, the connectivity matrix C is assumed to be non-symmetric. Moreover, since the examined networks are not fully coupled, the sequence of pulses contributing to E_{i} depends on the pool of neurons that are actually connected to neuron i. This is specified by the connectivity matrix C, as specified in the following formula

E_{i}(t) = \frac{1}{N} ∑_{n|t_{n}<t} C_{j(n),i}θ(t − t_{n})s(t − t_{n})  (3.2)

where θ(x) is the Heavyside function and the shape of the single pulse emitted at time t = 0 is given by s(t) = α^{2}t exp(−αt), where 1/α is the pulsewidth [28, 25]. Since also in the diluted case we will consider massively connected networks [35], where the average number of synaptic inputs per neuron varies proportionally to the system size, it is natural to scale the synaptic strength with N as done in Eq. (3.2).

The model is fully characterized by the two sets of equations (3.1) and (3.2). The dynamical system has a peculiar mathematical structure, with the field variable appearing as a memory kernel, which involves a summation over all past spiking events. We find much more convenient to turn the explicit equation (3.2) into the implicit differential equation

Ë_{i}(t) + 2αĖ_{i}(t) + α^{2}E_{i}(t) = \frac{α^{2}}{N} ∑_{n|t_{n}<t} C_{j(n),i}δ(t − t_{n}) .  (3.3)

As a result, the dynamics of the neural network model takes the more “canonical” form of a set of coupled ordinary differential equations (3.1,3.3), which can be analyzed with the standard methods of dynamical systems (see e.g., [28, 25, 33, 46]) .

3.2.1 Event-driven Map

The presence of δ-like pulses into the set of coupled differential equations (3.1) and (3.3) may still appear as an intrinsic technical difficulty for the estimation of the stability properties. Actually, the standard algorithms for the evaluation of Lyapunov exponents rely upon the integration of differentiable operators acting in tangent space (see [56]). However, one can easily get rid of this problem by transforming the differential equations into a discrete time event-driven mapping. This task can be accomplished by integrating Eq. (3.3) from time t_{n} to time t_{n+1} (t_{n} representing the time immediately after the emission of the n-ht pulse). An explicit mapping can be written by introducing the auxiliary variable P_{i} := αE_{i} + Ė_{i},

E_{i}(n + 1) = E_{i}(n)e^{−ατ(n)} + P_{i}(n)τ(n)e^{−ατ(n)}  (3.4)
P_{i}(n + 1) = P_{i}(n)e^{−ατ(n)} + C_{j(n+1),i} \frac{α^{2}}{N}  (3.5)
u_{i}(n + 1) = u_{i}(n)e^{−τ(n)} + a(1 − e^{−τ(n)}) + gH_{i}(n) .  (3.6)

Here τ(n) = t_{n+1} − t_{n} is the interspike time interval: it is determined by the largest membrane potential (identified by the label q(n)) reaching the threshold value u_{q} = 1 at time t_{n+1},

τ(n) = ln [\frac{a − u_{q}(n)}{a + gH_{q}(n) − 1}] ,  (3.7)

where

H_{i}(n) = \frac{e^{−τ(n)} − e^{−ατ(n)}}{α − 1} (E_{i}(n) + \frac{P_{i}(n)}{α − 1}) − \frac{τ(n)e^{−ατ(n)}}{(α − 1)} P_{i}(n) ;  (3.8)

for the parameter values considered in this paper (g > 0 and a > 1), H_{i}(n) > 0.

Altogether, the model now reads as a discrete time map of 3N − 1 variables, {E_{i}, P_{i}, u_{i}}: one degree of freedom, u_{q(n)} = 0, is eliminated as a result of having constructed the discrete-time dynamics with reference to a suitable Poincaré-section. At variance with the usual approach (see e.g. Ref. [57]), the evolution equation does neither involve δ-like temporal discontinuities, nor formally infinite sequences of past events: the map is a piecewise smooth dynamical system. Furthermore, the evolution of the model is performed by assuming that only one neuron at a time will reach the threshold. We have checked this assumption during the map evolution and verified that it is always fulfilled in our model with the chosen pulse shape and parameter values. The event-driven scheme here described should be adapted to take in account the case of more than one neuron firing together. In practice, the map evolution is performed in three steps: firstly the next firing neuron (denoted as q above) should be identified, this is done by evaluating the expected firing times for all the membrane potentials by employing Eq. (3.7) and by selecting the neuron q associated with the minimal τ value; secondly the map is evolved by employing Eqs. (3.4-3.6); finally, u_{q} is resetted to zero (for more details see Appendix A.1).

For the sake of simplicity, we assume that the entries of the connectivity matrix C_{j,i} are either 0 or 1 (the homogeneous fully-coupled case corresponding to C_{j,i} = 1 for all j’s and i’s). If the entries are chosen randomly, symmetries are lost and the only way to determine the asymptotic state for a finite N is by numerically simulating map (3.4-3.6). In what follows, we consider two different setups: the synaptic connections are randomly chosen and are constant in time (quenched disorder); each time a neuron fires, the neurons receiving the excitatory pulse are randomly chosen (annealed disorder).

3.2.2 Linear Stability Analysis

As usual, the stability of (3.4-3.6) can be analyzed by following the evolution of infinitesimal perturbations in the tangent space. The corresponding equations are obtained by linearizing (3.4-3.6) as follows,

δE_{i}(n + 1) = e^{−ατ(n)}δE_{i}(n) + τ(n)e^{−ατ(n)}δP_{i}(n)  (3.9)
− [αE_{i}(n)(ατ(n) − 1)P_{i}(n)] e^{−ατ(n)}δτ(n) ,
δP_{i}(n + 1) = e^{−ατ(n)} [δP_{i}(n) − αP_{i}(n)] δτ(n) ,  (3.10)
δu_{i}(n + 1) = e^{−τ(n)}[δu_{i}(n) − u_{i}(n)δτ(n)] + ae^{−τ}δτ(n) + gδH_{i}(n) (3.11)
i = 1, …, N − 1; δu_{q}(n + 1) ≡ 0 ,

where the condition δu_{q}(n + 1) ≡ 0 is due to the Poincaré map we are implementing to evolve the network dynamics.

An explicit expression of δτ(n) can be obtained by differentiating Eqs. (3.7,3.8)

δτ(n) = τ_{u}δu_{1}(n) + τ_{E}δE(n) + τ_{P} δP (n) ,  (3.12)

where

τ_{u} := \frac{∂τ}{∂u_{m}} , τ_{E} := \frac{∂τ}{∂E_{m}} , τ_{P} := \frac{∂τ}{∂P_{m}} .  (3.13)

Moreover, δH_{i}(n) is a short-cut notation for the linearization of expression (3.8), which in turn depends on δE_{i}(n), δP_{i}(n) and δτ(n). For more mathematical details see [33, 46, 26].

The degree of chaoticity of a given dynamical state is obtained by computing the Lyapunov spectrum, i.e. the set of 3N − 1 exponential growth rates λ_{i} along the independent directions in tangent space. The method here outlined to estimate the Lyapunov exponents originates from a methodology previously introduced for continuous time dynamics interrupted by discontinuous events occurring at discrete times [58, 59]. In particular, our approach differentiates from previous applications of the latter method to non-smooth flux associated to spiking neuronal models [57, 60], because we consider the evolution in the tangent space associated to a suitable Poincaré map connecting successive events.

3.3 Collective Dynamics

As discussed in Section 1.1 in the fully coupled homogeneous case the local fields E_{i}, P_{i} are independent of the index i and the number of equations reduces to N+1. Depending whether α is smaller or larger than a critical value α_{c}(g), the dynamics either converges to a so-called splay state, with constant E(t), or to a partially synchronized state, where E(t) and u_{i}(t) evolve periodically and quasi-periodically, respectively [28, 29].

This “mean-field” dynamics is expected to change in disordered networks. Given the neuron dependence of the fields E_{i} and P_{i}, we find it convenient to introduce the average variables,

Ē(n) = \frac{1}{N} ∑_{k=1}^{N} E_{k}(n) ;  P̄(n) = \frac{1}{N} ∑_{k=1}^{N} P_{k}(n) .  (3.14)

as a tool to characterize the resulting dynamical regimes and to compare with the homogeneous case. All the results reported in this paper have been obtained by fixing the fraction f of missing connections equal to f = 0.2; this amounts to p = 0.8. Numerical investigation for different choices of f (not reported here) show similar results. In analogy with the analysis of the homogeneous case, we have chosen α as the main control parameter (one can easily realize that choosing g would be equivalent).

In Fig. 3.1 we plot the maximum and the minimum values of Ē(n) for different values of α, for g = 0.4 and N = 1600 in the presence of quenched disorder. The bifurcation diagram is similar to the one observed in the globally coupled networks [28]. However, the splay state found for small α-values has been replaced by a fluctuating asynchronous state (AS), where the average field Ē(t) is only approximately constant (the difference between min Ē(n) and max Ē(n) is of the size of the symbols). The periodic collective state (PS) is analogously affected by small irregular fluctuations. This strong similarity Figure 3.1: Minimal and maximal values of Ē(n) for different α values for N = 1, 600 with quenched disorder and g = 0.4.

between globally coupled and the diluted networks is not surprising: for any finite value of f, upon increasing N, the differences among the fields E_{i}(t) should progressively disappear. In fact, in the limit N → ∞, we expect randomly diluted networks to behave as fully coupled ones, provided that variables and parameters are properly rescaled. More precisely, the average field Ē(t) in a network with a fraction f of missing links and a coupling constant g, is expected to be equivalent to the neural field E generated in a fully coupled network with coupling constant g(1 − f). This expectation is confirmed by the critical value ᾱ_{c} separating the two dynamical phases. As shown in Fig. 3.1, for g = 0.4 and f = 0.2, ᾱ_{c} ∼ 6.8, a value which coincides, within the statistical error, with the critical value found in a globally coupled network with g = 0.32 = 0.4 × (1 − 0.2)[33]. A further consequence of the fact that the dynamics of disordered networks reduces in the thermodynamic limit to that of homogeneous fully coupled ones is that the sample-to-sample fluctuations typical of quenched disorder should vanish as well. For this reason we have not performed averages of different realizations of the disorder.

A more detailed representation of the quenched dynamics above ᾱ_{c} is obtained by looking at the projection in the plane (Ē(n), P̄(n)). In Fig. 3.2 data sets are shown for α = 9 and increasing values of N: panels a and b correspond to the annealed and quenched case, respectively. This allows seeing that the two kinds of disorder yield indeed qualitatively similar results.

More precisely, the phase points cluster around closed curves, revealing a “noisy” periodic dynamics. In the annealed case, the fluctuations can be attributed to the stochasticity of the evolution rule. Surprisingly, they are even larger in the quenched case, which correspond to a deterministic dynamics, in which case they must be attributed to the presence of deterministic chaos.

Upon increasing N, the amplitude of the fluctuations decreases and the attractor shape converges to an asymptotic curve, which should correspond to that of a ho-Figure 3.2: Ē(n) versus P̄(n) for various system size: (a) annealed disorder and (b) quenched case here is shown also the annealed result for N = 100, 000. All data refers to α = 9.

mogeneous fully coupled network with a properly rescaled coupling strength g. This expectation is indeed confirmed in Fig. 3.3, where the attractor of a homogeneous network (with N = 1600 and g = 0.32) superposes to that of the annealed dynamics (with N = 10^{5} and g = 0.4). In the quenched case, the asymptotic shape is the same, but the convergence is slower. Since the representation of the attractor in the plane (gĒ, gP̄) seems to be the correct one to better appreciate the thermodynamic limit, in the following chapters we will perform a macroscopic characterization of the dynamical states of the network by using the average fields Ē, P̄ properly rescaled by the synaptic strength.

Figure 3.3: gĒ(n) versus gP̄(n). The (blue) curve represents the attractor of a globally coupled system with N = 1, 600 and g = 0.32, while the black dots refer to the annealed case with N = 100, 000 and g = 0.4. For both cases α = 9 and a = 1.3.

In order to investigate quantitatively the scaling behavior of the finite-size corrections, we have studied the N dependence of

ΔP̄ = ⟨P̄⟩(N) − ⟨P̄⟩(∞)  (3.15)

where the angular brackets denote the (time) average of P̄-value of all configurations falling within the lower gĒ(n)-window [0.36, 0.44]. Since the asymptotic value ⟨P̄⟩(∞) is independent of the setup, we have extrapolated it in the simpler context of a fully coupled network with g = 0.32 As a result, we find that ΔP̄ always converges to zero as power-law, N^{−β}, with β = 1 in the fully coupled network, β = 0.55 for quenched disorder, and β = 0.82 for annealed disorder (see Fig 3.4a). These latter values have to be considered as approximate estimates, affected both by statistical errors and finitesize corrections. A more accurate estimate is beyond our computational capabilities. However, for the present purpose, the relative differences are large enough to conclude that quenched disorder is characterized by a slower convergence.

Figure 3.4: (a) The quantity ΔP̄ (see Eq. 3.15) versus the network size; (black) triangles refer to the quenched case, while (red) squares to the annealed one: data fits are reported as dashed lines. The reported values have been obtained by averaging over 6 × 10^{6} to 3 × 10^{7} consecutive spikes. (b) Average standard deviation < σ > (see Eq. 3.16) vs N : circles and diamonds refer to the annealed and quenched case, respectively. The reported data have been averaged over M = 900 iterates. All data refers to f = 0.2, α = 9, g = 0.4 and a = 1.3.

As a further test of the collective dynamics, we have computed the standard deviation σ(n) of the instantaneous fields

σ(n) = (\frac{∑_{i=1}^{N} E_{i}^{2}(n)}{N} − Ē^{2}(n))^{1/2} ,  (3.16)

which is a measure of the degree of synchronization among the various fields. In Fig. 3.4b we plot the time average ⟨σ⟩ for annealed and quenched disorder: in both cases ⟨σ⟩ decreases with N as N^{−1/2}. These results confirm once more that in the limit N → ∞ the neural field dynamics converges to that of homogeneous networks, irrespectively of the disorder.

3.4 Lyapunov Analysis

In order to provide a more detailed characterization of the macroscopic as well as of the microscopic dynamics, in this section we analyze the Lyapunov spectra for the fully-coupled network and for its disordered variants.

3.4.1 Globally Coupled Network

In this case, the fields seen by the neurons are equal to one another and it is therefore sufficient to introduce a single pair of field variables E and P. The corresponding stability analysis for the splay state has been analytically carried out in Ref. [33], finding that the spectrum of Floquet exponents is composed of a band of values of order 1/N^{2}, which become of order O(1) at one band extreme, plus two isolated eigenvalues associated with the field dynamics. Therefore, it is convenient to start the numerical analysis by testing the effect of attaching a pair of field variables to each neuron. The results plotted in Fig. 3.5a refer to the Lyapunov spectrum of the splay state for α = 3 and g = 0.4: one can observe two bands and two isolated exponents. The first band is composed of N − 1 nearly vanishing exponents, while the two isolated exponents concern the dynamics of the coupling field. The second band, composed of 2(N − 1) exponents quantifies the transversal stability of the synchronization manifold E_{i} = E_{1}, P_{i} = P_{1} , ∀i, where, without loss of generality, we consider the field of the 1st neuron as the reference variable. In fact, if one introduces the suitable difference variables,

w_{i} = E_{i} − E_{1} z_{i} = P_{i} − P_{1} i = 2, ⋯, N  .  (3.17)

the corresponding tangent dynamics is ruled by the equations

δw_{i}(n + 1) = δw_{i}(n)e^{−ατ(n)} + δz_{i}(n)τ(n)e^{−ατ(n)}  (3.18)
δz_{i}(n + 1) = δz_{i}(n)e^{−ατ(n)} i = 2, ⋯, N .  (3.19)

This shows that the linearized dynamics of each pair of twin variables (δw_{i}, δz_{i}) is decoupled from the rest of the network. Its stability can be evaluated by solving the corresponding two-dimensional eigenvalue problem. Accordingly, we expect two bands of N − 1 eigenvalues each. However, since the two eigenvalues of Eq. (3.18) are both equal to −α, we do have a single band, as indeed observed in Fig. 3.5a. As a last check, we have verified the 1/N^{2} scaling of the first band of the splay state. The nice overlap of the three sets of exponents corresponding to N = 100, 200, and 400 with the analytical estimate [46] reported in Fig. 3.5c confirms the reliability of the numerical code.

The stability of CPOs^{1} arising when α > α_{c} = 8.34(1) [28] has not been investigated so far. However, the above arguments concerning the presence of a single degenerate band still hold, as confirmed by the spectrum plotted in Fig. 3.5b, which corresponds Figure 3.5: Lyapunov exponents λ_{k} versus k/N for the fully coupled network. Complete spectrum for N = 50 for α = 3 (a) and α = 9 (b). Rescaled Lyapunov eigenvalues of the first band λ_{k} × N^{2} for different system sizes, namely N = 50 (filled black circles), 100 (open red squares), 200 (green open triangles) and 400 (blue crosses), for α = 3 (c) and α = 9 (d). For α = 3 in (c) also the analytical expression (dashed black line) reported in Eq. (29) in [46] is shown. The data has been obtained by following the evolution in the tangent space of the event driven map for a number of consecutive spikes of the order of 10^{8} − −10^{9}, after discarding a transient composed by 50, 000 × N spikes.

1In th context of the fully synchronized network the CPO state corresponds to the partial synchronization discussed in Section 3.1.

to g = 0.4, and α = 9. This adds to the first band that is again located just below zero, with the exception of the first exponent which is exactly equal to zero, as a result of the quasi-periodic dynamics of the single neurons (the second zero Lyapunov exponent has been discarded while taking the Poincaré section).

Finally, in Fig. 3.5d we have plotted the Lyapunov exponents of the first band multiplied by N^{2}. The tendency to overlap of the spectra obtained for N = 50, 100, 200, and 400 indicate that we are again in the presence of a 1/N^{2} scaling as for the splay states, although finite-size corrections are more relevant in this case. Accordingly, we conjecture that the 1/N^{2} dependence is a general property, which depends on the shape of the velocity field (see [46]), rather than on the structure of the solution itself. This point will be examined in details in Chapter 2.

3.4.2 Diluted Network

Since the collective solutions observed in the globally coupled network are characterized by many weakly stable directions, it is reasonable to expect that generic perturbations Figure 3.6: (a) Maximal Lyapunov exponent λ_{1} as a function of the percentage f of broken links for N = 800 and α = 9. (b) Lyapunov spectrum {λ_{k}} versus k/N for α = 3 with networks of size N = 50 (filled black circles) and N = 100 (open red squares) and 20% of broken links. The data has been obtained by following the evolution of the event driven map and of the associated linearized equations ruling the evolution of the Lyapunov vectors for 3 − 4 × 10^{8} consecutive spikes, after discarding a transient composed by 300, 000 − 600, 000 spikes. The data refers to quenched disorder.

of the network dynamics yield a chaotic evolution. On the other hand, in the thermodynamic limit we expect the inhomogeneities induced by a random dilution to vanish. In fact, we have already seen (in the previous section) that a macroscopically periodic dynamics is eventually recovered.

First of all we have verified that the network dynamics is chaotic for finite N, as soon as f > 0. In Fig. 3.6a one can appreciate the growth of the maximal (positive) Lyapunov exponent with the fraction of missing links f (in the quenched case).

For α = 3, analogously to the globally coupled case, the spectrum is composed of two bands and two isolated eigenvalues (see Fig. 3.6b). On the one hand, the degeneracy among the exponents lying in the second band is lifted (as a consequence of the disordered network structure) and the second band acquires a finite width. On the other hand, the comparison between the spectrum obtained for N = 50 and N = 100, reveals that the band width shrinks with N around the value λ = −α. This result again confirms that the inhomogeneities of the disordered network vanish in the thermodynamic limit. A similar scenario occurs also in the CPO regime, as well as for annealed disorder.

Furthermore, in Fig. 3.7 (where we plot the first band of the Lyapunov spectrum for quenched and annealed disorder, in both dynamical phases), we see a variable number of positive exponents, indicating that finite disordered networks are typically chaotic. More precisely, panels a and b of Fig. 3.7, refer to quenched disorder. In both cases, we present the spectra resulting from three different realizations of the disorder. We see that sample-to-sample fluctuations are quite relevant in the second part of the band, while they affect less the largest exponents. Another qualitative observation Figure 3.7: (Color online) First band of the Lyapunov spectra {λ_{k}} versus k/N for g = 0.4 with 20% of broken links for quenched disorder with α = 3 (a) and 9 (b) and for annealed disorder with α = 3 (c) and 9 (d). For the quenched disorder only spectra corresponding to system size N = 50 are shown but for three different random network configurations, while in the annealed case spectra for network sizes N = 50 (filled black circles) and N = 100 (open red squares) are reported. The data has been obtained by following the evolution of the system and of its linearized copies for time lags similar to those reported in the caption of Fig. 3.6.

concerns the shape of the spectrum that seems to be smoother in the presence of CPOs. Unfortunately, it is almost impossible to perform a quantitative scaling analysis of the spectrum (given the need to average over different realizations and to consider yet larger network sizes). In the annealed case, Fig. 3.7c/d, there is no need to average over different realizations of the disorder and the Lyapunov spectra turn out to be smooth. However a scaling analysis is still beyond our computer capabilities: in both cases we report the spectra obtained for N = 50 and N = 100 which are far from exhibiting a clean scaling behavior (for this reason we do not even dare to formulate a conjecture).

Given the difficulties of dealing with the whole Lyapunov spectrum, we have concentrated our efforts on the maximum exponent λ_{1}, since we are thereby able to study larger networks. We limited ourselves to studying the CPO regime. As shown in Fig. 3.8, λ_{1} decreases with N as a power law: in the annealed case, λ_{1} ≈ 1/N, while in the quenched case λ_{1} ≈ 1/\sqrt{N}. This is at variance with the Kuramoto model, where in the “equivalent” quenched case, a 1/N behavior has been observed [61]. Accordingly, although annealed networks are more chaotic than the quenched ones over the numerically accessible network sizes, we conclude that the opposite will be true for yet larger networks.

Figure 3.8: Maximal Lyapunov exponent λ_{1} versus N: (a) annealed case; (b) quenched case. In the quenched case also a power-law fit is reported (red dashed line) with decay exponent β = 0.51 ± 0.01. The data has been obtained by following the evolution in the tangent space of the event driven map for 10^{8} − −10^{9} consecutive spikes, after discarding a transient of 200, 000 spikes. Both figures refer to 20 % of dilution, with parameters α = 9, g = 0.4 and a = 1.3.

3.5 Discussion

Our numerical analysis suggests that in the thermodynamic limit, a random, uncorrelated network behaves like a homogeneous globally coupled system with a rescaled value of the coupling constant to account for the different fraction of active links. This is because each neuron receives a random series of spikes: when the number of neurons increases, the number of spikes per unit time increases as well, while the relative fluctuations decrease and all neurons are affected by increasingly similar forcing fields (for a discussion on the thermodynamic limit in neural networks see also [35, 62]). While there are little doubts that the collective motion is independent of the presence of disorder, less compelling is the evidence that the same is true for the stability properties. This is because numerical simulations are computationally expensive and it is not possible (at least for us) to study the scaling behavior of the entire Lyapunov spectrum in disordered systems. This is doable in homogeneous networks, thanks also to a faster convergence to the thermodynamic limit, and we have indeed observed that the first band of the Lyapunov spectrum scales as 1/N^{2}, exactly as in the splay state, a case that is analytically known [46]. In disordered networks we have nevertheless been able to study the scaling behavior of the maximal Lyapunov exponent, finding that it is positive and scales as 1/N^{β}, with β close to 1/2 (resp. 1) in the quenched (resp. annealed) case. This means that deterministic chaos disappears in the thermodynamic limit. This observation is consistent with the fact that the collective motion is periodic: a periodically forced phase oscillator (such as a LIF neuron) cannot be chaotic. However, nothing seems to prevent the collective motion itself to be chaotic. Whether different connection topologies have to be invoked or yet unidentified constraints forbid this to occur, remains as an open question and will be clarified in the following chapters where different topologies will be considered. In particular collective chaos will be discuss in cap. 6.

The results presented in this Chapter have been published in Ref. [47].

Coherent periodic activity in Erdös-Renyi neural networks

4.1 Introduction

Neural collective oscillations have been observed in very many context in brain circuits, ranging from ubiquitous γ oscillations to θ rhythm in the hippocampus. The origin of these oscillations is commonly associated with the balance between excitation and inhibition in the network, while purely excitatory circuits are believed to lead to ”unstructured population bursts” [9]. However, recent ‘ex vivo” measurements performed on the rodent neocortex [52] and hippocampus [63] in the early stage of brain maturation reveal coherent activity patterns, termed as Giant Depolarizing Potentials. These collective oscillations emerge despite the fact that the GABA transmitter has essentially an excitatory effect on immature neurons [53]. Therefore, also in purely excitatory networks one can expect non trivial dynamics at a macroscopic level.

Numerical and theoretical studies of collective motions in networks of simple spiking neurons have been mainly devoted to balanced excitatory-inhibitory configurations, e.g see [40] and references therein. Only few studies focused on coherent periodic activity in fully coupled excitatory networks of leaky integrate-and-fire (LIF) neurons. These analysis revealed a regime characterized by a partial synchronization (PS) at the population level, while the single neurons perform quasi-periodic motions [28]. As shown in the previous Chapter the PS regime is quite robust to perturbations, since it survives to moderate levels of noise or dilution [29, 47].

Furthermore, in their recent study Bonifazi et al. [63] found that the functional connectivity of developing hippocampal networks is characterized by a truncated powerlaw distribution of the out-degrees with exponent γ = 1.1 − 1.3. This scaling has been shown to hold over one/two decades, thus not ensuring a scale-free distribution for the links over all scales, but surely indicating the presence of a large number of hub neurons, namely cells characterized by a high connectivity. At early developmental stages of the brain, GABAergic hub interneurons, performing complex excitatory/shunting inibitory actions [53], seem to be responsible for the orchestration of the coherent activity of hippocampal networks [63]. The relevance of hubs in rendering a neural circuit extremely hyperexcitable has been also demonstrated in simulation studies of a realistic model of the epileptic rat dentate gyrus, even in the absence of a scale-free topology [64].

Motivated by these studies, but without attempting to reproduce the experimental results, we focus on a very preliminary question: to which extent is the macroscopic neural dynamics influenced by the average degree of connectivity of the neurons ? Our specific aim is to analyze the key ingredients leading to the onset of coherent activity, as opposed to asynchronous dynamics, for different network size and topology.

Specifically, we consider the transition from an asynchronous regime to partial synchronization in the excitatory LIF pulse coupled neural networks introduced by Abbott and van Vreeswijk [25]. At variance with previous works, we study dynamical evolution on random Erdös-Renyi networks with an average connectivity growing (sub)-linearly with the network size N. In particular, we consider average connectivities scaling as N^{2−γ}, thus exhibiting the same system size dependence of connectivities associated to truncated power-law distributions with decay exponent 1 < γ < 2. In the limit γ → 1 the massively connected network, where the connectivity is proportional to N, is recovered; while for γ → 2 a sparse network, where the average probability to have a link between two neurons vanishes in the thermodynamic limit [35], is retrieved. The topology of Erdös-Renyi networks is modified by varying the parameter γ in the interval [1 : 2], in particular as far as γ ≤ 2 trees and cycles of any order are present in the network, while for γ → 1 complete subgraphs of increasing order appear in the system [65]. In the chapter we will study how and to which extent these topological modifications influence the macroscopic dynamics of the network, with particular emphasis on the transition from asynchronous to (partially) synchronous collective dynamics. In fact at a macroscopic level we observe two distinct dynamical phases: an asynchronous state corresponding to a desynchronized dynamics of the neurons (splay state) and a regime of partial synchronization associated with a coherent periodic activity of the network. At low connectivity the system is in an asynchronous state, while PS emerges above a certain critical average connectivity ⟨k⟩_{c}. For sufficiently large networks, ⟨k⟩_{c} saturates to a constant value suggesting that a minimal average connectivity is sufficient to observe coherent activity in systems of any size irrespectively of the kind of considered network: sparse or massively connected. However, this value depends on the nature of the synapses: reliable or unreliable. For unreliable synapses the critical value required to observe the onset of macroscopic behaviors is noticeably smaller than for reliable synaptic transmission. Due to the disorder present in the system, for finite number of neurons we have inhomogeneities in the neuronal behaviors, inducing a weak form of chaos, which vanishes in the thermodynamic limit. In such a limit the disordered systems exhibit regular (non chaotic) dynamics and their properties correspond to that of a homogeneous fully connected network for any γ-value. Apart for the peculiar exception of sparse networks, which remain intrinsically inhomogeneous and chaotic for any system size.

4.2 Model and Methods

4.2.1 The Model

We study a network of N LIF neurons with the membrane potential u_{i}(t) ∈ [0 : 1] of the neuron i evolving as described in Eq. (3.1). The synaptic current is still given by gE_{i}(t), but the field E_{i}(t), which represents the linear superposition of the pulses s(t) received by neuron i in the past, is now weighted by the effective number of in-degree connections, namely

E_{i}(t) = \frac{1}{k_{i}} ∑_{n|t_{n}<t} C_{j,i}Θ(t − t_{n})s(t − t_{n}) ,  (4.1)

where k_{i} is the number of pre-synaptic neurons connected to the neuron i (i.e. the in-degree of neuron i) and Θ(t) is the Heaviside function. The connectivity matrix C_{j,i} appearing in Eq. (4.1) has entry 1 (resp. 0) depending if the pre-synaptic firing neuron j is connected (resp. not connected) to neuron i and in general it is not symmetric.

Following van Vreeswijk [28] we assume, for the single pulse emitted at t = 0, the shape s(t) = α^{2}t exp(−αt). The explicit equation (4.1) can be thus rewritten as an implicit ordinary differential equation:

Ë_{i}(t) + 2αĖ_{i}(t) + α^{2}E_{i}(t) = \frac{α^{2}}{k_{i}} ∑_{n|t_{n}⟨t} C_{j,i}δ(t − t_{n}) .  (4.2)

The continuous time evolution of the network can be transformed in a discrete time event-driven map by integrating Eq. (4.2) from time t_{n} to time t_{n+1}, t_{n} being the time immediately after the n-th spike emission.

Analogously to the Eqs. (3.4-3.6) and introducing the explicit dependence on k_{i}, the event-driven map read as:

E_{i}(n + 1) = E_{i}(n)e^{−ατ(n)} + P_{i}(n)τ(n)e^{−ατ(n)}  (4.3)
P_{i}(n + 1) = P_{i}(n)e^{−ατ(n)} + C_{m,i} \frac{α^{2}}{k_{i}}  (4.4)
u_{i}(n + 1) = u_{i}(n)e^{−τ(n)} + a(1 − e^{−τ(n)}) + gH_{i}(n) .  (4.5)

where P_{i} ≡ αE_{i} +Ė_{i} is an auxiliary variable and τ(n) = t_{n+1} −t_{n} is the inter-spike time interval. This can be determined by solving the standard implicit relationship (3.7).

The explicit expression for the nonlinear function H_{i}(n) appearing in (4.5) is given by Eq. (3.8).

The evolution of the system is modeled, as in the previous chapter, with a discrete time map of 3N − 1 variables, {E_{i}, P_{i}, u_{i}}. In fact the construction of the event-driven map implies that the membrane potential of the neuron firing at time t_{n+1} is always equal to the threshold value, therefore this degree of freedom does no take part to the discrete time dynamics. At variance with the study reported in Chapter 3.1 the pulse amplitudes, appearing in Eqs. (4.1),(4.2) and (4.4), are normalized by the in-degree k_{i} of neuron i and not by the total number of neurons N. The model parameters were fixed as a = 1.3, g = 0.4 and α = 9, in order to ensure the emergence of a PS regime in the corresponding fully coupled network [28].

4.2.2 The Connectivity Matrix

In order to reveal the possible influence of network topology on the dynamics of the system, we consider Erdös-Renyi random network with an average connectivity growing (sub)-linearly with N [65]. The scaling of the average connectivity with the system size is controlled by a parameter γ and in the next Sections we will analyze how the dynamical properties of the network are influenced by the values of this parameter, i.e. by different ways of approaching the thermodynamic limit. In particular, we consider a directed Erdös-Renyi random graph, where the distribution of links k per neuron is well approximated by a Poisson distribution [65], namely

P (k) = e^{−⟨k⟩} \frac{⟨k⟩^{k}}{k!} .  (4.6)

According to Eq. (4.6) the degree distribution is completely defined once the value of ⟨k⟩ is given. To make a connection with recent experimental findings obtained for excitatory networks revealing power-law distributed connectivity [63], we choose

⟨k⟩ = \frac{p}{2 − γ} [N^{2−γ} − 1] .  (4.7)

The above expression corresponds to the average connectivity associated to a truncated power-law distribution P (k) = pk^{−γ}, once assumed that each neuron has at least one connection, and where p is a proportionality factor. We limit our analysis to 1 < γ < 2, since Bonifazi et al [63] have measured quite low values for the exponent γ, namely γ = 1.1 − 1.3, suggesting the existence of a large number of highly connected neurons, hubs, in the network.

As a matter of fact, by choosing for ⟨k⟩ the expression (4.7), the probability of existence of an unidirectional link connecting neuron j to i (i.e. the probability to have C_{j,i} = 1) is:

P r(N, γ) = \frac{⟨k⟩}{N} = \frac{p}{2 − γ} [N^{1−γ} − \frac{1}{N} ] .  (4.8)

In the limit γ → 1 the massively connected network is recovered [35], since the average connectivity ⟨k⟩ = p × (N − 1) is proportional to the system size and P r(N, 1) = p(1 − 1/N) is, apart finite-size corrections, constant and coincident with p. For γ < 2 (resp. γ > 2) the average number of synaptic inputs per neuron will grow (resp. decrease) with N, in the limiting case γ = 2 a sparse network will be essentially obtained ^{1} since ⟨k⟩ = p ln N will vary in a limited manner with respect to the system size. Indeed, by varying N by three orders of magnitude from 100 to 100,000 the value of ⟨k⟩ will modify from 3.7 to 9.2 with p = 0.8. In the following, we study networks of various sizes N, ranging from N = 100 to N = 200, 000, for different γ-values in the interval [1, 2]. The value p is usually fixed to 0.8, apart for γ = 1 (constant probability case) when p indeed coincides with the probability to have an unidirectional link. In such a case the dependence on p is examined in details in Section 3.

1A sparse network is usually defined as a network where the average connectivity remains constant in the thermodynamic limit, see [35].

Let us stress that in the present study the distributions of pre-synaptic (in-degree) and of post-synaptic (out-degree) connections are identical, and this is guaranteed by the above outlined procedure to determine unidirectional links for Erdös-Renyi networks.

In what follows, we consider two different ways to select the random connectivity matrix: in the first case the synaptic connections are randomly chosen before the simulation and they do not change in time (quenched disorder); while for the second procedure the post-synaptic neurons receiving the emitted pulse (i.e. the efferent connections) are randomly selected each time a neuron fires (annealed disorder). The latter choice can be justified from a physiological point of view by the fact that the synaptic transmission of signals is an unreliable process [66]. It should be noticed that in the annealed case, since the network modifies in time, the pulse amplitudes are normalized by the average in-degree ⟨k⟩ and not by k_{i} as in the quenched case.

4.2.3 Characterization of Macroscopic Attractors

In order to perform a macroscopic characterization of the dynamical states of the network we employ the average fields (3.14) rescaled by the synaptic strength as done in the previous chapter.

As a measure of the level of homogeneity among the neurons of the network, we consider the standard deviation σ(t) of the fields {E_{i}(t)} acting on the single neurons (3.16) rescaled by the synaptic strength, analogously to the rescaling of the average fields; for completely homogeneous systems, such as globally coupled networks, E_{i}(t) ≡ Ē and σ ≡ 0.

The degree of synchronization among the neurons is quantified by the modulus of the Kuramoto order parameter^{2} defined in Eq. (1.8)

R(t) = |⟨\frac{1}{N} ∑_{j=1}^{N} e^{iθ_{j}(t)}⟩| .  (4.9)

2In the original definition reported in [67], T_{q,n} = T_{j,n}, thus ensuring that the phase is bounded between 0 and 2π. However this definition requires to memorize all the spiking events within a certain time interval and then estimate the phases and R(t) over such time span. Our definition does not guaratee that the phase is strictly bounded within the [0 : 2π]-interval, but it allows to estimate R(t)

during the simulation run without memorizing any event. Indeed we have verified that, at least in the studied cases, no peculiar differences are observable by employing the two phase definitions and this is probably due to the fact that the dynamics of all the neurons in the networks are statistically equivalent.

A non zero-value of R represents an indication of partial synchronization among the neurons, while a vanishingly small R ∼ 1/\sqrt{N} is observable for asynchronous states in finite systems [31].

4.2.4 Lyapunov Analysis

The dynamical microscopic instabilities of a system can be characterized in terms of the maximal Lyapunov exponent λ_{1}: a positive λ_{1} being a measure of the degree of chaoticity of the considered system. In particular, the stability of (4.3,4.4,4.5) can be analyzed by following the evolution of infinitesimal perturbations in the tangent space. The corresponding equations are obtained by linearizing (4.3) and are given by (3.9,3.10,3.11). In addition to this an explicit expression of δτ(n) can be obtained by differentiating Eqs. (3.7,3.8), as reported in Eq. (3.12).

To estimate the maximal Lyapunov exponent λ_{1} we follow the usual procedure [56], in particular we measure the exponential growth rate of an infinitesimal perturbation by iterating the system (3.9). In order not to loose accuracy in the estimation the amplitude of the perturbation should be rescaled at regular intervals, in our case we have rescaled every 1, 000 spikes.

4.3 Phase Diagram

In this model two different macroscopic regimes can be observed: the asynchronous state and the partial synchronization. The asynchronous state is characterized by an incoherent dynamics of the neurons in the network leading to a spot-like macroscopic attractor in the (gĒ, gP̄) plane and an almost constant average field Ē [25], while to the coherent PS regime corresponds a closed curve attractor [29] and a periodic behavior of Ē in time as shown in Fig. 4.1a and b. The incoherent and coherent neural dynamics can be clearly appreciated, in the two regimes, also at microscopic level by examining the corresponding raster plots reported in Fig. 4.1c and 4.1d.

As a first aspect, we will investigate the occurrence of asynchronous and partially synchronized states for finite networks and different average connectivity ⟨k⟩. In order to distinguish the two regimes we have examined the macroscopic attractor shape, the extrema values of the average field Ē, and the synchronization parameter R as a function of N.

4.3.1 Erdös-Renyi networks with constant probability

Analogously to the study reported in Chapter 3 we consider Erdös-Renyi networks with unidirectional links chosen at random with a constant probability p for any network size N. This amounts to consider the limiting case γ → 1 and to have an average connectivity scaling linearly with the size, i.e. ⟨k⟩ = pN. Let us firstly examine the minima and maxima of Ē for a network of size N by varying the probability p Figure 4.1: Characterization of the asynchronous state (AS) and of the partially synchronized (PS) one in terms of macroscopic fields and single neuron dynamics. Panel a: macroscopic attractors in the (gĒ, gP̄) plane. Panel b: the average field gĒ as a function of time. Panel c and d: raster plots. The data refer to Erdös-Renyi networks with ⟨k⟩ = p × N, quenched disorder, and N = 1, 600, the black (resp. red) symbols correspond to the asynchronous state observable for p = 0.2 (resp. PS for p = 0.7).

between 0 and 1, as shown in Fig. 4.2a for N = 1, 600. At low probability one has an asynchronous state while the PS regime emerges only for sufficiently large p, both in the quenched and annealed case. This result is related to the fact that the presence of noise or disorder reduces the coherence needed to observe the PS regime [29]. As a matter of fact increasing p (i.e. diminishing the number of broken links in the network) the attractor size increases and finally reaches the fully coupled result. The degree of coherence can be measured in terms of the average synchronization indicator ⟨R⟩. As shown in Fig. 4.2b the system coherence steadily increases with p, except in the asynchronous regime where ⟨R⟩ ∼ 0, apart for finite size fluctuations.

Let us now report the phase diagram for the macroscopic activity of the network in the (N, ⟨k⟩) plane, for both annealed and quenched disorder. Increasing the average connectivity, keeping the system size fixed, leads to a transition from asynchronous to partially synchronized regimes (see Fig. 4.3). The transition occurs at a critical average connectivity ⟨k⟩_{c}, indicates by the asterisks connected by a black solid line in Fig. 4.3, which for low N increases steadily with N, but eventually saturates for N > 10, 000 to an asymptotically constant value which depends on the disorder realization: namely, ⟨k⟩_{as} = 725 ± 25 for the quenched case and ⟨k⟩_{as} = 225 ± 25 in the annealed one ^{3}.

3We have verified that the above reported scenario is not modified by considering fully homogeneous connectivity with P (k) = δ(k − k_{0}), the only difference concerns the asymptotic value ⟨k⟩_{as} at which the transition from asynchronous dynamics to partial synchronization occurs. As a matter of fact ⟨k⟩_{as} is slightly smaller than for Erdös-Renyi distributions due to the absence of fluctuations among the in-degree values in the homogeneous case.

Figure 4.2: a) Minima and maxima of the average field gĒ as a function of p, the circles joined by solid lines refer to the annealed disorder, while the stars connected by the dashed line to the quenched case. The dot-dashed (red) lines indicate the fully coupled results (corresponding to p ≡ 1). b) Synchronization indicator R averaged over time as a function of p for the quenched case. The data refer to Erdös-Renyi networks with constant probability (CP) and size N = 1, 600 and have been estimated, after discarding a transient of 4 × 10^{7} spikes, by averaging over a train of 1 − 2 × 10^{7} spikes.

Figure 4.3: Phase diagram for the macroscopic activity of the network in the (N, ⟨k⟩) plane: a) annealed disorder and b) quenched disorder. The (black) asterisks connected by the solid (black) line correspond to the transition values ⟨k⟩_{c} from asynchronous (AS) to partially synchronized (PS) regime estimated for Erdös-Renyi networks with constant probability. The other symbols refer to Erdös-Renyi with γ > 1: solid (resp. empty) symbols individuate asynchronous (resp. partially synchronized) states. In particular, (blue) squares refer to γ = 1.1, (red) circles to γ = 1.3, (green) diamonds to γ = 1.5, and (magenta) triangles to γ = 1.7. The reported data are relative to the state of the network after discarding transients ranging from 2 × 10^{7} spikes at the smaller sizes to 3 × 10^{8} spikes for the larger networks.

4.3.2 Erdös-Renyi networks with γ-dependent probability

To verify the generality of these results we investigate Erdös-Renyi networks with γ-dependent probability. In particular, we have estimated, for system sizes in the range 100 < N < 200, 000, the macroscopic attractors for various γ-values (namely, γ = 1.1, 1.3, 1.5 and 1.7), after discarding long transient periods. For small system sizes the network is in the asynchronous regime which is characterized by a spot-like attractor in the (gĒ, gP̄)-plane. For larger number of neurons, PS emerges in the system characterized by closed curve (macroscopic) attractors. Furthermore, similarly to the results reported in Ref. [47] and in Section 3.3, increasing N the curves tend to an asymptotic shape, corresponding to the fully coupled attractor, while fluctuations diminish. To exemplify this point various macroscopic attractors for γ = 1.3 and annealed disorder are reported in Fig. 4.4a.

As already reported for the constant probability networks (see Chapter 3), the systems with annealed disorder converge more rapidly with N towards the asymptotic fully coupled attractor with respect to the quenched case, as shown in Fig. 4.4b for γ = 1.1, 1.3, 1.5. Increasing γ we observe that the transition from the asynchronous to the partially synchronized regime occurs at larger and larger system size, both for annealed and quenched disorder. The results for γ = 1.7 are not shown in Fig. 4.4b, since for all the examined network sizes (up to N = 200, 000) Ē extrema coincide within the error bar, indicating that the system is in the asynchronous regime ^{4}.

Asynchronous (resp. partially synchronized) regimes are reported in the phase diagram displayed in Fig. 4.3 as filled (resp. empty) symbols for the investigated γ-values. The critical line ⟨k⟩_{c} = ⟨k⟩_{c}(N) (indicated by the solid black line in Fig. 4.3) denoting the transition from asynchronous to partially synchronized regime coincides with that determined previously for constant probability networks. We can thus safely affirm that the dynamical regimes of the Erdös-Renyi networks depend, at a macroscopic level, simply on the average connectivity, once the system size N is fixed. This could be expected from the fact that for Erdös-Renyi networks the distribution of links per neuron is completely determined by ⟨k⟩ (see Eq. (4.6)). However, the independence of ⟨k⟩_{c} from N at N > 10, 000 is unexpected and suggests that coherent behaviors, like PS regimes, can be observed in networks of any size for vanishingly small relative connectivity ⟨k⟩/N.

N_{c} = [\frac{p}{2 − γ} ⟨k⟩_{as} + 1]^{1/(2−γ)} ;  (4.10)

4A simple estimation of the minimal network size needed to observe PS for generic γ-distributions can be obtained by inverting Eq. (4.7) for the asymptotic values ⟨k⟩_{as}, namely for γ = 1.7 this leads to N_{c} = 131, 784, 000 (resp. N_{c} = 2, 740, 117) in the quenched (resp. annealed)

case. This confirms that coherent activity is unobservable for γ = 1.7 with our computational resources for the chosen coupling parameter, namely g = 0.4.

Figure 4.4: a) Macroscopic attractors in the (gĒ, gP̄) plane for a Erdös-Renyi networks with γ = 1.3 and annealed disorder, the curves from the interior to the exterior corresponds to increasing system sizes, from N = 1, 600 to 25,600. The most external (black) curve refers to a fully coupled network with N = 3, 200. b) Minima and maxima values of the average field gĒ as a function of N for various γ-values: namely, (black) circles γ = 1.1, (red) triangles γ = 1.3 and (blue) squares γ = 1.5. The empty (resp. filled) symbols refer to annealed (resp. quenched) disorder. The dot-dashed (magenta) lines indicate the fully coupled values.

4.4 Network Homogeneity

We observe that for any considered γ the fields (gE_{i}, gP_{i}) associated with the different neurons tend to synchronize for increasing N. Therefore in the thermodynamic limit (N → ∞) disordered networks tend to behave as fully coupled ones, where all the neurons are equivalent and a single field is sufficient to describe the macroscopic evolution of the system.

In order to quantify the level of homogeneity among the various neurons, we measured the standard deviation (Eq. (3.16)) relative to the fluctuations of the different local fields E_{i} with respect to the average field Ē. In Fig. 4.5a we plot the time average of the standard deviation rescaled by the synaptic strength, g⟨σ⟩, for annealed disorder and various γ values. We observe a power law decay ⟨σ⟩ ∝ N^{−β}, where the exponent β depends on the γ-parameter as β = 1 − γ/2 (see Fig. 4.5b). Furthermore, for the limiting case γ = 2.0 the decay of ⟨σ⟩ is consistent with a scaling 1/\sqrt{ln N} as displayed in the inset of Fig. 4.5a. Altogether, the reported dependencies suggest the following relationship to hold:

⟨σ⟩ ∝ \frac{1}{\sqrt{⟨k⟩}} ;  (4.11)

thus fields fluctuations are driven by the average in-degree value irrespectively of the total number of neurons. These results confirm once more that in the limit N → ∞ the neural field dynamics converges to that of homogeneous networks, for both quenched and annealed disorder.

Figure 4.5: (Color online) a) Average standard deviation g⟨σ⟩ versus the system size N for annealed disorder and various γ-values: γ = 1.1 (blue) squares, 1.3 (red) circles, 1.5 (green) diamonds and 1.7 (magenta) triangles. The dashed line represents best fits with a power-law N^{−β} to the reported data. The data in the inset (black asterisks) refers to a γ = 2.0, the dashed line is a guide for the eyes. b) Power-law exponents β for annealed (black circles) and quenched (red triangles) disorder in the network as a function of the parameter γ. The dashed (blue) line refers to the linear law β = 1−γ/2. The reported data have been estimated by averaging over trains made of 2 × 10^{6} − 10^{8} spikes, after discarding transients of 4 × 10^{5} − 4 × 10^{6} spikes.

The relationship among ⟨σ⟩ and the average connectivity reported in eq. (4.11) clearly indicates that for a sparse network, with constant average connectivity, ⟨σ⟩ will remain finite even in the thermodynamic limit. Therefore, we can conclude that a sparse network cannot be ever reduced to a fully coupled one by simply rescaling the synaptic coupling as done in Ref. [47], even for very large systems. We will analyze sparse networks in more details in the next Chapter.

4.5 Chaotic vs Regular Dynamics

As shown in Chapter 1, homogeneous fully connected pulse coupled networks exhibit regular dynamics [28]. In particular, for excitatory network and finite pulses the asynchronous state becomes a splay state characterized by all neurons spiking one after the other at regular intervals with the same frequency and by a constant mean field Ē, while the PS regime becomes perfectly periodic at a macroscopic level [25].

The introduction of disorder in the network leads to irregularity in the dynamics of the single neurons, which are reflected also at a macroscopic level. This kind of deterministic irregular behavior has been identified as weak chaos whenever the irregularity vanishes for sufficiently large system size [61]. The chaotic motion can be characterized in terms of the maximal Lyapunov exponent λ_{1}: regular orbits have non positive exponents, while chaotic dynamics are associated with λ_{1} > 0.

For finite size networks chaotic dynamics is observed both for annealed and quenched Figure 4.6: (Color online) Maximal Lyapunov exponents as a function of the system size N for various γ-values. The data have been obtained by discarding a transient of the order of 10^{8} − 10^{9} spikes and then by following the dynamics in the real and tangent space for an equivalent duration, moreover the data have been averaged over 3 to 5 different network realization with quenched disorder.

\begin{matrix}γ & δ \\ 1.1 & 0.26 ± 0.01 \\ 1.3 & 0.228 ± 0.007 \\ 1.5 & 0.174 ± 0.006 \\ 1.7 & 0.085 ± 0.009 \end{matrix}

Table 4.1: Power law exponents giving the decay of the maximal Lyapunov λ_{1} ∝ N^{−δ} for the data reported in Fig. 4.6a for quenched disorder.

disorder. However, for all γ-values examined in this work, λ_{1} tends to decrease for a sufficiently large number of neurons in the network. Therefore we can safely affirm that for any Erdös-Renyi network with average connectivity given by Eq. (4.7) the neuronal dynamics is weakly chaotic; i.e. the evolution will become completely regular for infinite networks. Numerical results for quenched disorder are reported in Fig. 4.6a for various γ-values. The maximal Lyapunov exponent exhibits clear power-law decays N^{−δ}, with δ decreasing and eventually vanishing for γ → 2 (see Table 4.5). These results generalize previous indications reported in Chapter 3 for a specific realization of diluted network.

Considering networks with annealed disorder it must be underlined that for sufficiently large number of neurons and γ < 2, the maximal Lyapunov reveals a tendency to decrease. However, clear scaling laws cannot be inferred for all the examined γ-values on the range of affordable system sizes, namely 50 ≤ N ≤ 12, 500. Larger network sizes are probably required with unreliable synapses to have clear scaling laws, due to the fact that the transition from asynchronous regime to partial synchronization occurs at critical system sizes within the investigated range. On the contrary, in the quenched case once γ is fixed an unique regime is observable for almost all the investigated sizes. In particular, the network is always in an asynchronous regime for γ = 1.7 and 1.5, while it is in partial synchronization for γ = 1.1 and 1.3 for (almost) all the considered number of neurons. These additional findings strengthen the above reported conclusions: finite size systems are weakly chaotic for any considered dynamical regime 5.

According to the results reported above, we expect that for sparse networks the maximal Lyapunov exponent will eventually saturate to some constant value, apart for possible logarithmic corrections. Thus, as we will show in the next Chapter, sparse networks remain chaotic for any system size, paralleling the behavior of the microscopic fluctuations shown in the previous Section. Therefore microscopic inhomogeneities and chaotic behavior appear as deeply related.

4.6 Discussion

Collective periodic oscillations in excitatory Erdös-Renyi networks can be observed only above a critical average in-degree. This latter quantity saturates to a constant value for networks with a sufficiently large number of neurons, thus suggesting that the key ingredient responsible for the emergence of collective behaviors is the number of pre-synaptic neurons both for massively connected networks as well as for sparsely connected ones. This result confirms and generalizes previous findings on the stability of the complete synchronized state for pulse-coupled Hindmarsh-Rose neurons [68]. Furthermore, our results indicate that the minimal network size required to observe collective oscillations diverges with the exponent γ ruling the scaling of the average connectivity with the number of neurons.

The presence of annealed disorder in the network (corresponding to unreliable synapses) favors the emergence of coherent activity with respect to the quenched case (associated with reliable synapses), since in the first case the asymptotic average connectivity required to observe collective oscillations is much smaller. This is probably due to the fact that in the annealed situation each neuron is on average subjected to the same train of stimuli, while with quenched disorder the dynamics of each neuron depends heavily on its neighbors. Furthermore, for sufficiently large networks the macroscopic behavior observable with reliable or unreliable synapses becomes identical. This seems to indicate that at the level of population dynamics the reliability or unreliability in the synaptic transmission can be irrelevant.

The average in-degree ⟨k⟩ also controls the fluctuations among different neurons, being the fluctuations proportional to the inverse of the square root of ⟨k⟩. Therefore, for Erdös-Renyi networks with average in-degree proportional to any positive power of N, the fluctuations will vanish in the limit N → ∞, leading to a homogeneous collective behavior analogous to that of fully connected networks. On the contrary, inhomogeneities among neurons will persist at any system size in sparse networks.

5We have carefully verified for various system sizes that the maximal Lyapunov exponents, averaged over a suitable time window, do not vary over time and does not depend on the duration of the chosen transitory.

Recent experimental results on the intact neuronal network of the barrel cortex of anesthetized rats seems to clearly suggest that the dynamics of this system is chaotic [69]. This result poses severe questions about the possibility of reliable neural coding. The evolution of our models is chaotic for any finite networks. However, in presence of coherent periodic activity, the chaoticity present in the system is not so strong to destroy the average collective motion. Thus the trial-to-trial variability, observable in the response of each single neuron and induced by chaos, does not prevent the possibility of the network to encode information. The information can be coded in some property associated to the global (almost regular) oscillations of the network which are robust to local fluctuations emerging at the neuronal level, and this collective coding can represent an alternative to the dilemma rate versus spike timing coding [69].

The level of chaos in the examined networks decreases with N and the dynamics becomes regular in the thermodynamic limit. However, the decrease of the maximal Lyapunov exponent with N slows down dramatically for γ → 2, suggesting an erratic asymptotic behavior for sparse networks. This point will be deeply analyzed in the next Chapter.

The results presented in this Chapter have been published in Ref. [70]

Collective dynamics and extensive chaos in sparse networks

5.1 Introduction

In this Chapter, we study the dynamical properties of “sparse networks” at a macroscopic and microscopic level, by investigating several different local dynamics. As discussed in Section 1.4 we can distinguish between two classes of systems [35]: massive networks, where K is proportional to the network size N; sparse (or strongly diluted) networks, where K ≪ N, and specifically K is independent of N as N → ∞. Regular lattice systems with short-range interactions belong to the latter class. While it is not surprising to observe the onset of a collective motion in massive networks, as shown in the previous chapters, it is less obvious to predict whether and when this can happen in sparse ones. As shown in Chapter 3 in a model of leaky integrate-and-fire (LIF) neurons, a finite connectivity can sustain a partially synchronized regime [71]. Here, we show that the emergence of a collective dynamics above a finite critical connectivity K_{c} is a general and robust property of sparse networks of oscillators. Since K_{c} turns out to be of the order of a few tens in all models we have investigated, macroscopic motion appears to be rather ubiquitous and possibly relevant in the context of neural dynamics. In our simulations, we have typically assumed that all nodes are characterized by the same connectivity K, but we have verified that the same scenario holds assuming a Poissonian degree distribution with average connectivity K, as in Erdös-Renyi graphs.

Finally, we analyze the microscopic dynamics, irrespective of the presence of the macroscopic phase, finding that it is always extensive (the number of unstable directions, as well as the power contained in the principal components, is proportional to the network size). This property is highly nontrivial, as the network dynamics is non additive (it cannot be approximated with the juxtaposition of almost independent substructures, see below). This is at variance with globally coupled systems, which exhibit a non-extensive component in the Lyapunov spectrum [72].

More specifically we study three classes of dynamical systems: (i) units that are chaotic by themselves (logistic maps - LM); (ii) units that may become chaotic as a result of a periodic forcing (Stuart-Landau oscillators - SL); (iii) phase-oscillators that cannot behave chaotically under any forcing (LIF neurons).

5.2 Coupled Maps

The dynamics on a network of N coupled logistic maps (LM) is defined as

x_{n+1}(i) = (1 − g)f(x_{n}(i)) + gh_{n}(i) ,  (5.1)

where x_{n}(i) represents the state of the ith node (i = 1, …, N) at time n, the logistic map f(x) = ax(1 − x) rules the internal dynamics, and g is the coupling strength. Finally,

h_{n}(i) = (1/K) ∑_{j=1}^{N} C_{ij}f(x_{n}(j))

denotes the local field, where C_{ij} is the connectivity matrix: C_{ij} = 1 if an incoming link from j to i is present, otherwise C_{ij} = 0. It is convenient to introduce the average field \overline{h}_{n} and its standard deviation σ_{h} (σ_{h}^{2} = ⟨\overline{h}_{n}^{2}⟩ − ⟨\overline{h}_{n}⟩^{2}) ^{1} .

In Fig. 5.1, σ_{h} is plotted versus the connectivity K for a = 3.9, g = 0.1 and increasing network sizes. For low connectivity, σ_{h} is quite small and decreases as 1/\sqrt{N} with the system size (see left inset), i.e. the deviation from zero is a finite-size effect. Above K_{c} ≃ 60, σ_{h} assumes finite values, independently of the system size, signaling the onset of a collective dynamics. In fact, the right inset in Fig. 5.1 reveals nontrivial collective oscillations for K = 500 and N = 20, 000 (we have verified that the thickness of the “curve” does not decrease upon increasing the network size - data not shown). The phase-portrait is analogous to that previously obtained in globally coupled maps [73]. This indicates that the evolution of a sparse network reduces, for K → ∞, to that of its corresponding mean-field version. What is new and a priori non obvious is that a finite and relatively small connectivity suffices to sustain a macroscopic motion.

As for the evolution of the single units, the most appropriate tool to investigate the microscopic dynamics is Lyapunov analysis. In Fig. 5.2a,b we can see that both below and above K_{c} the dynamics is characterized by extensive high-dimensional chaos [74, 75], since the spectra of the Lyapunov exponents (LE) {λ_{i}} collapse onto one another, when they are plotted versus the intensive variable i/N [76]. In the inset, one can appreciate that the convergence occurs also for the largest LE, at variance with the non-extensive behavior, recently detected in globally coupled networks [72].

In a sparse network, the field h_{n}(i) fluctuates with i, no matter how large the network is, since h_{n}(i) is the sum of a finite number of contributions. One way to characterize its variability is by determining the covariance exponents (CE), i.e. the eigenvalues μ_{i} of the matrix

S_{ij} = ⟨δh_{n}(i)δh_{n}(j)⟩ − ⟨δh_{n}(i)⟩⟨δh_{n}(j)⟩,

1Here and in the following, overlines and angular brackets denote averages over nodes and time, respectively.

Figure 5.1: Model LM. Standard deviation of the mean field, σ_{h} (averaged over 5 realizations of disorder) versus K for N = 5, 000 (red) circles, N = 10, 000 (green) squares, N = 20, 000 (black) triangles, and N = 40, 000 (blue) diamonds. The upper inset shows the rescaled σ_{h} at low K. In the lower inset the return map of \overline{h}_{n} for K = 500 and N = 20, 000. Symbols not connected by lines refer to K = 2N_{s} + k_{R} (N_{s} = 40 is the number of nearest neighbours and k_{R} of the random links) for N = 5, 000 (violet) crosses, N = 10, 000 (magenta) stars and N = 20, 000 (red) plus.

where δh_{n}(i) = h_{n}(i) − \overline{h}_{n}. In 1D spatial systems with periodic boundary conditions such an approach would correspond to determine the spatial Fourier spectrum. In this case, since there is no “wavelength” to refer to, it is natural to order the eigenvalues from the largest to the smallest one. The results for different network sizes are plotted in Fig. 5.2 versus i/N, (panel c and d refer to K values below and above K_{c}, respectively). The good data collapse confirms the extensivity of microscopic fluctuations. In both panels, μ_{i} has been rescaled by K, to emphasize extensivity in yet a different way; in fact, as the local field is the sum of K contributions, its variance is expected to be on the order of 1/K.

5.3 Stuart-Landau Oscillators

The second model we have analysed is a network of Stuart-Landau oscillators (SL),

ẇ_{i} = w_{i} − (1 + ic_{2})|w_{i}|^{2}w_{i} + g(1 + ic_{1})(W_{i} − w_{i}) ,  (5.2)

where w_{i} is a complex variable, and W_{i} = (1/K) ∑_{j=1}^{N} C_{ij}w_{j} represents the local field (C_{ij} is defined as before). Since the local variable is a complex number, it is convenient to introduce the global mean field W (t) ≡ |\overline{W} | (where | · | denotes the modulus operation), which essentially coincides with the Kuramoto order parameter (1.8) [31]. The data reported for ⟨W ⟩ in Fig. 5.3a reveals the discontinuous emergence of some form of synchronization (at least for our choice of the parameter values, c_{1} = −2, c_{2} = 3 and Figure 5.2: Model LM. LE spectra are reported in (a) for K = 10 and N = 100−200− 500 − 1, 000, and in (b) for K = 80 and N = 200 − 500 − 1000. CE spectra (rescaled by K) are reported in (c) for K = 10 and N = 200 − 400 − 800 − 1, 600, and in (d) for K = 100 and N = 800 − 1, 600 − 3, 200. In the insets a zoom of the largest values is shown. In this and in the following figures the arrow direction indicates data obtained for increasing system sizes.

g = 0.47^{2}). More precisely, there exists a finite parameter range (55 < K < 85), where, depending on the initial conditions ^{3}, ⟨W ⟩ may either vanish or take a finite value. A more precise characterization of the collective phase can be obtained by looking at the variance of the order parameter, namely σ_{w}^{2} = ⟨W^{2}⟩ − ⟨W ⟩^{2}. The data in Fig. 5.3b reveals that the discontinuous transition is accompanied by the birth of temporal fluctuations which increase with the connectivity. As shown in the inset of Fig. 5.3b, the global attractor exhibits an irregular dynamics. Moreover, it is remarkable that the attractor is qualitatively different from the one found in the mean-field version of the model [77], which is centered around W ≃ 0.35.

As for LM, we have verified that both below and above the transition region, the microscopic dynamics is chaotic. In both cases, there is a clear evidence of a convergence towards an asymptotic LE spectrum (the LE spectra for K = 10 are reported in Fig. 5.3c). Analogous conclusions can be drawn from the CE spectra (see Fig. 5.3d).

2Such values have been selected from [77], where the globally coupled version of the model was thoroughly investigated.

3The hysteresis has been detected by simulating a network where links are first slowly added and then removed.

Figure 5.3: Model SL. Average (a) and standard deviation (b) of W versus K for N = 10, 000 (red) squares, and N = 20, 000 (blue) triangles (averaged over 5 realizations of the disorder). The dashed lines in (a) show the region of bistability (see the text). The inset in (b) shows the macroscopic attractor for N = 10, 000 and K = 800. (c) LE spectra for K = 10 and N = 100 − 200 − 400 . (d) CE spectra (rescaled by K) for K = 10 and N = 40 − 80.

5.4 Leaky Integrate and Fire Neurons

Finally, we have considered LIF pulse-coupled neurons. They are among the most popular and yet simple models used in computational neuroscience, the field where understanding the onset of collective motion is likely to have the deepest impact. The evolution equation for the membrane potential u_{i}, writes as Eq. (3.1), where the local field E_{i} satisfies the equation (4.1) with k_{i} ≡ K , ∀i. Whenever the membrane potential u_{i} reaches the threshold u_{i} = 1, it is instantaneously reset to the value u_{i} = 0 and a so-called α-pulse is emitted towards the connected neurons (for more details see Chapter 2 [25]). In this case, we introduce the mean fields \overline{E}, \overline{P} as in Eq. (3.14) and the corresponding standard deviation σ_{E}. The model has been simulated for a = 1.3, excitatory synaptic strength g = 0.2 and inverse pulse-width α = 9. For such values, it is known that in the global coupling limit, there exists a periodic collective motion accompanied by a quasi-periodic microscopic dynamics [28]. For small K, the mean field is constant in the thermodynamic limit, revealing a so-called asynchronous state, while above a critical value K_{c} ≈ 9, it oscillates periodically, as seen in the inset of Fig. 5.4a. As in the previous systems, the Lyapunov analysis reveals an extensive behaviour, including the initial part (see the two insets), where at finitesize corrections of order 1/N are detected. This is to be contrasted with the initial non extensive “layer” observed in globally coupled systems [72]. Fully extensive behaviour was already observed for the Θ-neuron model in [78]. Finally, extensivity is confirmed by the orthogonal decomposition applied to the fluctuations of the local field E_{i} (see the CE spectra shown in Fig. 5.4c).

Figure 5.4: Model LIF. (a) Standard deviation of the mean field, σ_{E}, versus K for N = 1, 000 (black) circles, N = 5, 000 (red) squares, N = 10, 000 (green) triangles. The inset shows the macroscopic attractor for N = 5, 000 and K = 200 (\overline{P} = ˙ \overline{E} + α\overline{E}).

(b) LE spectra (in the lower inset a zoom of the largest values) for K = 20 and N = 240 − 480 − 960. In the upper inset the maximum Lyapunov exponent, λ_{1}, versus N is shown, the (red) line represents the nonlinear fit λ_{1} = 0.0894 − 2.3562/N and the (green) dashed line marks the asymptotic value. (c) CE spectra (rescaled by K) for K = 20 and N = 200 − 400 − 800 − 1600.

5.5 Discussion

By studying three different models, we find that a finite connectivity is able to sustain a nontrivial collective motion, as signalled by the appearance of finite temporal fluctuations of the mean field. This scenario emerges in LMs, whose entire Lyapunov spectrum is positive, as well as in LIF neurons that cannot behave chaotically as standalone devices, even when subject to an irregular forcing. The differences among the various models concern the nature of the transition (continuous in LMs and LIF neurons, discontinuous and hysteretic in the SL oscillators) and the complexity of the collective phase that is periodic in LIF neurons but certainly higher dimensional in the other two models. It is desirable to trace back analogies and differences to some general properties of the single systems. In this perspective, it is natural to assume that a finite connectivity acts as a mean field accompanied by an effective noise (due to the “statistical” differences among the fields seen by the different nodes). Accordingly, one expects that a small K corresponds to a large noise and is thereby unable to maintain a global coherence, as indeed observed. This picture is, however, rather qualitative, since the collective motion of SL oscillators markedly differs from that one generated in the mean-field limit, even when K ≈ 800. Altogether, it is remarkable that a few tens of connections consistently suffice to sustain a collective motion in such different environments.

In order to shed further light on the phenomenon, we have added a spatial structure, by organizing the nodes along a line and adding finite-distance interactions. When only the latter are present, no collective motion can arise, because of the low dimensionality of the lattice [79]. The mixed case corresponds to a “small-world” arrangement [80]; we have imposed that the connectivity is the sum of two contributions (K = 2N_{s} +k_{R}): a number k_{R} of random links and 2N_{S} links with the nearest nodes. The data reported in Fig. 5.1 for N_{s} = 40 shows that it is sufficient to add k_{R} = 20 links to establish a collective dynamics; in other words, in the presence of a lattice structure, a lower number of long-range connections may be necessary, although the overall connectivity is larger. This observation reveals that the network structure plays a nontrivial role in determining the number of links that can sustain macroscopic motion.

Next, we comment on the extensive nature of the microscopic motion, a property that is much less obvious than one could think. In fact, the existence of a limit Lyapunov spectrum in regular lattices is the natural consequence of the additivity of the dynamics [74]. Imagine to use a hyperplane P to divide a spatial domain of size N into two subdomains S_{1} and S_{2} of size N/2. The mutual interaction between S_{1} and S_{2} is negligible as it affects only the interfacial region around P (it is a “surface” effect). As a result, the entire system can be approximately seen as the juxtaposition of two independent subsets. In the case of sparse networks, it is not even obvious how to split them in two minimally-connected components S_{1}, S_{2}. This problem goes under the name of graph bipartitioning; it is known to be NP complete, and, more important, the solution involves a number of connections that is proportional to N itself, whenever K > 2 log 2 [81, 82]. Therefore, the “interface” cannot be likened to a “surface” and the evolution is necessarily non-additive. Accordingly, the extensive nature of the microscopic evolution is due to subtle properties, that have yet to be clarified.

The results presented in this Chapter have been published in Ref. [83].

Collective chaos and chimera states in a pulse-coupled network of networks

6.1 Introduction

A slightly more complex but meaningful setup is that of two symmetrically coupled populations of identical oscillators. This is the simplest instance of “network-of-networks” that is often invoked as a paradigm for neural systems. With reference to sinusoidal oscillators, this setup has revealed the onset of chimera states (one of the two populations is fully synchronized, while the oscillators of the other one are not synchronized at all [84]), as well as more complex macroscopic states with periodic [85] and quasi-periodic [86] collective oscillations. In the present Chapter we study the two-population setup with reference to LIF neurons characterized by an excitatory coupling. Upon varying the coupling strengths between and within the two populations we find symmetrybreaking transitions that lead to the onset of various chimera states (similar to those observed in [85, 86]) as well as to a new regime, where the two populations are both partially synchronized, but with a different degree. Symmetric collective states of increasing dynamical complexity are also observed. Very interesting is the parameter region where the collective motion is chaotic, as indicated by the finite-amplitude Lyapunov exponent (FALE) [87] and confirmed by the computation of the standard Lyapunov spectrum which reveals the existence of several positive exponents, indicating that the microscopic dynamics is high-dimensional.

6.2 The Model

We consider two fully coupled networks, each made of N LIF oscillators. Following 3.1, the membrane potential u_{j}^{(k)}(t) of the j − th oscillator (j = 1, …, N) of the kth population (k = 0, 1) evolves according to the differential equation

u̇_{j}^{(k)}(t) = a − u_{j}^{(k)}(t) + g_{s}E^{(k)}(t) + g_{c}E^{(1−k)}(t)  (6.1)

where a > 1 is the suprathreshold input current, while g_{s} > 0 and g_{c} > 0 gauge the self- and, resp., cross-coupling strength of the excitatory interaction. Whenever the membrane potential reaches the threshold u_{j}^{(k)} = 1, it is reset to u_{j}^{(k)} = 0, while a so-called α−pulse is sent and instantaneously received by all the neurons. The field E^{(k)}(t) represents the linear superposition of the pulses emitted within the population k in the past. It can be shown in Section 3.2 that E^{(k)}(t) satisfies the differential equation

Ë^{(k)}(t) + 2αĖ^{(k)}(t) + α^{2}E^{(k)}(t) = \frac{α^{2}}{N} ∑_{j,n} δ(t − t_{j,n}^{(k)}) ,  (6.2)

where t_{j,n}^{(k)} is the nth spiking time of the jth neuron within the population k, and the sum is restricted to times smaller than t. In the limit case g_{s} = g_{c} = g, the two populations can be seen as a single one made of 2N neurons with an effective coupling constant G = 2g.

The degree of synchronization within each population of neurons can be quantified by introducing the typical order parameter (1.8) used for phase oscillators, generalized to a two population setup

R^{(k)}(t) = |⟨e^{iθ_{j}^{(k)}(t)}⟩|,  (6.3)

where θ_{j}^{(k)} is the phase of the jth oscillator, that can be properly defined as a (suitably scaled) time variable [67], θ_{j}^{(k)}(t) = 2π(t − t_{j,n}^{(k)})/T_{q,n}, where t_{j,n} indicates the time of the last spike emitted by the j^{th} neuron, while T_{q,n} = t_{q,n+1} −t_{q,n}^{(k)} is the n-th inter-spike interval associated to the neuron q which was the last to fire in the population k. One can verify that this phase is bounded between 0 and 2π, as it should. It is interesting to see that the application of this definition to the PS regime described by van Vreeswijk [28] reveals that the order parameter fluctuates periodically with the same period as the macroscopic field. In other words PS differs from the regime observed in the Kuramoto model above the synchronization threshold, where the order parameter is constant in time [31].

6.3 Phase Diagram

The equations have been integrated by extending the event-driven approach described e.g. in Section 3.2. In practice, the (linear) equations of motion are solved analytically in between two consecutive spike-emissions, obtaining a suitable map. Since the ordering of the single-neuron potentials does not change within each population, the next firing event can be easily determined by comparing the neurons that are closest to threshold within each of two populations. In spite of the conceptual simplicity and the effectiveness of the code, one must be nevertheless careful in handling nearly Figure 6.1: Phase diagram in the (g_{c}, g_{s})-plane reporting the stability region of the observed various collective solutions. For the definitions of the different phases see the text.

singular cases, when many neurons almost cluster together. In order to avoid the spurious clustering, due to numerical roundoff, we have changed variables, introducing and monitoring the logarithm of the difference of the membrane potentials of two successive neurons. This requires some care in defining the right number of variables: given N potentials, one naturally has N − 1 differences that have to be complemented by a proper Nth variable (for more details see Appendix A.2).

The phase plane (g_{c}, g_{s}) shown in Fig. 6.1 has been obtained by studying the model (6.1,6.2) for a = 1.3 and α = 9. The diagram is semiquantitative in the sense that a much more detailed work would be needed to identify exactly the stability borders of the different regimes. Along the diagonal (g = g_{s} = g_{c}) the model reduces to that for a single population with coupling strength G = 2g. For our choice of a and α values, the system exhibits PS, since we are below the critical value G_{0} = 0.425 [33] above which the splay state is stable (the splay state is a regime characterized by a constant spiking rate and thereby a constant field, as discussed in Chapter 1). Below the diagonal, the evolution is still symmetric but fully synchronized (FS), i.e. all neurons of both populations fire together. More intriguing is the region above the diagonal, that is characterized by a spontaneous symmetry breaking: one population fully synchronizes, while the other is in a PS regime, i.e. we are in presence of a generalized chimera state, as shown in Fig. 6.3(e). This can be appreciated by looking at the synchronization parameter R^{(k)} of the two populations, one of which is equal to one, while the other oscillates periodically close to 0.8 (see Fig. 6.3(g)). By following [85], this state can be classified as a periodically breathing chimera. In this regime, the two populations are characterized by a microscopically periodic and quasi-periodic behaviour, respectively. In spite of this qualitative difference, the two (macroscopic) fields E^{(0)} and E^{(1)} are both periodic and phase locked (see Fig. 6.2(a)). This means that the neurons subject Figure 6.2: The macroscopic attractors displayed by reporting the fields E^{(0)} vs E^{(1)} for four different non chaotic phases, namely (a) PS+FS (g_{c} = 0.07, g_{s} = 0.1), (b) PS1-PS2 (g_{c} = 0.02, g_{s} = 0.17), (c) APS (g_{c} = 0.07, g_{s} = 0.35) and (d) TORUS (g_{c} = 0.07, g_{s} = 0.3) (for the exact definitions see the text). The red curves reported in the panel (d) are the Poincaré sections obtained by imposing that the sum E^{(0)} + E^{(1)} is maximal.

to two different linear combinations of E^{(0)} and E^{(1)} behave differently: a population locks with the forcing field, while the other one behaves quasi-periodically. Another even more interesting symmetry broken state can be observed for larger g_{s}-values and g_{c} < 0.055; in this case both populations exhibit PS, but their dynamics take place over two different attractors with two different degrees of synchronization (PS1-PS2 regime), as shown in Fig. 6.3(a). Like in the PS-FS regime, the two fields behave periodically (with the same period) and are phase locked, as it can be appreciated by looking at the closed curve E^{(0)} versus E^{(1)} in Fig. 6.2(b). However, at variance with PS-FS, here both populations exhibit quasi-periodic motions. In other words we are in presence of a different symmetry breaking, where two populations with distinct quasiperiodic motion spontaneously emerges. Pikovsky and Rosenblum reported recently the existence of a quasi-periodic breathing chimera where one population was FS while the other oscillates quasi-periodically [86]. However, in our case both populations undergo quasi-periodic motions. For yet larger g_{s} values the equivalence between the collective dynamics of the two population is restored with the only difference of phase shift between the two fields which oscillate in antiphase and this is why we term this regime Antiphase Partial Synchronization (APS). In the APS phase, for finite N, the instantaneous maximum Lyapunov exponent strongly fluctuates and we cannot rule out the possible existence of some form of weak chaos, analogous to the one discussed in Chapter 3 [47] for a model of diluted neural network, i.e. a chaotic behaviour that Figure 6.3: Macroscopic attractors displayed by reporting P ≡ E + αĖ vs E for a PS1-PS2 state (a), a chaotic phase (b), a PS-FS state (e) and a collective quasiperiodic motion (f). The time evolution of the corresponding order parameters R^{(0)} and R^{(1)} is also reported in (c), (d), (g), (h). The variables corresponding to population 0 (resp.

1) are shown in black (resp. red) in Figs (a)-(d). The variables corresponding to population 0 (resp. 1) are shown in blue (resp. maroon) in Figs (e)-(h). As regards the parameter values, (g_{c} = 0.02, g_{s} = 0.17) in (a),(c), (g_{c} = 0.08, g_{s} = 0.16) in (b),(d), (g_{c} = 0.08, g_{s} = 0.1) in (e),(g), (g_{c} = 0.04, g_{s} = 0.3) in (f),(h).

disappears in the thermodynamic limit. In a strip above the chaotic region (discussed below), one can observe collective quasiperiodic motion (see Fig. 6.3(f),(h)). This means that the quasiperiodic motion of the fields is accompanied by a dynamics of the single neurons along a torus T^{3}. An analogous regime has been previously reported in [88] in the context of a population of coupled Stuart-Landau oscillators. Here, we find it in a model where the single units are described by an single variable. Furthermore, we have characterized the motion on the macroscopic T^{2} attractor reported in Fig. 6.2, by estimating the winding numbers for various values of the coupling strengths. We find that the winding number is typically quite small – on the order of ∼ 10^{−2} – but independent of the system size, indicating that the torus survives in the thermodynamic limit. Finally, for yet larger g_{s}-values both populations converge towards a splay state. This is not surprising, as we already know that for the chosen α- and a-values, the splay state is stable in a single population of neurons for G > G_{0} ≡ 0.425.

6.4 Collective Chaos

In a limited region above the diagonal and for g_{c} > 0.055 the collective behaviour is irregular, this is clearly seen by observing the two macroscopic attractors, corresponding to the two populations, reported in Fig. 6.3b. Nonetheless, from Fig. 6.3d it appears Figure 6.4: Finite amplitude Lyapunov exponents λ_{F} versus the logarithm of the perturbation amplitude Δ for three different system sizes: namely, N = 400 (black circles), N = 800 (green squares) and N = 1, 600 (blue triangles). The amplitudes Δ have been estimated by considering the euclidean distance among the perturbed and unperturbed fields. The dashed (magenta) line indicates the maximal microscopic Lyapunov exponent λ_{1} for N = 1, 600 obtained by following the orbit over a time span containing 10^{8} spikes and after discarding a transient composed by 10^{6} spikes. λ_{F} has been estimated by averaging over 25, 000 − 50, 000 different trajectories. The results are obtained for (g_{c} = 0.08, g_{s} = 0.16).

that the associated order parameters behave almost periodically, but this periodicity is only apparent since it reflects only the larger time scale present in the system, the one associated to the modulation of the fields E^{(k)}, while the irregularity of the dynamics are more evident at smaller time scales. In particular, at least two other time scales are present: a scale O(1) associated to the firing period of each specific neuron and a scale O(1/N) corresponding to the interspike interval between two successive spike emissions in the network. In order to make a quantitative assessment of the chaoticity we have first studied the FALE λ_{F}. The FALE can be determined from the growth rate of a small finite perturbation at different amplitudes Δ of the perturbation itself (after averaging over different trajectories) [87]. This is done by randomly perturbing the coordinates (both fields and the membrane potentials of the two populations) of a generic configuration on the attractor. The results for g_{s} = 0.16, g_{c} = 0.08 and three different system sizes are plotted in Fig. 6.4. For small Δ values, λ_{F} grows with Δ, since the perturbation needs first to converge towards the most expanding direction, while the final drop is a manifestation of the saturation of the perturbation amplitude. It is the height of the intermediate plateau which measures the amplitude of the FALE. Since the height is independent of N, one can conjecture that the collective motion is chaotic and stays chaotic in the thermodynamic limit. Besides the data reported in Fig. 6.4, we have checked that the plateau height is not influenced by changing the amplitude of the initial perturbation and by performing tests up to N=6,400.

Figure 6.5: Lyapunov spectra λ_{i} versus i/N for three N-values: namely, N = 50 (filled black circles), N = 100 (empty red squares) and N = 200 (empty green triangles). In the inset the first three Lyapunov exponents are reported as a function of N: λ_{1} (blue circles), λ_{2} (turquoise squares) and λ_{3} (magenta triangles). The Lyapunov esponents have been obtained by following the dynamics in the real and tangent space for a time span containing 10^{8} − 5 × 10^{9} spikes, after discarding a transient period of 10^{6} spikes. The results are obtained for (g_{c} = 0.08, g_{s} = 0.16).

It is instructive to compare λ_{F} with the maximum λ_{1} of the standard Lyapunov spectrum. For small Δ values, the two indicators should coincide, but there is no reason for the agreement to persist at larger amplitudes. In fact, a second plateau has been detected in the context of globally coupled maps. The plateau occurring for larger Δ’s has been interpreted as an indication that the chaoticity of collective variables differs from that of the microscopic ones. Since, in the present case, we observe only a single plateau, it is crucial to verify its compatibility with λ_{1}. The horizontal dashed line in Fig. 6.4 corresponds to the Lyapunov exponent determined for N=1600 (for the dependence of λ_{1} on N, see below). The two indicators are consistent with each other and this means that collective variables are as chaotic as the microscopic ones.

Having established the existence of one unstable direction, the next question is to determine how many such directions are present. Unfortunately, the concept of FALE allows to determine just one exponent. As a consequence, we turn our attention to the usual Lyapunov spectrum, well aware that the “microscopic” exponents do not necessarily reproduce the chaoticity of the “macroscopic” variables.

In the present model, the Lyapunov spectrum {λ_{i}} is composed of 2N +3 exponents (i = 1, …, 2N +3). The phase space dimension is in fact equal to 2N +4 (N potentials plus 2 two-dimensional equations for the fields), but the direction corresponding to the zero-Lyapunov exponent is implicitly eliminated as a consequence of taking the Poincaré section [33], as already explained in Section 3.2.1. In Fig. 6.5 we have plotted the first part of the spectrum (with the exception of λ_{1}, for the clarity of presentation) with the usual normalization i/N of the x axis. The figure clearly indicates that the spectrum becomes increasingly flat and converges to zero. This can be understood in the following way. In the thermodynamic limit, the dynamics of globally coupled identical oscillators can be viewed as that of single oscillators forced by the same (selfconsistent) field [88]. As a result, in a first approximation, we expect all Lyapunov exponents to be equal to the conditional Lypunov exponent λ_{c} obtained by forcing a single LIF neuron (identical to the others) with the self-consistent fields (obtained by integrated the whole ensemble). We have found λ_{c} = 0 (within numerical accuracy). This justifies why the increasingly flat Lyapunov spectrum converges towards zero. Moreover, since a LIF neuron is described by a single variable, under no circumstance λ_{c} can be positive. On the other hand, λ_{c} can be negative, and we expect this to occur whenever the given population exhibits full synchronization. Since we know that no synchronization is observed in the chaotic regime, we can conclude that λ_{c} is not just very small, but it must be exactly equal to zero.

Having established that the Lyapunov spectra converge to zero, it is interesting to investigate its scaling behaviour. By comparing the spectra obtained for different system sizes, it turns out that they scale as 1/N^{β}, although we are unable to extrapolate the value of β from the analysis of the relatively small systems that we have simulated. By comparing the spectra obtained for N = 100 and 200 we can at most guess that β ≈ 1.5. This value is not far from β = 2, found analytically while studying the splaystate stability in single populations of LIF neurons, as shown in Chapter 2 [46] and numerically for the PS state, as shown in Chapter 3 [47].

Finally, we turn our attention to the first part of the spectrum, where the flatness hypothesis does not hold. More precisely, we investigate the N-dependence of the first three Lyapunov exponents which can be computed for larger lattices (up to N=1600). The maximal Lyapunov exponent appears to converge to a finite asymptotic value λ_{1} = 0.0195(3). On the other hand, the second and third exponents grow systematically with N, both becoming positive for N > 200, with no clear evidence of an eventual saturation. These results suggest that the microscopic chaos is highdimensional (there is no reason to believe that the number of positive exponents is just equal to three). However, we cannot tell whether the number of positive exponents is extensive (proportional to N) or sub-extensive.

6.5 Discussion

We have studied two symmetrically coupled populations of leaky integrate-and-fire neurons for different values of the self- (g_{s}) and cross-(g_{c}) coupling-strength. Some of the collective phenomena that we have identified are quite similar to those observed in the two-population setup in a Kuramoto-like model^{1}. This is not surprising, since it is known that, in the small coupling limit, an ensemble of LIF neurons behaves like the Kuramoto model [89], the only difference being that the coupling function is not purely sinusoidal. The onset of PS suggests that the correspondence can be extended to larger coupling strengths. However the observation that, in the Kuramoto setup, PS can be obtained only by invoking a more general kind of coupling [90] indicates a more complicate relationship between the two classes of models. In fact, here we have found new regimes, such as that where the two population follow different partially synchronized dynamics. A yet more intriguing phenomenon is the collective chaos that we have recognized as such from the computation of the finite-amplitude Lyapunov exponent. Altogether, a general question stands: To what extent pulse-coupled oscillators are equivalent to Kuramoto-like models? Studying the mutual relationships could lead to a relevant progress and allow for a cross fertilization in two distinct communities.

1By Kuramoto-like we mean sinusoidal oscillators with some sort of continuous-time coupling that depends on the phase-differences.

Another still open question is that of the degree of chaoticity of the collective dynamics. The fact that the finite-amplitude Lyapunov exponent coincides with the maximum of the standard Lyapunov spectrum suggests that the latter one is richer and possibly contains all elements necessary for a detailed description of the collective properties too. This was not a priori obvious. The computation of the FALE λ_{F} in globally coupled maps [91] reveals a different scenario, where λ_{F} differs from the maximum Lyapunov exponent. On the other hand, the most recent study of globally coupled Stuart-Landau oscillators [77] provides instead an example where, like here, the maximum exponent provides information on the collective motion.

The results reported in this Chapter have been published in Ref. [92].

Conclusion

In this dissertation we used simple mathematical models to describe neural network dynamics and to analyze the emergence of collective behaviour in these networks. In particular we investigated the role played by the topology in promoting coherent activity in excitatory pulse-coupled networks. The presence of collective oscillations assume greater importance in neuroscience since coherent oscillatory activities are prominent in the cortex of the awake brain during attention and are involved in higher level processes, such as sensory binding, storage of memories, and even consciousness. The origin of these oscillations has been commonly associated with the balance between excitation and inhibition in the network [9], while purely excitatory circuits are believed to lead only to fully synchronized regimes, therefore to a quite trivial macroscopic dynamics. Furthermore in neuroscience a completely synchronized firing is a symptom of neural disorders, like e.g. Parkinson disease. However, recent “ex vivo” measurements performed on the rodent neocortex [52] and hippocampus [63] in the early stage of brain maturation reveal coherent activity patterns, despite the fact that the synapses are completely excitatory at this stage while inhibitory synapses develop only later.

For the above reasons, numerical and analytical studies of collective motions in networks of simple spiking neurons have been mainly devoted to balanced excitatoryinhibitory configurations [40], while few studies focused on the emergence of coherent activity in purely excitatory networks. In particular pioneristic studies on coherent periodic activity in fully coupled excitatory networks of leaky integrate-and-fire neurons revealed a regime characterized by a partial synchronization at the population level, while the single neurons perform quasi-periodic motions [28]. Since in the neural cortex the presence of noise is an intrinsic property and real neural circuits are not expected to have a full connectivity, it is important to investigate whether the partially synchronized regime will survive to strong level of dilution and to the introduction of different topologies. On the other hand it has been already shown that the partially synchronized regime is quite robust to moderate levels of noise [29].

Therefore we analyzed the influence on the dynamics of these models of different topologies. In particular we considered an excitatory random network where neurons are connected as in a directed Erdös- Renyi graph with average connectivity scaling sublinearly with the number of neurons in the network [47, 70]. In these “massively connected” networks the dynamics of coherent collective states is similar to the one found in globally coupled networks. However, due to the disorder present in the system, for finite number of neurons we have inhomogeneities in the neuronal behaviors, inducing a weak form of chaos, which vanishes in the thermodynamic limit. In this limit, the disordered systems exhibit regular (non chaotic) dynamics thus recovering the properties of a homogeneous fully connected network. Similar results will be expected in different dynamical systems, once considered a massively coupled network.

On the other hand, the situation is completely different for a “sparse network” characterized by a constant connectivity, independent on the size of the network. In fact, by increasing the connectivity value, we found that a few tens of random connections are sufficient to sustain a nontrivial (and possibly irregular) collective dynamics [83]. In other words, collective motion is a rather generic and robust property and does not require an extremely high connectivity to be sustained. Moreover, the collective motion coexists with a microscopically chaotic dynamics that does not vanish in the thermodynamic limit and turns out to be extensive (the number of unstable directions is proportional to the network size). More specifically, we studied various classes of dynamical models on random sparse networks and in all cases, irrespective of the presence of the macroscopic phase, we found that the chaotic dynamics is always extensive. Extensive chaos has been already found in spatially extended system with nearest-neighbour coupling (diffusive coupling) induced by the additivity of the system. In our case this property is highly nontrivial, as the network dynamics is non additive and it cannot be approximated as the juxtaposition of almost independent sub-structures. This is one of the fundamental results found in this dissertation.

Furthermore we studied the dynamics of two symmetrically coupled populations of leaky integrate- and-fire neurons characterized by excitatory coupling: this is the simplest instance of “network-of- networks”, that is often invoked as a paradigm for neural system [92]. Even though the system is fully coupled, an inhomogeneity is present due to the fact that self-coupling among the neurons in the same population is different from the cross-coupling among the neurons in the other population. Upon varying both the self-coupling and the cross-coupling strengths we found various kinds of symmetric and symmetry broken collective states some of which are similar to states previously observed [85, 86], while other are completely new. In addition to this, since single neurons are symmetrically coupled, the symmetry breaking is completely spontaneous. This setup has revealed, for the first time in neural networks, the onset of “chimera states”, where one of the two populations is fully synchronized, while the oscillators of the other one are not synchronized at all. So far this kind of symmetry broken states has been observed only in oscillators’ networks [84, 86], while the existence of these states in neural networks have been just conjectured. Furthermore we have found macroscopic states of increasing complexity, characterized by collective periodic oscillations, quasi-periodic motion or even macroscopic chaos. Collective chaos has been observed so far in coupled map lattices or in oscillator networks but never in neural networks. The collective origin of the chaotic motion is revealed by the Finite Size Lyapunov Exponent analysis and confirmed by comparing the previous results with the standard Lyapunov analysis.

In addition to the numerical analysis previously described, we carried out a broad and prime analytical study about the stability of a largely investigated collective state, the so called “splay state” (or asynchronous state), in networks of N globally pulsecoupled phase-like models of neurons [34]. In a splay state, all neurons follow the same periodic dynamics and their phases are evenly shifted: they fire periodically with period T and two successive spike emissions occur at regular intervals T/N. Our model neurons are characterized by a membrane potential u that is continuously driven by the velocity field F (u), from the resetting value u = 0 towards the threshold u = 1. As threshold and resetting value can be identified with one another and thereby u interpreted as a phase, it will be customary to refer to the case F (1) ≠ F (0) as to that of a discontinuous velocity field. For the simple choice F (u) = a−u, the model reduces to the well known case of LIF neurons. In this case we considered both excitatory and inhibitory networks. The study of the stability requires determining the Floquet spectrum, i.e. the complex eigenvalues associated to a given periodic orbit of period T. By developing a perturbative approach that is valid for arbitrary coupling strength we were also able to determine the spectral shape and find it to be independent of the structure of the velocity field. In the large-N limit it is natural to consider 1/N as a smallness parameter and thereby to expand the evolution equations in powers of 1/N since two successive spikes are emitted in an interval T/N. Furthermore we proved that, in the thermodynamic limit, the Floquet spectrum scales as (1/N^{2}) for generic discontinuous velocity fields. For continuous fields, it has been numerically observed that the scaling of the spectrum is at least O(1/N^{4}) [46]. In other words the shape of the spectrum is universal, apart from a multiplicative factor that vanishes if and only if F (1) = F (0), i.e. for true phase rotators where u = 0 coincides with u = 1. The stability of the splay state can be inferred, for arbitrary coupling strength, from the sign of F (1) − F (0): in excitatory (inhibitory) networks, the state is stable whenever F (0) > F (1) (F (0) < F (1)). Finally, another obtained analytical result regards the dependence of the period T onto the size N which is of order o(1/N^{3}).

The results reported in this dissertation have lead to the publication of 5 research articles, namely Refs. [34, 47, 70, 83, 92].

Appendix A

Numerical algorithms

A.1 Simulation of Single Networks

The numerical simulation of the event-driven map (3.4,3.5,3.6) (resp. 4.3) discussed in Chapter 3 (resp. Chapter 4) requires an enormous effort in terms of time and storage memory. In order to improve the code performances and reduce the simulation time by a factor 6 − 7 we developed a new algorithm which consists in the following steps

1. We analyze the membrane potentials at a certain time and we order the potentials in a descending order: the first neuron to reach threshold will be the one with the highest potential, while the furthest to threshold is the last neuron in the queue.

2. We resolve equation (3.7) determining the τ_{i} putative firing times for the first M = N/10 neurons in the queue; the search for the smallest firing time τ is limited to this pool of neurons.

3. The event driven map of the entire network is evolved for I = M/4 spike events. At each event all neuron variables are evolved while the computation of the smallest firing time τ is carried out only among the pool of the first M neurons.

4. After I iterations all membrane potentials are reordered and labelled as in step 1.

5. The simulation is iterated by repeating the above 4 steps.

By calculating the minimum τ among only M neurons we could, in principle, overlook some relevant spike events and analyze false ones. Then it becomes important to avoid overtakings between neurons outside the pool and the ones in the M pool. In order to overcome this problem we refined the algorithm constraining it to work on a constant M-dimensional subset of neurons. Every time a neuron emits a spike, it is reset to zero; therefore the last to spike neuron is removed from the M pool and put at the end of the queue. This means that, iterating the above algorithm, we check a neuron less, among the pool of the first M neurons, each time the event driven map is evolved. In other words, if the map is evolved for I = M/4 spike events, as described in step 3, the computation of the smallest to reach threshold neuron τ is finally carried out among (M − I) neurons only. As soon as we consider highly diluted networks, this approximation turns out to be extremely poor and we should refine the algorithm in a way that it will work on a constant M-dimensional subset of neurons. Therefore each time a neuron has emitted a spike, being thus removed from the M pool and put at the end of the queue, we insert another neuron in the pool. The neuron we insert, is the first neuron out of the pool, i. e., the (M +k)th ordered neuron in the queue at the kth iteration. In this way we keep the number of checked neurons constant, thus ensuring the possibility to determine the shortest time to reach threshold τ always between the closest to threshold pool of neurons.

A.2 Simulation of Network of Networks

In the two symmetrically coupled populations setup discussed in Chapter 6 we find various kind of symmetric and symmetry broken states; in particular we find synchronized states. In such states, due to the fact that the membrane potentials lie in a very limited range of values, it is numerically difficult to determine the neuron nearest to threshold. Thus, to increase the numerical accuracy and to avoid spurious clustering due to numerical round off, we decided to modify variables. Instead of integrating the membrane potentials, we performed the integration of the logarithm of the difference of two successive neurons. This transformation is uniquely defined in globally coupled systems, since the order of the membrane potentials of the neurons is preserved in time. Therefore, it is possible to define an ordered list of the potentials and, on this basis, to define uniquely the “neighbours” of a neuron. Given a set of N membrane potentials {u_{j}^{(k)}}_{j=1,...,N}, with k = 0, 1 depending on the considered family, we introduce at a

generic time t the following N + 1 auxiliary variables:

ω_{1}^{(k)}(t) = ln [1 − u_{1}^{(k)}(t)]
ω_{j}^{(k)}(t) = ln [u_{j−1}^{(k)}(t) − u_{j}^{(k)}(t)]  j = 2, …, N
ω_{N+1}^{(k)}(t) = ln [u_{N}^{(k)}(t)]

where the threshold (resp. reset) value is u_{th} = 1 (resp. u_{R} = 0) and u_{1}^{(k)} is the next to threshold neuron.

Since we would like to define an event driven map for the two coupled families, it is necessary to find which neuron is going to fire and then evolve the membrane potentials of the two populations. The evolution of the two populations is different and it depends on the fact that the firing neuron belongs to the considered family or not. Let us schematize the algorithm in three steps:

1. As a first step we compare u_{1}^{(0)} with u_{1}^{(1)} to identify to which family the firing neuron belongs.

2. As a second step we check if the firing neuron belongs to a family which has already fired at the previous event or not. Depending on this, we have two possible alternatives: if the next and previous firing neurons belong to the same family we iterate the network as in point (a) below, otherwise as in point (b).

(a) Let us suppose that the firing population is the family (0). We evolve all the N +1 variables {ω_{j}^{(0)}}_{j=1,...,N+1} for the firing family, while, for the other

family, it is sufficient to evolve ω_{1}^{(1)} and ω_{N+1}^{(1)}. All the above mentioned variables are evolved for a lapse of time corresponding to the interval elapsed from the last firing event of the family (0).

(b) If the firing family is (1) and previously fired a neuron of family (0), the evolution is more complicated. The variables ω_{1}^{(k)} and ω_{N+1}^{(k)} of both families are integrated for the time interval elapsed from the last firing time of family (0). The N −1 variables {ω_{j}^{(1)}}_{j=2,...,N} should be instead evolved for a longer

time corresponding to the last interspike interval associated to family (1), because these variables have not been updated since the last firing of family (1).

3. The firing family is iterated in the comoving frame: this amounts to update the membrane potentials and to shift the index of all neurons by one unit. The membrane potentials of the other family are updated in the fixed reference frame.

4. The simulation is iterated by repeating the above three steps.

When we approach a linear stability analysis the previous algorithm is no more effective since now it is necessary to evolve all the variables at each time step in order to calculate the linearized equations in the tangent space. In this case we use the difference of the membrane potentials of two successive neurons instead of the logarithm. We still search for the first to fire neuron between the two populations and we treat differently the variables of the two populations at each time step depending on which neuron has emitted a spike previously. As in the previous case we employ different reference frames for the firing or not firing family (see point 3 above).

Appendix B

Linear Stability Analysis of the Splay State

B.1 Fixed-point expansion (general case)

The 1/N expansion of the exact expressions (2.15) for P̃ and Ẽ leads to

P̃ = \frac{α}{T^{(0)}} + [\frac{α}{2} − \frac{T^{(1)}}{T^{(0)}^{2}} ] \frac{α}{N} + [\frac{α^{2}T^{(0)}}{12} − \frac{T^{(2)}}{T^{(0)}^{2}} + \begin{matrix}T^{(1)} & 2 \\ T^{(0)} & 2 \end{matrix} ] \frac{α}{N^{2}}  (B.1)

+ [\frac{α^{2}T^{(1)}}{12} − \frac{T^{(3)}}{T^{(0)}^{2}} + 2\frac{T^{(1)}T^{(2)}}{T^{(0)}^{3}} − \begin{matrix}T^{(1)} & 3 \\ T^{(0)} & 4 \end{matrix} ] \frac{α}{N^{3}} + O (\frac{1}{N^{4}} ) ,

Ẽ = \frac{1}{T^{(0)}} − \frac{T^{(1)}}{T^{(0)}^{2}N} + [−\frac{α^{2}T^{(0)}}{12} − \frac{T^{(2)}}{T^{(0)}^{2}} + \begin{matrix}T^{(1)} & 2 \\ T^{(0)} & 3 \end{matrix} ] \frac{1}{N^{2}}  (B.2)

+ [−\frac{α^{2}}{12}T^{(1)} − \frac{T^{(3)}}{T^{(0)}^{2}} + 2\frac{T^{(1)}T^{(2)}}{T^{(0)}^{3}} − \begin{matrix}T^{(1)} & 3 \\ T^{(0)} & 4 \end{matrix} ] \frac{1}{N^{3}} + O (\frac{1}{N^{4}} ) ,

˙Ẽ = \frac{α^{2}}{2N} + \frac{α^{2}T^{(0)}}{6N^{2}} + \frac{α^{3}T^{(1)}}{6N^{3}} + O (\frac{1}{N^{4}} ) ,  (B.3)

where we have reported also the expansion of ˙ Ẽ that is necessary to pass from expression (2.13) to (2.18). Please notice that while the membrane potentials and the period are expanded up to O(1/N^{4}), as in (2.17) and (2.16), here we limit the expansion to O(1/N^{3}) terms, since the field variables appearing in the event-driven map are integrated over an interspike-interval (see (2.8).

To proceed further, we need also to introduce the expansions of the velocity field and of its derivatives,

F (ũ_{j}) = \overline{F}_{j} + \overline{F} ′_{j} \frac{ũ_{j}^{(1)}}{N} + \overline{F} ′_{j} \frac{ũ_{j}^{(2)}}{N^{2}} + \overline{F} ′_{j} \frac{ũ_{j}^{(3)}}{N^{3}} + \overline{F} ′′_{j} \frac{[ũ_{j}^{(1)}]^{2}}{2N^{2}}

+\overline{F} ′′_{j} \frac{ũ_{j}^{(1)}ũ_{j}^{(2)}}{N^{3}} + \overline{F} ′′′_{j} \frac{[ũ_{j}^{(1)}]^{3}}{6N^{3}} + O (\frac{1}{N^{4}} )

F′(ũ_{j}) = \overline{F} ′_{j} + \overline{F} ′′_{j} \frac{ũ_{j}^{(1)}}{N} + \overline{F} ′′_{j} \frac{ũ_{j}^{(2)}}{N^{2}} + \overline{F} ′′′_{j} \frac{[ũ_{j}^{(1)}]^{2}}{2N^{2}} + O (\frac{1}{N^{3}} )

F′′(ũ_{j}) = \overline{F} ′′_{j} + \overline{F} ′′′_{j} \frac{ũ_{j}^{(1)}}{N} + O (\frac{1}{N^{2}} ) ,

where the overline means that the function is computed in Ũ(x = j/N), which corresponds to the infinite N limit.

By replacing the membrane potentials, the period, the self-consistent fields and the velocity field with their expansions, the event-driven map (2.13) can be formally rewritten for the splay state as (2.18) with the introduction of the following auxiliary variables

Q^{(1)} = g + T^{(0)}\overline{F}_{j} ,  (B.4)
Q^{(2)} = T^{(1)}\overline{F}_{j} + [ũ_{j}^{(1)} + \frac{g}{2} + \frac{\overline{F}_{j}}{2} T^{(0)}] \overline{F} ′_{j}T^{(0)} ,  (B.5)

Q^{(3)} = [\overline{F} ′_{j}ũ_{j}^{(2)} + \frac{\overline{F} ′′_{j}}{2} [ũ_{j}^{(1)}]^{2} + g \frac{\overline{F} ′′_{j}}{2} ũ_{j}^{(1)} + \frac{g^{2}}{6} \overline{F} ′′_{j}  (B.6)

+ (\frac{2g}{3} \overline{F} ′′_{j}\overline{F}_{j} + \overline{F} _{j}^{′2}ũ_{j}^{(1)} + \overline{F} ′′_{j}\overline{F}_{j}ũ_{j}^{(1)} + \frac{g}{3}\overline{F} _{j}^{′2}) \frac{T^{(0)}}{2}

+ (\overline{F} ′′_{j}\overline{F}_{j} + \overline{F} _{j}^{′2}) \frac{T^{(0)}^{2}\overline{F}_{j}}{6} ] T^{(0)} + [ũ_{j}^{(1)} + \frac{g}{2} + \overline{F}_{j}T^{(0)}] T^{(1)}\overline{F} ′_{j} + T^{(2)}\overline{F}_{j} ,

Q^{(4)} = [\overline{F} ′_{j}ũ_{j}^{(3)} + (ũ_{j}^{(1)} + \frac{g}{2}) \overline{F} ′′_{j}ũ_{j}^{(2)} + (\frac{1}{6}[ũ_{j}^{(1)}]^{3} + \frac{g}{4}[ũ_{j}^{(1)}]^{2} + \frac{g^{2}}{6} ũ_{j}^{(1)} + \frac{g^{3}}{24}) \overline{F} ′′′_{j}] T^{(0)}

+ [(\overline{F} _{j}^{′2} + \overline{F}_{j}\overline{F} ′′_{j}) \frac{ũ_{j}^{(2)}}{2} + (3\overline{F} ′_{j}\overline{F} ′′_{j} + \overline{F}_{j}\overline{F} ′′′_{j}) \frac{[ũ_{j}^{(1)}]^{2}}{4}

+g (2\overline{F}′_{j}\overline{F} ′′_{j} + \overline{F}_{j}\overline{F} ′′′_{j}) \frac{ũ_{j}^{(1)}}{3} + \frac{g^{2}}{6} \overline{F} ′_{j}\overline{F} ′′_{j} + \frac{g^{2}}{8} \overline{F}_{j}\overline{F} ′′′_{j}] T^{(0)}^{2}

+ [(g + 4ũ_{j}^{(1)})\overline{F} _{j}^{′3} + (8g + 16ũ_{j}^{(1)})\overline{F}_{j}\overline{F} ′_{j}\overline{F} ′′_{j} + (3g + 4ũ_{j}^{(1)})\overline{F} _{j}^{2}\overline{F} ′′′_{j}] \frac{T^{(0)}^{3}}{24}

+ [\overline{F}_{j}\overline{F} _{j}^{′3} + 4\overline{F} _{j}^{2}\overline{F} ′_{j}\overline{F} ′′_{j} + \overline{F} _{j}^{3}\overline{F} ′′′_{j}] \frac{T^{(0)}^{4}}{24} + [\overline{F} ′_{j}ũ_{j}^{(2)} + \frac{\overline{F} ′′_{j}}{2} [ũ_{j}^{(1)}]^{2} + g \frac{\overline{F} ′′_{j}}{2} ũ_{j}^{(1)} + \frac{g^{2}}{6} \overline{F} ′′_{j}] T^{(1)}

+ [\frac{2g}{3} \overline{F} ′′_{j}\overline{F}_{j} + \overline{F} _{j}^{′2}ũ_{j}^{(1)} + \overline{F} ′′_{j}\overline{F}_{j}ũ_{j}^{(1)} + \frac{g}{3}\overline{F} _{j}^{′2}] T^{(0)}T^{(1)} + \frac{\overline{F}_{j}}{2} (\overline{F} ′′_{j}\overline{F}_{j} + \overline{F} _{j}^{′2}) T^{(0)}^{2}T^{(1)} + \overline{F}_{j}T^{(3)}

+ \{(ũ_{j}^{(1)} + \frac{g}{2}) T^{(2)} + \frac{\overline{F}_{j}}{2} [T^{(1)}^{2} + 2T^{(0)}T^{(2)}]\} \overline{F} ′_{j}.  (B.7)

B.2 Fixed-point expansion (LIF model)

In the case of the LIF neuron (see Eq. (2.2)), the fixed point of the event-driven map reads

u_{j−1} = e^{−τ}u_{j} + χ ,  (B.8)

where

χ = a(1 − e^{−τ}) + g \frac{e^{−τ} − e^{−ατ}}{α − 1} (E + \frac{P}{α − 1}) − g \frac{τ}{α − 1}e^{−ατ}P .  (B.9)

Its solution is

u_{j} = χ\frac{1 − e^{−NT +jτ}}{1 − e^{−τ}} .  (B.10)

By expanding Eq. (B.10) for j = 0 and for a generic j, one can derive perturbative expressions for the period T and the membrane potential, respectively. Let us start by substituting the expressions (2.16, B.1,B.2, B.3) in Eqs. (B.9). This leads to the expansion

χ = (a + \frac{g}{T} ) [τ − \frac{τ^{2}}{2} + \frac{τ^{3}}{6} − \frac{τ^{4}}{24}] + a\frac{τ^{5}}{120} + \frac{g}{T} \frac{τ^{5}}{720}K(α) + O(1/N^{4}) ,  (B.11)

where

K(α) = \frac{360α^{6} − 722α^{5} + 363α^{4} + 5α^{2} − 12α + 6}{(α − 1)^{2}} ,  (B.12)

accounts for the dependence on the field dynamics. Now, with the help of Eqs. (2.16,B.11) and expanding the exponential terms up to the fourth order, we obtain a closed equation for the interspike interval,

u_{0} = 1 = (a + \frac{g}{T^{(0)}} ) (1 − e^{−T^{(0)}}) + \frac{T^{(1)}}{N} ξ(T^{(0)}) + \frac{1}{N^{2}} [T^{(2)}ξ(T^{(0)}) + T^{(1)}W_{1}^{(2)}]

+\frac{1}{N^{3}} [T^{(3)}ξ(T^{(0)}) + T^{(1)}W_{1}^{(3)} + T^{(2)}W_{2}^{(3)}] + \frac{1}{N^{4}} [T^{(4)}ξ(T^{(0)}) + ζ(T^{(0)}) + T^{(1)}W_{1}^{(4)}

+T^{(2)}W_{2}^{(4)} + T^{(3)}W_{3}^{(4)}] ,  (B.13)

where

ζ(T^{(0)}) = −(1 − e^{−T^{(0)}})\frac{g}{120}T^{(0)}^{3}(1 − \frac{K(α))}{6}) ,

ξ(T^{(0)}) = (a + \frac{g}{T^{(0)}} ) e^{−T^{(0)}} − \frac{g}{T^{(0)}^{2}} (1 − e^{−T^{(0)}}) ,

while W_{k}^{(h)} identifies a term of order 1/N^{h} that is multiplied by T^{(k)}. Since, while proceeding from lower to higher-order terms, we find that T^{(k)} = 0 (for k < 4), it is not necessary to give the explicit expression of the W_{k}^{(h)} functions as they do not contribute at all.

One can equivalently expand u_{j}

u_{j} = u_{j}^{(0)} + \frac{u_{j}^{(1)}}{N} + \frac{u_{j}^{(2)}}{N^{2}} + \frac{u_{j}^{(3)}}{N^{3}} + \frac{u_{j}^{(4)}}{N^{4}} = (a + \frac{g}{T^{(0)}} ) (1 − e^{T^{(0)}( \frac{j}{N} −1)}) + \frac{T^{(1)}}{N} ς(T^{(0)})

+\frac{1}{N^{2}} [T^{(2)}ς(T^{(0)}) + T^{(1)}Z_{1}^{(2)}] + \frac{1}{N^{3}} [T^{(3)}ς(T^{(0)}) + T^{(1)}Z_{1}^{(3)} + T^{(2)}Z_{2}^{(3)}]

+\frac{1}{N^{4}} [ς(T^{(0)})T^{(4)} + σ(T^{(0)}) + T^{(1)}Z_{1}^{(4)} + T^{(2)}Z_{2}^{(4)} + T^{(3)}Z_{3}^{(4)}] ,  (B.14)

where

ς(T^{(0)}) = − (a + \frac{g}{T^{(0)}} ) e^{T^{(0)}( \frac{j}{N} −1)} (\frac{j}{N} − 1) − \frac{g}{T^{(0)}^{2}} [1 − e^{T^{(0)}( \frac{j}{N} −1)}] ,

σ(T^{(0)}) = \frac{g}{120}(1 − e^{T^{(0)}( \frac{j}{N} −1)})T^{(0)}^{3} (\frac{K}{6} − 1) ,

while we do not provide explicit expressions for Z_{k}^{(h)} as they turn out to be irrelevant. Now we are in the position to analyze the different orders.

B.2.1 Zeroth Order

By assembling the terms of order 1 in Eq. (B.13), we obtain

(a + \frac{g}{T^{(0)}} ) (1 − e^{−T^{(0)}}) = 1 .

This is an implicit definition of the asymptotic interspike time T^{(0)}

T^{(0)} = ln (\frac{aT^{(0)} + g}{T^{(0)}(a − 1) + g} ) .

Analogously, we can find an explicit equation for the membrane potential by assembling the terms of order 1 in Eq. (B.14)

u_{j}^{(0)} = (a + \frac{g}{T^{(0)}} ) [1 − e^{T^{(0)}(\frac{j}{N} −1)}] .

In the thermodynamic limit the solution for u_{j}^{(0)} becomes

U^{(0)}(x) = (a + \frac{g}{T^{(0)}} ) [1 − e^{T^{(0)}(x−1)}] ,

which coincides with Eq. (2.22) with F = a − U^{(0)}.

B.2.2 From first to third order

By separately assembling the terms of order 1/N^{h} (for h = 1, 2, 3) in Eq. (B.13), we obtain

T^{(h)}ξ(T^{(0)}) = 0 ,

which implies that T^{(h)} = 0 since ξ ≠ 0. Moreover, by assembling the terms of order 1/N^{h} in Eq. (B.14), we obtain

u_{j}^{(h)} = ς(T^{(0)})T^{(h)} , j = 1, ⋯, N ,

which thereby implies that first, second and third order corrections vanish also for the membrane potential.

B.2.3 Fourth Order

The order which reveals a different scenario is the fourth one. By assembling the terms of order 1/N^{4} in Eq. (B.13) we obtain

T^{(4)} = −\frac{ζ(T^{(0)})}{ξ(T^{(0)})} ,

whose explicit expression is reported in Eq. (2.5). By analogously assembling the terms of order 1/N^{4} in Eq. (B.14), we obtain

u_{j}^{(4)} = ς(T^{(0)})T^{(4)} + σ(T^{(0)}) ,

which becomes, in the thermodynamic limit,

U^{(4)}(x) = − (a + \frac{g}{T^{(0)}} ) T^{(4)}e^{T^{(0)}(x−1)}(x − 1) − g \frac{T^{(4)}}{T^{(0)}^{2}} [1 − e^{T^{(0)}(x−1)}]

+\frac{g}{120}(1 − e^{T^{(0)}(x−1)})T^{(0)}^{3} (\frac{K(α)}{6} − 1) .

B.3 Expansion in tangent space around the fixed point

B.3.1 Introduction

The auxiliary variables required to complete the definition of the tangent-space evolution rule (2.29,2.30,2.31,2.32) are as follows

Z_{1} = (1 − ατ + \frac{α^{2}}{2} τ^{2} − \frac{α^{3}}{6} τ^{3} + \frac{α^{4}}{24}τ^{4}) , Z_{2} = (−α + α^{2}τ − \frac{α^{3}}{2} τ^{2} + \frac{α^{4}}{6} τ^{3}) ,

Z_{3} = (τ − ατ^{2} + \frac{α^{2}}{2} τ^{3} − \frac{α^{3}}{6} τ^{4}) , Z_{4} = (1 − 2ατ + \frac{3}{2}α^{2}τ^{2} − \frac{2}{3}α^{3}τ^{3} + \frac{5}{24}α^{4}τ^{4}) ,

W_{1} = \{1 + \overline{F} ′_{j}τ + \frac{τ^{2}}{2} (\overline{F} ′′_{j}\overline{F}_{j} + \overline{F}′_{j}^{2}) + \frac{τ^{3}}{6} [\overline{F} ′′′_{j}\overline{F}_{j}^{2} + \overline{F} ′′_{j} (4\overline{F} ′_{j}\overline{F}_{j} + g ˙Ẽ) + \overline{F} _{j}^{2}\overline{F} ′_{j}]\} ,

W_{2} = [gτ + \frac{τ^{2}}{2} g\overline{F} ′_{j} + \frac{τ^{3}}{6} (2g\overline{F} ′′_{j}\overline{F}_{j} + g\overline{F} _{j}^{2})] ,

W_{3} = [g \frac{τ^{2}}{2} + \frac{τ^{3}}{6} (g\overline{F} ′_{j} − gα)] ,

W_{4} = \{\overline{F}_{j} + [\overline{F} ′_{j}\overline{F}_{j} + g ˙Ẽ] τ + [\overline{F} ′′_{j}\overline{F}_{j}^{2} + (\overline{F} ′_{j})^{2}\overline{F}_{j} + g\overline{F} ′_{j} ˙Ẽ − gα( ˙Ẽ + P̃)] \frac{τ^{2}}{2}

+ [(\overline{F} ′_{j})^{3}\overline{F}_{j} + 4\overline{F} ′_{j}\overline{F} ′′_{j}\overline{F}_{j}^{2} + \overline{F} ′′′_{j}\overline{F}_{j}^{3} + gα(2α − \overline{F} ′_{j})P̃] \frac{τ^{3}}{6} \} ,

W_{5} = \{gτ − g [\frac{1}{2}(\overline{F} ′_{1} + α) + \frac{1}{\overline{F}_{1}} g ˙Ẽ] τ^{2} + \frac{g}{\overline{F}_{1}} [gα( ˙Ẽ + \frac{P̃}{2} ) + \frac{\overline{F}_{1}}{6} (α^{2} − 5\overline{F′}_{1}^{2} − \overline{F} ′′_{1}\overline{F}_{1} + 2α\overline{F} ′_{1})

+\frac{1}{\overline{F}_{1}} (\overline{F} ′_{1}\overline{F}_{1} + g ˙Ẽ)^{2}] τ^{3}\} ,

W_{6} = \{\frac{g}{2}τ^{2} − \frac{g}{2}τ^{3} [\frac{2}{3}(α + \overline{F} ′_{1}) + \frac{1}{\overline{F}_{1}} g ˙Ẽ]\} ,

W_{7} = \{1 − \frac{g}{\overline{F}_{1}} τ ˙Ẽ + τ^{2} [\frac{gα}{2\overline{F}_{1}} ( ˙Ẽ + P̃) + \frac{1}{\overline{F}_{1}^{2}} (\overline{F} ′_{1}\overline{F}_{1} + g ˙Ẽ)^{2} − \overline{F′}_{1}^{2} − \frac{3}{2} \frac{g}{\overline{F}_{1}} \overline{F} ′_{1} ˙Ẽ]

−τ^{3} [\frac{1}{3\overline{F}_{1}} gα(α + \overline{F} ′_{1})P̃ + \frac{1}{\overline{F}_{1}^{2}} gα ˙Ẽ(\frac{1}{2}\overline{F} ′_{1}\overline{F}_{1} + gP̃) + \frac{1}{\overline{F}_{1}^{3}} g^{2} ˙Ẽ^{3} + \frac{1}{\overline{F}_{1}^{2}} g^{2} ˙Ẽ^{2}(α + \overline{F} ′_{1}) − \frac{2}{3}g\overline{F} ′′_{1} ˙Ẽ]\} .

where the dependence of τ on n has been dropped, since we are considering a linearization around the splay state.

In order to find the Floquet eigenvalue μ_{k}, one should substitute the Ansätze (2.33) into Eqs. (2.29,2.30). This allows to find explicit expressions for δP and δE as a function of μ_{k}, τ and δτ, namely

δP = −\frac{α^{2}}{μ_{k} − 1} [1 − \frac{ατ}{2} \frac{μ_{k} + 1}{μ_{k} − 1} + α^{2}τ^{2}M_{k} − \frac{α^{3}τ^{3}}{2} \frac{μ_{k}(μ_{k} + 1)}{(μ_{k} − 1)^{3}}] \frac{δτ}{T} ,  (B.15)

δE = −\frac{α}{μ_{k} − 1} [\frac{ατ}{2} \frac{(μ_{k} + 1)}{μ_{k} − 1} − 2α^{2}τ^{2}M_{k} − \frac{3α^{3}τ^{3}}{2} \frac{μ_{k}(μ_{k} + 1)}{(μ_{k} − 1)^{3}}] \frac{δτ}{T} ,  (B.16)

where we have introduced the shorthand notation

M_{k} = \frac{μ_{k}^{2} + 10μ_{k} + 1}{12(μ_{k} − 1)^{2}} .

By substituting δP and δE, as given by (B.15) and (B.16), into Eq. (2.32), we can express δτ directly in terms of δu_{1}

δτ = − \{\overline{F}_{1} + \frac{gα^{2}}{T} M_{k}τ^{2}  (B.17)

−\frac{gα^{2}}{T} [\frac{\overline{F} ′_{1}}{12} \frac{(μ_{k} + 5)}{(μ_{k} − 1)^{2}} + α\frac{(μ_{k} + 1)}{(μ_{k} − 1)^{3}} ] μ_{k}τ^{3}\} \frac{δu_{1}}{\overline{F}_{1}^{2}} ,

where we exploited the equality \overline{F}_{j} ≡ \overline{F}_{j} + \frac{g}{T} which follows from the fact that in the thermodynamic limit Ẽ = \frac{1}{T}.

By inserting the expressions in Eqs. (B.15, B.16, B.17) into Eq. (2.31), we find a single equation for the eigenvalues and eigenvectors,

μ_{k}δu_{j−1} = \{1 + \overline{F} ′_{j}τ + [\overline{F} ′′_{j}\overline{F}_{j} + \overline{F} _{j}^{′2}] \frac{τ^{2}}{2} + [\overline{F} ′′′_{j}\overline{F}_{j}^{2} + 4\overline{F} ′_{j}\overline{F} ′′_{j}\overline{F}_{j} + \overline{F} _{j}^{′3}] \frac{τ^{3}}{6} \} δu_{j}

− \{\overline{F}_{j} + \overline{F} ′_{j}\overline{F}_{j}τ + [\frac{\overline{F} ′′_{j}}{2} \overline{F}_{j}^{2} + \frac{\overline{F} _{j}^{′2}}{2} \overline{F}_{j} + \frac{gα^{2}}{T} M_{k} (\begin{matrix}\overline{F} & j \\ \overline{F} & 1 \end{matrix} − 1)] τ^{2}

+ [\frac{\overline{F} _{j}^{′3}}{6} \overline{F}_{j} + \frac{2}{3}\overline{F} ′_{j}\overline{F} ′′_{j}\overline{F}_{j}^{2} + \frac{\overline{F} ′′′_{j}}{6} \overline{F}_{j}^{3} \frac{gα^{2}}{T} \frac{5μ_{k} + 1}{12(μ_{k} − 1)^{2}} (\begin{matrix}\overline{F} & j \\ \overline{F} & 1 \end{matrix} \overline{F} ′_{1} − \overline{F} ′_{j})

+\frac{gα^{3}}{T} \frac{μ_{k}(μ_{k} + 1)}{(μ_{k} − 1)^{3}} (1 − \begin{matrix}\overline{F} & j \\ \overline{F} & 1 \end{matrix} ) + \frac{gα^{2}}{T} \begin{matrix}\overline{F} & j \\ \overline{F} & 1 \end{matrix} (\overline{F} ′_{j} − \overline{F} ′_{1})M_{k}] τ^{3}\} \frac{δu_{1}}{\overline{F}_{1}} .

By now substituting the μ_{k} expansion (2.34) and retaining the leading terms, we obtain Eq. (2.35).

B.3.2 N → ∞ limit

Once the continuous variables (2.37) have been introduced, it is necessary to estimate U(1/N) and δU(1/N), by expanding such variables around zero. By inserting the resulting expansion for U(1/N) into the expressions for \overline{F}_{1} and \overline{F}_{1}, we obtain, respectively

F (U(\frac{1}{N})) = F (U(0)) + \frac{F′(U(0))}{N} \frac{dU}{dx} |_{0} + \frac{F′(U(0))}{2N^{2}} \frac{d^{2}U}{dx^{2}} |_{0} + \frac{F′′(U(0))}{2} (\frac{1}{N} \frac{dU}{dx} |_{0})^{2}

+\frac{F′′(U(0))}{2N^{3}} \frac{dU}{dx} |_{0} \frac{d^{2}U}{dx^{2}} |_{0} + \frac{F′′′(U(0))}{6} (\frac{1}{N} \frac{dU}{dx} |_{0})^{3} + \frac{F′(U(0))}{6N^{3}} \frac{d^{3}U}{dx^{3}} |_{0} + O(\frac{1}{N^{4}} ) ,

\frac{1}{F(U(\frac{1}{N}))} ≡ \frac{1}{F (U(\frac{1}{N})) + \frac{g}{T}} = \frac{1}{F(U(0))} + \frac{1}{N} \frac{T F′(U(0))}{F(U(0))}

+\frac{1}{2N^{2}} \frac{[T F′(U(0))]^{2}}{F(U(0))} − \frac{1}{2N^{2}} F′′(U(0))T^{2} + \frac{T^{3}}{6N^{3}} [F′′′(U(0))F(U(0)) − 2F′(U(0))F′′(U(0))

−\frac{g}{T^{2}} (\frac{F′(U(0))}{F(U(0))})^{2} − 5\frac{[F′(U(0))]^{3}}{F(U(0))} + 6\frac{F′(U(0))}{F(U(0))}] + O(\frac{1}{N^{4}} ) .

An analogous procedure for δU(1/N) leads to

δU(1/N) = e^{iϕ_{k}} [Θ^{(0)}(\frac{1}{N}) + \frac{Θ^{(1)}(\frac{1}{N})}{N} + \frac{Θ^{(2)}(\frac{1}{N})}{N^{2}}] + Π^{(0)}(\frac{1}{N}) + \frac{Π^{(1)}(\frac{1}{N})}{N} + \frac{Π^{(2)}(\frac{1}{N})}{N^{2}} + O (\frac{1}{N^{3}} )

≡ C^{(0)} + \frac{C^{(1)}}{N} + \frac{C^{(2)}}{N^{2}} + O (\frac{1}{N^{3}} ) ,

where C^{(0)}, C^{(1)}, C^{(2)} are defined according to the following equations

C^{(0)} = e^{iϕ_{k}}Θ^{(0)}(0) + Π^{(0)}(0) ,

C^{(1)} = e^{iϕ_{k}} [\frac{dΘ^{(0)}}{dx} |_{0} + Θ^{(1)}(0)] + \frac{dΠ^{(0)}}{dx} |_{0} + Π^{(1)}(0) ,

C^{(2)} = e^{iϕ_{k}} [\frac{1}{2} \frac{d^{2}Θ^{(0)}}{dx^{2}} |_{0} + \frac{dΘ^{(1)}}{dx} |_{0} + Θ^{(2)}(0)] + \frac{1}{2} \frac{d^{2}Π^{(0)}}{dx^{2}} |_{0} + \frac{dΠ^{(1)}}{dx} |_{0} + Π^{(2)}(0) .

We now expand δU(1/N) up to the order O (\frac{1}{N^{3}} ), thus neglecting higher orders, because they contribute to the definition of Π variable and we need terms at least of order O (\frac{1}{N^{3}} ) (one order lower than needed to define Θ). By inserting the Ansatz (2.36) and the previous expansions in Eq. (2.35), we finally obtain a closed equation for the eigenvalues and eigenvectors,

e^{iϕ_{k}} \{Π^{(0)} + [Π^{(1)} − Π^{(0)}′ + Π^{(0)}Γ^{(1)}] \frac{1}{N} + [\frac{Π^{(0)}′′}{2} − Π^{(1)}′ + Π^{(2)} − Π^{(0)}′Γ^{(1)} + Π^{(1)}Γ^{(1)}

+Π^{(0)}Γ^{(2)}] \frac{1}{N^{2}} + [\frac{Π^{(1)}′′}{2} − \frac{Π^{(0)}′′′}{6} − Π^{(2)}′ + Π^{(3)} + \frac{Π^{(0)}′′Γ^{(1)}}{2} − Π^{(1)}′Γ^{(1)} + Π^{(2)}Γ^{(1)}

−Π^{(0)}′Γ^{(2)} + Π^{(1)}Γ^{(2)} + Π^{(0)}Γ^{(3)}] \frac{1}{N^{3}} \} − Π^{(0)} − [Π^{(1)} + Π^{(0)}A^{(1)}] \frac{1}{N} − [Π^{(2)} + Π^{(0)}A^{(2)}

+Π^{(1)}A^{(1)}] \frac{1}{N^{2}} − [Π^{(3)} + Π^{(0)}A^{(3)} + Π^{(1)}A^{(2)} + Π^{(2)}A^{(1)}] \frac{1}{N^{3}} + C^{(0)}B^{(0)} + [C^{(0)}B^{(1)}

+C^{(1)}B^{(0)}] \frac{1}{N} + [C^{(0)}B^{(2)} + C^{(1)}B^{(1)} + C^{(2)}B^{(0)}] \frac{1}{N^{2}} + \frac{B}{N^{3}} = e^{iϕ_{k}j} \{[Θ^{(0)}′ − Θ^{(0)}Γ^{(1)}

+Θ^{(0)}A^{(1)}] \frac{1}{N} + [(Θ^{(0)}A^{(2)} − \frac{Θ^{(0)}′′}{2} ) + Θ^{(1)}′ + Θ^{(0)}′Γ^{(1)} − Θ^{(1)}Γ^{(1)} − Θ^{(0)}Γ^{(2)}

+Θ^{(1)}A^{(1)}] \frac{1}{N^{2}} + [(\frac{−Θ^{(1)}′′}{2} + Θ^{(1)}A^{(2)}) + (\frac{Θ^{(0)}′′′}{6} + Θ^{(0)}A^{(3)}) + Θ^{(2)}′ − \frac{Θ^{(0)}′′Γ^{(1)}}{2}

+Θ^{(1)}′Γ^{(1)} − Θ^{(2)}Γ^{(1)} + Θ^{(0)}′Γ^{(2)} − Θ^{(1)}Γ^{(2)} − Θ^{(0)}Γ^{(3)} + Θ^{(2)}A^{(1)}] \frac{1}{N^{3}} \} ,  (B.18)

where we have introduced the shorthand notation B in order to characterize a term of order O (\frac{1}{N^{3}} ), whose explicit expression is not necessary, since it turns out to contribute to the definition of the Π variable, and it is therefore one order beyond what we need. Moreover, notice that the terms appearing within round brackets in the rhs of the above equation can be shown to be zero, due to exact algebric cancellations that emerge from the solution of the equation order by order. Finally,

A^{(1)}(U(x)) = T F′(U(x)), A^{(2)}(U(x)) = \frac{T^{2}}{2} \{F′′(U(x))F(U(x)) + [F′(U(x))]^{2}\} ,

A^{(3)}(U(x)) = \frac{T^{3}}{6} \{F′′′(U(x))[F(U(x))]^{2} + 4F′(U(x))F′′(U(x))F(U(x)) + [F′(U(x))]^{3}\} ,

B^{(0)}(U(x)) = \frac{F(U(x))}{F(U(0))}, B^{(1)}(U(x)) = [T F′(U(x)) + T F′(U(0))] \frac{F(U(x))}{F(U(0))} ,

B^{(2)}(U(x)) = T^{2} \{−\frac{g}{T} α^{2} \frac{e^{2iϕ_{k}} + 10e^{iϕ_{k}} + 1}{12(e^{iϕ_{k}} − 1)^{2}} \frac{F (U(0)) − F (U(x))}{[F(U(0))]^{2}}

+\frac{F′′(U(x))}{2} \frac{[F(U(x))]^{2}}{F(U(0))} + \frac{[F′(U(x))]^{2}}{2} \frac{F(U(x))}{F(U(0))}\}

+ \{T F′(U(x))T F′(U(0)) + \frac{[T F′(U(0))]^{2}}{2} \} \frac{F(U(x))}{F(U(0))}

−\frac{T^{2}}{2} F′′(U(0))F(U(x)) .

B.4 Comparison with previous results

Abbott and van Vreeswijk consider in [25] N identical oscillators with uniform all-to-all coupling in the large N limit. The state of an oscillator labeled by the index i is described by a state variable u_{i} ∈ [0, 1] which satisfies the following equation

\frac{du_{i}}{dt} = F (u_{i}) + G(u_{i})E(t) .  (B.19)

The function F characterizes the behavior of the oscillator in the absence of coupling and is an arbitrary positive-definite function. G determines the dependence of the coupling on u_{i} and E(t) is a dynamic variable characterizing the inputs coming from the other oscillators in the population. The sign of G determines whether the coupling is excitatory or inhibitory.

Equation B.19 determines the behavior of u_{i} in the range between zero and one. When u_{i} reaches the value one, it is immediately reset to zero. This resetting corresponds to the firing of a pulse and results in a contribution to the coupling variable E(t) which is incremented by an amount

E(t) → E(t) + \frac{α^{2}}{N} (t − t_{0})e^{α(t_{0}−t)},  (B.20)

where t_{0} is the time at which the oscillator i reaches u_{i} = 1.

The asynchronous state is characterized by a coupling function E(t) = E_{0} that is time independent corresponding to asynchronous pulsing. By assuming that the oscillators pulse at a constant, single-oscillator rate r then, from Eq. (B.20), they have E(t) = E_{0} = r. Thus E_{0} is the steady-state, single-oscillator pulsing or firing rate. Integrating Eq. (B.19) over one period 1/E_{0} Abbott and van Vreeswijk find that E_{0} must satisfy

T^{(0)} = \frac{1}{E_{0}} = ∫_{0}^{1} \frac{dx}{F (x) + E_{0}G(x)} .

This equation always has a solution (for F > 0) if F and G are continuous with bounded derivatives unless G > 0 and ∫_{0}^{1} dx/G(x) < 1. To analyze the asynchronous state it is convenient to make a change of variables, replacing u_{i} by

y_{i} = ∫_{0}^{u_{i}} \frac{E_{0}dx}{F (x) + E_{0}G(x)},

which varies between zero and one and satisfies the equation

\frac{dy_{i}}{dt} = E_{0} + Γ(y_{i})ε(t)

where ε(t) = E(t) − E_{0} and Γ(y) = \frac{E_{0}G(u)}{F (u)+E_{0}G(u)}. In the following we will consider G(u) = g = constant in order to better clarify the relationship between the results shown in Chapter 2 and in [25]. To describe a system of all-to-all coupled oscillators in absence of noise it is then necessary to introduce a density function ρ(y, t) and a flux J(y, t) in order to write a continuity equation in the range 0 < y < 1 with boundary conditions J(0, t) = J(1, t), expressing the fact that once an oscillator pulses or fires it is reset instantaneously to zero. The asynchronous solution described in terms of the variable y has ρ(y) = 1 and J(y) = E_{0}. To investigate whether the asynchronous state is stable, one should expand around this solution considering small fluctuations and look for solutions of the linearized equations. Following [25], the eigenvalues of the stability matrix are given by

E_{0}(e^{λ/E_{0}} − 1)(λ + α)^{2} = α^{2}λ ∫_{0}^{1} dyΓ(y)e^{λy/E_{0}}.  (B.21)

Equation (B.21) has an infinite number of solutions giving the full spectrum of eigenvalues of the infinite-dimensional stability matrix. Finding the eigenvalues is, in general, a difficult problem, but if the coupling function Γ is independent of y they can be found easily. If Γ(y) is a constant, there is a set of imaginary eigenvalues given (for integer k) by

λ = 2πikE_{0}, k ≠ 0.

Since these are pure imaginary , the asynchronous state with y-independent coupling is at best marginally stable. In addition to these imaginary eigenvalues, there are eigenvalues

λ = −α ± [α^{2}Γ]^{1/2}.

In order for the y-independent case to be marginally stable, we must have Γ < 1. A special case of independent coupling is no coupling at all, Γ = 0, in which case these eigenvalues become simply λ = −α. The eigenvalues obtained for zero coupling can be used to construct a weak-coupling perturbation expansion for the case of small Γ(y). Since the zero coupling case is marginally stable with imaginary eigenvalues, small coupling adds a real part to these eigenvalues determining the stability or instability of the asynchronous state. Since these imaginary eigenvalues always come in complex conjugate pairs, we will consider only the eigenvalues with k > 0. To first order in Γ, the eigenvalues that are purely imaginary at lowest order become

λ = 2πikE_{0} [1 + \frac{α^{2}(A_{k} + iB_{k})}{(2πikE_{0} + α)^{2}} ] ,

where A_{k} + iB_{k} = ∫_{0}^{1} dyΓ(y)e^{2πiky}. Stability of the asynchronous state requires that the real parts of these eigenvalues be negative or, equivalently

B_{k}(α^{2} − 4π^{2}k^{2}E_{0}^{2}) > 4παkE_{0}A_{k}.

For large k and non-pathological functions Γ(y), B_{k} ∝ 1/k and A_{k} ∝ 1/k^{2}. From this we see that, by neglecting the α^{2} terms in between the brackets, the real parts of the eigenvalues, in the large k limit, read as

Re(λ_{k}) = \frac{α^{2}B_{k}}{2πkE_{0}} .  (B.22)

In addition to this, for large k, the stability conditions requires B_{k} < 0, which, from standard theorems on Fourier transforms, implies under quite general conditions that Γ(1) > Γ(0). This is algebraically equivalent to the condition F (1) < F (0). The Eq. ( B.22) can be evaluated in the large k limit, thus recovering the Eq. ( 2.7) derived independently. In order to demonstrate this equivalence we compute

B_{k} = ∫_{0}^{1} dyΓ(y) sin(2πky) = ∫_{0}^{1} dy \frac{gE_{0} sin(2πky)}{F (y) + gE_{0}}.

By using Eq. (B.2) approximated to the zero order and approximating F (y) as a piecewise linear function, the previous equation reads as

B_{k} = ∑_{j=1}^{k} ∫ j/k_{\frac{j−1}{k}} dy \frac{g sin(2πky)}{T^{(0)}F (\frac{j−1}{k}) + g} [1 − \frac{T^{(0)}F′(y − \frac{j−1}{k})}{T^{(0)}F (\frac{j−1}{k}) + g}] = ∑_{j=1}^{k} \frac{gT^{(0)}F′(\frac{j−1}{k})}{[T^{(0)}F (\frac{j−1}{k}) + g]^{2}} \frac{1}{2πk^{2}} .

Since 1/k is small we can take the continuum limit and switch from the sum to an integral in the previous equation as follows

B_{k} ≃ ∫_{0}^{1} dx\frac{F′(x)}{[T^{(0)}F (x) + g]^{2}} \frac{gT^{(0)}}{2πk} = \frac{F (1) − F (0)}{(F (0) + \frac{g}{T^{(0)}} ) (F (1) + \frac{g}{T^{(0)}} )} \frac{g}{2πkT^{(0)}} . (B.23)

Finally, substituting Eq. (B.23) in Eq. (B.22) we obtain

Re(λ_{k}) = \frac{α^{2}}{(2πk)^{2}} \frac{F (1) − F (0)}{(F (0) + \frac{g}{T^{(0)}} ) (F (1) + \frac{g}{T^{(0)}} )} ,  (B.24)

which is exactly Eq. (2.7), since T^{(0)} = 1/E_{0}.

Bibliography

[1] L. Ingber, Generic mesoscopic neural networks based on statistical mechanics of neocortical interactions. Phys. Rev. A 45, R2183-2186 (1992).

[2] Y. Imry, Introduction to Mesoscopic Physics, New York, Oxford University Press (1997).

[3] W. J. Freeman, Neurodynamics: An Exploration in Mesoscopic Brain Dynamics, Berlin, Springer (2000).

[4] H. Haken, Mesoscopic levels in science - some comments, in Micro - Meso - Macro: Addressing Complex Systems Couplings, 1924, Eds. H. Liljenström, U. Svedin, London, World Scientific (2005).

[5] H. Haken, Synergetics, an Introduction: Nonequilibrium Phase Transitions and Self-Organization in Physics, Chemistry, and Biology, New York, Springer-Verlag (1983).

[6] D. Walgraef, Spatio-temporal pattern formation: With examples from physics, chemistry, and materials science, New York, Springer-Verlag (1997).

[7] B. Ermentrout, Neural networks as spatio-temporal pattern-forming systems , Rep. Prog. Phys. 61, 353-430 (1998).

[8] S. L. Bressler and V. Menon, Large-scale brain networks in cognition: emerging methods and principles, Trends in Cognitive Sciences 14, 277-290 (2010).

[9] G. Buzsáki, Rhythms of the Brain, Oxford University Press (2006).

[10] V. Gallese, L. Fadiga, L. Fogassi and G. Rizzolatti, Action recognition in the premotor cortex, Brain, a Journal of Neurology 119, 593-609 (1996).

[11] W. Singer, Putative functions of temporal correlations in neocortical processing, in Large-scale neuronal theories of the brain, 201-237, Eds. C. Koch and J. L. Davis, Cambridge, The MIT Press (1994).

[12] F. Rieke, D. Warland, R. de Ruyter van Steveninck and W. Bialek, Spikes, exploring the neural code, Cambridge, The MIT Press (1999).

[13] H. Liljenström, Autonomous learning with complex dynamics. Intl. J. Intelligent Systems 10, 119-153 (1995).

[14] H. Liljenström and M. E. Hasselmo ME, Cholinergic modulation of cortical oscillatory dynamics, J. Neurophysiol. 74, 288-297 (1995).

[15] X. Wu and H. Liljenström H, Regulating the nonlinear dynamics of olfactory cortex, Network: Computation in Neural Systems 5, 47-60 (1994).

[16] G. Laurent, Olfactory network dynamics and the coding of multidimensional signals, Nature Reviews Neuroscience 3, 884-895 (2002).

[17] G. Halnes, H. Liljenström and P. Århem, Density dependent neurodynamics, BioSystems 89, 126-134 (2007).

[18] H. Liljenström H and P. Århem, Investigating amplifying and controlling mechanisms for random events in neural systems, in Computational Neuroscience, 711-716, Ed. J. M. Bower, New York, Plenum Press. (1997).

[19] S. Basu and H. Liljenström, Spontaneously active cells induce state transitions in a model of olfactory cortex, Biosystems 63, 5769 (2001).

[20] K. Wiesenfeld and F. Moss, Stochastic resonance and the benefits of noise: from ice ages to crayfish and SQUIDs, Nature 373, 33-36 (1995).

[21] R. Benzi, A. Sutera and A. Vulpiani, The mechanism of stochastic resonance, J. Phys. A: Math. Gen. 14, L453-L457 (1981).

[22] B. Lindner, J. Garca-Ojalvo, A. Neiman and L. Schimansky-Geier, Effects of noise in excitable systems, Physics Reports 392:6, 321424 (2004).

[23] H. Liljenström, Global effects of fluctuations in neural information processing, Intl. J. Neural Systems 7, 497-505 (1996).

[24] D. A. Steyn-Ross and M. L. Steyn-Ross, Modeling Phase Transitions in the Brain, New York, Springer (2010).

[25] L. F. Abbott and C. van Vreeswijk, Asynchronous states in networks of pulse-coupled oscillators, Phys. Rev. E, 48, 1483 (1993).

[26] R. Zillmer, R. Livi, A. Politi, and A. Torcini, Desynchronization in diluted neural networks, Phys. Rev. E 74 036203 (2006)

[27] A. N. Burkitt, A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input, Biol Cybern 95, 119 (2006).

[28] C. van Vreeswijk, Partial synchronization in populations of pulse-coupled oscillators, Phys. Rev. E, 54, 5522 (1996).

[29] P. K. Mohanty and A. Politi, A new approach to partial synchronization in globally coupled rotators, J. Phys. A, 39, L415 (2006).

[30] R. Brette, Exact Simulation of Integrate-and-Fire Models with Exponential Currents, Neural Computation 19:10, 2604-2609 (2007).

[31] Y. Kuramoto, Chemical oscillations, waves, and turbulence, Berlin, Springer-Verlag, (1984).

[32] D. Z. Jin, Fast convergence of spike sequences to periodic patterns in recurrents networks, Phys. Rev. Lett. 89, 208102 (2002).

[33] R. Zillmer, R. Livi, A. Politi and A. Torcini, Stability of the splay state in pulse-coupled networks, Phys. Rev. E 76, 046102 (2007).

[34] S. Olmi, A. Politi and A. Torcini, Stability of the splay state in networks of pulse-coupled neurons, The Journal of Mathematical Neuroscience 2, 12 (2012).

[35] D. Golomb, D. Hansel and G. Mato, Mechanisms of Synchrony of neural Activity in large Networks, in Handbook of biological physics, 887-967, Eds. Gielen S and Moss F, Elsevier, Amsterdam (2001).

[36] R. E. Mirollo and S. H. Strogatz SH, Synchronization of pulse-coupled biological oscillators, SIAM J. Appl. Math., 50, 1645-1662 (1990).

[37] W. Gerstner and W. M. Kistler, Spiking Neuron Models: Single Neurons, Populations, Plasticity. Cambridge University Press (2002).

[38] P. Dayan and L. F. Abbott, Theoretical neuroscience, MIT press (2001).

[39] P. Goel and B. Ermentrout, Synchrony, stability, and firing patterns in pulse-coupled oscillators, Physica D, 163, 191216 (2002).

[40] N. Brunel, Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons, J. Comput. Neurosci., 8, 183 (2000).

[41] P. C. Bressloff, Mean-field theory of globally coupled integrate-and-fire neural oscillators with dynamic synapses, Phys. Rev. E, 60, 2160-2170 (1999).

[42] C. C. Chow and N. Kopell, Dynamics of Spiking Neurons with Electrical Coupling, Neural Comp., 12, 1643-1678 (2000).

[43] M. Dipoppa, M. Krupa, A. Torcini and B. S. Gutkin, Splay states in finite pulse-coupled networks of excitable neurons, SIAM J. Applied Dynamical Systems, 11, 864-894 (2012).

[44] N. Brunel and D. Hansel, How Noise Affects the Synchronization Properties of Recurrent Networks of Inhibitory Neurons, Neural Comp., 18, 1066-1110 (2006).

[45] A. Treves, Mean-field analysis of neuronal spike dynamics, Network: Computation in Neural Systems, 4, 259-284 (1993).

[46] M. Calamai, A. Politi and A. Torcini, Stability of splay states in globally coupled rotators, Phys. Rev. E 80, 036209 (2009).

[47] S. Olmi, R. Livi, A. Politi and A. Torcini, Collective oscillations in disordered neural networks, Phys. Rev. E, 81, 046119 (2010).

[48] S. Watanabe and S. H. Strogatz, Constants of motion for superconducting Josephson arrays, Physica D, 74, 197-253 (1994).

[49] A. Mauroy and R. Sepulchre, Clustering behaviors in networks of integrateand-fire oscillators, Chaos, 18, 037122 (2008).

[50] S. Olmi, A. Politi and A. Torcini, in preparation, 2012.

[51] M. A. Whittington and R. D. Traub, Interneuron Diversity series: Inhibitory interneurons and network oscillations in vitro , Trends in Neurosciences, 26, 676 (2003).

[52] C. Allene, A. Cattani, J. B. Ackman, P. Bonifazi, L. Aniksztejn, Y. Ben-Ari and R. Cossart, Sequential generation of two distinct synapse-driven network patterns in developing neocortex, The Journal of Neuroscience, 26, 12851-12863 (2008).

[53] Y. Ben-Ari, J. L. Gaiarsa, R. Tyzio, R. Khazipov, GABA: a pioneer transmitter that excites immature neurons and generates primitive oscillations, Phisiol. Rev., 87, 1215-1284 (2007).

[54] C. Koch, Biophysics of computation, Oxford University Press, New York (1999).

[55] W. Kinzel, On the stationary state of a network of inhibitory spiking neurons, J. Comput. Neurosci. 24, 105 (2008).

[56] I. Shimada and T. Nagashima, A Numerical Approach to Ergodic Problem of Dissipative Dynamical Systems, Prog. Theor. Phys. 61, 1605 (1979); G. Benettin, L. Galgani, A. Giorgilli and J.M. Strelcyn, Lyapunov Characteristic Exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part II: Numerical application, Meccanica, March 15 and 21 (1980).

[57] S. Coombes, Proceedings of Stochastic and Chaotic Dynamics in the Lakes, AIP Conference Proceedings, 502, 88 (1999)

[58] P. C. Müller, Calculation of Lyapunov Exponents for Dynamic Systems with Discontinuities, Chaos, Solitons & Fractals, 5, 1671-1681 (1995).

[59] C.H. Dellago, H. A. Posch, W. G. Hoover, Lyapunov instability in a system of hard disks in equilibrium and nonequilibrium steady states, Phys. Rev. E, 53, 1485 (1996).

[60] D. Zhou, Y. Sun, A. V. Rangan and D. Cai, Spectrum of Lyapunov exponents of non-smooth dynamical systems of integrate-and-fire type, J. Comput. Neurosci., 28, 229-245 (2010).

[61] O. V. Popovych, Y. L. Maistrenko and P. A. Tass, Phase chaos in coupled oscillators, Phys. Rev. E 71, 065201(R) (2005).

[62] T.P. Vogels, K. Rajan, and L. F. Abbott, Neural Network Dynamics, Annu. Rev. Neurosci. 28, 357 (2005).

[63] P. Bonifazi, M. Goldin, M. A. Picardo, I. Jorquera, A. Cattani, G. Bianconi, A. Represa, Y. Ben-Ari and R. Cossart, GABAergic Hub Neurons Orchestrate Synchrony in Developing Hippocampal Networks, Science, 326, 1419-1424 (2009).

[64] R. J. Morgan and I. Soltesz, Nonrandom connectivity of the epileptic dentate gyrus predicts a major role for neuronal hubs in seizures, Proceedings of the National Academy of Sciences USA, 105, 6179-6184 (2008).

[65] R. Albert and A. L. Barabási, Statistical mechanics of complex networks, Rev. Mod. Phys., 74, 4797 (2002).

[66] J. Friedrich, M. Kinzel, Dynamics of recurrent neural networks with delayed unreliable synapses: metastable clustering, J. Comput. Neurosci., 27, 65 (2009).

[67] A. T. Winfree, The Geometry of Biological Time, Springer Verlag, Berlin, (1980).

[68] I. Belykh, E. de Lange and M. Hasler, Synchronization of bursting neurons: what matters in the network topology, Phys. Rev. Lett., 94, 188101 (2005).

[69] M. London, A. Roth, L. Beeren, M. Häusser and E. Latham, Sensitivity to perturbations in vivo implies high noise and suggests rate coding in cortex, Nature, 446, 123-128 (2010).

[70] L. Tattini, S. Olmi and A. Torcini, Coherent periodic activity in excitatory Erdös-Reniy neural networks: the role of network connectivity , Chaos, 22, 023133 (2012).

[71] D. Golomb, D. Hansel, The Number of Synaptic Inputs and the Synchrony of Large, Sparse Neuronal Networks, Neural Comput., 12, 1095 (2000).

[72] K. A. Takeuchi, H. Chaté, F. Ginelli, A. Politi, and A. Torcini, Extensive and Subextensive Chaos in Globally Coupled Dynamical Systems, Phys. Rev. Lett., 107, 124101 (2011).

[73] T. Shibata, T. Chawanya, and K. Kaneko, Noiseless Collective Motion out of Noisy Chaos, Phys. Rev. Lett. 82, 4424, (1999).

[74] P. Grassberger, Information content and predictability of lumped and distributed dynamical systems, Physica Scripta 40, 346 (1989).

[75] M.R. Paul, M.I. Einarsson, P.F. Fischer, and M.C. Cross, Extensive chaos in Rayleigh-Bnard convection, Phys. Rev. E, 75, 045203(R) (2007).

[76] R Livi, A Politi, and S Ruffo, Distribution of characteristic exponents in the thermodynamic limit, J. Phys. A: Math. Gen., 19, 2033 (1986).

[77] K.A. Takeuchi, F. Ginelli, and H. Chaté, Lyapunov Analysis Captures the Collective Dynamics of Large Chaotic Systems, Phys. Rev. Lett., 103, 154103 (2009).

[78] M. Monteforte and F. Wolf, Dynamical Entropy Production in Spiking Neuron Networks in the Balanced State, Phys. Rev. Lett., 105, 268104 (2010).

[79] H. Chaté and P. Manneville, Evidence of Collective Behaviour in Cellular Automata, Europhys. Lett. 14, 409 (1991).

[80] D. J. Watts and S. H. Strogatz, Collective dynamics of ’small-world’ networks, Nature, 393, 440 (1998).

[81] W. Liao, Graph bipartitioning problem, Phys. Rev. Lett., 59, 1625 (1987).

[82] A. G. Percus, G. Istrate, B. Gonçalves, R. Z. Sumi and S. Boettcher, The peculiar phase structure of random graph bisection, J. Math. Phys., 49, 125219 (2008).

[83] S. Luccioli, S. Olmi, A. Politi and A. Torcini, Collective dynamics in sparse networks, Phys. Rev. Lett., 109, 138103 (2012).

[84] D. M. Abrams and S. H. Strogatz, Chimera States for Coupled Oscillators, Phys. Rev. Lett., 93, 174102 (2004).

[85] D. M. Abrams, R. E. Mirollo, S. H. Strogatz and D. A. Wiley, Solvable Model for Chimera States of Coupled Oscillators, Phys. Rev. Lett., 101, 084103 (2008).

[86] A. Pikovsky and M. Rosenblum, Partially Integrable Dynamics of Hierarchical Populations of Coupled Oscillators, Phys. Rev. Lett., 101, 264103 (2008).

[87] E. Aurell, G. Boffetta, A. Crisanti, G. Paladin and A. Vulpiani Growth of Noninfinitesimal Perturbations in Turbulence, Phys. Rev. Lett., 77, 1262 (1996).

[88] N. Nakagawa and Y. Kuramoto, Anomalous lyapunov spectrum in globally coupled oscillators, Physica D, 80, 307, (1995).

[89] D. Hansel, G. Mato and C. Meunier, Synchrony in Excitatory Neural Networks, Neural Comput., 7, 307 (1995).

[90] M. Rosenblum and A. Pikovsky, Self-Organized Quasiperiodicity in Oscillator Ensembles with Global Nonlinear Coupling, Phys. Rev. Lett., 98, 064101 (2007).

[91] M. Cencini, M. Falcioni, D. Vergni and A. Vulpiani, Macroscopic chaos in globally coupled maps, Physica D, 130, 58 (1999).

[92] S. Olmi, A. Politi and A. Torcini, Collective chaos in pulse-coupled neural networks, Europhysics Letters, 92, 60007 (2010).