Optimized brute-force algorithms for the bifurcation analysis of a spin-glass-like neural network model

TL;DR

This paper introduces optimized brute-force algorithms for semi-analytical bifurcation analysis of spin-glass-like neural networks with binary states, enabling detailed understanding of their multistability and oscillatory behavior.

Contribution

It presents novel algorithms for bifurcation analysis of discrete-time, binary neural networks, including an efficient method for sparse network connectivity.

Findings

Algorithms successfully identify bifurcation points and multistability regions.

Efficient analysis of sparse networks with arbitrary connectivity matrices.

Python implementation provided for practical application.

Abstract

Bifurcation theory is a powerful tool for studying how the dynamics of a neural network model depends on its underlying neurophysiological parameters. However, bifurcation theory has been developed mostly for smooth dynamical systems and for continuous-time non-smooth models, which prevents us from understanding the changes of dynamics in some widely used classes of artificial neural network models. This article is an attempt to fill this gap, through the introduction of algorithms that perform a semi-analytical bifurcation analysis of a spin-glass-like neural network model with binary firing rates and discrete-time evolution. Our approach is based on a numerical brute-force search of the stationary and oscillatory solutions of the spin-glass model, from which we derive analytical expressions of its bifurcation structure by means of the state-to-state transition probability matrix. The algorithms determine how the network parameters affect the degree of multistability, the emergence and the period of the neural oscillations, and the formation of symmetry-breaking in the neural populations. While this technique can be applied to networks with arbitrary (generally asymmetric) connectivity matrices, in particular we introduce a highly efficient algorithm for the bifurcation analysis of sparse networks. We also provide some examples of the obtained bifurcation diagrams and a Python implementation of the algorithms.

Figures

Tables

Table 1. Table 1: Synaptic connectivity matrices used to plot Fig. ( 5 ) . The non-zero entries of the matrix R α ​ β subscript 𝑅 𝛼 𝛽 R_{\alpha\beta} (see Eq. ( 9 )) are generated with probabilities p E ​ E = p I ​ E = 0.4 subscript 𝑝 𝐸 𝐸 subscript 𝑝 𝐼 𝐸 0.4 p_{EE}=p_{IE}=0.4 and p E ​ I = p I ​ I = 0.6 subscript 𝑝 𝐸 𝐼 subscript 𝑝 𝐼 𝐼 0.6 p_{EI}=p_{II}=0.6 , while the strengths of the synaptic connections are random integers generated from homogeneous and uniform distributions with J E ​ E min = − J I ​ I max = 80 superscript subscript 𝐽 𝐸 𝐸 min superscript subscript 𝐽 𝐼 𝐼 max 80 J_{EE}^{\mathrm{min}}=-J_{II}^{\mathrm{max}}=80 , J E ​ E max = − J I ​ I min = 100 superscript subscript 𝐽 𝐸 𝐸 max superscript subscript 𝐽 𝐼 𝐼 min 100 J_{EE}^{\mathrm{max}}=-J_{II}^{\mathrm{min}}=100 , J I ​ E min = − J E ​ I max = 30 superscript subscript 𝐽 𝐼 𝐸 min superscript subscript 𝐽 𝐸 𝐼 max 30 J_{IE}^{\mathrm{min}}=-J_{EI}^{\mathrm{max}}=30 and J I ​ E max = − J E ​ I min = 50 superscript subscript 𝐽 𝐼 𝐸 max superscript subscript 𝐽 𝐸 𝐼 min 50 J_{IE}^{\mathrm{max}}=-J_{EI}^{\mathrm{min}}=50 .

Equations

$$V_i\left(t+1\right)=\frac{1}{M_i}\sum _{j=0}^{N-1}{J}_{ij}\mathcal{H}\left(V_j\left(t\right)-{\theta }_j\right)+I_i,i=0,...,N-1.$$

$$\mathcal{H}\left(V-\theta \right)=\{\begin{array}{ll}0& \mathrm{if}V\le \theta \\ & \\ 1& \mathrm{otherwise},\end{array}$$

$$P\left(\mathbf{\nu }\left(t+1\right)|\mathbf{\nu }\left(t\right)\right)=\frac{1}{2^N}\prod _{j=0}^{N-1}\left[1+{\left(-1\right)}^{{\nu }_j\left(t+1\right)}\mathrm{sgn}\left({\theta }_j-\frac{1}{M_j}\sum _{k=0}^{N-1}{J}_{jk}{\nu }_k\left(t\right)-I_j\right)\right],$$

$$\mathrm{sgn}\left(x\right)=\{\begin{array}{ll}-1& \mathrm{if}x<0\\ & \\ 1& \mathrm{otherwise}.\end{array}$$

$$\frac{1}{2}\left[1+{\left(-1\right)}^{{\nu }_j}\mathrm{sgn}\left({\theta }_j-\frac{1}{M_j}\sum _{k=0}^{N-1}{J}_{jk}{\nu }_k-I_j\right)\right]=1,$$

$$\{\begin{array}{ll}{I}_j\le {\mathcal{I}}_j& \mathrm{if}{\nu }_j=0\\ & \\ I_j>{\mathcal{I}}_j& \mathrm{if}{\nu }_j=1.\end{array},{\mathcal{I}}_j\stackrel{\mathrm{def}}{=}{\theta }_j-\frac{1}{M_j}\sum _{k=0}^{N-1}{J}_{jk}{\nu }_k.$$

$$I\in \left({\Lambda }_I,{\Xi }_I\right],{\Lambda }_I\stackrel{\mathrm{def}}{=}\underset{j\in {\Gamma }_I:{\nu }_j=1}{\max }{\mathcal{I}}_j,{\Xi }_I\stackrel{\mathrm{def}}{=}\underset{j\in {\Gamma }_I:{\nu }_j=0}{\min }{\mathcal{I}}_j.$$

$$I\in \left({\Phi }_I,{\Psi }_I\right],{\Phi }_I\stackrel{\mathrm{def}}{=}\underset{t\in \mathcal{T}}{\max }\left(\underset{j\in {\Gamma }_I:{\nu }_j\left(t+1\right)=1}{\max }{\mathcal{I}}_j\left(t\right)\right),{\Psi }_I\stackrel{\mathrm{def}}{=}\underset{t\in \mathcal{T}}{\min }\left(\underset{j\in {\Gamma }_I:{\nu }_j\left(t+1\right)=0}{\min }{\mathcal{I}}_j\left(t\right)\right),$$

$$J=\left[\begin{array}[]{cc}\mathfrak{J}_{EE}&\mathfrak{J}_{EI}\\ \mathfrak{J}_{IE}&\mathfrak{J}_{II}\end{array}\right],\quad\boldsymbol{I}=\left[\begin{array}[]{c}\boldsymbol{I}_{E}\\ \boldsymbol{I}_{I}\end{array}\right].$$

$${\mathfrak{J}}_{\alpha \beta }=\{\begin{array}{ll}{J}_{\alpha \alpha }\left({\mathbb{I}}_{N_{\alpha }}-{\mathrm{Id}}_{N_{\alpha }}\right),& \mathrm{for}\alpha =\beta \\ & \\ J_{\alpha \beta }{\mathbb{I}}_{N_{\alpha },N_{\beta }},& \mathrm{for}\alpha \ne \beta \end{array},{\mathbf{I}}_{\alpha }=I_{\alpha }{1}_{N_{\alpha }}.$$

$$\mathfrak{J}_{\alpha\beta}=R_{\alpha\beta},\quad\quad\boldsymbol{I}_{\alpha}=\left[\begin{array}[]{c}\boldsymbol{0}_{N_{\alpha}-1}\\ I_{\alpha}\end{array}\right].$$

$$\left(k-1\right)!=$$

$$\sum _{k=0}^{\mathcal{N}}\left(\genfrac{}{}{0ex}{}{\mathcal{N}}{k}\right)x^k=$$

$$\sum _{k=2}^{\mathcal{N}}\left(\genfrac{}{}{0ex}{}{\mathcal{N}}{k}\right)\left(k-1\right)!=\sum _{k=2}^{\mathcal{N}}\left(\genfrac{}{}{0ex}{}{\mathcal{N}}{k}\right){\int }_0^{\infty }{x}^{k-1}{e}^{-x}dx={\int }_0^{\infty }\frac{\left(1+x\right){}^{\mathcal{N}}-\mathcal{N}x-1}{x}{e}^{-x}dx$$

$$\sim {\int }_0^{\infty }{x}^{\mathcal{N}-1}{e}^{-x}dx=\left(\mathcal{N}-1\right)!$$

Full text

**Optimized brute-force algorithms for the bifurcation analysis of a spin-glass-like neural network model

**

Diego Fasoli1,2,∗, Stefano Panzeri1

1 Laboratory of Neural Computation, Center for Neuroscience and Cognitive Systems @UniTn, Istituto Italiano di Tecnologia, 38068 Rovereto, Italy

2 Center for Brain and Cognition, Computational Neuroscience Group, Universitat Pompeu Fabra, 08002 Barcelona, Spain

$\ast$** Corresponding Author. E-mail: diego.fasoli@upf.edu**

Abstract

Bifurcation theory is a powerful tool for studying how the dynamics of a neural network model depends on its underlying neurophysiological parameters. However, bifurcation theory has been developed mostly for smooth dynamical systems and for continuous-time non-smooth models, which prevents us from understanding the changes of dynamics in some widely used classes of artificial neural network models. This article is an attempt to fill this gap, through the introduction of algorithms that perform a semi-analytical bifurcation analysis of a spin-glass-like neural network model with binary firing rates and discrete-time evolution. Our approach is based on a numerical brute-force search of the stationary and oscillatory solutions of the spin-glass model, from which we derive analytical expressions of its bifurcation structure by means of the state-to-state transition probability matrix. The algorithms determine how the network parameters affect the degree of multistability, the emergence and the period of the neural oscillations, and the formation of symmetry-breaking in the neural populations. While this technique can be applied to networks with arbitrary (generally asymmetric) connectivity matrices, in particular we introduce a highly efficient algorithm for the bifurcation analysis of sparse networks. We also provide some examples of the obtained bifurcation diagrams and a Python implementation of the algorithms.

1 Introduction

Neural complexity refers to the wide variety of dynamical behaviors that occur in neural networks [5, 27, 9]. This set of dynamical behaviors includes variations in the number of stable solutions of neuronal activity, the formation of neural oscillations, spontaneous symmetry-breaking, chaos and much more [16, 1]. Qualitative changes of neuronal activity, also known as bifurcations, are elicited by variations of the network parameters, such as the strength of the external input to the network, the strength of the synaptic connections between neurons, or other network characteristics.

Bifurcation theory is a standard mathematical formalism for studying neural complexity [20]. It allows the construction of a map of neuronal activity, known as bifurcation diagram, that links points or sets in the parameters space to their corresponding network dynamics. In the study of firing-rate network models, bifurcation theory has been applied mostly to graded (smooth) neural networks with analog firing rates (e.g. [4, 3, 5, 27, 15, 9]), proving itself as an effective tool for deepening our understanding of network dynamics. The bifurcation analysis of smooth models is based on differential analysis, in particular on the Jacobian matrix of the system. However, the Jacobian matrix is not defined everywhere for artificial neuronal models with discontinuous activation function, such as networks of binary neurons. For this reason, bifurcation theory of smooth dynamical systems cannot be applied directly to these models. On the other hand, while the bifurcation analysis of non-smooth dynamical systems has recently received increased attention, it has been developed only for continuous-time models, described by non-smooth differential equations or by differential inclusions [22, 2, 21, 24] (see [14] for an example of application to firing-rate models). Yet, many interesting problems in neuroscience involve the use of non-smooth discrete-time models [18], therefore this gap prevents us from understanding the changes of dynamics in these artificial systems.

If we define the states of a binary network of size $N$ as the collection of the rates at which its cells are firing, then the number of possible stationary states is $2^{N}$. Moreover, the number of possible oscillatory solutions is $\sim\left(2^{N}-1\right)!$ for large $N$, as we prove in Appx. (A). In particular, the bifurcation analysis of the model would unveil the actual solutions in these large sets of states, for any given set of network parameters. This would represent a tremendous simplification in our comprehension of neural dynamics, providing great insight into the operation of the network [23].

In this article, we develop brute-force algorithms for the bifurcation analysis of the spin-glass-like neural network model introduced in [10]. The network is composed of an arbitrary number of neurons with binary firing rates, connected through arbitrary (generally asymmetric) synaptic connections. The network dynamics is deterministic and evolves in discrete time steps. In [10] we proved that the bifurcation structure of the model can be studied semi-analytically, and in particular we were able to characterize the formation of multistability, neural oscillations and symmetry breaking in the stimuli space. While in [10] we derived the bifurcation diagrams of simple networks through hand calculations, in the present article we propose algorithms for automatic bifurcation analysis. Therefore our work complements standard numerical continuation softwares that are widely used for the bifurcation analysis of graded neuronal models, such as the MatCont Matlab toolbox [7] and XPPAUT [8].

Previous work focused on the emergence of complexity in ideal mean-field limits of spin-glass models, see e.g. [28, 6, 26]. On the contrary, in the present article we do not make use of any mean-field approximation. In particular, we consider exactly solvable finite-size networks to be used in real-world applications. In Sec. (2) we briefly review the model (SubSec. (2.1)) and we show how the analytical formula of the state-to-state transition probability matrix can be inverted for any network size $N$ to determine the portions of parameters space where multistability, neural oscillations and symmetry breaking occur (SubSec. (2.2)). The resulting formula can be written in any programming language. In particular, we propose a Python implementation, that the reader may find in the Supplementary Materials (scripts “Multistability_Diagram.py” and “Oscillation_Diagram.py”).

However, the inversion of the state-to-state transition probability matrix requires a priori knowledge of sets of candidate stationary and oscillatory solutions, as we showed in [10], and as we will explain in more detail in SubSec. (2.2). For this reason, in SubSec. (2.3) we introduce non-optimized algorithms that, through a brute-force searching procedure, generate the sets of all the potential stationary and oscillatory solutions (see the supplemental script “Non_Efficient_Algorithm.py” for a Python implementation). These algorithms can be applied to networks with any topology of the synaptic connections (dense or sparse), but their computational time increases exponentially with the network size.

Yet, the study of sparse systems is of particular interest in neuroscience. The average density of the synaptic connections, defined as the ratio between the actual and the maximum possible number of connections in the network, is approximately $10^{-6}-10^{-7}$ across the whole cortex, and it can increase up to $\sim 0.2-0.4$ in connection pathways linking cortical areas [19]. For this reason, in SubSec. (2.4) we propose an optimized brute-force algorithm, whose efficiency in generating sets of candidate stationary and oscillatory solutions increases with the network sparseness (see the supplemental script “Sparse_Efficient_Algorithm.py” for a Python implementation).

Then, in SubSec. (2.5) we introduce some widely-used examples of network topologies to be tested, while in Sec. (3) we show the corresponding bifurcation diagrams generated by our codes. To conclude, in Sec. (4) we discuss the importance and the biological implications of our results. In particular, we discuss how our work advances the comprehension of neural networks with respect to previous work (SubSec. (4.1)), new insights into the dynamics of binary network models revealed by our algorithms (SubSec. (4.2)), and future directions that need to be pursued to address the limitations of our work (SubSec. (4.3)).

2 Materials and Methods

In SubSec. (2.1) we describe the spin-glass-like neural network model whose dynamics we would like to investigate. Moreover, in SubSec. (2.2) we introduce a technique for plotting the bifurcation diagram of the model, provided a set of candidate stationary states and a set of candidate oscillatory solutions are known. Then, in SubSec. (2.3) we propose non-optimized brute-force algorithms for generating the candidate sets. While these algorithms can be applied to networks with any topology of the synaptic connections, their computational time increases exponentially with the network size. For this reason, in SubSec. (2.4) we introduce an efficient algorithm that takes advantage of the sparseness of biological networks for increasing its computational speed. To conclude, in SubSec. (2.5) we introduce two standard examples of network topologies, whose bifurcation structure will be derived in Sec. (3) through the use of our algorithms.

2.1 The Network Model

We study the bifurcation structure of the following spin-glass-like network model [10]:

[TABLE]

In Eq. (1), $N$ is the number of neurons in the network, while $V_{i}$ and $I_{i}$ are the membrane potential and the external stimulus of the $i$th neuron, respectively. $J_{ij}$ is the synaptic weight from the $j$th (presynaptic) to the $i$th (postsynaptic) neuron. The collection of the synaptic weights, for $i,j=0,\ldots,N-1$, defines the synaptic connectivity matrix $J$, which in this article is supposed to be arbitrary (generally asymmetric). In Eq. (1), $\mathscr{H}\left(\cdot\right)$ represents the Heaviside step function:

[TABLE]

where $\theta$ is the firing threshold. The firing rate of the $i$th neuron is the binary variable defined as $\nu_{i}\overset{\mathrm{def}}{=}\mathscr{H}\left(V_{i}-\theta_{i}\right)\in\left\{0,1\right\}$, so that $\nu_{i}=0$ if the neuron is not firing, and $\nu_{i}=1$ if it is firing at the maximum rate. Moreover, in Eq. (1) the parameter $M_{i}$ represents the number of presynaptic neurons that are directly connected to the $i$th (postsynaptic) neuron. Therefore $M_{i}$ is a normalization factor, that prevents the divergence of the sum $\sum_{j=0}^{N-1}J_{ij}\mathscr{H}\left(V_{j}\left(t\right)-\theta_{j}\right)$ for $N\gg 1$.

The state-to-state transition probability matrix $\mathcal{P}$ provides a convenient way to describe the dynamics of the firing rates $\nu_{i}$. By defining $\boldsymbol{\nu}\left(t\right)\overset{\mathrm{def}}{=}\nu_{0}\left(t\right)\nu_{1}\left(t\right)\ldots\nu_{N-1}\left(t\right)$ as the binary string obtained by concatenating the firing rates at time $t$, in [10] we proved that the entries of the $2^{N}\times 2^{N}$ matrix $\mathcal{P}$ are:

[TABLE]

where $\mathrm{sgn}\left(\cdot\right)$ is the sign function, defined as follows:

[TABLE]

$P\left(\boldsymbol{\nu}\left(t+1\right)|\boldsymbol{\nu}\left(t\right)\right)\in\left\{0,1\right\}$ represents the probability of the transition $\boldsymbol{\nu}\left(t\right)\rightarrow\boldsymbol{\nu}\left(t+1\right)$ to occur, for specific firing-rate states $\boldsymbol{\nu}\left(t\right)$ and $\boldsymbol{\nu}\left(t+1\right)$, given the network is in the state $\boldsymbol{\nu}\left(t\right)$ at the time instant $t$. For example, if at the time instant $t$ the network is in the state $\boldsymbol{\nu}\left(t\right)=00\cdots 0$ and we want to check if the neurons flip their firing rates at the next time step (so that $\boldsymbol{\nu}\left(t+1\right)=11\cdots 1$), then $P\left(\boldsymbol{\nu}\left(t+1\right)|\boldsymbol{\nu}\left(t\right)\right)=1$ if the transition $00\cdots 0\rightarrow 11\cdots 1$ occurs from $t$ to $t+1$, and $P\left(\boldsymbol{\nu}\left(t+1\right)|\boldsymbol{\nu}\left(t\right)\right)=0$ otherwise. In SubSec. (2.2) we will show how to use Eq. (3) for plotting the bifurcation diagram of the network.

2.2 Plotting the Bifurcation Diagrams

In this subsection we introduce a general way for plotting the bifurcation diagrams of networks with arbitrary topology of the synaptic connections. For a given network, its bifurcation diagram is composed of two panels, namely the multistability and the oscillation diagrams. These diagrams provide a complete picture of the relation between the stationary/oscillatory solutions of the network and the set of stimuli (see for example the left and right panels of Figs. (3) and (5)). In SubSecs. (2.2.1) and (2.2.2) we describe algorithms for calculating the multistability and the oscillation diagrams, which are implemented in the supplemental Python scripts “Multistability_Diagram.py” and “Oscillation_Diagram.py” respectively.

Finally, in the case of networks with homogeneous populations, we superimpose to these diagrams the regions of the stimuli space where the system undergoes spontaneous intra-population symmetry-breaking. In these regions, the stationary and oscillatory solutions calculated by our algorithms show non-homogeneous firing rates within one or more populations, despite the homogeneity of their neurophysiological parameters. More generally, given a neuron $i$ in population $\alpha$, our algorithms detect the formation of spontaneous symmetry-breaking in networks where all the terms $\frac{1}{M_{i}}\sum_{j\in\beta}J_{ij}$ and $I_{i}$ depend only on the populations $\alpha,\>\beta$, since this is the most general symmetry condition for the network, according to Eq. (1). Therefore, networks with homogeneous parameters $M_{i},\>J_{ij},\>I_{i}$ represent only a special case of this condition.

2.2.1 Multistability Diagram

The firing-rate state $\boldsymbol{\nu}$ is stationary if it satisfies the condition $P\left(\boldsymbol{\nu}|\boldsymbol{\nu}\right)=1$. From Eq. (3) we observe that this condition holds if:

[TABLE]

for $j=0,\ldots,N-1$. By inverting Eq. (4) with respect to the stimulus $I_{j}$, we get that the equation is satisfied whenever:

[TABLE]

Generally, if some neurons share the same stimulus, so that for example the neurons with indexes in the set $\Gamma_{I}$ receive the same external current $I$, then from Eq. (5) we get:

[TABLE]

Eq. (6) provides the ranges of the stimuli (if any) where the firing-rate state $\boldsymbol{\nu}$ is stationary. By calculating the ranges $\left(\Lambda_{I},\Xi_{I}\right]$ for every set $\Gamma_{I}$, we get a complete picture of the relation between the stationary states and the set of stimuli (see for example the left panels of Figs. (3) and (5)). If the ranges corresponding to $\mathcal{M}$ different states $\boldsymbol{\nu}$ overlap, the overlapping area has multistability degree $\mathcal{M}$. In other words, for combinations of stimuli lying in that area, the network has $\mathcal{M}$ stationary firing rates. Therefore Eq. (6) allows the analytical derivation of the multistability diagram, extending the analysis we performed in [10] to networks with arbitrary topology and size.

Note that Eq. (6) implies the calculation of the set $\left(\Lambda_{I},\Xi_{I}\right]$ for each (potentially) stationary firing-rate state $\boldsymbol{\nu}$. By defining $\mathscr{S}$ as the set of the candidate stationary states to be checked, we observe that $\boldsymbol{\nu}\in\mathscr{S}$ is actually stationary if $\Lambda_{I}<\Xi_{I}$. In SubSecs. (2.3) and (2.4) we will introduce two different algorithms for generating $\mathscr{S}$, which are called by the Python script “Multistability_Diagram.py”. Depending on the algorithm, the cardinality $\left|\mathscr{S}\right|$ may be different for a given network topology. This will affect the overall speed with which the multistability diagram is plotted.

2.2.2 Oscillation Diagram

In principle, the same approach discussed in the previous subsection for the multistability diagram of the network can be extended to the study of neural oscillations. A sequence $\mathcal{O}$, defined as $\boldsymbol{\nu}\left(0\right)\rightarrow\boldsymbol{\nu}\left(1\right)\rightarrow\ldots\rightarrow\boldsymbol{\nu}\left(\mathcal{T}\right)$ with $\boldsymbol{\nu}\left(0\right)=\boldsymbol{\nu}\left(\mathcal{T}\right)$, is an actual oscillatory solution (of period $\mathcal{T}\in\left\{2,\ldots,2^{N}\right\}$) if every transition $\boldsymbol{\nu}\left(t\right)\rightarrow\boldsymbol{\nu}\left(t+1\right)$ in $\mathcal{O}$ satisfies the condition $P\left(\boldsymbol{\nu}\left(t+1\right)|\boldsymbol{\nu}\left(t\right)\right)=1$. Similarly to SubSec. (2.2.1), if some neurons share the same stimulus, so that the neurons with indexes in the set $\Gamma_{I}$ receive the same external current $I$, then the sequence $\mathcal{O}$ is an oscillatory solution of the network dynamics if:

[TABLE]

where $\mathscr{T}\overset{\mathrm{def}}{=}\left\{0,\ldots,\mathcal{T}-1\right\}$ and $\mathcal{I}_{j}\left(t\right)\overset{\mathrm{def}}{=}\theta_{j}-\frac{1}{M_{j}}\sum_{k=0}^{N-1}J_{jk}\nu_{k}\left(t\right)$. Finally, by calculating the ranges $\left(\Phi_{I},\Psi_{I}\right]$ for every set $\Gamma_{I}$, we get a complete analytical picture of the relation between the oscillatory solutions and the set of stimuli (see for example the right panels of Figs. (3) and (5)).

Similarly to the case of the multistability diagram, we observe that Eq. (7) implies the calculation of the set $\left(\Phi_{I},\Psi_{I}\right]$ for each (potential) oscillatory solution. In SubSecs. (2.3) and (2.4) we will introduce two different algorithms for generating the set $\mathscr{O}$ of the candidate oscillatory solutions, which are called by the Python script “Oscillation_Diagram.py”.

2.3 Non-Efficient Algorithms for Generating the Stationary and Oscillatory

Solutions

In SubSec. (2.2) we described an algorithm for plotting the multistability and oscillation diagrams, provided some sets of candidate stationary and oscillatory solutions, $\mathscr{S}$ and $\mathscr{O}$ respectively, are known. In the present section we introduce non-optimized brute-force algorithms, implemented in the Python script “Non_Efficient_Algorithm.py”, that generate the sets $\mathscr{S}$ and $\mathscr{O}$ for networks with arbitrary topology.

2.3.1 Generation of the Set $\mathscr{S}$

The simplest way for generating the set $\mathscr{S}$ is to fill it with all the $2^{N}$ states of the firing rates, from $00\cdots 0$ to $11\cdots 1$. This approach allows the evaluation of the multistability diagram of small-size networks, since its computational time increases exponentially with $N$. We observe that only the states in $\mathscr{S}$ that satisfy the constraint $\Lambda_{I}<\Xi_{I}$ (see SubSec. (2.2)) are actually stationary, therefore the set $\mathscr{S}$ generated by this algorithm is usually oversized. In SubSec. (4.3) we will discuss a way to reduce $\left|\mathscr{S}\right|$ for specific network topologies.

2.3.2 Generation of the Set $\mathscr{O}$

Unfortunately, the brute-force approach described above for the generation of the set $\mathscr{S}$ is unfeasible for the study of neural oscillations, due to computational time. Indeed, the number of possible oscillations in a network of size $N$ is $n_{\mathcal{O}}=\sum_{k=2}^{2^{N}}\binom{2^{N}}{k}\left(k-1\right)!$, which grows as $\sim\left(2^{N}-1\right)!$ for $N\rightarrow\infty$ (see Appx. (A)). The best solution we found to reduce the computational time, given a network with arbitrary topology, is to obtain $\mathscr{O}$ by discretizing the stimuli space, and then by solving iteratively Eq. (1) for every combination of stimuli and for every initial condition of the firing rates. In other words, we discretized the stimuli space through a grid $G$ composed of $n_{G}$ points, each one representing a combination of stimuli. Then we solved iteratively the equation $V_{i}\left(t+1\right)=\frac{1}{M_{i}}\sum_{j=0}^{N-1}J_{ij}\nu_{j}\left(t\right)+I_{i}$, with $\nu_{j}\left(t\right)=\mathscr{H}\left(V_{j}\left(t\right)-\theta_{i}\right)$, for each of the $n_{G}$ combinations of stimuli of the grid $G$, and for each of the $2^{N}$ initial conditions $\boldsymbol{\nu}\left(0\right)$. If the firing-rate vector $\boldsymbol{\nu}\left(t\right)$ calculated through Eq. (1) oscillates according to a sequence $\mathcal{O}$ for a given initial condition $\boldsymbol{\nu}\left(0\right)$, then $\mathcal{O}\in\mathscr{O}$. Finally, we filtered the set $\mathcal{O}$ in order to remove duplicate oscillations, so that the set will contain every oscillation exactly once (for example, in a $2$-neurons circuit, the oscillations $01\rightarrow 10\rightarrow 11\rightarrow 01$ and $10\rightarrow 11\rightarrow 01\rightarrow 10$ are circularly identical, therefore one must be discarded). Through this approach, we can derive the actual oscillations in the set $\mathscr{O}$, by solving Eq. (1) $2^{N}n_{G}$ times. Generally $2^{N}n_{G},\,\left|\mathscr{O}\right|\ll n_{\mathcal{O}}$, so that if we calculate the range $\left(\Phi_{I},\Psi_{I}\right]$ through Eq. (7) for every oscillation in $\mathscr{O}$, the computational time required for deriving the oscillation diagram now increases as $2^{N}$ with the network size, rather than $\left(2^{N}-1\right)!$.

However, while being much faster, this algorithm does not guarantee that the resulting oscillation diagram is complete. Indeed, for some oscillations, the range $\left(\Phi_{I},\Psi_{I}\right]$ could be smaller than the grid resolution. For this reason, the parameter $n_{G}$ must be chosen accurately in order to avoid relevant information loss (which may occur when $n_{G}$ is too small for a given network) and an excessive computational load of the algorithm (which may occur when $n_{G}$ is too large for the computational power available).

To conclude, we observe that the oscillation diagram may be calculated in a purely numerical way, by solving Eq. (1) for every combination of stimuli on the grid, without making use of the analytical formula (7). However, by relying on Eq. (7), our semi-analytical approach allows the analytical derivation of the oscillation diagram, which would not be possible otherwise. Moreover, it is easy to verify that the use of Eq. (7) allows a reduction of the resolution of the grid $G$ compared to a purely numerical approach, resulting in a much faster derivation of the oscillation diagram.

2.4 An Efficient Algorithm for Generating the Stationary and Oscillatory

Solutions in Sparse Networks

We observe that the algorithms described in SubSec. (2.3) can be applied to networks with any topology of the synaptic connections (dense or sparse), but their computational time increases as $2^{N}$ with the network size. In this subsection we show that, for the specific case of sparse networks, it is possible to take advantage of the sparseness of the synaptic connections to define an efficient algorithm for the generation of the sets $\mathscr{S}$ and $\mathscr{O}$. This approach, which we implemented in the Python script “Sparse_Efficient_Algorithm.py”, may outperform of several orders of magnitude the non-optimized approaches introduced in SubSec. (2.3).

2.4.1 Generation of the Set $\mathscr{S}$

As we saw in SubSec. (2.2.1), in order to identify the stationary states of the network we need to check if the condition (4) is satisfied for each neuronal index $j=0,\ldots,N-1$, given a set of stimuli $I_{0},\ldots,I_{N-1}$. If the $j$th neuron receives inputs from all the other neurons in the network (i.e. if $J_{jk}\neq 0$ for $k=0,\ldots,N-1$), the condition (4) must be checked for all the $2^{N}$ binary states $\boldsymbol{\nu}$ of length $N$, from $00\cdots 0$ to $11\cdots 1$. On the contrary, in sparse networks some synaptic weights are equal to zero. If some neurons do not project a synaptic connection to the $j$th neuron (i.e. if $J_{jk}=0$ for some index $k$), the corresponding firing rates $\nu_{k}$ will not affect the sum $\sum_{k=0}^{N-1}J_{jk}\nu_{k}$. Therefore for sparse networks there is no need to check Eq. (4) for all the binary states of length $N$, unlike the case of dense networks. This observation allows us to introduce an efficient algorithm for finding the candidate stationary states of sparse neural networks, which is described below.

At step [math], we set $j=0$ and we call $P_{0}$ the set of neurons with indexes $k\neq 0$ that do not project a synaptic connection to the [math]th neuron (i.e. $k\in P_{0}$ if $J_{0k}=0$ and $k\neq 0$). Moreover, we call $Q_{0}$ the set of the remaining $N-\left|P_{0}\right|$ neurons. This is the set of neurons that can affect the condition (4) through their firing rates. In particular, the [math]th neuron belongs to $Q_{0}$ since it can affect the condition (4) through the term $\left(-1\right)^{\nu_{0}}$ even in the case when $J_{00}=0$ (therefore $\left|Q_{0}\right|=M_{0}$ if $J_{00}\neq 0$, and $\left|Q_{0}\right|=M_{0}+1$ if $J_{00}=0$, where $M_{0}$ is the incoming vertex degree of the [math]th neuron, see Eq. (1)). At this step, we only need to check Eq. (4) for all the binary states $s_{0}$ of length $\left|Q_{0}\right|$. If no state $s_{0}$ satisfies Eq. (4), the network has no stationary solution and the algorithm is stopped. Otherwise we call $\mathscr{S}$ the set of states $s_{0}$ that satisfy Eq. (4), and we switch to the neuronal index $j=1$, as described below.

At step $1$, we set $j=1$ and we call $P_{1}$ the set of neurons with indexes $k\neq 1$ that do not project to the $1$st neuron. Moreover we call $Q_{1}$ the set of the remaining $N-\left|P_{1}\right|$ neurons, and $R_{1}=Q_{1}\cap Q_{0}$. Then we generate all the binary states $s_{1}$ of length $\left|Q_{1}\right|-\left|R_{1}\right|$, and we use each of them to complete the states in $\mathscr{S}$, by creating new binary states of length $\left|Q_{0}\right|+\left|Q_{1}\right|-\left|R_{1}\right|$. For example, we suppose that at step [math] we got $Q_{0}=\left\{0,\>7,\>8\right\}$ and $\mathscr{S}=\left\{010,\>111\right\}$, while at step $1$ we got $Q_{1}=\left\{1,\>2\right\}$. Then the states $s_{1}=00,$ $01$, $10$, $11$ will be used to fill the states in $\mathscr{S}$ according to the index $k$, generating the new states $0\boldsymbol{0}\boldsymbol{0}10$, $0\boldsymbol{0}\boldsymbol{1}10$, $0\boldsymbol{1}\boldsymbol{0}10$, $0\boldsymbol{1}\boldsymbol{1}10$ and $1\boldsymbol{0}\boldsymbol{0}11$, $1\boldsymbol{0}\boldsymbol{1}11$, $1\boldsymbol{1}\boldsymbol{0}11$, $1\boldsymbol{1}\boldsymbol{1}11$ (the bits of the states $s_{1}$ are highlighted in bold). Then the algorithm checks the condition (4) for $j=1$ on the newly generated states (when $Q_{1}$ is empty, the algorithm checks the condition directly on the states in $\mathscr{S}$). If no state satisfies Eq. (4), then the network has no stationary solution, therefore the set $\mathscr{S}$ is cleared and the algorithm is stopped. Otherwise $\mathscr{S}$ is cleared and filled with the newly generated states that satisfy Eq. (4). Then we proceed to the neuronal index $j=2$.

In a similar way, for $2\leq j\leq N-1$, we define $P_{j}$ as the set of neurons with indexes $k\neq j$ which do not project to the $j$th neuron. Moreover we call $Q_{j}$ the set of the remaining $N-\left|P_{j}\right|$ neurons, and $R_{j}=Q_{j}\cap\left(\bigcup_{n=0}^{j-1}Q_{n}\right)$. The procedure described at step $1$ is repeated iteratively through the steps $2$, $3$, … by completing the states in $\mathscr{S}$ through the binary strings $s_{j}$ of length $\left|Q_{j}\right|-\left|R_{j}\right|$. If the algorithm does not stop, when the $j$th step has been completed the states in $\mathscr{S}$ have length $L_{j}=\sum_{n=0}^{j}\left|Q_{n}\right|-\sum_{n=1}^{j}\left|R_{n}\right|$. We observe that $L_{j}$ quickly tends to $N$ for increasing $j$. For example, if the network has a ring topology (i.e. if each neuron receives a connection only from the previous neuron and projects a connection only to the next one), we get $\left|Q_{n}\right|=2$ for $n=0,\ldots,N-1$, $\left|R_{n}\right|=1$ for $n=1,\ldots,N-2$ and $\left|R_{N-1}\right|=2$, so that $L_{j}\rightarrow N$ for $j\rightarrow N-1$.

At the end of the process, $\mathscr{S}$ will contain the actual stationary states of the network (if any), namely all the states $\boldsymbol{\nu}$ of length $N$ that satisfy the condition $P\left(\boldsymbol{\nu}|\boldsymbol{\nu}\right)=1$, if the stimuli $I_{0},\ldots,I_{N-1}$ are all known. The efficiency of the algorithm is directly proportional to the sparseness of the matrix $J$ and inversely proportional to the number of stationary states of the network. Moreover, further speed-up is achieved by sorting the neuronal indexes with ascending vertex degrees $M_{j}$, before running the algorithm. This is due to the fact that the algorithm slows down when both $\left|Q_{j}\right|-\left|R_{j}\right|$ and $\left|\mathscr{S}\right|$ are large. However, we observe that if at some step $\overline{j}$ a string of length $L_{\overline{j}}$ does not satisfy the condition (4), then during the next steps the algorithm will not check all the $2^{N-L_{\overline{j}}}$ binary states of length $N$ that contain that string. Since $L_{\overline{j}}\approx\sum_{n=0}^{\overline{j}}M_{n}-\sum_{n=1}^{\overline{j}}\left|R_{n}\right|$, if the vertex degrees $M_{n}$ are small $\forall n\leq\overline{j}$, then the string length $L_{\overline{j}}$ must be small. For this reason, the number of (non-stationary) states containing the string of length $L_{\overline{j}}$ is large ($2^{N-L_{\overline{j}}}$). Being non-stationary, this large number of states will not fill the set $\mathscr{S}$, which therefore remains small. When the algorithm will proceed by iterating over the indexes $j$ with large degree $M_{j}$ (which are computationally more expensive, since the algorithm has to iterate over all the binary states of length $\left|Q_{j}\right|-\left|R_{j}\right|\approx M_{j}-\left|R_{j}\right|$), the number of states in $\mathscr{S}$ will be small. This reduces the computational load in the slowest part of the algorithm, resulting in an overall speed-up. On the contrary, for the same reason the algorithm slows down when the neurons are sorted with descending $M_{j}$.

Fig. (1) shows a comparison between the computational time $T_{\mathrm{sparse}}$ required by our algorithm for calculating the stationary states of a given network topology for a fixed set of stimuli $I_{0},\ldots,I_{N-1}$, and the computational time $T_{\mathrm{non-opt}}$ required for performing the same calculation through the non-optimized brute-force approach introduced in SubSec. (2.3.1).

The figure shows that the speed gain $\frac{T_{\mathrm{non-opt}}}{T_{\mathrm{sparse}}}$ of the algorithm for sparse networks with respect to the non-optimized algorithm is generally much larger than $1$, and it increases with the network size and sparseness. The computational time of the non-optimized algorithm increases exponentially as $2^{N}$ regardless of the network sparseness. However, for small and highly dense networks this algorithm could outperform that for sparse networks ($\frac{T_{\mathrm{non-opt}}}{T_{\mathrm{sparse}}}<1$).

As we said, when the currents $I_{0},\ldots,I_{N-1}$ are all known, the algorithm for sparse networks generates a set $\mathscr{S}$ containing the actual stationary states of the network. However, since our purpose is to plot the multistability diagram, we observe that the values of the currents that define the parameter space of the diagram are generally not fixed (see e.g. Fig. (2), where the currents $I_{3}$ and $I_{4}$ are unspecified since they represent the bifurcation parameters of the network). For this reason, in order to derive the multistability diagram, the algorithm described above must be adapted for generating a set $\mathscr{S}$ whenever some inputs are unspecified. In what follows we propose two different solutions.

The first is to apply our algorithm for all the combinations of the bifurcation stimuli on a grid in the stimuli space, similarly to the method described in SubSec. (2.3.2) for the study of oscillations. This approach works for any set of stimuli (e.g. when all the neurons share the same stimulus or when each neuron receives a distinct external input), but it does not guarantee that the resulting multistability diagram is complete if the resolution of the grid is not high enough.

Unlike the first method described above, the second way to use our algorithm (which is the one we implemented in the Python script “Sparse_Efficient_Algorithm.py”) is specific to systems where the external currents that we vary during the bifurcation analysis are injected into a limited number of neurons. In other words, we suppose that the stimuli that span the bifurcation diagram are $I_{x},\ldots,I_{N-1}$, where $x$ is an integer close or equal to $N-1$ (see e.g. the network considered in Eq. (9), where $x=N-2$ after a proper rearrangement of the neural indexes). We also suppose that the remaining stimuli, i.e. $I_{0},\ldots,I_{x-1}$, are known and fixed, so that the algorithm for sparse networks can be applied for $j=0,\ldots,x-1$. The algorithm will generate a set $\mathscr{S}$ containing incomplete stationary states, namely binary strings of length $L_{x-1}$. Now, the $2^{y}$ binary states of length $y\overset{\mathrm{def}}{=}N-L_{x-1}$ from $00\cdots 0$ to $11\cdots 1$ must be used to complete the states in $\mathscr{S}$, creating $K=2^{y}\left|\mathscr{S}\right|$ binary strings of length $N$. For example, in the case of the $5$-neurons network shown in Fig. (2), we suppose we got $\mathscr{S}=\left\{1001,\>1101\right\}$ after the step $x-1$ (with $x=3$) has been completed, and that the states in $\mathscr{S}$ are the firing rates of the neurons with indexes $\left\{0,\>1,\>2,\>4\right\}$. Then $y=1$, therefore the complete binary strings are $\boldsymbol{\nu}=100\boldsymbol{0}1,\>100\boldsymbol{1}1,\>110\boldsymbol{0}1,\>110\boldsymbol{1}1$ (the bit corresponding to $y$ is highlighted in bold).

These are candidate stationary states, which can be used to derive the multistability diagram by inverting the relation $P\left(\boldsymbol{\nu}|\boldsymbol{\nu}\right)=1$ in the currents $I_{x},\ldots,I_{N-1}$. Note that the relation between $x$ and $y$ depends on the topology of the network. However, if $x$ is close or equal to $N-1$, then $y$ is generally small, because most of the binary digits in every stationary firing-rate state are determined, through our algorithm, by the known currents $I_{0},\ldots,I_{x-1}$. In this case, the overall efficiency of the algorithm will be high for large $N$, since $K$ does not diverge exponentially with the network size.

2.4.2 Generation of the Set $\mathscr{O}$

The algorithm introduced in SubSec. (2.4.1) can be naturally extended to find the oscillatory solutions of the network, as we describe below. In SubSec. (2.2.2) we observed that the sequence $\boldsymbol{\nu}\left(0\right)\rightarrow\boldsymbol{\nu}\left(1\right)\rightarrow\ldots\rightarrow\boldsymbol{\nu}\left(\mathcal{T}\right)$, with $\boldsymbol{\nu}\left(0\right)=\boldsymbol{\nu}\left(\mathcal{T}\right)$, is an oscillatory solution if every transition $\boldsymbol{\nu}\left(t\right)\rightarrow\boldsymbol{\nu}\left(t+1\right)$ satisfies the condition $P\left(\boldsymbol{\nu}\left(t+1\right)|\boldsymbol{\nu}\left(t\right)\right)=1$. For this reason, the initially empty set $\mathscr{O}$ of the oscillatory solutions with fixed period $\mathcal{T}$ is populated as follows. First, at the $j$th step of the algorithm (for $j=0,\ldots,N-1$), we need to generate all the binary states of length $\left|Q_{j}\right|-\left|R_{j}\right|$ (where $Q_{j}$ and $R_{j}$ are defined as in SubSec. (2.4.1), and $R_{0}\overset{\mathrm{def}}{=}\left\{\right\}$) for each time instant $t=0,\ldots,\mathcal{T}-1$. In other words, the algorithm creates all the $2^{\mathcal{T}\left(\left|Q_{j}\right|-\left|R_{j}\right|\right)}$ $\mathcal{T}$-tuples of states of length $\left|Q_{j}\right|-\left|R_{j}\right|$. Then, similarly to SubSec. (2.4.1), we have to use these states for completing those in the set $\mathscr{O}$, and to check if the states in $\mathscr{O}$ satisfy the condition $\frac{1}{2}\left[1+\left(-1\right)^{\nu_{j}\left(t+1\right)}\mathrm{sgn}\left(\theta_{j}-\frac{1}{M_{j}}\sum_{k=0}^{N-1}J_{jk}\nu_{k}\left(t\right)-I_{j}\right)\right]=1$ for $t=0,\ldots,\mathcal{T}-1$. At the end of the $j$th step, the set $\mathscr{O}$, if not empty, is cleared and filled with the newly generated states that satisfy the condition, before proceeding to the $\left(j+1\right)$th step of the algorithm. The algorithm is stopped only if $\mathscr{O}$ is empty or after the $\left(N-1\right)$th step is completed. Then, similarly to the method of SubSec. (2.3.2), we discard the duplicate oscillations in the set $\mathscr{O}$ . At the end of the process, if the currents $I_{0},\ldots,I_{N-1}$ are all known, the set will contain exactly once all the actual oscillatory solutions of the given network topology with fixed period $\mathcal{T}$ . Whenever some external stimuli are unspecified, the algorithm is adapted as we discussed in SubSec. (2.4.1) for the stationary states.

We observe that this approach, which we implemented in the script “Sparse_Efficient_Algorithm.py”, detects only the oscillations with a fixed period $\mathcal{T}$. Therefore this procedure must be repeated for $\mathcal{T}=2,\ldots,\mathcal{T}_{\mathrm{max}}$, where the maximum period $\mathcal{T}_{\mathrm{max}}$ is a user-defined parameter, that should be chosen according to the computational power available. Generally, for $\mathcal{T}_{\mathrm{max}}<2^{N}$ the algorithm is not guaranteed to find all the oscillatory solutions. However, setting $\mathcal{T}_{\mathrm{max}}$ close to or larger than $Z\overset{\mathrm{def}}{=}\frac{N}{\max_{j}\left(M_{j}-\left|R_{j}\right|\right)}$ results in an overall (at least exponential with $N$) slowing down of the algorithm, due to the generation of the $2^{\mathcal{T}\left(\left|Q_{j}\right|-\left|R_{j}\right|\right)}\approx 2^{\mathcal{T}\left(M_{j}-\left|R_{j}\right|\right)}$ $\mathcal{T}$-tuples. Therefore, similarly to SubSec. (2.3.2), the study of oscillations proves to be computationally expensive, requiring a compromise between precision and computational speed, which in this case is generally achieved by setting $2\ll\mathcal{T}_{\mathrm{max}}\ll Z$ for large $Z$. We observe that $Z$ increases with the network sparseness for a fixed network size $N$, therefore in sparse networks the algorithm will detect more oscillation periods in a fixed amount of time, compared to dense networks. Unfortunately, it is not possible to know a priori the effective maximum period of the oscillations generated by a network. However, in the sparse topologies we analyzed, this was usually small, and oscillations with period equal to or larger than the network size appear much more rarely (see e.g. the right panels of Fig. (5), where random sparse networks of size $N=4,\,6,\,8$ showed only oscillations with period $2$ and $4$). Typically, oscillation periods larger than $N$ occur in systems with built-on-purpose, and usually very regular, topologies.

To conclude, we observe that more generally this technique can be used to calculate efficiently the whole state-to-state transition probability matrix $\mathcal{P}$. According to Eq. (3), this is defined as the $2^{N}\times 2^{N}$ binary matrix of the conditional probabilities $P\left(\cdot|\cdot\right)$ of the firing rates. More precisely, $\mathcal{P}_{ij}=P\left(\boldsymbol{b}_{N}\left(i\right)|\boldsymbol{b}_{N}\left(j\right)\right)\in\left\{0,1\right\}$, for $i,\>j=0,\ldots,2^{N}-1$, where the firing rate $\boldsymbol{b}_{N}\left(k\right)$ is the binary representation of length $N$ of the index $k$ (e.g., if $N=5$, then $\boldsymbol{b}_{N}\left(14\right)=01110$). Our algorithm provides all the firing-rate pairs $\left(\boldsymbol{b}_{N}\left(i\right),\boldsymbol{b}_{N}\left(j\right)\right)$ such that $P\left(\boldsymbol{b}_{N}\left(i\right)|\boldsymbol{b}_{N}\left(j\right)\right)=1$, from which we can calculate the corresponding (decimal) indexes $\left(i,j\right)$. In other words, the algorithm calculates $\mathcal{P}$ efficiently and in a convenient sparse-matrix notation, since it provides only the coordinates $\left(i,j\right)$ of the non-zero entries of the matrix. The matrix $\mathcal{P}$ describes all the possible transitions between the firing rates, which are not restricted to stationary and oscillatory solutions only. Eventually, the stationary and oscillatory solutions of the network may be calculated from $\mathcal{P}$ through a cycle-finding algorithm (in particular, the stationary solutions may be considered as oscillations with period $\mathcal{T}=1$). To the best of our knowledge, the fastest cycle-finding algorithm was introduced by Johnson [17], and is implemented in the function simple_cycles of the Python library NetworkX. However, this approach proved to be even slower than the non-optimized algorithm of SubSec. (2.3), therefore we will not consider it here any further.

2.5 Examples of Network Topologies

In this section we report two standard examples of network topologies, whose bifurcation diagrams will be studied in Sec. (3): fully-connected networks (SubSec. (2.5.1)) and sparse networks with random synaptic weights (SubSec. (2.5.1)). In these examples, we suppose that the networks are composed of one excitatory ($E$) and one inhibitory ($I$) population. Two-population networks are commonly considered a good approximation of a single neural mass [12], however our analysis may be extended to systems composed of an arbitrary number of populations, if desired. We define $N_{E}$ ($N_{I}$) to be the size of the excitatory (inhibitory) population, with $N=N_{E}+N_{I}$. Note that $N_{E,I}$ can be arbitrary, but for illustrative purposes in this section we consider the case $N_{E}=N_{I}$ (rather than the $N_{E}\backslash N_{I}=4$ ratio experimentally observed in biological systems [25]), because we found that the network complexity increases with the size of the inhibitory population. Interestingly, the same phenomenon was found to occur also in multi-population networks with graded activation function [9]. Moreover, without further loss of generality, we index the neurons of the excitatory population as $i=0,\ldots,N_{E}-1$ and the inhibitory neurons as $i=N_{E},\ldots,N-1$, so that the synaptic connectivity matrix $J$ and the stimulus vector $\boldsymbol{I}$ can be written as follows:

[TABLE]

$\mathfrak{J}_{\alpha\beta}$, for $\alpha,\beta=E,I$, is a $N_{\alpha}\times N_{\beta}$ matrix that describes the synaptic connections from the population $\beta$ to the population $\alpha$. The entries of the matrices $\mathfrak{J}_{EE}$ and $\mathfrak{J}_{IE}$ must be non-negative, while those of the matrices $\mathfrak{J}_{II}$ and $\mathfrak{J}_{EI}$ must be non-positive, since the populations $E$ and $I$ are composed of excitatory and inhibitory neurons, respectively. Moreover, self-connections are not present in biological networks, so that the main diagonals of the matrices $\mathfrak{J}_{EE}$ and $\mathfrak{J}_{II}$ should be set to zero (even though, more generally, our algorithms could be applied also to networks with self-connections). In a similar way, $\boldsymbol{I}_{\alpha}$ represents the collection of stimuli to the population $\alpha$. Note that $\mathfrak{J}_{\alpha\beta}$ and $\boldsymbol{I}_{\alpha}$ depend on the structure of the specific network we study, as we will show below. Finally, we will call $\theta_{\alpha}$ the (homogeneous) firing threshold of all the neurons in the population $\alpha$.

2.5.1 Fully-Connected Networks

Our first example is a fully-connected network with homogeneous intra-population inputs, whose parameters are:

[TABLE]

The parameters $J_{\alpha\beta}$ describe the strength of the homogeneous synaptic connections from the population $\beta$ to the population $\alpha$, while $I_{\alpha}$ represents the stimulus current to each neuron of the population $\alpha$. Moreover, $\mathbb{I}_{N_{\alpha},N_{\beta}}$ is the $N_{\alpha}\times N_{\beta}$ all-ones matrix (here we use the simplified notation $\mathbb{I}_{N_{\alpha}}\overset{\mathrm{def}}{=}\mathbb{I}_{N_{\alpha},N_{\alpha}}$), while $\mathrm{Id}_{N_{\alpha}}$ is the $N_{\alpha}\times N_{\alpha}$ identity matrix and $\boldsymbol{1}_{N_{\alpha}}$ is the $N_{\alpha}\times 1$ all-ones vector.

2.5.2 Sparse Random Networks

The second example we consider is a sparse network with random synaptic weights, whose parameters are:

[TABLE]

$R_{\alpha\beta}$ is a $N_{\alpha}\times N_{\beta}$ sparse random matrix, whose entry $\left[R_{\alpha\beta}\right]_{ij}$, either for $\alpha=\beta$ and $i\neq j$, or for $\alpha\neq\beta$ and $\forall i,\>j$, is equal to a non-zero random number with probability $p_{\alpha\beta}$, while it is equal to zero with probability $1-p_{\alpha\beta}$. For example, we suppose that the non-zero entries of $R_{\alpha\beta}$ are generated from a homogeneous and uniform (i.e. rectangular) distribution with support $\left[J_{\alpha\beta}^{\mathrm{min}},J_{\alpha\beta}^{\mathrm{max}}\right]$, where the parameters $J_{\alpha\beta}^{\mathrm{min}}$ and $J_{\alpha\beta}^{\mathrm{max}}$ describe the minimum and maximum strength of the synaptic connections, respectively. Moreover, in Eq. (9), $I_{\alpha}$ represents the stimulus current to one neuron of the population $\alpha$, while $\boldsymbol{0}_{N_{\alpha}-1}$ is the $\left(N_{\alpha}-1\right)\times 1$ all-zeros vector (so that the remaining neurons of the population do not receive any external input).

3 Results

In this section we report the bifurcation diagrams generated by our algorithms for the network topologies described in SubSec. (2.5). In Fig. (3) we show the (codimension two) bifurcation diagrams of the fully-connected network (see SubSec. (2.5.1)) in the $I_{E}-I_{I}$ plane, for $N=4,\,6,\,8$. The left and right panels have been obtained from Eqs. (6) and (7) respectively, for $\Gamma_{I_{E}}=\left\{0,\ldots,N_{E}-1\right\}$ and $\Gamma_{I_{I}}=\left\{N_{E},\ldots,N-1\right\}$.

In particular, we observe that the case $N=6$ is the same considered in [10], and whose bifurcation diagram was derived through hand calculations. Due to the homogeneity of the synaptic connections, the bifurcation diagrams show symmetric structures, characterized (also for $N>8$, results not shown) by an increase of the maximum degree of multistability with the network size and by low-period oscillations. Interestingly, while spontaneous symmetry-breaking may occur in both the populations during neural oscillations, only the inhibitory one may undergo the formation of heterogeneous activity during the stationary states. As we discuss in SubSec. (4.3), for any homogeneous multi-population network it is generally possible to take advantage of this phenomenon, in order to speed up the non-optimized algorithm of SubSec. (2.3.1).

Moreover, in Fig. (4) we show some examples of changes of dynamics that occur in the state-to-state transition probability matrix of the fully connected network for $N=4$. This figure shows that the graph of the matrix, as given by Eq. (3), changes its structure as a function of the stimuli. In particular, the graph undergoes changes in the degree of multistability and in the period of the oscillations, which are more conveniently described by Fig. (3) (top panels).

To conclude, Fig. (5) shows the bifurcation diagrams of the sparse random networks (see SubSec. (2.5.2)) in the $I_{E}-I_{I}$ plane, which have been obtained for $N=4,\,6,\,8$ and $\Gamma_{I_{E}}=\left\{N_{E}-1\right\},\,\Gamma_{I_{I}}=\left\{N-1\right\}$. As expected, the randomness of the synaptic connections is reflected in the irregular structure of both the multistability and the oscillation diagrams, as opposed to the regular structure of the fully-connected networks shown in Fig. (3). In the random networks, high multistability degrees and high-period oscillations can occur only by chance, depending on the randomness of the synaptic weights and on the sparseness of the synaptic topology.

4 Discussion

We studied how the dynamics of the spin-glass-like neural network model (1) depends on its underlying parameters. The network is composed of a finite and arbitrary number of neurons with binary firing rates, that interact through arbitrary (generally asymmetric) synaptic connections. While in artificial neural networks with graded (smooth) activation functions this study is performed by means of standard bifurcation theory [20], the discontinuity of the activation function of the spin-glass model at the firing threshold prevents us from using this powerful tool.

Due to the discrete (binary) nature of the network, the number of possible stationary states and neural oscillations is finite. For this reason, in SubSec. (2.3) we developed brute-force algorithms for studying the bifurcations of the model. By taking advantage of the state-to-state transition probability matrix, these algorithms find the actual stationary and oscillatory solutions of the network, and build the corresponding bifurcation diagram in the parameters space. In particular, in Sec. (3) we provided examples that show how the network dynamics depends on the external stimuli, the network size and the topology of the synaptic connections. Similarly to the case of graded firing-rate network models [4, 3, 5, 27, 15, 9], this analysis revealed a complex bifurcation structure, encompassing several changes in the degree of multistability of the network, oscillations with stimulus-dependent frequency, and various forms of spontaneous symmetry-breaking of the neural activity.

The computational time of these algorithms increases as $2^{N}$ with the network size, regardless of the network topology. In particular, the study of neural oscillations proved very challenging due to combinatorial explosion, causing us to find a compromise between speed and precision. The best solution we found was to discretize the parameter space of the oscillation diagram and to run a searching algorithm for every combination of the parameters on the grid. While this solution allowed us to perform the study in exponential time, it may generate incomplete oscillation diagrams, due to the finite grid resolution.

Since biological networks are sparse [19], in SubSec. (2.4) we introduced an efficient algorithm that takes advantage of the limited number of synaptic connections. This algorithm may outperform of several orders of magnitude the non-optimized approaches introduced in SubSec. (2.3). While the computational time of the non-optimized algorithms increases as $2^{N}$ regardless of the network topology, the speed of the efficient algorithm increases with the network sparseness and is inversely proportional to the number of stationary/oscillatory solutions. For these reasons, the algorithms introduced in SubSec. (2.4) prove to be particularly convenient when applied to large sparse networks. However, for highly dense networks their speed decreases significantly, therefore in this case the algorithms of SubSec. (2.3) should be used instead.

While the study of multistability did not reveal any complication, the study of neural oscillations proved again to be computationally demanding, also for sparse networks. Our algorithm does not rely on the discretization of the parameter space, therefore it can detect areas at any resolution. Nevertheless, it outperforms the non-optimized algorithm of SubSec. (2.3) only in detecting oscillations with period $\mathcal{T}\ll Z\overset{\mathrm{def}}{=}\frac{N}{\max_{j}\left(M_{j}-\left|R_{j}\right|\right)}$. Since $Z$ increases with the network sparseness for a fixed network size $N$, in sparse networks the algorithm will detect more oscillation periods in a fixed amount of time, compared to dense networks. While generally it is not possible to know a priori the effective maximum period of the oscillations generated by a network (an exception is represented by symmetric networks, which can sustain only oscillations with period $\mathcal{T}=2$, see [11]), in the sparse random topologies we analyzed, this period was usually small (see e.g. the right panels of Fig. (5), where we showed that the networks corresponding to the topologies of Tab. (1) undergo only oscillations with period $2$ and $4$). Oscillations with period equal to or larger than the network size occur much more rarely by chance, usually in very small regions of the parameter space (a formal explanation of this phenomenon is reported in SubSec. (4.2)). Nevertheless, special (usually very regular) topologies may show large oscillation periods, therefore a complete analysis of these networks is beyond the capability of our algorithm for large $N$. Moreover, while for most of the sparse random topologies we analyzed the algorithm proved to be very fast, for some networks the script “Sparse_Efficient_Algorithm.py” generated large sets of candidate oscillatory solutions (only few of which were actual neural oscillations), resulting in a considerable slowing down of the bifurcation analysis. This phenomenon occurs in those networks where the activity of the neurons that receive the fixed external stimuli does not provide enough information for determining which oscillations are allowed and which are not. In other words, depending on the topology of the synaptic connections, in these special networks the number of oscillatory solutions was strongly influenced by both the fixed stimuli and the free stimuli that define the parameter space of the oscillation diagram. Interestingly, this phenomenon occurred despite the free currents were injected only into a limited number of neurons by hypothesis (see SubSec. (2.4)). In these cases, a grid implementation of the algorithm (similar to that used in the script “Non_Efficient_Algorithm.py”) could be considerably faster, since at every point of the grid in the parameter space the algorithm for sparse networks would generate only actual oscillatory solutions (see SubSec. (2.4)).

In the following we discuss the advances of our results with respect to previous work (SubSec. (4.1)), the implications of our work to better understand neural network dynamics (SubSec. (4.2)), and future directions that need to be pursued to address the limitations of our results (SubSec. (4.3)).

4.1 Progress with Respect to Previous Work on Bifurcation Analysis

Previous studies on the complexity of dynamics in firing-rate network models focused either on ideal mean-field limits of spin-glass-like models with discontinuous firing rates (e.g. [28, 6, 26]), or on systems with graded (smooth) firing rates (e.g. [4, 3, 5, 27, 15, 9]). On the contrary, in this article we studied finite-size networks with discontinuous activation functions. The analysis of these models is typically more complex, since we can rely neither on mean-field approximations, nor on the powerful methods developed in bifurcation theory for smooth systems. Moreover, our model evolves in discrete-time steps according to the recurrence relation (1). This also prevents us from using the methods recently developed for the bifurcation analysis of non-smooth dynamical systems, which can be applied only to continuous-time models described by non-smooth differential equations or by differential inclusions [22, 2, 21, 24]. This proves that the model we considered in this article is particularly hard to tackle with standard mathematical techniques.

However, in [10] we showed that it is possible to derive the bifurcation structure of this model by combining a brute-force search of the stationary and oscillatory solutions of the network, with the analytical formula of the state-to-state transition probability matrix. Our proof of concept was limited to the specific case of a fully-connected network, whose bifurcation structure was calculated through tedious hand calculations. In the present article we developed fast algorithms that perform this analysis automatically, for any topology of the synaptic connections. Therefore our work complements standard numerical continuation softwares such as the MatCont Matlab toolbox [7] and XPPAUT [8], that are widely used in neuroscience for the bifurcation analysis of graded neuronal models (see e.g. [13, 29, 30, 9]).

4.2 New Insights into the Dynamics of Discrete Network Models

In his pioneering work [23], Little underlined the importance of identifying the states that characterize the long-term behavior of a neural network, which would provide a great simplification in our comprehension of the network dynamics. In particular, given a deterministic network model, this problem is equivalent to identifying the actual stationary and oscillatory solutions for $t\rightarrow\infty$. While the actual stationary states must be identified out of the set of $2^{N}$ possible solutions in a network of size $N$, the problem of selecting the actual oscillations is even more formidable, since in large asymmetric networks the number of possible oscillatory solutions is $\sim\left(2^{N}-1\right)!$ (see Appx. (A)).

In our work, we developed algorithms that identify the actual stationary and oscillatory solutions of Eq. (1). The algorithms are particularly efficient when applied to sparse networks, providing great insight into the operation of the network. Moreover, the algorithms make use of this knowledge to deepen further our understanding of the network dynamics, by calculating the corresponding bifurcation structure. Interestingly, our algorithms revealed the formation of complex bifurcation diagrams also in small networks. The diagrams strongly depend on the network size, $N$, and on the topology of the synaptic connections, $J$. However, a detailed analysis of the network dynamics as a function of $N$ and $J$ is computationally very expensive (due to the large number of possible topologies for each fixed network size) and beyond the purpose of this article. Nevertheless, our study revealed that the spin-glass-like network model (1) can undergo three different kinds of bifurcations.

The first is a change in the degree of multistability of the network. This phenomenon occurs through the formation and destruction of stationary solutions, similarly to the limit-point (also known as saddle-node) bifurcations that occur in graded models [20].

The second kind of bifurcation is the formation of neural oscillations. If this phenomenon occurs through the annihilation of a stationary state (see e.g. the formation of the period-$4$ oscillation in the middle panels of Fig. (3)), then it may be interpreted as the discrete counterpart of the Andronov-Hopf bifurcation of graded models [20]. Otherwise, the number of oscillatory solutions may change with no variation in the number of stationary states (see e.g. the transition between the blue and red areas in the top-right panel of Fig. (3)), similarly to the limit point of cycles bifurcation of graded models [20]. It is interesting to observe that usually large-period oscillations occur more rarely and in smaller regions of the parameter space (for example, compare the green area of the period-$4$ oscillation in the middle-right panel of Fig. (3), with the blue and red areas of the period-$3$ and period-$2$ oscillations). The reason of this phenomenon can be easily deduced from Eq. (7): the larger the period of the oscillation, the more probable is that the function $\Phi_{I}=\underset{t\in\mathscr{T}}{\max}\left(\cdot\right)$ (respectively, $\Psi_{I}=\underset{t\in\mathscr{T}}{\min}\left(\cdot\right)$) will be large (respectively, small). In turn, this implies that the range $\left(\Phi_{I},\Psi_{I}\right]$ where the large-period oscillation occurs will be narrower, and in particular if $\Phi_{I}\geq\Psi_{I}$ the oscillation will not occur at all.

The third kind of bifurcation is characterized by spontaneous symmetry-breaking. This phenomenon is similar to the branching-point bifurcations [20] that occur in homogeneous multi-population graded models (see e.g. [9]). Here we observe the formation of heterogeneous intra-population neural activity, even if Eq. (1) does not contain any term that breaks explicitly the network symmetry. In particular, the symmetry may be broken at the stationary states or during oscillations, in both excitatory and inhibitory neural populations, depending on the network topology and parameters.

To conclude, we observe that our algorithms also calculate the number of stationary and oscillatory solutions, as well as the number of areas with distinct degrees of multistability, the number of oscillations with a given period, etc. We propose that these measures could be used to quantify rigorously the complexity of the neural dynamics in relation to the network size.

4.3 Future Directions

Neural networks with discrete activation functions may be considered as caricatures of the corresponding graded models. By discretizing the activation function, we gain the possibility to solve exactly the network equations [10], but currently it is not known to which extent this approximation oversimplifies the dynamical behavior of the system. For this reason, in future work we will investigate a potential reduction of the underlying neural complexity of the network by comparing the two classes of models.

Moreover, note that while in this article we introduced new algorithms for the bifurcation analysis of Eq. (1), we did not perform an exhaustive analysis of the relationship between the neural dynamics and the network parameters. In Sec. (3) we showed some relevant examples of the codimension two bifurcation diagrams for specific networks sizes and topologies, while in future work we will systematically investigate, through a large scale analysis, to which extent the network parameters correlate with the neural complexity.

To conclude, we observe that while we introduced fast algorithms for the bifurcation analysis of sparse networks, efficient solutions for dense networks are still missing. Ad hoc solutions can be developed for specific network topologies. For example, in the case of multi-population fully-connected networks with homogeneous weights, such as that considered in Eq. (8), it is possible to prove that the symmetry of the excitatory populations is never broken when the network is in a stationary state. The same phenomenon occurs in homogeneous multi-population networks with graded activation function [9]. Accordingly, the method described in SubSec. (2.3.1) can be easily modified in order to include in the set $\mathscr{S}$ only the states with homogeneous excitatory firing rates. In the case of the two-populations network described by Eq. (8), this approach would generate a set with cardinality $\left|\mathscr{S}\right|=2^{1+N_{I}}$, rather than $\left|\mathscr{S}\right|=2^{N_{E}+N_{I}}$, reducing considerably the computational time for $N_{E}\gg 1$. However, this idea can be applied only to specific network topologies, and a general method for speeding up the bifurcation analysis of dense networks with arbitrary topology is still missing. This represents another challenge that needs to be addressed in future work, and that will complete the tools at our disposal for studying the dynamical behavior of binary neuronal network models.

Appendix A: Asymptotic Estimation of the Number of Possible Oscillatory

Solutions in Large Networks

In this appendix we want to determine how the number of possible oscillatory solutions in an asymmetric network of size $N$ grows for $N\rightarrow\infty$. Given any finite $N$, we define $\mathcal{N}\overset{\mathrm{def}}{=}2^{N}$ and we observe that $\binom{\mathcal{N}}{k}$ is the total number of sets (known as $k$-combinations) containing $k$ distinct neural states among which the network may oscillate. For example, for $k=3$, the oscillations $0\rightarrow 1\rightarrow 2\rightarrow 0$ and $0\rightarrow 1\rightarrow 3\rightarrow 0$ occur among the states in the $3$-combinations $\left\{0,1,2\right\}$ and $\left\{0,1,3\right\}$, respectively. Now we observe that any given set of $k$ states gives rise to $\left(k-1\right)!$ distinct oscillations. For example, $0\rightarrow 1\rightarrow 2\rightarrow 0$ and $0\rightarrow 2\rightarrow 1\rightarrow 0$ are two distinct oscillations generated by the set $\left\{0,1,2\right\}$, while the oscillations $0\rightarrow 1\rightarrow 2\rightarrow 0$ and $1\rightarrow 2\rightarrow 0\rightarrow 1$ are (circularly) identical. Therefore, by summing over all the oscillation periods $k$ from $2$ to $\mathcal{N}$ (the case $k=1$ corresponds to the stationary states and therefore is not considered), we get that the total number of distinct oscillations is $\sum_{k=2}^{\mathcal{N}}\binom{\mathcal{N}}{k}\left(k-1\right)!$ (see also [17]). Then, by observing that:

[TABLE]

we finally get the following asymptotic expansion:

[TABLE]

in the limit $\mathcal{N}\rightarrow\infty.\;\;\blacksquare$

Acknowledgments

We thank Anna Cattani for useful feedback. This research was supported by the Autonomous Province of Trento, Call “Grandi Progetti 2012,” project “Characterizing and improving brain mechanisms of attention - ATTEND”.

The funders had no role in study design, data collection and analysis, decision to publish, interpretation of results, or preparation of the manuscript.