Optimized brute-force algorithms for the bifurcation analysis of a binary neural network model

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 of neural networks has been developed mostly for mean-field limits of infinite-size spin-glass models, for finite-size dynamical systems whose units have a graded, continuous output, and for models with discrete-output neurons that evolve in continuous time. To allow progress on understanding the dynamics of some widely used classes of neural network models with discrete units and finite size, which could not be studied thoroughly with the previous methodology, here we introduced algorithms that perform a semianalytical bifurcation analysis of a finite-size firing-rate neural network model with binary firing rates and discrete-time evolution. In particular, we focus on the case of small networks composed of tens of neurons, to which existing statistical methods are not applicable. Our approach is based on a numerical brute-force search of the stationary and oscillatory solutions of the model, from which we derive analytical expressions of its bifurcation structure by means of the state-to-state transition probability matrix. Our algorithms determine how the network parameters affect the degree of multistability, the emergence and the period of the neural oscillations, and the formation of spontaneous 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.

Figure 1

Speed test of the algorithm for sparse networks. This figure shows the improvement of performance achieved by the algorithm for sparse networks described in Sec. 2d1, compared to the general nonoptimized algorithm of Sec. 2c1. We tested the algorithms by calculating the stationary firing-rate vectors of a sparse circulant network with topology $J=\overline{J}\mathrm{circ}\left(0,1,\cdots ,1,0,\cdots ,0\right)$. Each neuron has incoming vertex degree $M\in \left\{0,\cdots ,N-1\right\}$ ($M=0$ for a fully disconnected network, and $M=N-1$ for a fully connected one without loops), so that the density of the synaptic connections is proportional to $M$, according to the formula $\frac{\#\mathrm{synapses}}{N^2-N}=\frac{M}{N-1}\in \left[0,1\right]$. Moreover, we set the overall synaptic strength to $\overline{J}=10$ and the external stimuli to ${\mathfrak{I}}_0=\cdots ={\mathfrak{I}}_{N-1}=0$. For these values of the parameters, the network has two stationary firing-rate vectors, $\forall M,N$ considered in this figure. (a) Computational times $T_{\mathrm{sparse}}$ and $T_{\mathrm{non}-\mathrm{opt}}$ (see text) required for calculating the stationary vectors by means of an Intel^® Core™ i5-5300U CPU clocked at 2.30 GHz with 16 GB RAM. In particular, we observe that $T_{\mathrm{sparse}}\ll T_{\mathrm{non}-\mathrm{opt}}$ if $M$ is sufficiently smaller that $N-1$, while for small, highly dense networks the general nonoptimized algorithm may outperform the algorithm for sparse networks. Moreover, note that $T_{\mathrm{non}-\mathrm{opt}}$ does not depend on $M$. (b) Ratio $T_{\mathrm{non}-\mathrm{opt}}/T_{\mathrm{sparse}}$, namely the speed gain of the algorithm for sparse networks with respect to the nonoptimized one, as a function of $N$ and $M$.

Figure 2

An example of the progression of the efficient algorithm for a specific connectivity matrix. This figure shows the steps performed by the algorithm introduced in Sec. 2d1 during the derivation of the set $\mathcal{S}$, for a sparse network composed of five neurons. (a) Directed graph of the connectivity matrix. Each node represents a neuron in the network, while a black arrow from a node $j$ to a node $i$ represents a synaptic connection with weight $J_{ij}$. Each node $i$ receives an external input current ${\mathfrak{I}}_i$: In the specific case considered in this figure, the inputs to the neurons 0, 1, 2 are known (i.e., ${\mathfrak{I}}_0=-1, {\mathfrak{I}}_1=0, {\mathfrak{I}}_2=3$), so that the multistability diagram is calculated with respect to the unspecified inputs ${\mathfrak{I}}_3=I_3$ and ${\mathfrak{I}}_4=I_4$. (b) Progression of the efficient algorithm for sparse networks during the evaluation of the set of the candidate stationary vectors, $\mathcal{S}$.

Figure 3

Examples of bifurcation diagrams for small fully connected networks. This figure reports the bifurcation diagrams obtained for the network structure of Eq. (9), with $N=4$, 6, 8, and $J_{EE}=-J_{II}=80, J_{IE}=-J_{EI}=70, {\theta }_E={\theta }_I=1$. [(a)–(c)] Multistability diagrams. These panels show how the degree of multistability of the network, namely the number of stationary solutions, depends on the external currents $I_{E,I}$. Each color represents a different degree of multistability (e.g., blue = tristability). [(d)–(f)] Oscillation diagrams. These panels show how the number of oscillatory solutions and their period are affected by the stimuli. The notation $x:y$ reveals the presence of $y$ different oscillations with period $x$ in a given region of the diagram. Note that, unlike the multistability diagrams, generally there is no correspondence among the colors of the three oscillation diagrams. The script “Oscillation_Diagram.py” assigns the user-defined colors to each region of the oscillation diagrams regardless of the specific pairs $x:y$ that form that region (e.g., $\mathrm{red}=2:2$ for $N=4$, while $\mathrm{red}=3:1$ for $N=8$). Since the total number of possible regions of the oscillation diagrams is very large and the actual regions are not known a priori, this prevents the user from defining too many colors, most of which typically remain unused.

Figure 4

Examples of bifurcations in the state-to-state transition probability matrix. This figure shows the changes of dynamics that occur in the fully connected network of Eq. (9) for $N=4$, when varying the input current $I_E$ for $I_I=-30$. The nodes in these graphs are the (decimal representations of the) possible $2^N=16$ firing-rate vectors of the network (e.g., the node 3 corresponds to the firing rate $\mathbf{\nu }=0011$), while the arrows describe the transitions between the firing rates according to the state-to-state transition probability matrix [see Eq. (4)]. The figure highlights in red all the stationary and oscillatory solutions of the network (compare with the areas crossed by the dashed line in Figs. 3 and 3, when moving from left to right).

Figure 5

Examples of bifurcation diagrams for small sparse networks with random weights. This figure reports the bifurcation diagrams obtained for the network structure of Eq. (10), with $N=4$, 6, 8, the synaptic connectivity matrices in Table 1, and ${\theta }_E={\theta }_I=1$. Similar to Fig. 3, panels (a)–(c) and (d)–(f) show the multistability and oscillation diagrams, respectively, in the stimuli space for each network size.

Figure 6

Examples of bifurcation diagrams for sparse networks of tens of neurons. This figure shows the bifurcation diagrams obtained for two different realizations (left and right panels) of the network structure of Eq. (10), with $N=50$ and $p_{EE}=p_{EI}=p_{IE}=p_{II}=0.03$, while all the other parameters are set as in Fig. 5. (a) Synaptic connectivity matrices of the network realizations (the matrices are also reported in the supplemental files “Sparse_J_1.csv” and “Sparse_J_2.csv”) [58]. [(b), (c)] Multistability and oscillation diagrams obtained for ${\mathcal{T}}_{\max }=2$; see text.