## Abstract

Despite their biological plausibility, neural network models with asymmetric weights are rarely solved analytically, and closed-form solutions are available only in some limiting cases or in some mean-field approximations. We found exact analytical solutions of an asymmetric spin model of neural networks with arbitrary size without resorting to any approximation, and we comprehensively studied its dynamical and statistical properties. The network had discrete time evolution equations and binary firing rates, and it could be driven by noise with any distribution. We found analytical expressions of the conditional and stationary joint probability distributions of the membrane potentials and the firing rates. By manipulating the conditional probability distribution of the firing rates, we extend to stochastic networks the associating learning rule previously introduced by Personnaz and coworkers. The new learning rule allowed the safe storage, under the presence of noise, of point and cyclic attractors, with useful implications for content-addressable memories. Furthermore, we studied the bifurcation structure of the network dynamics in the zero-noise limit. We analytically derived examples of the codimension 1 and codimension 2 bifurcation diagrams of the network, which describe how the neuronal dynamics changes with the external stimuli. This showed that the network may undergo transitions among multistable regimes, oscillatory behavior elicited by asymmetric synaptic connections, and various forms of spontaneous symmetry breaking. We also calculated analytically groupwise correlations of neural activity in the network in the stationary regime. This revealed neuronal regimes where, statistically, the membrane potentials and the firing rates are either synchronous or asynchronous. Our results are valid for networks with any number of neurons, although our equations can be realistically solved only for small networks. For completeness, we also derived the network equations in the thermodynamic limit of infinite network size and we analytically studied their local bifurcations. All the analytical results were extensively validated by numerical simulations.

## 1 Introduction

In networks of neurons in real biological brains, the vast majority of synapses are asymmetric, meaning that the synapse from a neuron to another has a different strength from the synapse in the reverse direction (DeFelipe, Marco, Busturia, & Merchán-Pérez, 1999). In many cases, a synapse may be present between one neuron and another along one direction but not in the other direction. Because of this, the most biologically relevant mathematical models of neural networks are those with asymmetric synaptic weights. These models are widespread and have been intensively studied numerically (see, e.g., Bray & Moore, 1984; Mattia & Del Giudice, 2000; Amit & Mongillo, 2003; Markram, 2006; Morrison, Straube, Plesser, & Diesmann, 2007; Brette et al., 2007; Paik, Kumar, & Glaser, 2009; Yucesoy, Katzgraber, & Machta, 2012; Baños et al., 2012; Sanz Leon et al., 2013; Wang, Machta, & Katzgraber, 2015; Zhu, Ochoa, & Katzgraber, 2015). Yet neural network models with asymmetric synaptic connections that can be mathematically solved either exactly in closed form or through approximations are rare and typically are made of model neurons with discrete, spin-like outputs (see, e.g., Derrida, 1980; Crisanti & Sompolinsky, 1988; Roudi & Hertz, 2011; Sakellariou, Roudi, Mezard, & Hertz, 2012; Dahmen, Bos, & Helias, 2016). Developing closed-form analytical solutions of such neural network models is important because asymmetric networks typically exhibit a rich set of oscillatory behavior that does not occur in models with symmetric connections. In symmetric networks, only a limited set of oscillations can occur, and only if the network is synchronously updated (Goles-Chacc, Fogelman-Soulie, & Pellegrin, 1985) or if time delays are present (Marcus & Westervelt, 1989). Closed-form solutions of asymmetric models would thus make it possible to thoroughly study the formation and the functional role of a wide range of neural oscillations. In particular, a full analytical understanding of asymmetric neural network models would help us to better understand how the emergence of neural oscillations depends on the network parameters in systems composed of excitatory and inhibitory populations, to understand how oscillations can be stored in the neural tissue, and to study how oscillations influence the integration of information among neurons.

The added difficulty in studying asymmetric neural networks is the impossibility of applying the powerful methods of equilibrium statistical mechanics because no energy function exists for these systems. Due to their complexity, in general, models with spin-like discrete neurons admit exact analytical solutions only for specific network topologies and distributions of the synaptic weights, and in some limiting cases or mean-field approaches, such as the replica trick (Mézard, Parisi, & Virasoro, 1986), the limit of an infinite number of spin components (Hertz, Grinstein, & Solla, 1987), or the thermodynamic limit of infinite network size (Crisanti & Sompolinsky, 1987; Derrida et al., 1987). Yet many biologically important neural networks have a finite or even small size, such as microcolumns in the primate cortex (Mountcastle, 1997) or the nervous system of some invertebrates (Williams & Herrup, 1988). Thus, it is desirable for methods to mathematically study networks with asymmetric connections to be able to analytically solve the dynamics of networks of arbitrary size, including small sizes.

Here, we make progress on these issues by showing how a spin-like (discrete output) model of neural networks with asymmetric connections can be solved exactly without resorting to a limiting case or a mean-field approximation. In particular, our solutions are valid for networks composed of an arbitrary number of neurons, do not require additional constraints (such as the synaptic dilution assumption adopted in Derrida, Gardner, & Zippelius, 1987), and can be derived in the presence of noisy inputs with any statistics.

Our approach is based on calculating exact analytical solutions for the conditional probability distributions of the membrane potentials and the firing rates, as well as for the joint probability distributions in the stationary regime (see section 3.1). The ability to study the network analytically enabled us to define a new learning rule that safely stores stationary and oscillating solutions even in the presence of noise (see section 3.2), study the qualitative changes in network dynamics as a function of the input parameters by means of bifurcation theory (see section 3.3), and compute experimentally observable quantities such as correlations among the activity of different neurons (see section 3.4). While all the results are valid for networks composed of an arbitrary number of neurons, for completeness, we also derived the mean-field equations in the thermodynamic limit of infinite network size, and we found exact analytical expressions of the local codimension 1 bifurcations (see section 3.5).

## 2 Materials and Methods

### 2.1 The Stochastic Firing-Rate Network Model

In this section, we describe the neural network model we study.

The activity of each neuron in the network is influenced (through synapses that may have asymmetric weights) by the activity of other neurons in the network, as well by an external input that has both a deterministic and a noisy component. For clarity and simplicity, in the main text, we study and solve a model in which neurons are driven by independent noise sources with gaussian distribution. In the supplemental information, we show how to solve the model with an arbitrary noise distribution.

To conclude, we observe that in the special case when $^J$ is symmetric with zeros on the main diagonal ($^Jii=0$) and $Ii=\sigma iB=0\u2200i$, equation 2.5 describes the synchronous version of the network model introduced by Hopfield (1982). Note, however, that when the connections are asymmetric, our model will exhibit important differences in dynamics and behavior with respect to the Hopfield model, as we show in section 3.

### 2.2 Numerical Simulations

To further validate our analytical results, we compared them with numerical simulations (see section 3, where all the results in the figures are expressed in simulation units).

The joint probability distributions of the membrane potentials and the firing rates in the stationary regime have been obtained through a Monte Carlo method by iteratively solving equation 2.1 in the temporal interval $t=0,100$. Then the probability distributions have been evaluated from the collected data at $t=100$ through a kernel density estimator. We assume that at $t=100$, the probability distributions have already reached a stationary regime, which is confirmed by the good agreement between the analytical and numerical results. The network statistics have been calculated numerically in a similar way, by means of the sample variance and sample (distance) correlation over the Monte Carlo repetitions.

To conclude, the bifurcation diagrams of the network in the thermodynamic limit have been obtained numerically by the MatCont Matlab toolbox (Dhooge, Govaerts, & Kuznetsov, 2003).

## 3 Results

In this section, we study the dynamical and statistical properties of the network defined above. In section 3.1, we report the exact solutions of the conditional probability distributions of the membrane potentials and the firing rates, as well as the joint probability distributions in the stationary regime. In section 3.2, we invert the formula of the conditional probability distribution of the firing rates, which allows us to define a new learning rule to safely store stationary and oscillatory patterns of neural activity even in the presence of noise. In section 3.3, we study analytically the bifurcations of the network dynamics in the zero-noise limit $\sigma B\u21920$, in terms of its codimension 1 and codimension 2 bifurcation diagrams. In section 3.4, we derive exact analytical expressions for the groupwise cross-correlations of the membrane potentials and the firing rates for any noise intensity. Finally, in section 3.5, we report the mean-field equations of the network, which are exact in the thermodynamic limit $N\u2192\u221e$, and we analytically study its local bifurcations. All the analytical results are extensively validated by numerical simulations.

### 3.1 Conditional and Joint Probability Distributions

$J=040-3660-361040-4032-404080040-85260-560-843664-44480,I=-10-220,\sigma B=21123$ |

$J=040-3660-361040-4032-404080040-85260-560-843664-44480,I=-10-220,\sigma B=21123$ |

We also observe that while equations 3.1 to 3.3 represent probability density functions (pdfs), namely, probability distributions of continuous random variables (the membrane potentials), equations 3.4 to 3.6 represent probability mass functions (pmfs), namely, probability distributions of discrete random variables (the firing rates). In particular, $P\nu |\nu '$ is known as a state-to-state transition probability matrix in the context of the theory of Markov processes (Feller, 1971). Examples of $P\nu $ for $N=5$ and $N=2$ are shown in Figures 2D and 3D, which report the comparison between the analytical and numerical probability distribution of the firing rates over the $25=32$ and $22=4$ states of the network, respectively.

It is important to observe that while the synchronous version of the Hopfield network (which has symmetric connections) can sustain only oscillations with period $T=2$ (known as two-cycles; Goles-Chacc et al., 1985), Figure 3 (see also Table 1) shows that the asymmetric network can undergo oscillations with larger period. In this figure, the oscillations are perturbed by the presence of the noisy terms $\sigma iBBit$, so that the network switches randomly between an oscillation with $T=3$, $00\u219210\u219201\u219200$, and an oscillation with $T=4$, $00\u219210\u219211\u219201\u219200$. In particular, Figure 3A shows an example of oscillatory temporal evolution of the stochastic membrane potentials, while Figure 3B represents the corresponding dynamics in the phase space of the network. Then, Figure 3C describes the formation of neural oscillations in terms of the state-to-state transition probability matrix.

### 3.2 A New Learning Rule for Storing Point and Cyclic Attractors

At time $t$, the state of the neural network is described by the vector of the firing rates, $\nu t=[\nu 0t,\u2026,\nu N-1t]T$, which represents the activity pattern of the system at time $t$. In the context of content-addressable memories, one aims to determine a synaptic connectivity matrix $J$ that stores one or more desired sequences of activity patterns. The way such a matrix is built defines a learning rule for storing these patterns.

In particular, we observe that oscillations of period $T$ correspond to the special case $Li=T$ with $\nu it=\nu it+T$. For $T=1$, this condition allows us to store stationary states. Examples of $M$-stable systems obtained through this method are shown in Figure 4 (see also Table 3) for $M=2,3,4$, while examples of oscillatory dynamics with period $T=2,3,4$ are shown in Figure 5 (see also Table 4).

$J=0-3993-39-840-84171123-660-66-1554-150,J=0-1234848-51001711713781230-25513854-840,J=09393-1711710-339171123-2550123-8454540$ |

$I=-212-3,\sigma B=2431$ |

$J=0-3993-39-840-84171123-660-66-1554-150,J=0-1234848-51001711713781230-25513854-840,J=09393-1711710-339171123-2550123-8454540$ |

$I=-212-3,\sigma B=2431$ |

Note: The synaptic connectivity matrices $J$ have been obtained from equation 3.7 for $Kji,ni=10\u2200i,j,ni$.

$J=0306330-840-843398442042-1584-150,J=0-93939317103-339123-25503-303840,J=0-8493-84-1710171-1714242042-42-42540$ |

$J=0306330-840-843398442042-1584-150,J=0-93939317103-339123-25503-303840,J=0-8493-84-1710171-1714242042-42-42540$ |

Note: The synaptic connectivity matrices $J$ have been obtained from equation 3.7 for $Kji,ni=10\u2200i,j,ni$.

Learning rule equation 3.7 can be easily extended by including further constraints. For example, it is possible to relax the full-connectivity assumption and set a given portion of the synaptic weights to zero. This allows us to define a learning rule for sparse neural networks. The storage capacity is expected to depend on the network topology and the noise intensity, but a detailed investigation is beyond the purpose of this letter.

To conclude, we observe that in order to store the $D$ sequences of patterns, according to equation 3.7, we need to invert $N$ linear systems of size $L\xd7N-1$. Therefore, the learning rule has polynomial complexity in the network size $N$.

### 3.3 Bifurcations in the Deterministic Network

Bifurcation analysis is a mathematical technique for investigating the qualitative change in the neuronal dynamics produced by varying model parameters. Therefore, it represents a fundamental tool for performing a systematic analysis of the complexity of the neuronal activity patterns. In particular, in this section, we study how the dynamics of the network depends on the external stimuli $Ii$. Moreover, we perform our analysis in the zero-noise limit $\sigma B\u21920$, as is common practice in bifurcation theory (see, e.g., Borisyuk & Kirillov, 1992; Grimbert & Faugeras, 2006; Haschke & Steil, 2005; Fasoli, Cattani, & Panzeri, 2016). We observe that local bifurcation analysis in graded (i.e. smooth) models is performed through the calculation of the eigenvalues of the Jacobian matrix (Kuznetsov, 1998). However, the Jacobian matrix of our model is not defined at the discontinuity of the activation function, equation 2.2, thus preventing the use of the powerful methods of bifurcation analysis developed for graded models. The bifurcation analysis of nonsmooth dynamical systems has been developed only for continuous-time models, described by nonsmooth differential equations or differential inclusions (Leine, Van Campen, & Van De Vrande, 2000; Awrejcewicz & Lamarque, 2003; Leine & Van Campen, 2006; Makarenkov & Lamb, 2012). For this reason, we need to define an alternative approach for studying bifurcations in our binary network model.

This approach requires prior knowledge of the stationary and oscillatory solutions of the network. One possibility is to determine these solutions numerically, for example, by solving, for $\sigma iB=0$, equation 2.1 iteratively for all the $2N$ initial conditions of the firing rates and for all the combinations of the currents $I0,\u2026,IN-1$ on a sufficiently dense discretization of the stimulus space. Another possibility is to perform a numerical calculation of the conditional probability distribution (similarly to Figures 4 and 5), from which the stationary and oscillatory solutions can be detected through a search of the simple cycles of length $L$ of the $2N\xd72N$ binary matrix $P\nu |\nu '$. In this approach, the stationary states correspond to loops of the matrix $P\nu |\nu '$ (i.e., to cycles of length $L=1$). In a similar way, oscillations correspond to simple cycles of length $L=T>1$, where $T$ represents the period of the oscillation. More efficient techniques are considered in Fasoli and Panzeri (2017).

The bifurcation diagrams of the network strongly depend on its connectivity matrix $J$. However, a detailed analysis of the relation between the bifurcation structure and the network topology is beyond the purpose of this letter. For an example, we apply our method to a fully connected network composed of $NE=3$ excitatory and $NI=3$ inhibitory neurons, even though this technique can be easily employed for calculating the bifurcation diagrams of networks with any size and topology. We suppose that each excitatory (resp. inhibitory) neuron receives an external stimulus $IE$ (resp. $II$), and we derive the codimension 2 bifurcation diagram of the network in the $IE-II$ plane. The remaining parameters of the network are reported in Table 5.

$\theta =1\vdots 1,J=JEEJEIJIEJII,J\alpha \alpha =J\alpha \alpha 011101110,J\alpha \beta =J\alpha \beta 111111111for\alpha \u2260\beta $ |

$JEE=-JII=80,JIE=-JEI=70$ |

$\theta =1\vdots 1,J=JEEJEIJIEJII,J\alpha \alpha =J\alpha \alpha 011101110,J\alpha \beta =J\alpha \beta 111111111for\alpha \u2260\beta $ |

$JEE=-JII=80,JIE=-JEI=70$ |

By representing the firing rates of the excitatory neurons through the three top entries of the vector $\nu $, we found numerically that the network admits the stationary states 0 to 7 and 56 to 63 (in decimal representation) for particular combinations of $IE,I$. Thus, for example, the state $\nu =[0,0,0,1,0,1]T$ (i.e., the state 5 in decimal representation), which is characterized by two active inhibitory neurons (while the remaining neurons in the network are not firing), is a stationary state for some values of the stimuli. Moreover, we found numerically that the network undergoes the oscillations $0\u21927\u21920$, $56\u219263\u219256$, $0\u219256\u219263\u21920$, $0\u219263\u21927\u21920$, and $0\u219256\u219263\u21927\u21920$. Now, by inverting the conditions provided by the conditional probability distribution 3.8, we get the portions of the $IE-II$ plane where each stationary state and oscillation occurs. The details of the analytical calculations are shown in sections S4.1 and S4.2 of the supplemental information for the stationary and oscillatory solutions, respectively. The resulting codimension 2 bifurcation diagram is shown in Figure 6 (see also Table 5).

^{7}$IE,I$ is calculated analytically in section S4 of the supplemental information. In particular, we observe that in the excitatory population, the stationary firing rates (the three top entries of the binary representation of the states 0 to 7 and 56 to 63) and the corresponding membrane potentials are homogeneous ($\nu i=\nu EIE,II$ and $Vi=VEIE,II$ for $i=0,1,2$), while in the inhibitory population, they are generally heterogeneous ($\nu NE+i=\nu I,iIE,II$ and $VNE+i=VI,iIE,II$ for $i=0,1,2$). This is an example of spontaneous symmetry breaking that occurs in the inhibitory population (see the shaded areas in Figure 6A). On the contrary, according to the numerical solutions, there is no spontaneous symmetry breaking during the oscillatory dynamics in this example (see Figure 6B); therefore, in this case, the firing rates are homogeneous in both the neural populations. Equation 3.9 is plotted in Figures 7A and 7B (see also Table 5). Figures 7C to 7F show the variations of the state-to-state transition probability matrix caused by the bifurcation points of the model. In particular, the panels show that for $II=-20$ and for an increasing value of $IE$, the network exhibits monostability (C), quadristability (D), tristability and an oscillation with period $T=3$ (E), and tristability and an oscillation with period $T=2$ (F). Compare Figures 7C to 7F with the four ranges of the stimulus $IE$ delimited by the vertical dashed lines in Figures 7A and 7B (for a derivation of the boundaries of these ranges, see section S4 in the supplemental information).

### 3.4 Groupwise Cross-Correlations

The study of correlations among neurons is a topic of central importance to systems neuroscience. Pairwise and groupwise cross-neuron correlations are key to understanding the information encoding capabilities of neural populations (Abbott & Dayan, 1999; Pola, Thiele, Hoffmann, & Panzeri, 2003; Pillow et al., 2008; Cohen & Maunsell, 2009; Moreno-Bote et al., 2014) and to making inferences about how neurons exchange and integrate information (Singer, 1993; Tononi, Sporns, & Edelman, 1994; David, Cosmelli, & Friston, 2004; Rogers, Morgan, Newton, & Gore, 2007; Friston, 2011).

$J=0-1545-12015090-28590-1050-3165-60750,I=I1111,\sigma =1322$ |

$J=0-1545-12015090-28590-1050-3165-60750,I=I1111,\sigma =1322$ |

Note: The values of the stimulus $I$ are specified in the figures.

Examples of cross-correlations for $n>2$ are shown in Figure S3 of the supplemental information. Our analysis shows that the external stimulus $I$ dynamically switches the network between synchronous (i.e., highly correlated) and asynchronous (i.e., uncorrelated) states. The conditions under which these states may occur are discussed in section S5.1.1, S5.1.2, and S5.2.1 of the supplemental information. In general, we did not observe a fixed relation between firing rate and membrane potential correlations, and low (resp. high) correlations between the membrane potentials did not necessarily correspond to low (resp. high) correlations between the firing rates. In Figure 9, we consider the network reported in Table 6 for $I=-3$, and we show an example of state characterized by highly correlated membrane potentials and low correlated firing rates. The membrane potentials of neurons 0 and 1 fluctuate randomly around the same mean, which for $I=-3$ is close to the firing threshold $\theta =1$ (represented by the dashed line in the figure). The sources of noise have different standard deviations ($\sigma 0B=1$, $\sigma 1B=3$), so that the membrane potentials fluctuate with standard deviations $\sigma 1V>\sigma 0V$. For this reason, the membrane potential of neuron 1 is more likely to cross the firing threshold than the potential of neuron 0. Therefore, the firing rate of neuron 1 switches often and randomly between the rate states 0 (neuron not firing) and 1 (neuron firing at the maximum rate). On the contrary, neuron 0 does not fire very often; therefore, the correlation between the firing rates is small even though the corresponding membrane potentials are highly correlated.

### 3.5 Mean-Field Limit

$JEE=-JII=80,$ | $\theta E=\theta I=1$ |

$JIE=-JEI=70,$ | $\sigma EB=1,\sigma IB=2$ |

$JEE=-JII=80,$ | $\theta E=\theta I=1$ |

$JIE=-JEI=70,$ | $\sigma EB=1,\sigma IB=2$ |

Unlike the original network equations, 2.1, the activation functions in the mean-field equations 3.16 are continuous and differentiable everywhere for $\sigma B>0$. Networks with real-valued activation functions such as that described by equation 3.16 (and the corresponding activity-based versions, which can be obtained through the change of variables, equation 2.3) are known in neural modeling as discrete-time recurrent neural networks (Jin, Nikiforuk, & Gupta, 1994; Tsoi & Back, 1997; Haschke & Steil, 2005). Thus, when noise is present in every population, the bifurcation structure of equation 3.16 can be studied through the bifurcation theory of graded systems (Kuznetsov, 1998). For an example, we focus on the case of a network composed of two neural populations, one excitatory and one inhibitory (even though the bifurcation structure may be studied for every $P$). From now on, it is convenient to change the notation slightly and consider $\alpha =E,I$ rather than $\alpha =0,1$. Beyond global bifurcations, which generally cannot be studied analytically, the mean-field network may undergo limit point, period doubling, and Neimark-Sacker local bifurcations. Limit points correspond to the formation of network bistability. We observe that in this model, bistability coexists with hysteresis phenomena, which were proposed as a candidate neurophysiological mechanism of working memory (Wang, 2001). Period doubling and Neimark-Sacker bifurcations correspond to the formation of neural oscillations with period $T=2$ and $T>2$, respectively.

## 4 Discussion

### 4.1 Comparison with Previous Analytical Solutions of Networks with Asymmetric Synaptic Weights

We studied a synchronously updated firing-rate neural network model with asymmetric synaptic weights and discrete-time evolution that can be exactly solved for any network size and topology. Network models with asymmetric synaptic weights are important because real neural networks in the brain have asymmetric synapses. Exact analytical study, though, has been limited to a small number of studies valid in some limiting cases or in some mean-field approaches, such as the replica trick (Mézard et al., 1986), the limit of an infinite number of spin components (Hertz et al., 1987), or the thermodynamic limit of infinite network size (Crisanti & Sompolinsky, 1987; Derrida et al., 1987).

The main mathematical difficulty in studying asymmetric neural networks is the impossibility to apply the powerful methods of equilibrium statistical mechanics because no energy function exists for these networks. Our formalism allows exact solutions of such network models without resorting to a limiting case or a mean-field approximation. In particular, our results are exact for any network size. We used the stochastic recurrence relation, equation 2.1, rather than Little's (1974) definition of temperature (in particular, compare our equation 3.6 with Little's equation 4). In this respect, our work follows the approach described in Wang, Pichler, and Ross (1990), but with some important differences. Wang et al. (1990) specifically considered a diluted network with gaussian noise and studied the dynamical evolution of the overlap between the state of the network and the stored point attractors in the thermodynamic limit. In our work, we considered a network with an arbitrary (generally nondiluted) synaptic connectivity matrix and arbitrary noise statistics, and in section S2 of the supplemental information, we derived exact solutions for any network size $N$ without further assumptions.

In this letter (see section S1 of the supplemental information), we showed the relationship between our approach and the mathematical techniques developed for Markovian systems (Coolen, 2001a, 2001b). In particular, we showed that the original voltage-based equation 2.1 can be transformed into the activity-based equation 2.5, which can be solved with standard methods from the theory of Markovian processes. This approach requires the normalized connectivity matrix $^J$ (defined in equation 2.4) to be invertible. Our approach generalizes this previous work in that it does not rely on this transformation, and therefore it can be applied to networks with any matrix $^J$.

In particular, through this technique, we derived exact solutions for the conditional probability distributions of the membrane potentials and the firing rates, as well as for the joint probability distributions in the stationary regime. Due to the asymmetry of the synaptic weights, the network we studied can undergo oscillations with period $T\u22652$, while synchronous Hopfield networks (which have symmetric connections) can sustain only oscillations with period $T=2$, known as two-cycles (Goles-Chacc et al., 1985). Going from activity-based to voltage-based equations through equation 2.3 requires the hypothesis of synchronous updates of the neural activity. For this reason, we cannot extend the formalism in a straightforward way to asynchronous dynamics.

Compared to small-size graded networks (Fasoli et al., 2015, 2017), where the impossibility to using statistical methods restricts the derivation of (approximate) analytical formulas of the joint probability distributions and of the cross-correlations to simple network topologies, here we derived exact solutions that are valid for any connectivity matrix $J$.

### 4.2 Implications of Analytical Solvability for Learning Rules

The analytical solution of the conditional probability distribution of the firing rates allowed us to define a new learning rule to store point and periodic attractors. Point attractors correspond to stable states, while periodic attractors represent oscillatory solutions of the network activity. The synchronously updated Hopfield network (which can be obtained from equation 2.1 in the special case of symmetric connections) can store only point attractors and two-cycles, while our learning rule for asymmetric networks can be used to store any oscillatory solution or temporal pattern that is compatible with the network topology and size. From a biological point of view, learning specific temporal patterns of neural activity is a general problem involved in several cognitive and sensory-motor tasks (Grossberg & Pearson, 2008), and therefore it represents an essential ingredient in advancing our understanding of cortical activity.

Learning rules for deterministic networks are important theoretical tools for investigating pattern storage in idealized, noise-free models (Hagan et al., 1996). However, cortical activity is often characterized by the presence of noise (Deco, Rolls, & Romo, 2009). This underlines the importance of developing new learning rules to safely store activity patterns in a noisy medium such as the cortical tissue. The rule that we introduced belongs to this class of learning algorithms for stochastic networks. Therefore, our rule ensures that a stored sequence of activity patterns will not be disrupted by the noise, and it does so by setting the membrane potentials far away from the firing thresholds during their temporal evolution. This distance is proportional to the noise intensity and ensures that undesired random switching of the firing rates will be very unlikely (compare with Figure 3B, where the membrane potentials fluctuate near the firing threshold, causing the random switching of the firing rates described by Figure 3C).

For simplicity, in section 3.2, we showed the learning rule for a fully connected network driven by gaussian noise; however, our derivation can be easily extended to networks with any topology of the synaptic connections and driven by noise with any distribution. The storage capacity depends on the topology of the network (random or regular), the sparseness of the synaptic connections, and the intensity and distribution of the noise, but a detailed investigation is beyond the purpose of this letter.

Several learning rules have been proposed for storing activity patterns in networks of model neurons with discrete, spin-like output. One of the first studies was proposed by Amari (1972), who analyzed the effects of random spin flips on Hebbian learning of activity patterns. Sompolinsky and Kanter (1986) introduced an infinite-size network where the response of the neurons is not instantaneous but has a dynamic memory that slows the transition between patterns. This ad hoc assumption proves sufficient for generating arbitrarily long pattern sequences by preventing the network state to become progressively mixed among all the activity patterns to be stored. A similar model was introduced independently by Kleinfeld (1986). A different approach was followed by Buhmann and Schulten (1987), who proved that noise and sufficiently asymmetric connections can induce random transitions between patterns. Interestingly, our learning rule does not rely on any of these assumptions, since it can also be used for storing pattern sequences if the neurons have instantaneous response and the network is deterministic. Our result can be considered an extension of the associating learning rule introduced by Personnaz, Guyon, and Dreyfus (1986) which is also known as the pseudoinverse learning rule (Zhang, Dangelmayr, & Oprea, 2013). While Personnaz and coworkers derived the learning rule by making implicit use of the state-to-state transition probability matrix of their deterministic network model (see equation 10 in Personnaz et al., 1986), we generalized this approach by taking advantage of the transition probability matrix described by equation 3.4 to define a new learning rule for storing arbitrary sequences of activity patterns in stochastic networks.

### 4.3 Bifurcations in Neural Activity

We performed an analytical study of the bifurcations to understand how the formation of stable and oscillatory solutions is affected by a variation of the external stimuli to the network. The method we proposed can be applied to networks with any size and topology, but for an example, we considered the case of a network composed of six fully connected neurons. As is common practice, we performed the bifurcation analysis in the zero-noise limit. We derived analytical expressions for the codimension 1 and codimension 2 bifurcation diagrams, showing how the external stimuli affect the neuronal dynamics. It is important to observe that the bifurcation analysis of graded (i.e., smooth) networks is based on differential analysis—in particular, on the Jacobian matrix of the system. However, the Jacobian matrix of the model we studied in this letter is not defined everywhere due to the discontinuity of the activation function. For this reason, bifurcation theory of smooth dynamical systems cannot be applied directly to the binary model we discussed. Although the bifurcation analysis of nonsmooth dynamical systems has recently received increased attention, it has been developed only for continuous-time models, described by nonsmooth differential equations or differential inclusions (Leine et al., 2000; Awrejcewicz & Lamarque, 2003; Leine & Van Campen, 2006; Makarenkov & Lamb, 2012). For this reason, in our work we followed another approach. By combining the analytical formula of the conditional probability distribution of the firing rates with a numerical brute force search of the stationary and oscillatory solutions, we determined the combinations of the external stimuli for which the network undergoes multistability, oscillations, or symmetry breaking. In particular, for a fully connected network composed of one excitatory and one inhibitory population in a stationary state, we observed the formation of spontaneous symmetry breaking only in the inhibitory population. Intuitively, this phenomenon is similar to the magnetic spin ordering (i.e., neighboring spins pointing in opposite directions) observed in antiferromagnetic materials in a weak external field below their Néel temperature (Néel, 1952). In antiferromagnetic materials, the interaction weights between spins are negative, as in the inhibitory neural population of our network model. Moreover, when the network is noise free (as in section 3.3), it can be interpreted as a spin model at zero temperature. Thus, the noise-free network we considered is below its Néel temperature, which intuitively explains the formation of an “antiferromagnetic” phase characterized by heterogeneous activity in the inhibitory population. In a similar way, the absence of symmetry breaking observed in the excitatory population is similar to spin alignment observed in ferromagnetic materials (which are characterized by positive interaction weights) below their Curie temperature (Néel, 1952). However it is important to bear in mind that a first difference between the spin systems and our network model is that our model has two distinct populations, whose neurons may interact in more complex ways than purely ferromagnetically or antiferromagnetically. Moreover, while magnetic materials are usually described by models where spins are arranged on a lattice, here we considered a neural network model with fully connected topology.

### 4.4 Analytical Estimates of Correlations of Neural Activity in Asymmetric Networks

Thanks to the analytical solvability of the model, we derived exact expressions of the groupwise correlation structure of noisy networks at any order (beyond simple pairwise correlations often considered in both theoretical and empirical studies of neural correlations) in the stationary regime for both the membrane potentials and the firing rates. In the case of the time-continuous graded version of the network model, Fasoli et al. (2015, 2017) found analytical (approximate) solutions of the correlation structure through a perturbation analysis of the neural equations in the small-noise limit. A consequence of this approximation is that the correlations between the membrane potentials and those between the firing rates have the same mathematical expression in the graded model. On the contrary, in this letter, we derived exact expressions of the correlations for any noise intensity. Due to the discontinuous activation function, equation 2.2, the two correlations structures are never identical, even in the small-noise limit. Moreover, similarly to the case of graded networks (Fasoli et al., 2017), we found that the external stimuli can dynamically switch the neural activity from asynchronous (i.e., uncorrelated) to synchronous (i.e., highly correlated) states, with two important differences. The first is that low (resp. high) correlations between the membrane potentials do not necessarily correspond to low (resp. high) correlations between the firing rates. The second is that while in graded networks synchronous states may occur through critical slowing down (Fasoli et al., 2017), the discrete network considered here relies on different mechanisms for generating highly correlated activity that we have only partially covered.

### 4.5 Behavior of Asymmetric Networks in the Large-Population Limit

For completeness, and to link our results to previous work on asymmetric models, we derived the mean-field equations of the network in the thermodynamic limit of large network size (see equation 3.16). Due to the limitations of Sznitman's mean-field theory, we derived these equations only for sufficiently dense multipopulation networks driven by independent sources of noise. These hypotheses allowed obtaining a dimensionally reduced set of equations to describe neural activity, which proved very convenient in studying the dynamical behavior of the multipopulation network in the thermodynamic limit.

The presence of a noisy component in the external stimulus smoothed out the activation function of the network in the thermodynamic limit. Unlike the case of finite-size networks, this allowed us to study the local bifurcations of the infinite-size network for any nonzero noise intensity by applying standard techniques developed for smooth dynamical systems (Kuznetsov, 1998; Haschke & Steil, 2005; Fasoli et al., 2016). In particular, we derived exact analytical expressions for the local codimension 1 bifurcations of the mean-field model in terms of the external stimuli. This method can be applied to networks composed of an arbitrary number of populations, but for an example, we considered the simple case of two populations. We also derived analytically part of the codimension 2 bifurcation diagram of the network and we compared it with the bifurcation structure obtained numerically by the MatCont Matlab toolbox (Dhooge et al., 2003).

Equation 3.16 can be considered the discrete-time version of the mean-field equations derived in Touboul et al. (2012) for a continuous-time stochastic network composed of an arbitrary number of neural populations. While Touboul et al. (2012) were able to derive simple mean-field equations for a graded network with error-function activation functions, in our work we considered a network with binary firing rates (see equation 2.2). Despite these differences, through the bifurcation analysis tools that we developed in the this letter, it is possible to prove that our network model exhibits similar phenomena to those that Touboul et al. (2012) studied, such as noise-induced formation and destruction of stationary and oscillatory solutions (results not shown).

It is important to observe that the use of Sznitman's mean-field theory is limited to networks driven by independent sources of noise and with sufficiently dense synaptic connectivity. In Fasoli et al. (2015), the authors proved that in a graded firing-rate network model, propagation of chaos does not occur in the thermodynamic limit whenever the sources of noise are correlated or the synaptic connections are too sparse. The same results can be proven for the discrete network described by equation 2.1. The study of infinite-size neural networks driven by correlated sources of noise typically requires more advanced techniques such as large deviations theory, which is discussed in Faugeras and MacLaurin (2015).

We also observe that in Sznitman's mean-field theory, the neurons are identically distributed within each population. For this reason, the theory does not describe intrapopulation symmetry breaking of neural activity, as we already observed in Fasoli et al. (2016) for a graded firing-rate network model. This proves that the dimensionally reduced set of mean-field equations generally oversimplifies the dynamical properties of the neural network; therefore, it may neglect important neurophysiological processes that can be captured only by finite-size network models.

### 4.6 Computational Complexity and Limits of Applicability of the Arbitrary-Network-Size Theory

Although we set out to develop a formalism that can be applied to arbitrary network sizes, the complexity of the various analyses makes them applicable in practice to selected regimes. In the following, we discuss the practical domain of applicability of different aspects of our techniques.

The analytical calculations of the stationary distributions and the cross-correlation structure require an exponentially increasing number of calculations in the network size $N$. These formulas become intractable for networks of a few tens of neurons. However, our theoretical framework provides new insights into the behavior of cross-correlations that are expected to hold for any network size $N$ and therefore are not restricted to small-size circuits. Indeed, several properties of cross-correlations between the membrane potentials or the firing rates can be shown to hold for any $N$ (see section S5 of the supplemental information). Therefore, our analytical solutions complement the knowledge that can be obtained from numerical simulations and provide new insights into the behavior of arbitrary-sized neural networks. In contrast, our learning rule (see section 3.2) has polynomial complexity; therefore, its applicability is not restricted to small-size networks only. Finally, efficient algorithms for the bifurcation analysis of noiseless networks of spin-like neurons in principle may be applicable also to large networks, provided the density of their synaptic connections is small enough (Fasoli & Panzeri, 2017).

### 4.7 Concluding Remarks

To conclude, we observe that solvable finite-size network models are invaluable theoretical tools for studying the brain at its multiple scales of spatial organization. Studying how the complexity of neuronal activity changes with increasing network size is of fundamental importance for unveiling the emergent properties of the brain. In this letter, we made an effort in this direction, trying to fill the gap in the current neuroscientific literature. Asymmetric synaptic connections, which are widely considered as mathematically challenging for a rigorous understanding of the network's dynamics, increase the biological plausibility of the model and allow a more complete description of neural oscillations. In future work, we will rigorously determine how the two main assumptions of the model, the discrete-time evolution and the binary firing rates, affect its capability to describe realistic neuronal activity.

## Acknowledgments

We thank Davide Corti for contributing to the initial stages of this work. 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.

## References

## Author notes

D.F. and S.P. designed the study. D.F. performed the study through analytical calculations and numerical simulations. A.C. performed the numerical study of the bifurcations of the mean-field network model. D.F., S.P., and A.C. wrote the paper.