article: Mean-field approximations are a powerful tool for studying large neural networks. However, they do not describe well the behavior of networks composed of a small number of neurons. In this case, major differences between the mean-field approximation and the real behavior of the network can arise. Yet, many interesting problems in neuroscience involve the study of mesoscopic networks composed of a few tens of neurons. Nonetheless, mathematical methods that correctly describe networks of small size are still rare, and this prevents us to make progress in understanding neural dynamics at these intermediate scales. Here we develop a novel systematic analysis of the dynamics of arbitrarily small networks composed of homogeneous populations of excitatory and inhibitory firing-rate neurons. We study the local bifurcations of their neural activity with an approach that is largely analytically tractable, and we numerically determine the global bifurcations. We find that for strong inhibition these networks give rise to very complex dynamics, caused by the formation of multiple branching solutions of the neural dynamics equations that emerge through spontaneous symmetry-breaking. This qualitative change of the neural dynamics is a finite-size effect of the network, that reveals qualitative and previously unexplored differences between mesoscopic cortical circuits and their mean-field approximation. The most important consequence of spontaneous symmetry-breaking is the ability of mesoscopic networks to regulate their degree of functional heterogeneity, which is thought to help reducing the detrimental effect of noise correlations on cortical information processing.
  The brain is a complex system organized at multiple spatial scales, and the concerted interactions between these multiple scales of organization are probably crucial for the emergence of its computational power [24]. At the coarsest level, the brain can be divided into macroscopic circuits spanning several areas and containing millions of neurons. These macroscopic circuits are believed to accomplish very complex functions, ranging from motor control to cognition. At the finest level of organization, basic computations are performed at the single cell level by neurons, which can be seen as elementary computational units [20]. However, there is an intermediate, mesoscopic level of organization between the macroscopic and microscopic one: neurons are organized in microcircuits, whose size can vary from several thousands of cells as in cortical columns, to a few tens of cells as in micro-columns [45–47]. This mesoscopic level of investigation has received considerable attention in recent years both from the theoretical [25–29] and experimental [21,22] point of view, and it is often seen as a middle ground that is fundamental to link single neuron activity to behavior [23]. For this reason, finding an appropriate mathematical description of the brain at the mesoscopic scale is of fundamental importance for unveiling its emergent properties. At the mesoscopic scale, the brain is often described as a collection of neural masses, i. e. homogeneous neuronal populations within a cortical column [30]. Usually, these groups of neurons are described by the so-called neural-mass models [31]. A typical example is the well-known Jansen-Rit model [32–34], which describes local cortical circuits as a population of excitatory pyramidal cells, which receive inhibitory and excitatory feedback from local interneurons, and excitatory input from other regions such as the thalamus. This class of models may be studied using mean-field theory [35]. Mean-field theory is a mathematical tool that approximates the behavior of large networks [36–38], and it is useful to study the mass activity of few thousands of neurons, which constitutes the upper limit of mesoscopic descriptions of neural circuits and can adapt well to describe for example measures based on LFP/MEG/EEG recordings. This approximation becomes exact in the limit of networks with an infinite number of neurons, the so-called thermodynamic limit. For finite-size networks, however, the mean-field theory provides only an approximation of the real behavior of the system, and therefore may neglect important phenomena, such as qualitative differences in the transitions between static regimes and chaos [39], or in the degree or nature of correlations among neurons [40]. Clearly, these macroscopic differences in the dynamical and statistical behavior of finite and infinite-size networks may have important consequences on the information processing capability of the system, as potential oversimplifications in the mean-field approximation may hide important neural processes that are fundamental for the comprehension of neural circuits of finite size. To go beyond the mean field limit, several mathematical techniques have been developed to quantify finite-size effects. These finite-size methods (such as the linear noise approximation [41], the density functional approach [42], large-deviations theory [43] and path-integral methods [44]) typically can only be applied to mesoscopic circuits composed of a finite but large number of neurons. On the other side, most of the literature about small neural circuits deals with networks composed of 2–4 neurons (e. g. [1–3]), while methods for the analysis of networks made of a few tens of neurons, which represent the lower bound of the mesoscopic scale, are still limited. The purpose of our work is to make progress in the mathematical methodology to study the dynamics of such networks. The analysis of the dynamics of small neural networks was pioneered by Beer, who studied the bifurcations of networks of arbitrary size with highly symmetric assumptions on the strength of the synaptic weights [48], and through asymptotic approximations of the bifurcation manifolds [4]. In our article we extend his analysis to a more biologically plausible network of arbitrary size by deriving exact expressions of the bifurcation manifolds, and with less rigid constraints on the synaptic weights. In more detail, we consider a deterministic finite-size firing-rate network model with excitatory and inhibitory populations composed of an arbitrary number of homogeneous neurons in each population. Moreover, the model is characterized by homogeneous and arbitrarily strong weights between the populations. Then, we perform a numerical analysis of the global bifurcations that emerge by varying the external input currents (i. e. the stimuli) to the network and the strength of inhibition, and we introduce a mathematical theory that describes local bifurcations analytically. We find qualitative differences with the mean-field approximation when the system has strong inhibitory synaptic weights. In this case, through a phenomenon of spontaneous symmetry-breaking, the neural network undergoes a special bifurcation known as pitchfork or branching point [49], from which multiple solutions of the neural equations emerge. On the new branches, new bifurcations can occur, enriching considerably the complexity of the bifurcation diagram of the neural network. This dynamics is not revealed by the mean-field approximation. This article is organized as follows. In Materials and Methods we describe the neural model we use. We then explain intuitively in Results the formation of new branches of solutions through spontaneous symmetry-breaking, depending on the strength of inhibition. This is followed by a numerical and analytical study of the bifurcations of the network in weak and strong-inhibition regimes for different sizes of the inhibitory population. In Discussion we conclude by examining the biological implications of our results and by comparing our approach with previous work on homogeneous networks.
  Here we describe the assumptions we made about the model of finite-size neural circuits, whose structure is schematized in Fig 1, and whose dynamics we would like to investigate mathematically. These assumptions represent the best compromise we could find between biological detail, biological plausibility and mathematical tractability. We perform a numerical and analytical study considering the case of two neural populations of excitatory and inhibitory neurons respectively. The populations contain an arbitrary finite number of neurons, which are connected with each other through non-symmetric synaptic connections with arbitrarily strong weights. In order to make the network analytically tractable, we make some simplifying assumptions. In particular, we assume (as it is often made when considering local cortical circuits [30,31]) that the neurons in each population have homogeneous properties and parameters, that all neurons in the network are connected to each other (fully-connected topology), and that the axonal delays are negligible. Moreover, we study the deterministic neural equations, as it is common practice in the analysis of bifurcations. Therefore we do not consider typical sources of randomness such as stochastic fluctuations in the firing rates and random distributions of the strength of the synaptic weights. However, we note that some of these assumptions can be overcome by extensions of this formalism (see Discussion). In more detail, we consider a widely used rate model to describe the dynamics of single neurons [30,35,38,40,48,53]: d V i t d t = - 1 τ i V i t + 1 M i ∑ j = 0 N - 1 J i j A j V j t + I i, i = 0, . . . , N - 1, (1) where N ≥ 4 represents the number of neurons in the network. The function Vi (t) is the membrane potential of the ith neuron, while τi is its membrane time constant. Mi, the total number of connections to the ith neuron, is a normalization factor that has been introduced to avoid the explosion of the total synaptic input ∑ j = 0 N - 1 J i j A j (V j (t) ) in Eq (1) in the thermodynamic limit N → ∞. The synaptic connectivity matrix J specifies the strength of the synaptic connections in the network. Jij is the synaptic weight from the jth to the ith neuron, and for simplicity it is assumed to be deterministic and constant in time, for every pair of neurons. A j (·) is the activation function of the jth neuron and converts its membrane potential into the corresponding firing rate ν j = A (V j). Typically, piecewise-linear activation functions are required in order to obtain analytical results [11,55,56]. However, here we will show how to obtain explicit expressions for the equilibrium points and the local bifurcations of the network when A (V) is a smooth (i. e. infinitely differentiable) function. In more detail, we will consider S-shaped (i. e. sigmoidal) activation functions, since they are biologically realistic [6,54]. Particularly convenient from the mathematical point of view is the so-called algebraic activation function, which is defined as follows: Aj (V) =νjmax2[ 1+Λj2 (V−VjT) 1+Λj24 (V−VjT) 2 ]. (2) Here νmax is the maximum firing rate, Λ determines the slope of the activation function when νmax is fixed, and VT is the horizontal shift, therefore it can be interpreted as a firing threshold for the membrane potentials. In Eq (1), Ii are external input currents, namely the stimuli to the network, which here are supposed to be constant in time, in order to perform the bifurcation analysis of the network with Ii as bifurcation parameters. In order to make our analysis analytically tractable, we suppose that all the parameters of the system are homogeneous for each population of neurons considered. Thus, these parameters will be indexed only at the population level. By defining NE (NI) to be the size of the excitatory (inhibitory) population (with NE, I ≥ 2), and by indexing the neurons of the excitatory population as i = 0, …, NE − 1 and the inhibitory neurons as i = NE, …, N − 1 (with N = NE + NI), the synaptic connectivity matrix is structured as follows: J=[ JEEJEIJIEJII ], Jαβ={ Jαα (INα−IdNα), JαβINα, Nβ, for α=β for α≠β where the block J α β is a Nα × Nβ matrix that represents the connections from the population β to the population α, with α, β ∈ {E, I}. Moreover, I N α, N β is the Nα × Nβ all-ones matrix (here we use the simplified notation I N α = def I N α, N α), while IdNα is the Nα × Nα identity matrix. Since approximately 80% of neocortical neurons are excitatory and the remaining 20% are inhibitory [58], we set N E N I = 4. This is a typical choice in many theoretical works (e. g. [9,10]). Then we set the diagonals of J E E and J I I to zero because of the absence of self-connections in biological networks. JEE, JII, JEI, JIE are free parameters that describe the strength of the interactions between and within the neural populations. Clearly we have JEE, JIE > 0, and JII, JEI < 0, which also means that ME = MI = N − 1. It is important to observe that, compared to [21] where the authors divided the synaptic weights by the population size, we have chosen a different normalization. In both cases, the total synaptic input ∑ j = 0 N - 1 J i j A j (V j (t) ) is convergent in the thermodynamic limit, but the two normalizations behave in different ways when the network is mixed with both finite-size and infinite-size populations. According to our normalization, if one population has infinite size and the other has a finite size, then the contribution of the finite-size population to the total synaptic input of a neuron is negligible compared to the one of the infinite-size population. The normalization introduced by [38] leads instead to the infinite-size and finite-size populations both providing finite contributions to the total synaptic input. However, if required, it is possible to switch from our normalization to the other one by performing the following substitution: J α β → N - 1 N β J α β in all the analytical results reported in Results and in S1 Text. Given the notations described above, the external input currents are defined by two vectors, IE and II, representing the inputs to the excitatory and inhibitory populations: I α = I α 1 N α where 1 N α = def I N α, 1 is the Nα × 1 all-ones vector. The same subdivision between excitatory and inhibitory populations is performed for the other parameters of the network, namely τ, νmax, Λ, VT. To summarize, the equations of our neural system (1) can be written explicitly as follows:
  In mathematics, the branching-point bifurcations are described by the so-called equivariant bifurcation theory [50], namely the study of bifurcations in symmetric systems. The latter being rather technical, here we prefer to follow a more intuitive approach to the problem. First of all, we observe that according to bifurcation theory, local bifurcations are calculated by means of the eigenvalues of the Jacobian matrix of the network evaluated at the equilibrium points. Therefore, we set d V i (t) d t = 0 ∀ i in Eq (3) to compute the stationary solutions of the system. For the subsequent analysis, it is useful to introduce the following parameter: ψ = def τ I J I I ν I max Λ I 4 N - 1 (4) and we will show that, as long as inhibition is weak (ψ < 1), the equilibrium points in each population are always homogeneous: μ = μ E, …, μ E ︷ N E - times, μ I, …, μ I ︷ N I - times. (5) Specifically, μE and μI are the solutions of the following system of algebraic non-linear equations, obtained from Eq (3): F μ E, μ I = def - 1 τ E μ E + N E - 1 N - 1 J E E A E μ E + N I N - 1 J E I A I μ I + I E = 0 G μ E, μ I = def - 1 τ I μ I + N E N - 1 J I E A E μ E + N I - 1 N - 1 J I I A I μ I + I I = 0. (6) The curves defined by Eqs. F (x, y) = 0 and G (x, y) = 0 ∀ (x, y) ∈ ℝ 2 are the so-called nullclines of the network [17,18]. Fig 2 (top) shows an example obtained for JII = −10, while the remaining parameters are chosen as in Table 1. For completeness, in S1 Text we show how to get approximate analytical solutions for μE, I, even though our bifurcation analysis can be performed without knowing them explicitly. From Eqs (3) + (5), we obtain that the Jacobian matrix J of the linearized system evaluated at the stationary solution on the primary branch Eq (5) is: J=[ jEEjEIjIEjII ], jαβ={ −1ταIdNα+JααN−1Aα′ (μα) (INα−IdNα), for α=βJαβN−1Aβ′ (μβ) INα, Nβ, for α≠β (7) From this we can prove (see S1 Text) that the eigenvalues of J are: λ 0,1 R = Y + Z ± (Y − Z) 2 + 4 x 2, λ E = − [ 1 τ E + J E E N − 1 A E ′ (μ E) ], λ I = − [ 1 τ I + J I I N − 1 A I ′ (μ I) ], (8) where: x = N E N 1 (N − 1) 2 J E I J I E A E ′ (μ E) A I ′ (μ I), Y = − 1 τ E + N E − 1 N − 1 J E E A E ′ (μ E), Z = − 1 τ I + N I − 1 N − 1 J I I A I ′ (μ I). (9) The inequalities ψ < 1 and λI ≥ 0 imply: A I ′ μ I ≥ N - 1 τ I J I I > ν I max Λ I 4 while, according to Eq (2), the derivative of the activation function cannot exceed A I ′ (V I T) = ν I max Λ I 4. Therefore λI is always negative in the weak-inhibition regime. According to bifurcation theory [49], the system undergoes special bifurcations when one of its eigenvalues is equal to zero. In particular, whenever λI = 0, the network may undergo a so-called branching-point bifurcation. Since λI is always negative in the weak-inhibition regime, branching-point bifurcations occur only in the strong-inhibition regime (ψ ≥ 1). This special bifurcation corresponds to the formation of heterogeneous membrane potentials in the inhibitory population. In other words, strong inhibition may break the symmetry of the system. Intuitively, this can be understood as in Fig 3. In the weak-inhibition regime there is only one valley or basin in the “energy landscape” of the network. However, strong inhibition may lead to the formation of multiple valleys, and a small perturbation determines the valley the inhibitory potential will converge to. For this reason, for strong inhibition multiple new branches of the equilibrium points may emerge, which are described by the following stationary solutions: μ = μ E, …, μ E ︷ N E - times, μ I, 0, …, μ I, N I - 1 (10) where μE and μI, i are the solutions of the following system of algebraic non-linear equations, obtained from Eq (3) in the stationary regime: { −1τEμE+NE−1N−1JEEAE (μE) +JEIN−1∑j=0NI−1AI (μI, j) +IE=0−1τIμI, i+NEN−1JIEAE (μE) +JIIN−1∑j=0j≠iNI−1AI (μI, j) +II=0for i=0, …, NI−1. (11) For example, for NI = 2, Eq (11) can be written more explicitly as follows: F μ E, μ I, 0, μ I, 1 = def - 1 τ E μ E + N E - 1 N - 1 J E E A E μ E + J E I N - 1 A I μ I, 0 + A I μ I, 1 + I E = 0 G μ E, μ I, 0, μ I, 1 = def - 1 τ I μ I, 0 + N E N - 1 J I E A E μ E + J I I N - 1 A I μ I, 1 + I I = 0 H μ E, μ I, 0, μ I, 1 = def - 1 τ I μ I, 1 + N E N - 1 J I E A E μ E + J I I N - 1 A I μ I, 0 + I I = 0. (12) The surfaces F (x, y, z) = 0, G (x, y, z) = 0, H (x, y, z) = 0, ∀ (x, y, z) ∈ ℝ 3 are a higher-dimensional extension of the nullclines F (x, y) = 0 and G (x, y) = 0 ∀ (x, y) ∈ ℝ 2 that we encountered in the weak-inhibition regime. These surfaces have been called nullsurfaces in [48]. Fig 4 (top) shows an example obtained for JII = −100, while the remaining parameters are chosen as in Table 1. For the sake of clarity, here we treat in detail only the case NI = 2, and we will discuss briefly some results in the strong-inhibition regime for NI > 2. For NI = 2, from Eq (12) the Jacobian matrix on the secondary branches Eq (10) is: J = j 00 j 01 j 02 j 10 j 11 j 12 j 20 j 21 j 22 (13) where: J00=−1τEIdNE+JEEN−1AE' (μE) (INE−IdNE), J01=JEIN−1AI' (μI, 0) 1NE, J02=JEIN−1AI' (μI, 1) 1NE, J10=JIEN−1AE' (μE) 1NEt, J11=1τI, J12=JIIN−1AI' (μI, 1), J20=JIEN−1AE' (μE) 1NEt, J21=JIIN−1AI' (μI, 0), J22=1τI (here t is the transpose operator applied to a matrix), as derived in more detail in S1 Text. Intuitively, on the primary branch the Jacobian matrix was a 2 × 2 block matrix (see Eq (7) ), because we had only one excitatory and one inhibitory membrane potential (μE and μI respectively), while on the secondary branches it is a 3 × 3 block matrix since now the two inhibitory neurons have different potentials (μI, 0 and μI, 1). However, the Jacobian matrix on the new branches can be calculated for a network with an arbitrary number of inhibitory neurons (see S1 Text for more details). Finally, we observe that for any NI it is possible to find a relation between any pair (i, j) of inhibitory membrane potentials in the strong-inhibition regime, which proves very useful when we calculate analytically the local bifurcations of the system on the secondary branches. Indeed, from the second equation of the system (11), we get: - 1 τ I μ I, i + J I I N - 1 A I μ I, j = - 1 τ I μ I, j + J I I N - 1 A I μ I, i. (14) After some algebra Eq (14) can be converted into the following fourth-order polynomial equation: a ^ μ I, j 4 + b ^ μ I, j 3 + c ^ μ I, j 2 + d ^ μ I, j + e ^ = 0 (15) whose coefficients depend on μI, i as follows: a^=ΛI24τI2b^=−ΛI22τI (ϕ^+VITτI) c^= ΛI24[ ϕ^2+ (VITτI) 2+4τIVITϕ^ ]+1τI2−ξ^d^=−ΛI22ϕ^VIT (VITτI+ϕ^) −2τIϕ^+2ξ^VITe^= (ΛI2ϕ^VIT) 2+ϕ^2−ξ^ (VIT) 2ϕ^= 1τIμI, i+JIIN−1AI (μI, i) −νImaxJII2 (N−1) ξ^= (νImaxΛIJII4 (N−1) ) 2. Eq (15) can be solved analytically, providing an explicit expression of μI, j as a function of μI, i, which will be used later to evaluate the local bifurcations in the strong-inhibition regime. In Fig 5 we plot this relation by removing the spurious portions of the solutions of Eq (15) that do not satisfy Eq (14). This figure proves the formation of three branches of solutions of the stationary membrane potentials at the branching-point bifurcations of the network. Even though in the article we focus on the case NE = 8 and NI = 2, the analytical results derived in this subsection are valid for arbitrary NE, I. As we said, we want to understand how the network’s dynamics changes when we vary the external input currents IE, I and the strength of the synaptic weight JII. In this subsection we show the results that we obtain when we vary these parameters one by one, because this allows us to introduce the concepts of codimension one and codimension two bifurcation diagrams. For NI = 2, we first fix II to −10 and the inhibitory strength JII to −10. As illustrated in Fig 2 (top), Eq (6) admits multiple solutions depending on the specific value chosen for the current IE. Thus, we vary IE continuously while keeping II and JII fixed. In this way we can plot how the solutions μE, I of the system (6) change as a function of IE. As we anticipated in the previous subsection, the curve that we obtain (Fig 2 (bottom) ) is called the primary branch of the network. This figure represents the codimension one bifurcation diagram of the network from which we individuate two special points, called local saddle-node bifurcations (shortened to LP hereafter). They identify the value of IE for which the number of solutions of Eq (6) changes (notice the correspondence with the top panels of Fig 2). These are the first examples of (local) bifurcation points that the neural network exhibits, and that lead to the formation of hysteresis (see left panel in Fig 6). Briefly, hysteresis was suggested to describe the reverberation and persistence of neural activity sustained endogenously after an external input is removed, which is thought to be crucial for phenomena such as working memory [62–64]. It is important to observe that even if reverberation requires bistability, the latter can be present without hysteresis, but very often they coexist, as in our model. We then continuously vary both IE and II, while keeping JII fixed, and we plot the solutions of Eq (6) as a function of both the external currents, so that μE, I = μE, I (IE, II) now define three dimensional manifolds (see Fig 7 (left) for μE). On these manifolds, the special points LP depend also on II, therefore they form a set of points called saddle-node curves (see the blue curves in Fig 7, left). For visual convenience, these curves are projected on the IE−II plane in Fig 7 (right), defining the so-called codimension two bifurcation diagram of the network. This diagram will be the main instrument of the analysis presented hereafter. However, from Fig 7 (left) and the bottom panels of Fig 2, we can see that the saddle-nodes are not the only bifurcations we have in our network. For some values of the pair IE−II, the system can undergo also a local Andronov-Hopf bifurcation, or H for short. These bifurcations are represented by the red curves in Fig 7, and correspond to the emergence of neural oscillations, which are thought to play a key role in many cognitive processes [65]. Hereafter, we list all the bifurcations the system undergoes, dividing them in groups depending on the codimension of the bifurcation, which is defined as the number of parameters (in our case IE, I) that must be varied for the bifurcation to occur. Although only few of them are represented in Fig 2 (bottom), the complete set of codimension one bifurcations our system undergoes is: local saddle-node bifurcations (LP), for which new equilibria arise, or collide and annihilate each other; local Andronov-Hopf bifurcations (H), where stable or unstable self-sustained oscillations, described by limit cycles, arise or disappear; global homoclinic bifurcations, where limit cycles vanish in a special equilibrium point (i. e. a neutral saddle, see the next subsection, or an LP bifurcation point), giving rise to an orbit with infinite period; global limit point of cycles bifurcations (LPC), at which new limit cycles arise, or collide and annihilate each other. The codimension one diagrams collecting all these bifurcations are shown in S1 Text. Moreover, on the curves defined by these bifurcations, and that are obtained by varying both IE and II, the following codimension two bifurcations appear (see Fig 7, right): local cusp bifurcations (CP), on the LP curves; local generalized Hopf bifurcations (GH), which divide subcritical H curves from supercritical ones; local Bogdanov-Takens bifurcations (BT), which represent the contact point between the LP, H (ending here) and homoclinic curves; global cusp bifurcation of cycles (CPC), on the limit point of cycles curves; global saddle-node on invariant circle (SNIC), where an LP bifurcation occurs simultaneously with a homoclinic bifurcation. It is worth remarking that the bifurcation diagram in Fig 7 (right), obtained from the voltage-based model (3) in the weak-inhibition regime, is qualitatively similar to that of the activity-based (mean-field) Wilson-Cowan model (see Fig. 2. 12 in [66]). This strong similarity results from the fact that under the hypotheses of stationary inputs, invertible connectivity matrix J and homogeneous membrane time constants, the two models are mathematically equivalent (apart from a rescaling by τ of the external input). Indeed, it is possible to prove that the membrane potentials Vi (t) and the neural activities Ai (t) are formally related by the linear expression V i (t) = 1 M i ∑ j = 0 N - 1 J i j A j (t) + I i [66]. However, in the next subsection we will show that things may change significantly in the strong-inhibition regime if we take into account finite-size effects. Interestingly, Fig 7 presents two of the so-called catastrophe manifolds [67], one of which is shown in the right panel of Fig 6. This figure emphasizes the ability of the model to describe three different behaviors: leaky integrator, perfect integrator and switch. This triad represents the main ingredient for describing a mechanism which was proposed to explain interval timing by neural integration [68,69]. According to the theory, by changing some parameters of the network, it is possible to modify the shape of the equilibrium curve on the catastrophe manifold generating a ramping activity that can explain Weber’s law of time perception [70]. This phenomenon can easily occur in our model, where the shape of the equilibrium curve can be changed dynamically by varying the input currents. Now we prove some of our previous results on local bifurcations from the analytical point of view (global bifurcations are often harder to study analytically, therefore they are not considered here). We base our analysis on [5], where the authors developed a technique for deriving analytical expressions of local codimension one bifurcations. Furthermore, for the first time we will extend their method to the analysis of spontaneous symmetry-breaking in systems with arbitrary size and to some local bifurcations with codimension larger than one. We start by observing that, according to bifurcation theory [49,61], LP is defined by the condition that one real eigenvalue of the Jacobian matrix becomes equal to zero. Therefore, from Eq (8), we conclude that for our network this bifurcation occurs whenever λ 0 R = 0 or λ 1 R = 0, because λE is always negative, while λI is always negative in the weak-inhibition regime. For example, if Y + Z < 0, the condition λ 0 R = 0 is equivalent to X = Y Z which, according to Eq (9), provides: A I ′ (μ I) = − 1 τ E τ I + 1 τ I N E − 1 N − 1 J E E A E ′ (v) − 1 τ E N I − 1 N − 1 J I I + 1 (N − 1) 2 [ (N E − 1) (N I − 1) J E E J I I − N E N I J E I J I E ] A E ′ (v) (16) where we have defined the parameter v = def μ E. Now we invert A I ′ (μ I) (more details are provided in S1 Text), obtaining: μI (v) =VIT±2ΛI (νImaxΛI4AI′ (μI) ) 23−1 (17) and from Eq (6) we get: I E v = 1 τ E v - N E - 1 N - 1 J E E A E v - N I N - 1 J E I A I μ I v I I v = 1 τ I μ I v - N E N - 1 J I E A E v - N I - 1 N - 1 J I I A I μ I v. (18) These are parametric equations in the parameter v ∈ (v a, v b), where: vb, a=VET±2ΛE (NE−1N−1JEEνEmaxΛE4τE) 23−1 (19) and they define analytically the blue curves in Fig 7 (right) (the same result is obtained if Y + Z > 0, in which case the condition λ 1 R = 0 is equivalent to X = Y Z). As we said, this is not sufficient to prove that Eqs (16) – (19) describe LP curves, since we should check also the corresponding non-degeneracy conditions. Nevertheless, we observe a perfect agreement between these analytical curves and those obtained numerically by Cl_MatCont and XPPAUT, therefore for simplicity we do not check the remaining conditions and we leave them to the most technical readers. We adopt the same approach for the remaining bifurcations we are about to describe. Now we focus on the H bifurcations. According to [49,61], they appear whenever the network has complex-conjugate purely imaginary eigenvalues. Since λE, I are always real, this condition can be satisfied only by λ 0,1 R, by setting Y + Z = 0 and (Y − Z) 2 + 4 X < 0. In particular, from the equation Y + Z = 0 we get: A I ′ μ I = N - 1 N I - 1 J I I 1 τ E + 1 τ I - N E - 1 N - 1 J E E A E ′ v, (20) where v = def μ E as before. Following the same procedure introduced before for the LP curves, we obtain a set of parametric equations for the pairs IE−II that generate the H curves, with parameter v ∈ [ v f, v d ] ∪ [ v c, v e ], where: vc, d= VET±2ΛE (νEmaxΛE (NE−1) JEE4 (N−1) 11τE+1τI−νImaxΛI (NI−1) JII4 (N−1) ) 23−1ve, f=VET±2ΛE (νEmaxΛE4z) 23−1z=−b−b2−4ac2aa= (NE−1N−1JEE) 2−NENI (NE−1) (N−1) 2 (NI−1) JEEJEIJIEJIIb=−2τENE−1N−1JEE+NENI (N−1) (NI−1) JEIJIEJII (1τE+1τI) c=1τE2. (21) As described above, BT represents the point where the LP and H curves meet each other, and identifies also the end of the H curve. From the condition λ 0 R = 0 or λ 1 R = 0 that defines the LP curves, and the condition λ 0,1 R = ± ι ω (where ι represents the imaginary unit) that defines the H curves, we get λ 0 R = λ 1 R = 0. This is the condition that defines analytically the BT points, or equivalently v BT = v e, f as given by Eq (21), from which the coordinates of the BT points in the IE−II plane can be easily obtained through Eq (6). The remaining local bifurcations (i. e. CP and GH) are analytically intractable. To conclude this subsection, we describe briefly the effect of the variation of the remaining synaptic weights on the codimension two bifurcation diagram, considering the weights in Table 1 as reference point. Given that their variation does not generate interesting phenomena such as the branching-point bifurcation which is obtained by varying JII, we describe the effects of these parameters succinctly. The reader may verify the following results by the supplemental Python code (see S1 File). For JEE ≫ 10, the two LP curves become larger and larger on the IE axis (i. e. the distance between their vertical asymptotes increases). Moreover, the curves get closer and closer to each other, by shifting on the IE axis, until they intersect and their oblique parts (i. e. those between the BT points) overlap. If we increase JEE further, the LP curves split again in two disjoint parts, each one presenting two BT and two CP bifurcations (thus the total number of CP points increases from two to four). Between each pair of BT points (on the same LP curve) there is an H curve. These curves are very close to the corresponding LP curves, and if we increase JEE further they disappear, together with the BT bifurcations. Thus for very large JEE we get only two disconnected LP curves, or in other terms, for very strong excitation the oscillatory activity cannot be sustained anymore. Also on the opposite side, namely for weak excitation (i. e. JEE → 0), the H curves disappear. More precisely, according to Eq (21), these curves vanish whenever the following condition is satisfied: ν E max Λ E 4 z ≤ 1 or b 2 - 4 a c ≤ 0. Moreover, for weak excitation the width on the IE axis of the LP curves decreases, i. e. the distance between their vertical asymptotes becomes smaller and smaller, until the asymptotes collapse on each other and the LP curves disappear as well. In particular the curves vanish whenever the following condition is satisfied (see Eq (19) ): N E - 1 N - 1 J E E ν E max Λ E 4 τ E ≤ 1. For ∣JEI∣ ≫ 70, the width of the two LP curves remains almost constant, while the distance between them (and therefore also the length of the H curves) increases continuously. On the other side, for ∣JEI∣ → 0, the two LP curves get closer and closer to each other, until they intersect and then their oblique parts between the BT points overlap. If we decrease ∣JEI∣ further, similarly to the case with large JEE, the two LP curves split in two disjoint parts, each one presenting two BT and two CP bifurcations. For even smaller values of ∣JEI∣, the BT, CP and H bifurcations disappear, while the LP curves disappear for ∣JEI∣ = 0. For JIE ≫ 70, the LP curves are stretched vertically and shifted downwards along the II axis. Clearly, in the opposite direction (i. e. for JIE → 0), they are compressed, and while the IE coordinates of the vertical asymptotes remain almost unchanged with JIE, the two CP points get closer and closer to each other. At the same time, the two H curves tend to overlap. At some value of JIE, the two CP bifurcations touch each other and disappear, thus the two LP curves become tangent. If we further decrease JIE, the two LP curves split again, one over the other, and the BT and H bifurcations disappear. All the phenomena that we have just described are qualitatively similar for different values of JII, thus they occur also in the strong-inhibition regime for the primary branch. In the strong-inhibition regime (in particular here we consider the cases JII = −34 and JII = −100), most of the features of the weak-inhibition bifurcation diagram are preserved. However, besides the bifurcations explained in the previous subsection, from Figs 8,9 and 10, we can see that the system undergoes also the following codimension one bifurcations: local branching-point bifurcations (BP), at which two or more equilibrium curves emanate or collapse; local torus bifurcations (TR), at which the network exhibits quasi-periodic oscillations; this is a local bifurcation of the Poincaré map of a limit cycle of the network [49], whose change of stability determines the formation of solutions of Eq (3) containing two incommensurable frequencies; and the following codimension two bifurcations: local zero-Hopf (neutral saddle) bifurcations (ZH), at which the H curves and the BP curves intersect each other. This point identifies a change in stability of the H bifurcation curve. In particular, the BP bifurcations lead our model to show multiple branches of stationary solutions for suitable current values. This is a finite-size effect due to the finite number of neurons in each population, which leads to a richer set of eigenvalues than that obtained by using methods based on the reduction of the number of equations, such as the mean-field approximation (see S1 Text for more details). In order to thoroughly investigate the bifurcations the system undergoes in presence of strong inhibition, we start by analyzing the codimension one bifurcation diagram for JII = −34. In particular, the diagram in Fig 8 is obtained by setting II = −10. It turns out that, in addition to the primary equilibrium curve (black line), new branches of stationary states (violet lines) emanate and collapse in two BP bifurcations. These secondary branches hold supercritical H bifurcations that give rise to stable limit cycles. Instead, on the primary branch, we find two LP bifurcations and a subcritical H bifurcation, whose unstable limit cycles vanish in a homoclinic orbit. We also observe that the inhibitory neurons, as well as the excitatory neurons, for NI = 2 have the same bifurcation diagram (see Fig 8). However, differently from the excitatory case, this does not mean that the inhibitory membrane potentials are homogeneous. Indeed, when an inhibitory neuron is on the upper secondary branch (see the violet curve above the primary equilibrium curve in Fig 8, right), the other one is on the lower secondary branch, thus they are heterogeneous. For JII = −100, secondary branches of equilibrium points are still present, see Fig 9. Together with a subcritical H bifurcation, they unveil also LP bifurcations. In particular, the former generates unstable limit cycles that become stable after having crossed the limit point of cycles bifurcation (dark green loop). For increasing values of IE, the stable limit cycles collapse into the unstable limit cycles originated from the subcritical H bifurcation belonging to the primary equilibrium curve. Before collapsing, the stable limit cycles undergo a torus bifurcation (gray loop). By varying also II, we obtain the codimension two bifurcation diagrams displayed in Fig 10. It is worth noting that the branching points the system undergoes generate two bifurcation curves (light green dot-dashed lines) that pass through the whole IE−II domain. The presence of these curves is the most relevant difference with the weak-inhibition regime and the classic (mean-field) Wilson-Cowan model (compare with Fig. 2. 12 in [66]). Furthermore, the LP and H bifurcations that belong to the secondary branches give rise to further bifurcation curves (purple and light blue lines, respectively) in the IE−II domain, as shown in Fig 10. Interestingly, since the BP bifurcation increases the dimension of the network from 2 (for λI < 0, see Eq (6) ) to 3 (for λI > 0, see Eq (12) ), the network can exhibit more complex dynamics, such as quasi-periodic motions originated from the torus bifurcations. The biological importance of quasi-periodic oscillations in neural communication was discussed in [71]. Now we want to study the local bifurcations from the analytical point of view. On the primary branch, the LP, H and BT bifurcations have the same expressions obtained in the weak-inhibition regime, therefore they are valid for arbitrary NE, I, as well as the formulas that we derive below for the BP and ZH bifurcations. On the contrary, unlike the weak-inhibition regime, the formulas of the LP, H and BT bifurcations on the secondary branches in the strong-inhibition regime are valid only for NI = 2 (while NE can still be arbitrary), even though our formalism can be extended to the case NI > 2. We start by considering the BP bifurcations, which are defined by the condition λI = 0, as we saw before. From Eq (8) this condition implies: A I ′ μ I BP = N - 1 τ I J I I, (22) thus the solutions of this equation are: μ I BP = V I T ± 2 Λ I ψ 2 3 - 1. (23) Eq (23) shows that μI (BP) is defined only for ψ ≥ 1, and therefore that the BP bifurcations occur only in the strong-inhibition regime. Now, from the second equation of the system (6) (we can also use Eq (12), since for λI = 0 they are equivalent) we get: μ E BP = V E T ± 2 Λ E 1 1 2 N - 1 ν E max N E J I E 1 τ I μ I BP - N I - 1 N - 1 J I I A I μ I BP - I I - 1 2 - 1 (24) while from the first equation of Eq (6) we get: I E = 1 τ E μ E BP - N E - 1 N - 1 J E E A E μ E BP - N I N - 1 J E I A I μ I BP (25) where μI (BP) and μE (BP) are given by Eqs (23) and (24) respectively. Since μE (BP) depends on II, Eq (25) defines two explicit functions I E = F ± (I I), that provide the curves on which we have a BP bifurcation (see the light green lines in Fig 10 for JII = −34 and JII = −100. More details can be found in S1 Text). It is important to observe that for NI > 2 the network may undergo special branching points that are not determined by the condition λI = 0, rather by the fact that one of the eigenvalues of the reduced Jacobian matrix introduced in S1 Text tends to zero. These branching points can be studied analytically through our approach, but since the complexity of the corresponding eigenvalues strongly depends on NI and on the degree of heterogeneity of the inhibitory population, we do not analyze them in detail. The points where the H and BP curves meet each other define the ZH bifurcations. From this definition, we see that they can be calculated analytically from the conditions λ 0,1 R = ± ι ω and λI = 0, from which in turn we get: AE′ (μE (ZH) ) =N−1 (NE−1) JEE (1τE+NIτI) AI′ (μI (ZH) ) =N−1τI| JII | and therefore: μE± (ZH) =VET±2ΛE (νEmaxΛE (NE−1) JEE4 (N−1) (1τE+NIτI) ) 23−1μI± (ZH) =VIT±2ΛI (νImaxΛIτI| JII |4 (N−1) ) 23−1. As usual, if we substitute these expressions of the membrane potentials into Eqs (6) or (12), we obtain the coordinates of the ZH points in the IE−II plane. On the secondary branches that are generated by the branching points, new bifurcations can occur (in the case NI = 2, see for example the LP and H bifurcations in Figs (8) and (9), and the corresponding light blue and purple curves in Fig 10), also new branching points (for NI > 2), from which tertiary branches emerge, and so on. To study them, according to bifurcation theory, we need the Jacobian matrix of the network on the secondary (tertiary, and so on) branches, as we will explain more clearly in the next subsection. Here we focus again on the case NI = 2, therefore we can determine the local bifurcations on the secondary branches by means of Eq (13). Now we start with the LP bifurcations. We know that they are defined by the condition that one of the eigenvalues of Eq (13) is equal to zero. From it, as explained in S1 Text, we obtain that: A E ′ μ E = b ´ a ´ (26) where: a´=1τI2NE−1N−1JEE+1τINE (N−1) 2JEIJIE[ AI′ (μI, 0) +AI′ (μI, 1) ] +1 (N−1) 3[ 2NEJEIJIEJII− (NE−1) JEEJII2 ]AI′ (μI, 0) AI′ (μI, 1) b´=1τE[ 1τI2− (JIIN−1) 2AI′ (μI, 0) AI′ (μI, 1) ]. Therefore if we invert Eq (26) and we use the solution of Eq (15), we obtain the expression of μE as a function of μI, 0. If we replace the solutions μE and μI, 1 in the system (12), we get parametric equations for IE, I as a function of a single parameter, which is now defined as v = def μ I, 0. These equations are an analytical description of the light blue curves shown in Fig 10 (right) for JII = −100. Similarly, for the H bifurcations we obtain the following condition: A E ′ μ E ± = - b ` ± b ` 2 - 4 a ` c ` 2 a ` (27) where: a`=NE−1N−1JEE[ 2τINE−1N−1JEE+NE (N−1) 2JEIJIE (AI′ (μI, 0) +AI′ (μI, 1) ) ]b`=2NE (N−1) 3JEIJIEJIIAI′ (μI, 0) AI′ (μI, 1) − (1τE+1τI) [ 4τINE−1N−1JEE+NE (N−1) 2JEIJIE (AI′ (μI, 0) +AI′ (μI, 1) ) ]c`=2τI[ (1τE+1τI) 2− (JIIN−1) 2AI′ (μI, 0) AI′ (μI, 1) ], therefore again it is possible to describe these bifurcations analytically, obtaining the same results we got numerically in Fig 10 for JII = −34 and JII = −100 (see the purple curves in both the panels). However, unlike the primary branch, our theory does not allow us to calculate the range of the parameter v = μ I, 0 on the secondary branches, since the resulting equations that define the range are analytically intractable. In the same way, it is not possible to calculate explicitly the coordinates of the new BT bifurcations, where the LP and H curves that emanate from the secondary branches meet each other. Therefore they can be evaluated only through analytical approximations (which are beyond the purpose of the article), or through a numerical approach (as in S1 File). For an arbitrary NI, if spontaneous symmetry-breaking has not occurred, and if the initial conditions are homogeneous, then all the neurons in each population are described by the same equation. Therefore the network can be interpreted as a two-dimensional system, and for the Poincaré-Bendixson theorem it cannot show chaotic behavior. Nonetheless, when a symmetry-breaking occurs, the inhibitory neurons split in subgroups, and the network becomes at least three-dimensional depending on the number of inhibitory subpopulations arisen. In this case, it is possible to check numerically the presence of chaos by computing the Lyapunov exponents [62]. In particular, when the symmetry-breaking occurred in the simplest case NI = 2, we excluded the presence of chaos since the Lyapunov exponents were always negative. However, the emergence of chaos is still possible for NI > 2, as we will discuss in the next subsection. We conclude by describing briefly the effect of the variation of the remaining synaptic weights. When we discussed the weak-inhibition regime, we showed that the LP and H curves on the primary branch change in a qualitatively similar way for different values of JII, and the same conclusions apply to the strong-inhibition regime. On the other side, now we want to analyze the behavior of the BP curves and of the bifurcations on the secondary branches. For JII = −100 and increasing JEE, the most notable phenomenon is the overlap between the oblique parts of the BP and the LP curves of the primary branch. The latter finally collapse on each other and split in two disjoint parts, as in the case JII = −10 discussed for the weak-inhibition regime, while the bifurcations on the secondary branches do not show any interesting variation. Furthermore, when JEE → 0, we observe first of all the disappearance of the ZH bifurcations. This occurs because the H curves on both the primary and the secondary branches do not meet the BP curve anymore. If we further decrease JEE, the two CP bifurcations on the LP curve of the secondary branches get closer and closer until they annihilate each other and the curve disappears. This phenomenon implies also the disappearance of the BT bifurcations on the secondary branches. For smaller values of JEE, the H curves on both the primary and secondary branches disappear, and finally also the LP curve on the primary branch (see the weak-inhibition regime) and the BP curve. To conclude, for large ∣JEI∣ or large JIE, the LP curve on the secondary branches disappears again through the annihilation of its CP points (as explained above), while on the other side, when at least one of the two synaptic weights is small, we do not observe any interesting variation of the bifurcations on the secondary branches. The same analysis can be performed on networks with an arbitrary number NI of inhibitory neurons. When λI in Eq (8) goes to zero, we observed the formation of secondary branches from the primary one. This means that the inhibitory neurons split into two subsets [52,78], and all the neurons in a given subset behave identically. For this reason we can reinterpret the system as a network with an excitatory population with NE neurons, and two inhibitory populations. Furthermore, when we change the current IE (while keeping λI > 0) with II fixed, at some point one of the eigenvalues of the Jacobian matrix of the system (11) (see S1 Text for their analytical calculation, supposing that the sizes of the clusters are known) may go to zero, generating a new branching point on the secondary branches. In this case we observe the formation of tertiary branches, and so on. Whenever a group of K inhibitory neurons splits into two clusters with sizes P and Q = K−P and membrane potentials μ I (P) and μ I (Q) respectively, the Jacobian matrix of the network has eigenvalues λ I (P) = − [ 1 τ I + J I I N − 1 A I ′ (μ I (P) ) ] and λ I (Q) = − [ 1 τ I + J I I N − 1 A I ′ (μ I (Q) ) ] with multiplicity Q−1 and P−1 respectively. Immediately beyond the BP bifurcation that originated the two clusters, one of λ I (P) and λ I (Q) is positive. Thus if P and/or Q is larger than one, the new branches are unstable near the bifurcation point. These eigenvalues do not affect the stability of the new branches only in the case P = Q = 1 (having multiplicity equal to zero), which is compatible with the fact that for NI = 2 the secondary branches are stable near the branching point (see Figs 8 and 9). In the general case P = Q = 1 the stability of the new branches depends on the eigenvalues of the reduced Jacobian matrix introduced in S1 Text. These eigenvalues are strongly affected by the value of the parameter NI and by the degree of heterogeneity of the inhibitory population. For this reason it is not easy to draw general conclusions about the stability of the new branches in the case P = Q = 1. In Fig 11 we show an example of formation of secondary branches in the case NI = 4, obtained as usual for the values of the parameters reported in Table 1. In this case K = 4, P = 3 and Q = 1, therefore the secondary branches are unstable near the branching points. It is apparent there is an important difference compared to the case NI = 2. For a network with only two inhibitory neurons, Eq (12) implies that they both have the same codimension one bifurcation diagram (see the right panels of Figs 8 and 9). This is just a special case, because in general, for NI > 2, the inhibitory neurons have different codimension one bifurcation diagrams. Although Fig 11 is useful to get an idea of how a bifurcation diagram for NI = 4 would look like, it is important to underline that this diagram could not be entirely thorough. Indeed, since in this case λI in Eq (8) has multiplicity larger than one, the Cl_MatCont toolbox may fail in providing us the complete picture [78]. Nonetheless, in principle it is possible to check the completeness of the diagram since, for an arbitrary NI, local bifurcations can still be calculated analytically through our approach. In fact, from Eqs (14) + (15), it is possible to express NI−1 inhibitory membrane potentials as functions of the remaining one, which can be used as a parameter for the parametric equations in the codimension two bifurcation diagram. To conclude, we observe that for NI > 2 new kinds of bifurcations may appear, which do not occur in the case NI = 2. For example, for NI = 3, if the network has three different inhibitory potentials, the characteristic equation of the Jacobian matrix has the form p (λ) = (λ − λ E) N E − 1 p R (λ), where p R (λ) (the characteristic polynomial of the reduced Jacobian matrix introduced in S1 Text) is a fourth order polynomial. This means that in principle, for some values of the parameters, the network may have two pairs of purely imaginary eigenvalues. This condition corresponds to the formation of a double-Hopf bifurcation, which in turn may imply a local birth of chaos (see for example [72,73]).
summary: The mesoscopic level of brain organization, describing the organization and dynamics of small circuits of neurons including from few tens to few thousands, has recently received considerable experimental attention. It is useful for describing small neural systems of invertebrates, and in mammalian neural systems it is often seen as a middle ground that is fundamental to link single neuron activity to complex functions and behavior. However, and somewhat counter-intuitively, the behavior of neural networks of small and intermediate size can be much more difficult to study mathematically than that of large networks, and appropriate mathematical methods to study the dynamics of such networks have not been developed yet. Here we consider a model of a network of firing-rate neurons with arbitrary finite size, and we study its local bifurcations using an analytical approach. This analysis, complemented by numerical studies for both the local and global bifurcations, shows the emergence of strong and previously unexplored finite-size effects that are particularly hard to detect in large networks. This study advances the tools available for the comprehension of finite-size neural circuits, going beyond the insights provided by the mean-field approximation and the current techniques for the quantification of finite-size effects.
keywords: medicine and health sciences
  neural networks
  nervous system
  bifurcation theory
  membrane potential
  condensed matter physics
  electrophysiology
  neuroscience
  systems science
  mathematics
  algebra
  network analysis
  computer and information sciences
  animal cells
  neural pathways
  mesoscopic physics
  physics
  cellular neuroscience
  neuroanatomy
  cell biology
  linear algebra
  anatomy
  physiology
  neurons
  biology and life sciences
  cellular types
  physical sciences
  eigenvalues
title: The Complexity of Dynamics in Small Neural Circuits