Mathematical studies of the dynamics of finite-size binary neural networks: A review of recent progress

1.   Introduction

Understanding the dynamics of networks of neurons, and of how such networks represent, process and exchange information by means of the temporal evolution of their activity, is one of the central problems in neuroscience. Real networks of neurons are highly complex both in terms of structure and physiology. Introducing details of this complexity greatly complicates the tractability of the models. Thus, a mathematical model of neural networks needs to be carefully designed, by finding a compromise between the elements of biological complexity and plausibility that are introduced, and the analytical tractability of the resulting model ^[1].

A wide set of mathematical models have been proposed to investigate the behavior of biological neural networks ^[2,3,4]. Typically, these models attempt to simplify as much as possible the original system they describe, without losing the properties that give rise to the most interesting emergent phenomena observed in biological systems. Binary neural network models ^[5,6,7,8,9] represent one of the most successful examples in finding a good compromise between keeping simplicity to enhance tractability, while yet achieving a rich dynamics with a set of complex emergent network properties.

Binary models describe the dynamical properties of networks composed of threshold units, which integrate their inputs to produce a binary output, namely a high (respectively low) output firing rate when their membrane potential does (respectively does not) exceed a given threshold (see section (2) for more details). Among the wide set of neural network models proposed by computational neuroscientists, threshold units represent one of the most convenient tools for studying the dynamical and statistical properties of neural circuits. The relative ease with which these models can be investigated analytically, is a consequence of their thresholding activation function, which can be considered as the simplest, piecewise-constant, approximation of the non-linear (and typically sigmoidal shaped ^[10,11]) graded input-output relationship of biological neurons. Despite their simplicity, as shown both by classic work ^{[5,6,7,12,13]}, as well as by our work reviewed here ^[14,15,16], the jump discontinuity of their activation function at the threshold is sufficient to endow binary networks with a complex set of useful emergent dynamical properties and non-linear phenomena, such as attractor dynamics ^[17], formation of patterns and oscillatory waves ^[18], chaos ^[19], and information processing capabilities ^[20], which are reminiscent of neuronal activity in biological networks.

The importance of binary network models is further strengthened by their close relationship with spin networks studied in physics ^{[21,22,23,24]}. The temporal evolution of a binary network in the zero-noise limit is isomorphic to the dynamics of kinetic models of spin networks at absolute temperature ^[25]. This allowed computational neuroscientists to study the behavior of large-size binary networks, by applying the powerful techniques of statistical mechanics already developed for spin models (see e.g. ^[26]).

Sizes of brains and of specialized neural networks within brains change considerably across animal species, ranging from few tens of neurons in invertebrates such as rotifers and nematodes, to billions of neurons in cetaceans and primates ^[27]. Network size changes also across levels of spatial organizations, ranging from microscopic and mesoscopic levels of organization in cortical micro-columns and columns (including from few tens ^[28] to few tens of thousands of neurons ^[29,30]), to several orders of magnitude more in large macroscopic networks, such as the resting state networks in the human brain, that involve many brain areas ^[31,32]. For this reason, it is important to be able to study mathematically the dynamics of binary neural network models (or of any network model, see e.g. ^[33,34,35]), for a wide range of different network sizes.

Large-scale networks composed of several thousands of neurons or more, are typically studied by taking advantage of the powerful techniques of statistical mechanics, such as the law of large numbers and the central limit theorem, under the assumption of random connectivity matrices and/or noise-driven dynamics (see e.g. ^{[26,13,36,14,16]}). These theories, such as those developed in physics for spin models, typically approximate the interaction of all the other neurons to a given neuron with a mean field, namely an effective interaction which is obtained by averaging the interactions of the other neurons. This allows one to dimensionally reduce the model, by transforming the set of equations of a large network into a single-neuron equation. Therefore the mean-field theories represent a powerful tool for gaining insight into the behavior of large networks, at a relatively low cost.

However, statistical mechanics does not apply to small-scale networks containing only few tens of neurons, which therefore prove much more difficult to study mathematically. Importantly, several studies, see e.g. ^{[37,36,38,35,14,16]}, have shown that the dynamics of neural network models in the large-size limit can be qualitatively, not only quantitatively, very different from that of the same models with small or finite size. For this reason, the computational investigation of networks composed of a finite number of neurons requires the development of new specialized analytical and numerical techniques. Some notable examples of the techniques that have been recently introduced for studying finite-size networks include linear response theory ^{[39,40,41,42]} and mesoscopic population equations ^[43] for networks of integrate-and-fire neurons, the cumulant expansion of binary networks ^[44], and much more (see Subsection (6.2) for a deeper discussion).

Recently, we proposed an alternative theory to describe the dynamics of finite-size neural networks composed of binary threshold units ^[14,15,16], which we critically review in this paper. In section (2) we introduce the binary network model that we analyze in this review, while in section (3) we focus specifically on the zero-noise limit of small networks, and we characterize mathematically the bifurcation points of the network dynamics in single network realizations. In section (4) we describe the techniques that we developed for studying small networks when external sources of noise are added to the neural equations. In particular, in Subsection (4.1) we describe the dualism between neural activity and membrane potentials, and we derive a complete description of the probabilistic behavior of the network in the long-time limit. In Subsections (4.2) and (4.3), we introduce, under some assumptions on the nature of the noise sources, exact analytical expressions of the cross-neuron correlations, and a learning rule for storing patterns and sequences of neural activity. Then, in section (5), we extend the study of bifurcations of section (3) to networks with quenched disorder across multiple network realizations. To conclude, in section (6) we discuss the advantages and weaknesses of our techniques, also in relation to previous analytical work on finite-size networks. Moreover, we highlight key challenges that remain open, as well as future directions for further progress in the mathematical study of binary networks, and possible implications of our results for neuroscience.

2.   The binary network model

In this review we assume that neural activity evolves in discrete time steps, and that the threshold units are synchronously updated. These assumptions are often used in studying network dynamics, see e.g. ^[5,6,12,2]. A threshold unit, here referred to as artificial neuron, is a logic gate or a mathematical function that mimics the working mechanisms of a biological neuron. Typically the unit receives several inputs, which can be loosely interpreted as postsynaptic potentials at neural dendrites, and sums them to produce a binary or digital output (also known as activation or neural activity). Usually, each input of a threshold unit is multiplied by a so called synaptic weight, which represents the strength of connections between pairs of neurons. Moreover, the sum of the weighted inputs is passed through a piecewise-constant (or Heaviside) thresholding function, also known as activation function or transfer function. If the sum of the weighted inputs exceeds a threshold, the output is set to one, and the artificial neuron is said to fire at that rate. On the contrary, if the sum is below the threshold, the output is set to zero, and the neuron is quiescent. For this reason, the binary output of a threshold unit can be loosely interpreted as the firing rate of the postsynaptic neuron, namely as the number of spikes per second of its action potential, which propagates along its axon toward other neurons in the network.

The basic set of equations defining the dynamics of the discrete-time binary network is:

which describes the temporal evolution of the neural activity $A_{i}$ of the $i$th neuron, from the time instant $t$ to the time instant $t+1$. In Eq. (2.1), $N$ represents the size (namely the number of threshold units) of the network. The matrix $J = \left[J_{i, j}\right]_{i, j = 0, \cdots, N-1}$ is the (generally asymmetric) synaptic connectivity matrix of the network, whose entries $J_{i, j}$ are time-independent and represent the strength or weight of the synaptic connection from the $j$th (presynaptic) neuron to the $i$th (postsynaptic) neuron. Moreover, $\mathfrak{I}_{i}+\sum_{j = 0}^{N-1}\sigma_{i, j}^{\mathcal{N}}\mathcal{N}_{j}\left(t\right)$ represents the total external input current (i.e. the stimulus) to the $i$th neuron. In more detail, $\mathfrak{I}_{i}$ is the time-independent deterministic component of the stimulus, while its stochastic component is the sum of $N$ random variables $\sigma_{i, j}^{\mathcal{N}}\mathcal{N}_{j}\left(t\right)$, each one having zero mean and standard deviation $\sigma_{i, j}^{\mathcal{N}}$. The vector $\boldsymbol{\mathcal{N}}\overset{\mathrm{def}}{ = }\left[\mathcal{N}_{0}, \; \cdots, \; \mathcal{N}_{N-1}\right]^{T}$ represents a collection of stochastic variables with unit standard deviation, whose joint probability distribution $p_{\boldsymbol{\mathcal{N}}}$ is arbitrary. Then, in Eq. (2.1), $\mathscr{H}\left(\cdot\right)$ is the Heaviside activation function with threshold $\theta$, which is defined as follows:

It is important to observe that, unlike the classic Hopfield network ^[7], which is symmetric and asynchronously updated, a Lyapunov function for synchronous networks with asymmetric synaptic weights, like ours, is generally not known. For this reason, the analytical investigation of the network dynamics determined by Eq. (2.1) proves much more challenging. In Sections (3)-(5), we will review the techniques, that we developed in ^[14,15,16], for investigating the dynamical and probabilistic properties of Eq. (2.1) in small networks.

3.   Analysis of bifurcations in deterministic networks

An important problem in the theory of binary networks is represented by the study of the qualitative changes in the dynamics of their neuronal activity, which typically are elicited by variations in the external stimuli. These changes of dynamics are named in mathematical terms as bifurcations ^[45]. Seminal work in physics focused on the study of bifurcations in infinite-size spin networks, see e.g. ^[21,46,23]. On the other hand, the theory of bifurcations of non-smooth dynamical systems composed of a finite number of units, including those with discontinuous functions such as binary networks, has been developed mostly for continuous-time models, see e.g. ^{[47,48,49,50,51]}, and for piecewise-smooth continuous maps ^[52,53]. Discontinuous maps are formally described by the following recurrence relation:

where the function $f$ is discontinuous across the boundaries of the regions $S_{1}, \cdots, S_{n}$ of the phase space, while $p$ is a bifurcation parameter. Despite the important role that discontinuous maps play in computational neuroscience, the development of new techniques for studying their bifurcation structure has received much less attention ^[54].

In ^[14,15] we tackled the problem of deriving the bifurcations in the dynamics of the neural activity for finite-size binary networks with arbitrary connectivity matrix, which are described by the discontinuous map Eq. (2.1). As is common practice, we performed the bifurcation analysis in the zero-noise limit, i.e. for $\sigma_{i, j}^{\mathcal{N}}\rightarrow0\; \forall i, j$ (see Eq. (2.1)). In particular, we studied how the dynamics of neural activity switches between stationary states and neural oscillations, when varying the external stimuli to the network. Because of the discrete nature of the neural activity, there exists only a finite number of stationary and oscillatory solutions to Eq. (2.1). This allowed us to introduce a combinatorial brute-force approach for studying the bifurcation structure of binary networks, which we describe briefly below.

We introduce the vector $\boldsymbol{A}\overset{\mathrm{def}}{ = }\left[A_{0}, \; \cdots, \; A_{N-1}\right]^{T}$ containing the activities of all $N$ neurons, and the sequence $\mathcal{S}\left(0, \mathcal{T}\right)$ of activity vectors $\boldsymbol{A}^{\left(0\right)}\rightarrow\boldsymbol{A}^{\left(1\right)}\rightarrow\cdots\rightarrow\boldsymbol{A}^{\left(\mathcal{T}\right)}$, for some $1\leq\mathcal{T}\leq2^{N}$. Given a network with $\mathfrak{P}$ distinct input currents $I_{0}, \cdots, I_{\mathfrak{P}-1}$, we also define $\Gamma_{I_{\alpha}}$ to be the set of neurons that share the same external current $I_{\alpha}$ (namely $\Gamma_{I_{\alpha}}\overset{\mathrm{def}}{ = }\left\{ i\in\left\{ 0, \cdots, N-1\right\} :\; \mathfrak{I}_{i} = I_{\alpha}\right\}$), and $\Gamma_{I_{\alpha}, x}^{\left(j\right)}\overset{\mathrm{def}}{ = }\left\{ i\in\Gamma_{I_{\alpha}}:\; A_{i}^{\left(j\right)} = x\right\}$ for $x\in\left\{ 0, 1\right\}$. Then, in ^[15] we proved that the sequence $\mathcal{S}\left(0, \mathcal{T}\right)$ is a solution of Eq. (2.1) in the time range $\left[0, \mathcal{T}\right]$ (i.e. $\boldsymbol{A}\left(t = j\right) = \boldsymbol{A}^{\left(j\right)}$ for $j = 0, \cdots, \mathcal{T}$), for every combination of stimuli $\left(I_{0}, \cdots, I_{\mathfrak{P}-1}\right)\in\mathfrak{V} = \mathcal{V}_{0}\times\cdots\times\mathcal{V}_{\mathfrak{P}-1}$, where:

In Eq. (3.1), $\mathcal{V}_{\alpha}$ for $\alpha = 0, \cdots, \mathfrak{P}-1$ represent the edges of the $\mathfrak{P}$-dimensional hyperrectangle $\mathfrak{V}$ within which the sequence $\mathcal{S}\left(0, \mathcal{T}\right)$ is stable. $\mathcal{S}\left(0, \mathcal{T}\right)$ loses its stability, namely it turns into another sequence, at the boundaries $\Lambda_{\alpha}$ and $\Xi_{\alpha}$, which therefore represent the coordinates of the bifurcation points of the neural activity sequence. It is important to observe that $\Lambda_{\alpha}$ and $\Xi_{\alpha}$ are explicitly calculated from Eq. (3.1) in terms of the connectivity matrix and the firing thresholds. For this reason, by using our approach, there is no need to perform any fine (and computationally expensive) binning of the stimulus space in order to detect the bifurcation points of the sequence $\mathcal{S}\left(0, \mathcal{T}\right)$.

In this review, we focus specifically on the subset of sequences that satisfy the additional constraint $\boldsymbol{A}^{\left(0\right)} = \boldsymbol{A}^{\left(\mathcal{T}\right)}$: these sequences represent the candidate oscillatory solutions with period $\mathcal{T}$ of Eq. (2.1). We also observe that, in the special case $\mathcal{T} = 1$, we obtain the set of candidate stationary solutions of the network equations. For this reason, the bifurcation diagram of a binary network can be decomposed into two panels, the oscillation and the multistability diagrams. These diagrams describe the relationship between the oscillatory/stationary solutions of Eq. (2.1), and the set of stimuli. In other words, these diagrams display the fragmentation of the stimulus space into areas where several oscillatory solutions occur, and/or where the network is (multi)stable.

It is important to note that only the sequences whose hyperrectangles $\mathfrak{V}$ have positive hypervolumes (i.e. the sequences that satisfy the condition $\Lambda_{\alpha} < \Xi_{\alpha}$, for every $\alpha$ and $j$ such that $\Gamma_{I_{\alpha}, 0}^{\left(j+1\right)}, \; \Gamma_{I_{\alpha}, 1}^{\left(j+1\right)}\neq\emptyset$) are solutions of Eq. (2.1), for some combinations of stimuli. On the contrary, if the hypervolume of $\mathfrak{V}$ is zero, the corresponding neural sequence is never a solution of Eq. (2.1). Unfortunately, the sequences with positive hypervolumes are not known a priori, therefore they must be found through a brute-force searching procedure. Because of the combinatorial explosion of the number of possible sequences for increasing $N$, typically brute-force algorithms have at least exponential complexity with respect to the network size. Therefore they can be applied only to small networks (typically $N < 30$), regardless of the density of their synaptic connections. However, real cortical circuits are typically very sparse. At macroscopic scales of spatial organization, the average density of the synaptic connections (i.e. the ratio between the actual and the maximum possible number of connections among the neurons within the considered network), is approximately $10^{-6}-10^{-7}$ across the whole cortex, and it can increase up to $0.2-0.4$ in connection pathways linking cortical areas ^[55]. Moreover, mesoscopic structures such as cortical columns, contain thousands of neurons and millions of synaptic connections. In ^[56], the authors simulated a model of cortical column consisting of about $80,000$ neurons and $0.3$ billion synapses, so that the average density of the synaptic connections of the model is about $0.047$. Another estimation of the sparseness of the neural circuits within a cortical column was reported in ^[57], where the authors showed that more than $80\%$ of the synapses on each neuron in a cortical column originate from cells outside the local volume. To account for network sparseness, in ^[15] we developed an efficient algorithm specifically designed for networks with a low density of the synaptic connections. This efficient algorithm takes advantage of the information provided by the absence of the synaptic connections among the threshold units to speed up the detection of the oscillatory and stationary solutions of Eq. (2.1). In other words, the sparse-efficient algorithm avoids checking the sequences of neural activity vectors that are not compatible with the topology of the synaptic connections, resulting in a much faster calculation of the bifurcation structure of the network model. For example, in ^[15] we showed that, for networks composed of $50$ neurons and about $70$ synapses (network sparseness $\sim0.028$), the algorithm derived the complete multistability diagram in less than $2$ seconds on laptop, by identifying the stationary states of the network among a set of about $10^{15}$ possible states. Moreover, the algorithm detected the whole set of network oscillations with period $\mathcal{T} = 2$ in $7$ minutes, among a set of about $10^{29}$ possible oscillations with that period. The interested reader is referred to ^[15] for a detailed discussion of the algorithm.

In Figure (1) we show an example of bifurcation diagram, that we obtained from Eq. (3.1), in the specific case of the network parameters reported in Table (1). In this example, we consider a network with heterogeneous random synaptic weights, which is composed of $3$ excitatory neurons and $2$ inhibitory neurons. Moreover, in Figure (2), we show two examples of state-to-state transitions, obtained by solving Eq. (2.1) in the zero-noise limit, for all the $2^{N}$ initial conditions of the network dynamics (i.e. from $\boldsymbol{A}\left(t = 0\right) = \left[0, \; \cdots, \; 0\right]^{T}$ to $\boldsymbol{A}\left(t = 0\right) = \left[1, \; \cdots, \; 1\right]^{T}$).

Figure 1.  An example of bifurcation diagram. This figure shows the bifurcation diagram of the binary network, obtained for the parameters in Table (1). We supposed that the excitatory neurons with indexes $i = 0$ and $i = 1$ receive an arbitrary external stimulus $\mathfrak{I}_{0} = \mathfrak{I}_{1} = I_{E}$, which represents the first bifurcation parameter, while the excitatory neuron with index $i = 2$ receives a fixed stimulus $\mathfrak{I}_{2} = 10$. Moreover, we assumed that the inhibitory neuron with index $i = 3$ receives an arbitrary stimulus $\mathfrak{I}_{3} = I_{I}$, which represents the second bifurcation parameter, while the inhibitory neuron with index $i = 4$ receives a fixed stimulus $\mathfrak{I}_{4} = 5$. Then, we plotted the multistability and oscillation diagrams in the $I_{E}-I_{I}$ plane, according to Eq. (3.1). A) Multistability diagram. Each color represents a different degree of multistability (white = astable, red = monostable, green = bistable, blue = tristable). B) Oscillation diagram. Each color represents a different set of oscillatory solutions of Eq. (2.1) (the notation $x:y$ reveals the formation of $y$ distinct oscillations with period $\mathcal{T} = x$). For example, for every combination of stimuli $\left(I_{E}, I_{I}\right)$ that lies in the yellow area, Eq. (2.1) has $2$ oscillatory solutions with period $\mathcal{T} = 2$, while for every combination in the green areas, the equation has an oscillatory solution with period $\mathcal{T} = 3$. Note that, for other connectivity matrices $J$, oscillations with distinct periods may coexist in the same area.

DownLoad: Full-Size Img PowerPoint

Table 1.  Network parameters 1. This table reports the values of the network parameters that we used for plotting Figures (1)-(4).

$J=\left[\begin{array}{ccccc} 0 & 17 & 17 & -43 & -6\\ 25 & 0 & 15 & -3 & -32\\ 10 & 1 & 0 & -10 & -7\\ 50 & 29 & 6 & 0 & -15\\ 7 & 28 & 5 & -95 & 0 \end{array}\right], \quad\theta_{0}=\cdots=\theta_{4}=1$

---

 | Show Table

DownLoad: CSV

Figure 2.  Examples of state-to-state transitions of a binary network. This figure shows the allowed transitions between states of neural activity, obtained for the network parameters in Table (1) and $\mathfrak{I}_{2} = 10$, $\mathfrak{I}_{4} = 5$. The nodes in the graphs represent the $2^{N}$ states of the neural activity vector $\boldsymbol{A}$, namely they are numbered taking the activity states as binary words (e.g. the node $26$ corresponds to the state $\boldsymbol{A} = [1, \; 1, \; 0, \; 1, \; 0]^{T}$), while the arrows represent the allowed transitions between these states. A) State-to-state transitions, obtained for $I_{E} = -10.5$ and $I_{I} = 6$. We highlighted in red $3$ stationary states (i.e. the nodes $2$, $5$ and $14$) and $2$ oscillations of period $\mathcal{T} = 2$ (i.e. $0\rightarrow7\rightarrow0$ and $6\rightarrow10\rightarrow6$). Note that, as expected, the point in the stimulus plane with coordinates $\left(I_{E}, I_{I}\right) = \left(-10.5, 6\right)$ lies in the blue area of the multistability diagram (see Figure (1), panel A), which corresponds to tristability, and in the yellow area of the oscillation diagram (see Figure (1), panel B), which corresponds to the formation of $2$ oscillations with period $\mathcal{T} = 2$. B) State-to-state transitions, obtained for $I_{E} = -5$ and $I_{I} = -55$. We highlighted in red an oscillations of period $\mathcal{T} = 5$ (i.e. $6\rightarrow8\rightarrow21\rightarrow29\rightarrow31\rightarrow6$). Note that the point in the stimulus plane with coordinates $\left(I_{E}, I_{I}\right) = \left(-5, -55\right)$ lies in the white area of the multistability diagram, where Eq. (2.1) has no stationary solutions, and in the cyan area of the oscillation diagram, which corresponds to the formation of an oscillation with period $\mathcal{T} = 5$.

DownLoad: Full-Size Img PowerPoint

Figure 3.  Probability distributions in stochastic networks. This figure shows the probability mass function $P\left(\boldsymbol{A}, t\right)$ and the marginal probability density function $p\left(V_{i}, t\right)$ of the $0$th neuron, obtained for the network parameters in Table (1), $\boldsymbol{\mathfrak{I}} = \left[0, \; -3, \; 2, \; -4, \; 0\right]^{T}$, and $\Psi = \mathrm{diag}\left(2, \; 3, \; 2, \; 3, \; 3\right)$. A) - B) Probability distributions obtained for correlated normally-distributed noise sources, with uniform cross-correlation (namely $\mathrm{Corr}\left(\mathcal{N}_{i}\left(t\right), \mathcal{N}_{j}\left(t\right)\right) = 0.8$ $\forall i, \, j$ such that $i\neq j$, and $\mathrm{Corr}\left(\mathcal{N}_{i}\left(t\right), \mathcal{N}_{j}\left(s\right)\right) = 0$ $\forall t, \, s$ such that $t\neq s$). C) - D) Probability distributions obtained for independent noise sources with Laplace distributions, namely $p_{\boldsymbol{\mathcal{N}}}\left(\boldsymbol{x}\right) = \prod_{i = 0}^{N-1}p_{\mathcal{N}_{i}}\left(x_{i}\right)$, where $p_{\mathcal{N}_{i}}\left(x_{i}\right) = \frac{\sqrt{2}}{2}\exp\left(-\sqrt{2}\left|x_{i}\right|\right)$. The red bars in panels A and C are calculated analytically from Eq. (4.6), while the red curves in panels B and D are calculated according to Eq. (4.8). The blue bars in panels A, C, and the blue dots in panels B, D are calculated numerically through a Monte Carlo method, namely by solving iteratively Eqs. (2.1) and (4.2) $10,000$ times in the time range $\left[0,100\right]$, and then by calculating the probability distributions across the repetitions at $t = 100$.

DownLoad: Full-Size Img PowerPoint

Figure 4.  Cross-neuron correlations. This figure shows the dependence of the cross-neuron correlation of the binary network on the external stimulus, in the specific case of independent normally-distributed noise sources. Correlation is calculated between the neurons with indexes $i = 0$ and $i = 4$, for the values of the parameters reported in Table (1), $\boldsymbol{\mathfrak{I}} = I\left[1, \; 1, \; 1, \; 1, \; 1\right]^{T}$, $I\in\left[-14, 14\right]$, and $\Psi = \mathrm{diag}\left(2, \; 3, \; 2, \; 3, \; 3\right)$. A) Correlation between neural activities. The red curve is calculated analytically from Eq. (4.9). B) Correlation between membrane potentials. The red curve is calculated according to Eq. (4.10). In both panels, the blue dots are calculated numerically through a Monte Carlo method over $10,000$ repetitions, collected at $t = 100$. Moreover, the vertical dashed lines correspond to the stimuli $I = -12$ and $I = 8.5$. We chose these values of the stimuli to show that low (respectively high) correlations between neural activities do not necessarily correspond to low (respectively high) correlations between the membrane potentials ($\mathrm{Corr}\left(A_{0}, A_{4}\right)\approx0.99$ and $\mathrm{Corr}\left(V_{0}, V_{4}\right)\approx0.02$ for $I = -12$, while $\mathrm{Corr}\left(A_{0}, A_{4}\right)\approx0.06$ and $\mathrm{Corr}\left(V_{0}, V_{4}\right)\approx0.65$ for $I = 8.5$).

DownLoad: Full-Size Img PowerPoint

4.   Stochastic networks

4.1.   Dualism between neural activity and membrane potentials

For networks of synchronously-updated binary neurons, such as those considered in this article, it is possible to derive explicit dynamical equations of the membrane potentials, and to show the existence of a dualism between them and the neural activity states. The membrane potentials and the activity states are intrinsically related, but have distinct dynamical aspects. On the other hand, networks of asynchronously-updated binary neurons (see e.g. ^[7]), as well as maximum-entropy models and information theoretic approaches (see e.g. ^[58,59]), do not model membrane dynamics, therefore they are not considered in this review.

As we explained in section (2), the term $\sum_{j = 0}^{N-1}J_{i, j}A_{j}\left(t\right)$ in Eq. (2.1) can be loosely interpreted as the weighted sum of postsynaptic potentials at neural dendrites. Therefore this term, plus the eventual external stimulus to the $i$th neuron, can be interpreted as the total membrane potential $V_{i}$ of that neuron, namely:

Interestingly, we show that it is possible to exactly derive the set of equations satisfied by the membrane potentials, and that these equations provide a complementary description of the network dynamics with respect to Eq. (2.1). Under the change of variables Eq. (4.1), we observe that Eq. (2.1) can be transformed into the following set of equations for the membrane potentials:

Note also that Eqs. (2.1) and (4.1) imply $A_{i}\left(t\right) = \mathscr{H}\left(V_{i}\left(t\right)-\theta_{i}\right)\in\left\{ 0, 1\right\}$, and that the binary output of a threshold unit can be loosely interpreted as the firing rate of that neuron: $A_{i}\left(t\right) = 0$ if the $i$th neuron is not firing at time $t$, and $A_{i}\left(t\right) = 1$ if it is firing at unit rate.

Since Eq. (2.1) does not depend on the variables $V_{i}\left(t\right)$, it can be solved without knowing the behavior of the membrane potentials. Solving the equations that govern the microscopic dynamics of the membrane potentials in rate models is typically difficult, and several authors opted either for calculating population-averaged potentials, or for deriving approximate solutions in the large-network-size limit (see e.g. ^[60,36]). For these reasons, the exact calculation of the microscopic dynamics of $V_{i}\left(t\right)$ from Eq. (4.2) has always been neglected in the literature. On the contrary, in ^[14] we introduced a technique for solving Eq. (4.2), and we found the exact joint probability distribution of the microscopic membrane potentials in the long-time limit, which is reported later in this subsection.

We observe that the neural activities are discrete random variables, therefore they are described by probability mass functions (pmfs). We introduce the vector $\boldsymbol{A} = \left[A_{0}, \; \cdots, \; A_{N-1}\right]^{T}$ containing the activities of all $N$ neurons at time $t+1$, and $\boldsymbol{A}'$ the vector of the activities of all neurons at time $t$. We also define:

$\Psi$ represents the matrix of the standard deviations of the noise sources, and it is invertible by hypothesis. $\boldsymbol{\theta}$ and $\boldsymbol{\mathfrak{I}}$ represent, respectively, the vectors of the firing thresholds and of the deterministic external stimuli, while $\pmb{\mathscr{H}}$ is the element-wise Heaviside activation function. Moreover, we introduce the matrix $\mathcal{C} = \left[\mathcal{C}_{i, j}\right]_{i, j = 0, \cdots, 2^{N}-2}$, such that:

where $\delta_{i, j}$ is the Kronecker delta, and:

where we remind that $p_{\boldsymbol{\mathcal{N}}}$ is the joint probability distribution of the stochastic variables $\mathcal{N}_{0}, \; \cdots, \; \mathcal{N}_{N-1}$. In Eq. (4.4), $\pmb{\mathscr{B}}_{i}^{\left(N\right)}$ is the $N\times1$ vector whose entries are the digits of the binary representation of the index $i$ (e.g. $\pmb{\mathscr{B}}_{26}^{\left(5\right)} = \left[1, \; 1, \; 0, \; 1, \; 0\right]^{T}$). Moreover, the set $\mathscr{V}_{i}^{\left(N\right)}$ is defined as follows:

where $\mathscr{D}\left(\boldsymbol{\nu}\right)$ is the decimal representation of the binary vector $\boldsymbol{\nu}$. For example, for $N = 2$, we get:

(note for example that, for any $\left(x_{0}, x_{1}\right)\in\mathscr{V}_{3}^{\left(2\right)}$, we get $\left[\begin{array}{c} \mathscr{H}\left(x_{0}-\theta_{0}\right)\\ \mathscr{H}\left(x_{1}-\theta_{1}\right) \end{array}\right] = \left[\begin{array}{c} 1\\ 1 \end{array}\right]$, whose decimal representation is $3$). Finally, we define:

Then, in ^[14] we proved that the conditional probability distribution of $\boldsymbol{A}$ given $\boldsymbol{A}'$, and the stationary joint distribution of $\boldsymbol{A}$ in the limit $t\rightarrow+\infty$, are:

Note that, according to the theory of Markov processes (see e.g. ^[8,9]), $\underset{t\rightarrow+\infty}{\lim}P\left(\boldsymbol{A}, t\right)$ can be alternatively obtained as the eigenvector to the eigenvalue $1$ of the transition matrix $P\left(\boldsymbol{A}, t+1|\boldsymbol{A}', t\right)$. However, the Markov's formulation would obscure the dualism between the neural activity and membrane potentials, as we will see in what follows. For this reason, Eq. (4.6) is a more convenient description of the joint probability distribution of the neural activity.

On the other hand, the membrane potentials $V_{i}\left(t\right)$ are continuous variables, therefore they are described by probability density functions (pdfs). By introducing the vector $\boldsymbol{V}\overset{\mathrm{def}}{ = }\left[V_{0}, \; \cdots, \; V_{N-1}\right]^{T}$, which contains the membrane potentials of all $N$ neurons at time $t+1$, and the vector $\boldsymbol{V}'$ of the membrane potentials of all neurons at time $t$, the conditional probability distribution of $\boldsymbol{V}$ given $\boldsymbol{V}'$, and the stationary joint distribution of $\boldsymbol{V}$ in the limit $t\rightarrow+\infty$, are ^[14]:

Note that Eq. (4.8) can be obtained in two different, non-trivial ways: either from Eqs. (4.1) and (4.6), or by solving Eq. (4.2) in the long-time limit. In this work we followed the latter approach, whose details are reported in the supplemental information of ^[14]. There we also showed that $H_{j} = \underset{t\rightarrow+\infty}{\lim}\int_{\mathscr{V}_{j}^{\left(N\right)}}p\left(\boldsymbol{V}, t\right)d\boldsymbol{V}$, therefore the coefficient $H_{j}$ in Eqs. (4.6) and (4.8) can be finally interpreted as the probability that the vector $\boldsymbol{A}\left(t\right)$ equals the binary representation of the index $j$ in the long-time limit (for example, in a network of size $N = 5$, we get $H_{26} = P\left(\boldsymbol{A}\left(t\right) = \left[1, \; 1, \; 0, \; 1, \; 0\right]^{T}\right)$ for $t\rightarrow+\infty$).

By comparing Eqs. (4.5) and (4.6) with, respectively, Eqs. (4.7) and (4.8), we observe that there exists a close relationship between the neural activities $\boldsymbol{A}$ and the membrane potentials $\boldsymbol{V}$, which we further investigate in Subsection (4.2). It is also important to note that, generally, the integrals in Eqs. (4.4)-(4.6) can be calculated only numerically or through analytical approximations. However, in the specific case when the matrix $\Psi$ is diagonal, and the stochastic variables $\mathcal{N}_{i}\left(t\right)$ are independent, exact analytical solutions can be found in terms of the cumulative distribution functions of the noise sources (see ^[14]).

In Figure (3) we show an example of the joint probability distributions $P\left(\boldsymbol{A}, t\right)$ and $p\left(\boldsymbol{V}, t\right)$ in the large-time limit, obtained for the network parameters that we reported in Table (1). In this figure we considered two distinct distributions of the noise sources (namely correlated normally-distributed variables $\mathcal{N}_{i}\left(t\right)$, as well as independent sources with Laplace distributions), and we showed how they differently shape the probability distributions of the neural activity and of the membrane potentials.

4.2.   Cross-neuron correlations

The joint probability distributions $P\left(\boldsymbol{A}, t\right)$ and $p\left(\boldsymbol{V}, t\right)$, that we reported in the previous section, provide a complete probabilistic description of the network in the limit $t\rightarrow+\infty$. In particular, these distributions can be used for calculating cross-neuron correlations, which represents a powerful tool for quantifying the exchange of information between neurons. In ^[14], we derived exact analytical expressions of the Pearson correlation coefficient for $t\rightarrow+\infty$, in the case when the matrix $\Psi$ is diagonal (i.e. $\Psi = \mathrm{diag}\left(\sigma_{0}^{\mathcal{N}}, \; \cdots, \; \sigma_{N-1}^{\mathcal{N}}\right)$, where $\sigma_{i}^{\mathcal{N}}\overset{\mathrm{def}}{ = }\sigma_{i, i}^{\mathcal{N}}\; \forall i$), while the noise sources are independent and normally distributed. By applying Eq. (4.6), we found that the pairwise correlation between the neural activities of two neurons with indexes $i$ and $j$ (such that $i\neq j$) is:

where $\mathscr{B}_{n, m}^{\left(N\right)}\overset{\mathrm{def}}{ = }\left[\pmb{\mathscr{B}}_{n}^{\left(N\right)}\right]_{m}$ is the $m$th entry of the $N\times1$ vector $\pmb{\mathscr{B}}_{n}^{\left(N\right)}$. In a similar way, from Eq. (4.8) we derived the following formula for the pairwise correlation between the membrane potentials:

Eqs. (4.9) and (4.10) show that the cross-neuron correlations are generally very complex functions of the network parameters. However, these formulas strongly simplify in special limiting cases, such as the weak-noise or the strong-stimuli limits, which are not discussed here. The interested reader is referred to the supplemental information of ^[14] for a rigorous derivation of these simplified formulas.

In Figure (4) we plotted some examples of cross-correlations, obtained for the the network parameters that we reported in Table (1). This figure shows that variations of the external stimuli $\boldsymbol{\mathfrak{I}}$ switch the binary network between synchronous (i.e. highly correlated) and asynchronous (i.e. uncorrelated) states. Moreover, we observe that low (respectively high) correlations between neural activities do not necessarily correspond to low (respectively high) correlations between the membrane potentials. In other words, the linear relationship between the neural activity and the membrane potentials, as given by Eq. (4.1), is not reflected by the correlation structure of these variables. This result proves that, despite the similarity of the corresponding equations (which we already observed in Subsection (4.1), by comparing Eqs. (4.5) and (4.6) with, respectively, Eq. (4.7) and (4.8)) and their linear relationship, neural activity and the membrane potentials represent two distinct aspects of binary networks.

The interested reader is referred to ^[14] for a detailed description of the conditions under which synchronous and asynchronous states occur in the network, and for the extension of Eqs. (4.9) and (4.10) to encompass higher-order (i.e. groupwise) correlations among an arbitrary number of neurons.

4.3.   Pattern storage and retrieval in presence of noise

In this section we consider the problem of storing $\mathcal{D}$ sequences of neural activity vectors $\boldsymbol{A}^{\left(i, 0\right)}\rightarrow\boldsymbol{A}^{\left(i, 1\right)}\rightarrow\cdots\rightarrow\boldsymbol{A}^{\left(i, \mathcal{T}_{i}\right)}$, for some $\mathcal{T}_{i}$ (which represents the number of transitions between activity vectors in the $i$th sequence), and $i = 0, \cdots, \mathcal{D}-1$. In the context of content-addressable memories, one aims to determine a synaptic connectivity matrix $J$ that stores these sequences in the binary network, so that each sequence can be retrieved by initializing the network state to $\boldsymbol{A}\left(t = 0\right) = \boldsymbol{A}^{\left(i, 0\right)}$, even in the presence of noise. Any method for calculating such a connectivity matrix is typically called learning rule.

Each transition $\boldsymbol{A}^{\left(i, n_{i}\right)}\rightarrow\boldsymbol{A}^{\left(i, n_{i}+1\right)}$ in the sequences is noise-resistant whenever $P\left(\boldsymbol{A}^{\left(i, n_{i}+1\right)}, t_{n_{i}}+1|\boldsymbol{A}^{\left(i, n_{i}\right)}, t_{n_{i}}\right)\approx1$, since under this condition the probability that the state $\boldsymbol{A}^{\left(i, n_{i}\right)}$ switches to a state other than $\boldsymbol{A}^{\left(i, n_{i}+1\right)}$ at the time instant $t_{n_{i}}+1$ is negligible. Therefore, according to Eq. (4.5), the sequences of neural activity can be stored in the network by solving the following set of equations:

with respect to the connectivity matrix $J$. Generally, these equations can be solved only numerically or through analytical approximations. However, in the specific case when the matrix $\Psi$ is diagonal and the stochastic variables $\mathcal{N}_{i}\left(t\right)$ are independent, exact analytical solutions can be found.

In ^[14], we considered the case of independent normally-distributed noise sources, and we found that Eq. (4.11) reduces to:

(where $\mathrm{erf}\left(\cdot\right)$ is the error function), or in other words:

for $n_{i} = 0, \cdots, \mathcal{T}_{i}-1$, $i = 0, \cdots, \mathcal{D}-1$ and for any sufficiently large and positive constant $K_{j}^{\left(i, n_{i}\right)}$. If the network is fully-connected without self-connections (so that $J_{i, i} = 0$), Eq. (4.12) can be interpreted as the following $N$ systems of linear algebraic equations:

In Eq. (4.13), $\boldsymbol{J}^{\left(j\right)}$ is the $\left(N-1\right)\times1$ vector with entries $J_{j, k}$ for $k\neq j$. Moreover, if we define $\mathscr{T}\overset{\mathrm{def}}{ = }\sum_{i = 0}^{\mathcal{D}-1}\mathcal{T}_{i}$, $\mathfrak{T}_{i}\overset{\mathrm{def}}{ = }\sum_{k = 1}^{i}\mathcal{T}_{k-1}$ (for $i > 0$) and $\mathfrak{T}_{0}\overset{\mathrm{def}}{ = }0$, then $\boldsymbol{u}^{\left(j\right)}$ is a $\mathscr{T}\times1$ vector with entries:

Moreover, in Eq. (4.13), $\Omega^{\left(j\right)}$ is the $\mathscr{T}\times\left(N-1\right)$ matrix obtained by removing the $j$th column of the following matrix:

In particular, we observe that whenever $\boldsymbol{A}^{\left(i, 0\right)} = \boldsymbol{A}^{\left(i, \mathcal{T}_{i}\right)}$, the $i$th neural sequence is an oscillatory solution of Eq. (2.1) with period $\mathcal{T}_{i}$, so that if the matrix $J$ is calculated by solving Eq. (4.13) and the network is initialized to any state of the oscillation, the network will cycle repeatedly through the same set of states. Moreover, in the special case $\mathcal{T}_{i} = 1$, the neural sequence represents a stationary solution of Eq. (2.1).

In Figure (5), we show some examples of storage of stationary patterns and oscillatory sequences with $\mathscr{T} = 3$. This figure is obtained for the the network parameters reported in Table (2), and proves that the learning rule Eq. (4.13) can be used to store safely sequences of neural activity also in very noisy networks. To conclude, we observe that generally the entries of the matrices $J$ that store these patterns in very noisy networks have large absolute values (see Table (2)). This is a consequence of the fact that $J_{ij}\propto\sigma^{\mathcal{N}}$ for $\sigma^{\mathcal{N}}\gg1$, according to Eqs. (4.13) and (4.14). Large entries of the connectivity matrix ensure that the membrane potentials fluctuate, under noise perturbations, far away from the firing thresholds, so that undesired transitions between the activity states are unlikely. This mechanism protects the sequences of activity patterns from the disruptive effects of noise fluctuations, allowing their storage also in highly noisy media such as the cortical tissue.

Figure 5.  Storage of activity patterns and neural sequences. This figure shows examples of three distinct sets of activity patterns and neural sequences (highlighted in red), stored in a stochastic binary network driven by independent and normally distributed noise sources, for $N = 4$, $\boldsymbol{\mathfrak{I}} = \left[\begin{array}{cccc} 0, & 0, & 0, & 0\end{array}\right]^{T}$, and $\Psi = \sigma^{\mathcal{N}}\mathrm{diag}\left(1, \; 1, \; 1, \; 1\right)$. The color gradation of the blue arrows is proportional to the magnitude of $P\left(\boldsymbol{A}^{\left(i, n_{i}+1\right)}, t_{n_{i}}+1|\boldsymbol{A}^{\left(i, n_{i}\right)}, t_{n_{i}}\right)$ (see Eq. (4.5)), so that the arrows are white for every pair of states $\left(\boldsymbol{A}^{\left(i, n_{i}\right)}, \boldsymbol{A}^{\left(i, n_{i}+1\right)}\right)$ such that $P\left(\boldsymbol{A}^{\left(i, n_{i}+1\right)}, t_{n_{i}}+1|\boldsymbol{A}^{\left(i, n_{i}\right)}, t_{n_{i}}\right)\approx0$, while they are blue if the conditional probability is close to $1$. The entries of the matrix $J$ are calculated from Eq. (4.13) for $\sigma^{\mathcal{N}} = 50$ (see the highlighted central panels) and $K_{j}^{\left(i, n_{i}\right)} = 10\; \forall i, \:j, \:n_{i}$, in order to store $3$ stationary states (panel B), an oscillation with period $\mathcal{T} = 3$ (panel E), and a stationary state plus an oscillation with period $\mathcal{T} = 2$ (panel H). The matrices $\Omega$ and the corresponding solutions $J$ are reported in Table (2). Then the stability of the stored patterns under noise perturbations is tested for $\sigma^{\mathcal{N}} = 0.1$ (panels A, D, G) and $\sigma^{\mathcal{N}} = 1000$ (panels C, F, I), without changing the matrix $J$. Note that the stationary states and the oscillations are noise-resistant, despite the presence of strong noise sources ($\sigma^{\mathcal{N}} = 50$). Moreover, in the case of weak noise sources ($\sigma^{\mathcal{N}} = 0.1$), new stable patterns may appear, while for extreme noise intensity ($\sigma^{\mathcal{N}} = 1000$) the stored patterns become unstable, and the network is fully random.

DownLoad: Full-Size Img PowerPoint

Table 2.  Network parameters 2. This table reports the matrices $\Omega$ and the corresponding matrices $J$, that we used for plotting Figure (5).

Panel A	Panel B	Panel C
$\Omega=\left[\begin{array}{cccc} 1 & 1 & 0 & 1\\ 0 & 1 & 1 & 1\\ 1 & 0 & 1 & 1 \end{array}\right]$	$\Omega=\left[\begin{array}{cccc} 1 & 0 & 1 & 1\\ 1 & 1 & 0 & 1\\ 1 & 0 & 1 & 0 \end{array}\right]$	$\Omega=\left[\begin{array}{cccc} 0 & 1 & 0 & 1\\ 1 & 1 & 1 & 0\\ 1 & 0 & 0 & 1 \end{array}\right]$
$J\approx\left[\begin{array}{cccc} 0 & -1414 & -1414 & 2122\\ -1414 & 0 & -1414 & 2122\\ -1414 & -1414 & 0 & 2122\\ 354 & 354 & 354 & 0 \end{array}\right]$	$J\approx\left[\begin{array}{cccc} 0 & 708 & 708 & 0\\ -2120 & 0 & 1414 & 1414\\ 708 & 1414 & 0 & -1414\\ 1 & -707 & 707 & 0 \end{array}\right]$	$J\approx\left[\begin{array}{cccc} 0 & 0 & -706 & 708\\ -1414 & 0 & 2122 & 708\\ -1060 & 354 & 0 & 354\\ 708 & -706 & 706 & 0 \end{array}\right]$

---

 | Show Table

DownLoad: CSV

5.   Networks with quenched disorder

The results reported in Sections (3) and (4) are valid for binary neural networks with arbitrary topology of the synaptic connections (which does not evolve over time). For this reason, they can be applied to networks with regular connectivity matrices $J$, as well as to random networks with frozen synaptic weights (see e.g. the network parameters reported in Table (1)). In other words, the connectivity matrix of random networks can be interpreted as a single realization of the synaptic wiring among neurons, generated according to some known probability distribution $p_{J}$. These models are said to present quenched disorder ^[21,22,61].

Each realization of the connectivity matrix, generated according to the distribution $p_{J}$, usually produces a distinct matrix $J$, which in turn gives rise to distinct dynamical properties of the neural activity. In particular, each realization typically produces distinct bifurcation diagrams. For this reason, in order to obtain statistically representative results, one needs to average the coordinates of the bifurcation points over the variability of the matrix $J$. More generally, one would be interested in determining the probability distribution of the bifurcation points over the matrix $J$.

In ^[16] we derived semi-analytical expressions of these probability distributions. For simplicity, we focused on the bifurcation points of the stationary states in the zero-noise limit $\sigma_{i, j}^{\mathcal{N}}\rightarrow0\; \forall i, j$ (the possibility to extend our approach to the study of neural oscillations is discussed briefly in Subsection (6.3.4)). In ^[16] we supposed that the entries $J_{i, j}$ of the connectivity matrix can be decomposed as the product of a synaptic weight, $W_{i, j}$ (which represents the random strength of interaction of the $j$th neuron on the $i$th neuron), with another random variable, $T_{i, j}$ (which represents either the presence, for $T_{i, j} = 1$, or the absence, for $T_{i, j} = 0$, of a synaptic connections from the $j$th to the $i$th neuron). The random variable $W_{i, j}$ is supposed to be continuous, and distributed according to some distribution $p_{W_{i, j}}$. On the other hand, the variable $T_{i, j}$ is discrete, and such that $T_{i, j} = 1$ with probability $\mathcal{P}_{i, j}\in\left[0, 1\right]$, while $T_{i, j} = 0$ with probability $1-\mathcal{P}_{i, j}$. We also supposed that the variables $\left\{ W_{i, j}, T_{i, j}\right\} _{i, j = 0, \cdots, N-1}$ are statistically independent. Because of these assumptions, the random variable $\mathcal{I}_{i}$ (as given by Eq. (3.1); note that the superscript $\left(j\right)$ can be omitted in the case of stationary states studied in this section) is distributed as follows:

where $\mathbb{P}\left(R\right)$ represents the power set of $R\overset{\mathrm{def}}{ = }\left\{ i\in\left\{ 0, \cdots, N-1\right\} :\; A_{i} = 1\right\}$. Moreover, we call $F_{\mathcal{I}_{i}}$ the cumulative distribution function of $\mathcal{I}_{i}$. Then, the coordinates of the bifurcation points, $\Lambda_{\alpha}$ and $\Xi_{\alpha}$ (see again Eq. (3.1)), are distributed as follows:

for $X\in\left\{ \Lambda_{\alpha}, \Xi_{\alpha}\right\}$. In Eq. (5.1), $\delta\left(\cdot\right)$ is the Dirac delta function, $p_{X^{c}}$ is the component of $p_{X}$ that describes the statistical behavior of the continuous values of $X$, and $F_{X}$ is the cumulative distribution function of $X$. Since, given the activity vector $\boldsymbol{A}$, the coordinates $\Lambda_{\alpha}$ and $\Xi_{\alpha}$ are, respectively, the maximum and minimum of conditionally independent variables $\mathcal{I}_{i}$ (see Eq. (3.1)), they must be distributed according to order statistics ^{[62,63,64,65]}. By calling $\mathrm{per}\left(\cdot\right)$ the matrix permanent (which is an analog of a matrix determinant, but with all the signs in the Laplace expansion by cofactors taken as positive), in ^[16] we proved that $X^{c}$ is distributed as follows:

where $\gamma_{\alpha, u}\overset{\mathrm{def}}{ = }\left|\Gamma_{I_{\alpha}, u}\right|$, while $\left[\boldsymbol{a}_{\alpha, 1}\left(\boldsymbol{\theta}-\boldsymbol{x}\right), \; \boldsymbol{F}_{\alpha, 1}^{\left(\gamma_{\alpha, 1}-1\right)}\left(x\right)\right]$ and $\left[\boldsymbol{a}_{\alpha, 0}\left(\boldsymbol{\theta}-\boldsymbol{x}\right), \; \mathbb{I}_{\gamma_{\alpha, 0}, \gamma_{\alpha, 0}-1}-\boldsymbol{F}_{\alpha, 0}^{\left(\gamma_{\alpha, 0}-1\right)}\left(x\right)\right]$ are $\gamma_{\alpha, 1}\times\gamma_{\alpha, 1}$ and $\gamma_{\alpha, 0}\times\gamma_{\alpha, 0}$ matrices respectively, $\boldsymbol{a}_{\alpha, u}\left(\boldsymbol{\theta}-\boldsymbol{x}\right)\overset{\mathrm{def}}{ = }\left[a_{i}\left(\theta_{i}-x\right)\right]_{i\in\Gamma_{I_{\alpha}, u}}$ and $\boldsymbol{F}_{\alpha, u}\left(x\right)\overset{\mathrm{def}}{ = }\left[F_{\mathcal{I}_{i}}\left(x\right)\right]_{i\in\Gamma_{I_{\alpha}, u}}$ are $\gamma_{\alpha, u}\times1$ column vectors, $\boldsymbol{F}_{\alpha, u}^{\left(v\right)}\left(x\right)\overset{\mathrm{def}}{ = }\underset{v-\mathrm{times}}{\left[\underbrace{\boldsymbol{F}_{\alpha, u}\left(x\right), \; \cdots, \; \boldsymbol{F}_{\alpha, u}\left(x\right)}\right]}$ is a $\gamma_{\alpha, u}\times v$ matrix, and $\mathbb{I}_{\gamma_{\alpha, 0}, \gamma_{\alpha, 0}-1}$ is the $\gamma_{\alpha, 0}\times\left(\gamma_{\alpha, 0}-1\right)$ all-ones matrix. Moreover, in Eq. (5.1), $\left\{ x_{q}\right\} _{q\in D}$ represents the set of the discrete values of $X$, at which the cumulative distribution function $F_{X}$, namely:

is (possibly) discontinuous. Note that $D = \Gamma_{I_{\alpha}, 1}$ and $D = \Gamma_{I_{\alpha}, 0}$, for $\Lambda_{\alpha}$ and $\Xi_{\alpha}$ respectively, while $x_{q} = \theta_{q}$.

The mean multistability diagram of the network is the plot of the bifurcation points $\Lambda_{\alpha}$ and $\Xi_{\alpha}$, averaged over the realizations of the synaptic connectivity matrix $J$. In other words, the mean bifurcation points $\left\langle \Lambda_{\alpha}\right\rangle$ and $\left\langle \Xi_{\alpha}\right\rangle$ (where the brackets $\left\langle \cdot\right\rangle$ represent the mean over the realizations) correspond to the values of the stimulus $I_{\alpha}$ at which a given neural activity state $\boldsymbol{A}$ loses its stability on average, turning into another stationary state or an oscillation. In ^[16] we proved that:

where the functions $p_{X^{c}}\left(x\right)$ and $F_{X}\left(x\right)$ are given, respectively, by Eqs. (5.2) and (5.3).

The probability that a given activity state $\boldsymbol{A}$ is stationary for a fixed combination of stimuli $\widehat{\boldsymbol{I}} = \left[\widehat{I}_{0}, \; \cdots, \; \widehat{I}_{\mathfrak{P}-1}\right]^{T}$, corresponds to the probability that $\widehat{\boldsymbol{I}}\in\mathfrak{V}$, where the coordinates of the hyperrectangle $\mathfrak{V}$ (i.e. the bifurcation points where $\boldsymbol{A}$ becomes unstable) are random variables, whose density functions are calculated from Eq. (3.1). In ^[16] we proved that $P\left(\widehat{\boldsymbol{I}}\in\mathfrak{V}\right)$ can be calculated from the cumulative distribution of the bifurcation points as follows:

Moreover, the probability to observe the state $\boldsymbol{A}$ in the whole multistability diagram of a single realization of the matrix $J$ (i.e. the probability that $\boldsymbol{A}$ is stationary, regardless of the specific combination of stimuli), corresponds to the probability that the hyperrectangle $\mathfrak{V}$ has positive hypervolume $\mathrm{vol}\left(\mathfrak{V}\right)$. In ^[16] we proved that, according to Eq. (3.1), this probability has the following expression:

Eqs. (5.1)-(5.6) provide a complete description of the statistical properties of the stationary states in networks with quenched disorder. It is also important to note that these equations are semi-analytical, since they are expressed in terms of 1D integrals containing the distribution $p_{W_{i, j}}$. These integrals may be calculated exactly for some $p_{W_{i, j}}$, for example in the case of normally-distributed weights. However, for simplicity, in this review and in ^[16], they are calculated through numerical integration schemes, because fully-analytical expressions may be very cumbersome.

In Figures (6) and (7) we show an example of these results for a specific distribution of the connectivity matrix. We consider a network composed of $3$ excitatory neurons (with indexes $i = 0, 1, 2$), and $2$ inhibitory neurons ($i = 3, 4$). The excitatory and inhibitory neurons receive, respectively, external stimuli $I_{E}$ and $I_{I}$. We also assume that the synaptic weights $W_{i, j}$ are distributed according to the following powerlaw distribution:

Figure 6.  Probability distribution of the synaptic connections. This figure reports the probability distribution of the network topology and of the synaptic weights, given the values of the network parameters reported in Table (3). A) Probability distribution of the network topology, namely the graph of the matrix $\left[\mathcal{P}_{i, j}\right]_{i, j = 0, \cdots, N-1}$. The $N = 5$ nodes in the graph represent the neurons in the network, while the arrow from the $j$th to the $i$th neuron represents the probability to observe the $i\leftarrow j$ synaptic connection in a single realization of the matrix $J$ (note that the color gradation is proportional to $\mathcal{P}_{i, j}$). B) Examples of the powerlaw probability distributions of the synaptic weights $W_{i, j}$, see Eq. (5.7).

DownLoad: Full-Size Img PowerPoint

Table 3.  Network parameters 3. This table contains the values of the parameters that we used for plotting Figures (6) and (7) The symbol $\times$ in the matrices $\mathfrak{S}$ and $\mathfrak{W}$ means that the probability distributions of the stationary states and of the bifurcation points are not affected by those parameters, since the corresponding synaptic connections are absent ($\mathcal{P}_{i, j} = 0$).

$\mathcal{P}=\left[\begin{array}{ccccc} 0.2 & 0.1 & 0 & 1 & 0\\ 0 & 0.2 & 0.7 & 0.3 & 0\\ 0.5 & 0 & 0 & 0.8 & 0.2\\ 0 & 0 & 0.3 & 0.8 & 1\\ 1 & 0 & 0.8 & 0.4 & 0 \end{array}\right],$	$\boldsymbol{\theta}=\left[\begin{array}{c} 0\\ 3\\ 2\\ 1\\ 2 \end{array}\right]$
$\mathfrak{S}=\left[\begin{array}{ccccc} 5 & 5 & \times & -1 & \times\\ \times & 3 & 1 & -2 & \times\\ 3 & \times & \times & -3 & -4\\ \times & \times & 2 & -1 & -1\\ 5 & \times & 1 & -3 & \times \end{array}\right],$	$\mathfrak{W}=\left[\begin{array}{ccccc} 1 & 2 & \times & 1 & \times\\ \times & 1 & 1 & 2 & \times\\ 2 & \times & \times & 1 & 3\\ \times & \times & 2 & 1 & 1\\ 1 & \times & 1 & 2 & \times \end{array}\right]$

---

 | Show Table

DownLoad: CSV

Figure 7.  Stationary behavior of a binary network with quenched disorder. This figure reports the probabilistic properties of the stationary states of a binary network with quenched disorder. The matrix $J$ is generated randomly from the powerlaw distribution Eq. (5.7), and for the network parameters in Table (3) (see also Figure (6)). A) Cumulative distribution function $F_{\Lambda_{E}}\left(x\right)$ of the activity state $\boldsymbol{A} = \left[1, \; 0, \; 1, \; 1, \; 0\right]{}^{T}$. The red curve is derived semi-analytically from Eq. (5.3), while the blue dots are calculated numerically through a Monte Carlo method over $10,000$ repetitions of the synaptic connectivity matrix. B) Mean multistability diagram of the network, obtained semi-analytically from Eq. (5.4). C) -D) Occurrence probability of the stationary states, obtained for the fixed stimuli $\left[ \begin{array}{l} {I_E}\\ {I_I} \end{array} \right] = {\boldsymbol{\hat I}} = \left[ \begin{array}{l} 1\\ 0 \end{array} \right]$ (panel C), and regardless of the stimuli (panel D). The red bars in panel C (respectively, panel D) are derived semi-analytically from Eq. (5.5) (respectively, Eq. (5.6)), while the blue bars are calculated numerically through a Monte Carlo method (see ^[16] for more details).

DownLoad: Full-Size Img PowerPoint

where $\mathfrak{S}$ and $\mathfrak{W}$ represent, respectively, the horizontal shift and the width of the support of the distribution. To conclude, in Table (3) we reported the values of the parameters $\mathcal{P}$, $\boldsymbol{\theta}$, $\mathfrak{S}$ and $\mathfrak{W}$ that we chose for this network. Figure (6) reports the graph of the matrix $\left[\mathcal{P}_{i, j}\right]_{i, j = 0, \cdots, N-1}$ and some examples of the powerlaw distribution of the synaptic weights $W_{i, j}$ (see Eq. (5.7)). Moreover, in Figure (7) we show the mean multistability diagram of the network, as well as the occurrence probability of the activity states for fixed stimuli (i.e. $I_{E} = 1$ and $I_{I} = 0$) and regardless of the stimuli (see, respectively, Eqs. (5.4), (5.5) and (5.6)).

6.   Discussion

New mathematical techniques for analytically investigating finite-size and small-size neural network models are invaluable theoretical tools for studying the brain at its multiple scales of spatial organization, that complement the already existing mean-field approaches. Studying how the complexity and dynamics of neuronal network activity change with the network size is of fundamental importance to understand why networks in the brain appear organized at multiple spatial scales. In this article, we reviewed the effort we made in this direction, trying to fill the gap in the current neuro-mathematical literature. In the following, we discuss strengths and weaknesses of the approach we developed so far, the implications of our work for specific issues related to bifurcation dynamics and learning, and for future progress in the understanding of the function of networks in the brain.

6.1.   Advantages and weakness of our approach

6.1.1.   Bifurcation analysis of binary networks

An effective tool for studying spin networks in physics is represented by their energy (Hamiltonian) function. In order to study the low-temperature physical properties of the network at the thermodynamic equilibrium, one is often interested in finding out the global (and possibly degenerate) minimum energy state of the network. This is known as the ground state of the system, and it can be calculated by minimizing the energy function over the space of all possible spin configurations at absolute temperature. This is an optimization problem, which, for networks on non-planar or three- or higher-dimensional lattices, has been proven to be NP-hard ^[66].

In the study of discrete-time binary neural networks, energy (or, more generally, Lyapunov) functions typically are known only for asynchronously updated neurons with symmetric synaptic connections ^[7]. For neural networks with asymmetric connectivity matrices and/or synchronous update, like the one we considered in this review, the search for the ground state(s) turns into the more general problem of determining the long-time (non-equilibrium) states in the zero-noise limit. It is important to observe that this problem is even more formidable than the search for the ground states, due to the intractable number of oscillatory sequences that can eventually be observed during the network dynamics. We are not aware of any algorithm that performs efficiently this highly demanding combinatorial analysis in synchronously-updated networks with asymmetric connections. The algorithm that we introduced in ^[15] represents an attempt to tackle this problem, in the specific case of networks with sparse synaptic connections.

Once the set of all the possible stationary and oscillatory solutions has been evaluated, Eq. (3.1) provides a fast, analytical way to calculate the bifurcation diagram of the network. On the other hand, the bifurcation analysis of networks composed of neurons with graded output typically requires numerical continuation techniques ^[45], which do not provide any analytical intuition of the mechanisms underlying the changes of dynamics.

6.1.2.   Probability distributions and cross-neuron correlations

The analytical results reported in Subsections (4.1) and (4.2) provide a complete description of the probabilistic behavior of the neural activity and of the membrane potentials in the long-time regime of single network realizations. These results provide new qualitative insights into the mechanisms underlying stochastic neuronal dynamics, which hold for any network size, and therefore are not restricted to small networks only.

However, the main drawback of Eqs. (4.6) and (4.8) (and, as a consequence, also of Eqs. (4.9) and (4.10)), is represented by the quantitative evaluation of the probability distributions and of the cross correlations. This requires the calculation of a number of coefficients $H_{i}$ that increases exponentially with the network size, and therefore proves intractable already for networks composed of a few tens of neurons. A more efficient estimation of these quantities for large networks can be performed numerically through Monte Carlo methods.

An alternative approach we developed to deal with single realizations of large-size networks, is the derivation of their Sznitman's mean-field equations. In ^[14] we focused on the specific case of fully-connected multi-population networks driven by independent sources of noise, and we found that, under this assumption, the neurons become independent in the thermodynamic limit. This phenomenon simplifies considerably the complexity of the network equations, and allowed us to perform a detailed bifurcation analysis of the model.

In order to provide a complete description of the stationary behavior of networks with quenched disorder, the calculation of the probability distributions of the bifurcation points across network realizations, that we reported in section (5) (see Eqs. (5.1)-(5.4)), as well as the calculation of the occurrence probability of the stationary states (Eqs. (5.5) and (5.6)), must be performed for every stationary state of the model. Unfortunately, these states are not known a priory, therefore the calculation of the probability distributions must be repeated for all the $2^{N}$ combinations of the neural activity states. A possible solution to this problem, in the specific case of sparse networks, is represented by the sparse-efficient algorithm that we introduced in ^[15]. This algorithm allows a fast evaluation of the stationary solutions of the network, so that the calculation of the probability distributions can be performed only on the actual stationary states detected by the sparse-efficient algorithm.

Similarly to the case of single realizations of large networks that we discussed above, an alternative approach for studying large networks with quenched disorder is represented by the derivation of asymptotic formulas. In ^[16] we focused specifically on the case of arbitrarily sparse networks composed of statistically-homogeneous populations, which are often considered to be a good approximation of biological circuits, see e.g. ^[67,61,68]. In particular, by applying the Fisher-Tippett-Gnedenko theorem, we found that the coordinates of the bifurcation points $\Lambda_{\alpha}$ and $\Xi_{\alpha}$ in large networks are distributed according to the Gumbel's (double exponential) probability distribution.

To conclude, another disadvantage of Eqs. (5.1)-(5.6) is represented by the numerical calculation of the matrix permanent, which is computationally demanding. The fastest known technique for calculating the permanent of arbitrary matrices is the Balasubramanian-Bax-Franklin-Glynn (BBFG) formula ^{[69,70,71,72]}, which has complexity $O\left(2^{N-1}N^{2}\right)$. In order to alleviate this computational bottleneck, in ^[16] we derived a closed-form analytical expression of the permanent of uniform block matrices, which proved much faster than the BBFG formula. This solution allowed us to speed up considerably the calculation of the bifurcation structure of statistically-homogeneous multi-population networks with quenched disorder.

6.1.3.   Learning rule

The algorithm we introduced in Subsection (4.3) for storing some desired sequences of neural activity, was obtained in ^[14] by manipulating the conditional probability distribution of the neural activity (see Eq. (4.5)). This distribution does not depend on the coefficients $H_{i}$, therefore the learning rule has polynomial complexity. For this reason, its applicability is not restricted to small networks only.

Another interesting property of this learning rule is represented by the possibility to store sequences of neural activity also in noisy networks. Typically, noise can break a neural sequence if the stochastic fluctuations are sufficiently strong. However, our algorithm is designed for being noise-resistant, namely the probability of breaking the sequence under the influence of noise can be made arbitrarily small.

In the literature, several learning rules have been proposed for networks of binary neurons. A mechanism for storing and retrieving static patterns of neural activity in networks with symmetric connectivity was proposed by Hopfield ^[7]. The storage of sequences of temporally evolving patterns was investigated by Sompolinsky, Kanter and Kleinfeld ^[73,74], for networks with non-instantaneous synaptic transmission between neurons. In ^[75], Dehaene et al introduced an alternative approach based on temporally evolving synapses. Then, Buhmann and Schulten ^[76] showed that asymmetric connectivity and noise are sufficient conditions for storing and retrieving temporal sequences, without further assumptions on the biophysical properties of the synaptic connections.

It is important to observe that the learning rule introduced in Subsection (4.3) does not require transmission delays, time-dependent synaptic strengths, or the presence of noise. Only asymmetric synaptic connections are required. This is compatible with experimental observation, in that the vast majority of synapses in real biological networks are asymmetric ^[77]. Our result can be considered as an extension to stochastic networks of the associating learning rule for deterministic models, introduced by Personnaz et al in ^[78].

6.2.   Comparison with previous analytical solutions of finite-size networks

Unlike previous mathematical techniques introduced for studying finite-size effects in neural network models, our analytical results are exact for arbitrary network size (also small ones), in that they do no rely on any mathematical approximation. In contrast to linear response theories ^{[39,40,41,42]}, our approach works far away from stationary states and can be applied also to the study of neural oscillations. Compared to the method introduced in ^[43], which represents a mean-field description for the mesoscopic dynamics of interacting populations of integrate-and-fire neurons, our technique does not rely on any mean-field approximation and it can be applied to networks of arbitrary size. Moreover, while the cumulant expansion introduced in ^[44] for networks of binary neurons yields an approximate analytical description of correlations between the activities of individual pairs of units, which goes beyond averaged pairwise correlations, our approach calculates exact group-wise correlations among an arbitrary number of neurons (see ^[14]). An alternative technique for the calculation of the correlation structure in finite-size networks of binary neurons was introduced in ^[79]. This approach assumes that the fluctuations of neuronal activity are sufficiently small and it can be applied only to asynchronous and irregular network states. In ^[80] the authors describe neuronal dynamics by means of a spike-response model with escape noise, and they perform a linear-noise approximation of the equations which allows them to investigate the theoretical properties of finite-size activity fluctuations in the asynchronous regime. We observe that, unlike ^[79,80], our approach can describe both asynchronous and synchronous states with arbitrary levels of correlations between neurons. In ^[81] the authors consider a network composed of phase neurons, and they develop a perturbation theory in the inverse system size $N^{-1}$. Consequently, this approach is expected to be inaccurate for small network sizes, while our technique provides exact results also for arbitrarily small networks. The approach introduced in ^[82] applies to large but finite networks of sparsely connected rate-based neurons. Unlike this technique, our results are not restricted to sparse networks only, and they provide exact formulas of neural activity fluctuations also in highly dense networks.

To conclude, it is important to observe that while our approach to the study of finite-size effects in recurrent neural networks presents, as we discussed above, several advantages compared to previous work, our network equations describe binary firing rates that evolve in discrete-time steps. On the other hand, most of the techniques cited above deal with continuous-time equations of spiking neurons, which are a more realistic description of biological networks. Moreover, our derivation of the exact solutions of the binary equations comes with a cost, namely the exponential increase of complexity of our results with the network size. For this reason, our technique has a higher computational load compared to other methods, as we discussed in more detail in Subsection (6.1.2).

6.3.   Open problems and future directions for mathematical developments

6.3.1.   Dense networks

While in ^[15] we proposed an efficient solution for performing the bifurcation analysis of binary networks with sparse connectivity, fast algorithms for dense networks proved more difficult to develop. In ^[15] we showed that, in the specific case of homogeneous networks with regular topologies, it is possible to take advantage of the symmetries of the network equations to speed up the calculation of the bifurcation diagram. This observation allowed us to introduce an algorithm which runs in linear time with respect to the network size. However, the development of efficient algorithms for dense networks with arbitrary topology of the synaptic connections represents a much more difficult challenge, and it still remains an open problem to be addressed in future work.

6.3.2.   Large-size networks

Another open problem is represented by the bifurcation analysis of large and medium-size networks. Since computer's processing power increases over time, the size of the networks that can be studied through combinatorial approaches, such as the one we introduced in ^[15], is expected to increase accordingly. For this reason, we believe that these methods are key to the future development of computational neuroscience and of the physics of complex systems. Beyond brute-force processing power, other techniques can be developed for accelerating combinatorial algorithms; in particular, the algorithm that we developed in ^[15] lends itself to be parallelized over several processors. It also is important to note that our algorithm calculates exact bifurcation diagrams, while the calculation of approximate bifurcation diagrams, through a heuristic search of the oscillatory and stationary solutions of the network equations, would prove much faster.

6.3.3.   Noisy networks

In section (3) we calculated the bifurcation points of the network in the zero-noise limit. An interesting extension of our formalism, that deserves further investigation, is the study of bifurcations in noisy networks. Note that in Subsection (4.3) we studied the conditions under which the stored sequences of activity patterns are disrupted by noise. For this reason, it is natural to speculate that a theory of bifurcations in noisy networks of binary units must be closely related to the theory we developed in that subsection.

6.3.4.   Neural oscillations

In ^[83] the authors introduced a technique for quantifying the average number of oscillations in asymmetric fully-connected binary networks with quenched disorder. Their approach can be applied to networks in the large-size limit, provided the systems have no autapses (i.e. self-connections) and no external input currents. It is important to observe that in principle the techniques developed in ^[16] for studying the stationary states in networks with quenched disorder can be extended to neural oscillations, by applying order statistics twice to Eq. (3.1). This approach would provide a theoretical estimation of the average number of oscillations which, unlike that derived in ^[83], is not restricted to large fully-connected networks.

6.3.5.   Correlated synaptic weights

To conclude, an important open problem in the study of networks with quenched disorder, is represented by models with correlated synaptic weights. In section (5), we calculated the probability distribution of the bifurcation points of the stationary states, through the results derived in ^{[62,63,64,65]} for the order statistics of a set of independent random variables. On the other hand, it is known that synaptic weights in real cortical circuits are correlated, as a consequence, for example, of synaptic plasticity mechanisms. However, a generalization of order statistics to sets of arbitrarily correlated random variables is still out of reach. More generally, the analytical investigation of networks with correlated synaptic weights has proven a formidable problem in mathematical neuroscience, that has challenged also the theories of large-size systems ^[84,85]. Because of its biological relevance, the presence of synaptic correlations represents an important ingredient in the study of neural network dynamics, which needs to be addressed in future research for increasing the biological plausibility of the models.

6.4.   Possible implications of this mathematical progress for neuroscience

Being able to understand analytically the dynamics of finite-size networks from the circuit's equations is potentially important for improving our understanding of how neural circuits work and of how and why their function is impaired in certain neural disorders. The pattern and power of the external inputs, the pattern of anatomical connectivity, the synaptic strength and the relative firing rate of excitatory and inhibitory neurons, are all key elements in determining the functional organization and output of a circuit. Yet, it is not known how exactly these factors combine to produce brain functions and to cause dysfunctions. Mathematical work to understanding neural networks of arbitrary size could be useful to address two questions relevant to these issues.

The first question regards the relationship between anatomy, functional coupling and population coding in local neural circuits. Our own work (e.g. ^{[86,87,88,89]}), and that of many others ^{[90,91,92,93,94,95]}, has shown that both the dynamics of individual neurons and of populations, and the functional coupling between cells, is crucial for shaping information in population codes and for behavior. Functional coupling between cells may arise both because of anatomical connectivity between neurons, but also because of other factors such as common inputs. Thus, the relationship between functional coupling, circuit's anatomy, and population-level information coding has remained largely unaddressed. However, recent advances in experimental techniques allow the simultaneous functional imaging of activity of several neurons in mice during sensory or cognitive tasks, as well as the post-mortem measure by Electron-Microscopy of the anatomical connectivity of the same set of neurons that were functionally imaged in vivo ^[96]. The mathematical work reviewed here develops a set of tools that could be used to complement the measure of anatomy and physiology from the same circuits, and help bridging the gap between these two measures. Our tools, when coupled with modern experimental techniques such as those described above, could be used, in particular, to understand what is the consequence of specific patterns of recurrent anatomical connectivity on population coding and circuit dynamics, and then to test these theoretical relationships on real data.

A second possible direction of relevance for neuroscience of our work is to use these tools to understand the neural origin of certain brain disorders. For example, it is thought that Autism Spectrum Disorders (ASD) result, at least in part, from abnormal changes in the functional organization and dynamics of neural circuits ^{[97,98,99,100]}. However, although many changes in parameters such as the strength of synaptic connections and/or the change in firing properties of certain classes of neurons have been observed in ASD, it is still unclear how different elementary changes in neural parameters combine to change the circuit's function. Being able to test directly, and understand mathematically at a deep level, how the changes of such basic neural properties affect the circuit's dynamics at different spatial scales, including scales that involve finite-size neural networks, could be useful to understand the origin and consequences of aberrant circuit function in ASD conditions, as well as in other neural disorders.

Acknowledgments

This research was in part supported by the Simons Foundation (SFARI, Explorer Grant No. 602849) and by the BRAIN Initiative (Grant No. R01 NS108410). The funders had no role in study design, data collection and analysis, decision to publish, interpretation of results, or preparation of the manuscript.

Conflict of interest

The authors declare no conflict of interest.