Numerical Procedure and Convergence

All simulations were run in an event-based fashion following Ref. [51, 54, 58], where the exact map (2.3) was iterated from spike to spike in the phi representation of the theta neuron model with homogeneous coupling strengths and homogeneous external currents for all neurons in each population. To calculate the next spike time in the network it was sufficient to find the neuron with the largest phase in either population and then calculate its next spike time

t_s+1=t_s+min

i

π−φ_i(t_s) ω_i

, (2.48)

since the external currents of the neurons were identical in either population. In the case of two populations, the minimum of the two calculated next spike times would be the next spike time in the network. Then, all neurons’ phases were evolved until the next spike time using Eq. (2.35c) and the spike receiving neurons updated with Eq. (2.36c). These three steps compose one iteration and give numerically exact spike times and phases of the neurons.

With the exact phases of the neurons before spike reception, the single spike Jacobians (2.37) were evaluated using Eq. (2.38a). These were used to numerically calculate all Lyapunov ex-ponents in the standard procedure [68, 69]. After a warmup of the network dynamics, of typ-ically 100 spikes per neuron on average, we started with a random N-dimensional orthonor-mal system that was evolved in each iteration with the single spike Jacobian. After some itera-tions of about N/K spikes, the evolved vectors were reorthonormalized with the Gram-Schmidt-orthonormalization procedure, yielding the norms of the orthogonalized vectorsn_i(t_s)and the or-thonormal system to be used in the next iteration. After a short warmup of the oror-thonormal system of about one spike per neuron, these norms were used to calculate the N Lyapunov exponents λ_i=lim_p→∞t¹

p∑_s=1^p logn_i(t_s).

2.5 Networks of Theta neurons

0.01 0.1 1 10

t (s) -80

-40 0 40

λi (s-1 )

0.01 0.1 1 10

t (s) -80

-40 0 40

λi (s-1 )

λ₁ λ₂₀₀ λ₆₀₀ λ₁₂₀₀ λ₁₈₀₀ λ₂₀₀₀ C

B A

Figure 2.3–Convergence of Lyapunov spectra versus time in inhibitory networks. (logarithmic time scale) (A) Convergence of Lyapunov spectrum for one initial condition, (B) grey lines: some Lyapunov exponents for ten different initial phases, straight color lines: averages, dotted color lines: averages ± double standard errors, (C) as in (B) but for different network realizations in each run (parameters:N=2000, ν¯ =1 Hz,K=100,J₀=1,τm=10 ms).

An example code for MATLAB^®illustrating the principle steps to calculate the Lyapunov spec-tra is provided Appendix B. All full calculations were performed in custom code written in C++

with double precision. The GNU Scientific Library (GSL) was used for the random number gener-ator (Mersenne-Twister), the Automatically Tuned Linear Algebra Software (ATLAS) for matrix multiplications in the Gram–Schmidt procedure and the Message Passing Interface (MPI) for the parallel implementation of the simulations. The sparseness of the networks was efficiently used for the storage of the coupling matrices, the updates of the postsynaptic neurons and the matrix mul-tiplications of the orthonormal system with the sparse single spike Jacobians. For the reorthonor-malization, we chose a parallel recursive blocked version of the Gram–Schmidt procedure [91].

The Lyapunov spectra converged rather quickly over time to their asymptotic shape (Fig. 2.3).

One should note that the non-converged Lyapunov exponents are meaningless (they do not reflect the local or finite-time Lyapunov exponents). The converged Lyapunov exponents capture the asymptotic network dynamics. Figure 2.3A displays the convergence towards the full Lyapunov spectrum on logarithmic time scale. This calculation was repeated for different initial phases.

Figure 2.3B shows the results of ten such runs for six of the Lyapunov exponents (grey lines), together with their averagesλ_i=₁₀¹ ∑¹⁰_η=1λ_i,r (straight color lines) and confidence intervals (dotted color lines) of the double standard error 24λ_i=2

q1

10∑¹⁰_η=1(λi,r−λ_i)². Figure 2.3C shows the results of ten runs with different initial phases and different network realizations. The Lyapunov spectrum was independent of the initial phases as well as network realizations. Generally, all calculations of the Lyapunov spectra were repeated ten times with different initial phases and network realizations. Numerical errors were smaller than the symbol sizes in the presented figures.

The convergence of the largest Lyapunov exponents for different parameter sets in inhibitory and excitatory-inhibitory networks is depicted in Fig. 2.4 and Fig. 2.5, respectively. The evo-lution of the largest Lyapunov exponent was plotted versus the logarithm of the average number of spikes per neuron S (the number of all spikes in the network divided by the number of neu-rons). S=1 means that every neuron has spiked approximately once. A characteristic scale is S_D=log(N)/log(K), corresponding to the diameter of the random graph (the largest number of neurons on the shortest paths between any pair of neurons). In excitatory-inhibitory networks, the neurons have K inhibitory presynaptic neurons andK excitatory presynaptic neurons. Therefore we usedS_D=log(N)/log(2K)in the excitatory-inhibitory networks. S_Dis the approximate num-ber of spikes after which a perturbation of any neuron has influenced all others. This quantity is

0

Figure 2.4–Convergence of largest Lyapunov exponent in inhibitory networks versus average number of spikes per neuronS≈νt.¯ (logarithmic scale) (A)-(C) Different network sizesN(K=100, ¯ν=1 Hz), (D)-(F) different mean indegrees K (N =10 000, ¯ν=1 Hz), (G)-(I) different average firing rates ¯ν (N= 10 000, K=100), (grey lines: ten runs with different network realizations, black straight lines: averages of these ten runs, dotted lines: averages±double standard errors,S_D=log(N)/log(K), other parameters:

J₀=1,τ_m=10 ms).

Figure 2.5–Convergence of largest Lyapunov exponent in excitatory-inhibitory networks versus aver-age number of spikes per neuronS≈νt.¯ (logarithmic time scale) (A)-(C) Different excitatory-inhibitory feedback loop activationε (K=100, ¯ν=1 Hz), (D)-(F) different mean indegreesK (ε =0.3, ¯ν=1 Hz), (G)-(I) different average firing rates ¯ν (ε =0.3,K=100), (grey lines: ten runs with different network re-alizations, black straight lines: averages of these ten runs, dotted lines: averages±double standard errors, SD=log(N)/log(2K), other parameters:NE=8000,NI=2000,J₀=1,τm=10 ms).