Introduction

Firing-rate fluctuations and irregular spiking are ubiquitous in the neocortex [1, 2]. Furthermore this temporal variability is often observed to be correlated across spatial scales ranging from local cortical circuits to the entire brain: in local cortical circuits both in membrane potential fluctuations [3] and on the level of spiking [4–7], in the coherency measured in brain-wide EEG signals [8, 9], and in the global signal observed across all voxels in fMRI measurements [10–12].

A class of theoretical models has successfully accounted for temporal variability via internally generated chaotic dynamics of recurrent networks, whether through excitation-inhibition balance in spiking models [13, 14] or the more abstract models of rate chaos in randomly connected networks [15]. Yet a key emergent feature of these models is the decorrelation of neural activity such that the macroscopic, population activity remains nearly constant in time. Population-wide coherence or synchrony can be generated in a variety of ways for example by introducing spatial modes with self-excitation, but this comes at a cost of drowning out the chaotic fluctuations and yielding fixed points [16]. Indeed a major challenge to theorists has been to produce network models which generate spatially coherent, temporally irregular fluctuations which can account for broad spatial correlations observed in experiments.

Two recent studies have shown that excitation-inhibition balance networks can generate spatially modulated correlations [17, 18]. In both of these studies the correlations are driven by common input from an external source, and the average correlation across the network remains small. It remains an open question whether a network can internally generate correlated fluctuations that are coherent across the entire network.

The chaotic dynamics of a network of randomly connected firing-rate neurons has been well-studied [15, 19]. In such a network each individual neuron’s firing rate is given by a non-linear function of its input, which is in turn a weighted sum of the firing rates of all other neurons in the network. The network exhibits a phase transition from a fixed point to chaotic activity in which the randomness of the weights reverberates uncorrelated fluctuations throughout the network. Typically in this chaotic regime pairwise correlations are small and no coherent fluctuations emerge. Here we extend this model by adding a low-rank structured component of connectivity to the random connections. The structured connectivity sums the fluctuations along one spatial mode and projects them along a second, orthogonal mode yielding coherent fluctuations without drowning out the individual neuron variability which continues to drive the chaotic dynamics. A previous work studied a specific example of this structure and focused primarily on analyzing the non-chaotic regime [20]. Here we focus on the chaotic regime and show that this form of structured connectivity together with random connections provides a basic mechanism for internally generating coherent fluctuations.

Results

We study a network of N neurons in which the connectivity between neurons has two components: a random component, J, and a rank-1 structured component, , an outer product of a pair of orthogonal vectors both of which have elements of O(1) and norm , with strength parameter J₁. We restrict the elements of ξ to be binary, ξ_i = ±1, which will be important for some of the results to come, and we will comment on when this restriction can be relaxed. We can think of the row vector, ν^T, as an “output mode” performing a read-out of the network activity, and the column vector, ξ, as a corresponding “input mode” along which the output mode activity is fed back to the network (Fig 1A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Random network with structured connectivity generates coherent chaos.

(A) Network schematic showing neurons connected via random matrix J and rank-one structured connectivity. The structured component is represented schematically as a loop with drive through the output mode, ν, and feedback through the input mode, ξ. In our model these two vectors are orthogonal. The standard deviation of the random component is given by and the strength of the structured component is . (B) Sample network dynamics without structured component, i.e. J₁ = 0. Colored traces show a random collection of ten neural activity patterns, ϕ_j, black trace shows the coherent mode activity, , which exhibits only miniscule fluctuations. (C) Sample dynamics for J₁ = 1. Coherent mode displays substantial fluctuations. (D) Coherence, χ (definition in text), as a function of the strength of structured connectivity component, J₁ for small values of J₁. Simulation and theory (valid in the weak structured connectivity regime—J₁ ≪ g) shown for both g = 1.5 and g = 2. Bars show standard deviation over 60 realization of the random connectivity. (E) Passive coherent chaos. With weak structured connectivity fluctuations of the coherent mode follow the fluctuations of the independent residual components. Normalized autocorrelation of the coherent component of the current, , in red circles. Average normalized autocorrelation of the residuals, q_δ(τ), in blue ‘x’s. Both are averaged over 60 realizations of the random connectivity with J₁ = 0.1. Prediction from theory in black. N = 4000 and g = 2 in all panels unless stated otherwise.

https://doi.org/10.1371/journal.pcbi.1006309.g001

The random component of the connectivity is given by the matrix J consisting of identically distributed independent Gaussian elements with mean 0 and variance , where g is an O(1) parameter.

The state of each neuron is defined by its synaptic current, h_i(t), with its firing rate given by ϕ_i ≡ ϕ(h_i(t)), where ϕ is a sigmoidal function. For some later results it will be necessary to assume that ϕ has a first derivative that is an even function. We therefore assume here for concreteness ϕ(h) = tanh (h) unless otherwise noted, and we comment on when this assumption can be relaxed.

The dynamics of the synaptic current vector, h, is given by (1) We have scaled the strength of the structured connectivity such that when the scaling-parameter J₁ ∼ O(1) the contribution of the structure to individual synapses is of the same order of magnitude as the typical random connection.

We will be particularly interested in the coherent activity and the coherent current, i.e. the spatial overlap of both the firing rate and the synaptic current with the input mode. These are defined respectively as (2)

We define the measure of spatial coherence as the fraction of the total power of the network current h_i that is shared along the input mode, ξ: (3) where 〈〉 represents average over time. This is a useful measure as it varies from 0 to 1, and will yield for entirely independent, uncoupled fluctuations, and χ = 1 for complete synchrony along ξ.

Without the structured component, i.e. J₁ = 0, the network exhibits a phase transition at g = 1 from a zero fixed point to chaos [15]. In the chaotic state the randomness of the connectivity decorrelates the input current and yields an asynchronous state in which neurons fluctuate with negligible correlations such that both and are nearly constant in time, and (Fig 1B).

Setting J₁ = 1, we observe significant correlations along the input mode, such that the coherent mode activity, , fluctuates significantly. For this network with N = 4000 we find the coherence χ ≈ 0.4, which is about 25 times larger than is observed in the asynchronous state. (Fig 1C).

In order to analyze this system we decompose the dynamics of the synaptic currents into coherent component, , and the vector of residuals currents, . We also decompose the activity into its coherent component, , and vector of residual activity, . Because of the orthogonality between input and output mode, the output mode ignores the coherent component, , and projects only the activity of the residuals: ν^T ϕ = ν^T δ ϕ.

By projecting the dynamic equations (Eq 1) onto the input mode, ξ, on the one hand, and onto its orthogonal complement on the other, we can write decomposed dynamics for the coherent component, , and the residual synaptic current vector, δh: (4) (5) where the effective connectivity matrix in the residual dynamics (discussed more fully in Methods) is (6) which projects the output of J into the (N − 1) dimensional subspace orthogonal to ξ. This guarantees that the decomposed dynamics satisfy the constraint ξ^T δh = 0. For most of what follows this constraint can be ignored as it contributes only to each synaptic current. Nevertheless it plays a role later and we introduce it here for completeness. In Eq 4 we have ignored the small projection of the random connectivity along the coherent mode () and we discuss this the impact of making this approximation in S1 Appendix. Note that even in this approximation, the nonlinearity of ϕ in these equations couples the coherent and residual degrees of freedom.

In order to attempt to solve the approximate system we could assume that fluctuates according to some known random process and then consider the dynamics of the δh_i with the firing rate of individual units as given by . However, for general J₁ we are unable to analytically close the loop and self-consistently compute the statistics of .

Random connectivity with “row balance”

We observe that as J₁ increases to values near g, the network displays significant variability in the dynamics from realization to realization. The coherent mode autocorrelation function, , for example, is no longer self-averaging (Fig 2A). As we increase system size, N, we find that the realization-to-realization variability in saturates to a finite value (Fig 2B). Moreover, as J₁ increases we observe realization-dependent transitions out of chaos to either fixed points or limit cycles (Fig 2C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Row balance preserves chaos, and increases coherence.

(A) Full coherent mode autocorrelation, of 60 individual realizations. Thick black line shows average over realizations. Network with random connectivity J without row balance exhibits significant difference between realizations. (B) Coherent mode variance, , as a function of network size for 300 individual realizations. Gray line shows average and gray region with black boundary shows one standard deviation over realizations. Without row balance the standard deviation (over realizations) saturates to a finite value as the network size increases, indicating that the variability between realizations is not a finite-size effect. (C) Without row balance a network with moderate structured connectivity (J₁ = 2.5) exhibits a fixed point. (D)-(F) Same as (A)-(C) respectively, but network has “row balance” random connectivity, . (D) Individual realizations of Δ(τ) are all very close to the average. (E) With row balance the standard deviation of Δ(0) over realizations shrinks with N, suggesting that the variability between realizations is a finite-size effect. (F) Same realization of J as in (C), but with row balance. Chaos is preserved. (G)-(H) Networks without row balance in blue, with row balance in red (G) Fraction of realizations (out of 30 realizations) that lead to chaotic dynamics, as a function of structural connectivity, J₁. Row balance keeps nearly all realizations chaotic. (H) Coherence, χ, as a function of J₁ computed for the realizations from (E). Row balance increases coherence. Error bars display standard deviation. Full line shows all realizations, dashed line displays average coherence restricted to the chaotic realizations. J₁ = 0.2, g = 2, and N = 4000 for all panels unless otherwise noted.

https://doi.org/10.1371/journal.pcbi.1006309.g002

The reason for this realization dependence is that as J₁ increases and the fluctuations in the coherent mode grow, feedback is generated through the interaction between the random connectivity, J, and the input and output modes, ξ and ν. First of all, J maps the coherent activity back along the input mode with a small realization-dependent component which we ignored in Eq 4 driving feedback directly to the coherent current, . Secondly, J maps the coherent activity along the the residuals in a realization-dependent direction biasing the residual fluctuations δh. This direction in turn will have realization-dependent component along the output mode, ν, and therefore the coherent activity may additionally drive feedback by pushing the residual fluctuations along the output mode. See S1 Appendix for more details.

To suppress the strong realization dependence of the dynamics, we refine the random connectivity matrix by defining (11) This subtracts from each row of J its weighted average along the input mode. This “row balance” subtraction has been previously observed to remove realization-dependent outliers from the eigenspectrum of the full connectivity matrix [21, 22].

We find here that row balance suppresses realization-to-realization variability in the nonlinear chaotic dynamics, for example reducing the variability in the autocorrelation of the coherent mode, (Fig 2D). We observe that with row balance this variability drops as a function of increasing system size, suggesting (although not proving) that the dynamics are now self-averaging in the limit of large N, at least for these values of J₁ (Fig 2E).

The impact of row balance on the chaotic dynamics can be understood by noting that the resulting connectivity matrix, , now has a null-space, and the input mode, ξ, lies within it (). The result of row balance then is to ensure that the random connectivity matrix filters out any coherent activity fluctuations, : (12) This prevents the coherent mode activity from driving feedback to the dynamics of the coherent current (S1 Appendix).

Interestingly, we find that row balance allows chaotic fluctuations to persist for larger values of J₁, whereas without row balance a substantial fraction of realizations exhibit fixed points or limit cycles for J₁ > g (Fig 2F and 2G). Furthermore, the chaotic dynamics are more coherent with row balance than without (Fig 2H). As the structured connectivity is made stronger row balance appears to enable the dynamics to grow increasingly coherent even as chaotic fluctuations persist.

Row balance yields self-tuned coherent chaos

When we increase the strength of structured connectivity, J₁, we find that with row balance the network yields chaotic dynamics which are strikingly coherent and display switching-like macroscopic activity (Fig 3A). In contrast to the case of weak structured connectivity, the coherent mode dynamics are no longer passively driven by the fluctuations in the residual synaptic currents. This is evidenced by the normalized autocorrelation of the coherent mode, , which is no longer close to the normalized residual autocorrelation, q_δ (τ), but rather has qualitatively different dynamics including longer time-scales (Fig 3B as compared to Fig 1E).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Strong structured connectivity with row balance generates high coherence even as chaos persists.

(A) Sample activity of 10 randomly chosen neurons, ϕ_j, and coherent mode activity, , in black. Strong structured connectivity with row balance subtraction to the random component of connectivity yields chaotic activity that is highly coherent with switching-like behavior. (B) Normalized autocorrelation of coherent mode, , in red. Average normalized autocorrelation of the residuals, q_δ (τ), in blue. Shaded regions show standard deviation over 25 initial conditions of the same connectivity. Strong structured connectivity yields coherent mode dynamics that are qualitatively different from those of the residuals. J₁ = 15.8 for both (A) and (B), compare to Fig 1C and 1E, respectively. (C) Coherence, χ, as a function of J₁ is independent of network size. Coherence appears to approach 1 as J₁ is increased demonstrating that chaos persists even as fluctuations in the residuals shrink (See also S4 Fig). Coherence is averaged over 30 realizations of the connectivity for each N, excluding the few fixed points and limit cycles that occur for larger J₁ (2 out of 30 or less for the largest values of N). (D) Largest Lyapunov exponent as a function of J₁. Thick line shows average over 10 realizations, small dots show values for individual realizations, and shaded region is standard deviation. All but a small fraction of realizations are chaotic, even in the region where χ > 0.9. N = 4000 and g = 2 in all panels unless noted otherwise.

https://doi.org/10.1371/journal.pcbi.1006309.g003

We find that the coherence, χ, increases steadily as a function of the strength of structured connectivity, J₁, and notably it is independent of system size (Fig 3C).

To check whether this highly coherent state is chaotic, we calculate the largest Lyapunov exponent and verify that the dynamics are indeed chaotic for a vast majority of realizations even as the fluctuations are highly coherent (Fig 3D).

We now examine the qualitative changes in the chaotic state of the network with row balance as J₁ increases. We observe that for J₁ ≲ 1 the fluctuations are unimodal with an approximate Gaussian shape. The temporal fluctuations are dominated by a single time constant as with J₁ = 0 (Fig 4A and 4B). On the other hand, for larger values of J₁ the fluctuations deviate dramatically from Gaussian and instead become sharply bimodal (Fig 4A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Self-tuned coherent chaos.

(A)-(B) Comparison between weak structured connectivity (Left: J₁ = 0.8) and stronger structured connectivity (Right: J₁ = 15.8) Both with N = 4000. (A) Histogram of values of the coherent mode current, . Mild structured connectivity yields coherent fluctuations with a peak at zero and a distribution that appears not far from Gaussian. For stronger structured connectivity the histogram is clearly non-Gaussian and highly bimodal. (B) Top: Sample activity of 10 randomly chosen neurons, ϕ_i (t) and coherent mode activity, . Middle: Speed of network during same epoch of activity, (definition in text). Bottom: Instantaneous population variance of the residual currents δh_i (t). For mild structured connectivity, fluctuates around zero (top), speed is roughly constant throughout the trial (middle), residual currents maintain large variance throughout (bottom). On the other hand, for stronger structured connectivity, there is state-switching between bouts of high and low coherent-mode activity (top), these same bouts are associated with vanishing speed (middle), and with small residual currents (bottom). Gray shaded regions show epochs of speed lower than 0.18, which was the lowest instantaneous speed achieved without structured connectivity. (C) Statistical mode (most frequent value) of as a function of J₁. The results indicate a crossover to self-tuned coherent chaos, defined by the bimodal peaks of reaching a constant value. The crossover occurs very rapidly and independently of N. Dashed line shows . (D) The statistical mode of as a function of g with fixed J₁ = 10 and N = 8000. Dashed line shows . In all other panels g = 2.

https://doi.org/10.1371/journal.pcbi.1006309.g004

Furthermore, we find that with larger values of J₁ the network exhibits intermittent switching between two different values of (Fig 4B, top right). We observe that the dynamics at both of these values of are qualitatively distinct as reflected by the speed of the dynamics and the level of coherence. We define the speed of the dynamics as the norm of the vector of first time-derivatives per neuron: . Similar measures have been used to find fixed points and slow dynamics in highly non-linear dynamics [23]. We observe that the two distinct values of are both associated with vanishing speed (Fig 4B, middle right). Additionally they are associated with very high levels of coherence, i.e. small residuals as quantified by the population variance, (Fig 4B, bottom right). Evidently, in the network with larger values of J₁ there are two distinct states with slow dynamics and high levels of coherence, and the network switches rapidly between these two. In contrast, with mild structured connectivity both the speed and the variance of the residuals remain roughly constant throughout the trial (Fig 4B, left).

To gain insight into the emergence of the switching dynamics, we examine the most frequent value of as a function of J₁ and we find that there is a rapid crossover from a state with unimodal fluctuations around zero, to a state with bimodal peaks and that the most frequent value of saturates quickly with J₁ and then remains constant (Fig 4C). What is the nature of this regime? And what determines the two values of that come to dominate the dynamics?

Because the bouts of slow dynamics are associated with small residuals we linearize the dynamics around and assume that δh_i ≪ 1. This gives (13) Note that we have made use of the fact that ξ_i = ±1 and that ϕ is an odd function. Because ξ is in the null-space of the row-balanced connectivity (, as discussed above) the residual dynamics (Eq 5) become (14)

We observe that in these linearized dynamics the coherent mode current, , plays the role of a dynamic gain via the slope of the transfer function, ϕ′.

In the linearized dynamics we turn to the eigenvectors, u^(k), of and decompose the residual dynamics according to δh = ∑_k c_k u^(k). Given an instantaneous value of , the independent dynamics of the eigenmodes are (15) where λ_k is the kth eigenvalue of . (We show in S1 Appendix that has the same eigenvalues as ).

By the well-known circular law of random matrices, the leading eigenvalue, λ₁, has real part approximately equal to g. If is small then and there are many modes that diverge exponentially. If is large then and then all the modes decay exponentially. However, there are two critical values of which yield marginal and therefore slow dynamics for the leading mode, c₁. These are the values, , for which the slope of the transfer function is equal to : (16) If and the residuals are small then the time constant of fluctuations in the leading eigenmode, , are very long.

Indeed we find that the most frequent value of as a function of g fits the curve h^c(g) very well (Fig 4D).

We conclude that the switching between two states each with slow dynamics and a high level of coherence observed in Fig 4B reflects a distinct dynamical regime of self-tuned coherent chaos. In this regime the coherent mode self-adjusts to a critical value so that the dynamics of the small residuals are near-marginal, giving rise to slow dynamics. The above linearized dynamics (Eq 14) are not exact and therefore non-linear interactions eventually destabilize the system and precipitate a state-switch. Nevertheless, the linearized dynamics dominate the dynamics of the small residuals during the bouts of high coherence.

As observed in Fig 4C, when we increase J₁, the most frequent value of rapidly increases until it saturates at a value very near to . Moreover, we find that this crossover to the regime of self-tuned coherent chaos occurs at moderate values of J₁ (on the order of g), independently of network size.

Notice the crucial role of row balance in facilitating self-tuned coherent chaos: row balance filters out the direct contribution of the coherent mode activity to the dynamics of the residuals and enables the coherent mode to act as a dynamic gain. The coherent mode then self-adjusts to cancel the leading eigenvalue of and yield bouts of slow, highly coherent dynamics.

We note that we can loosen the constraint on the symmetry of the transfer function and allow any smooth sigmoidal transfer function if we restrict the input mode to be uniform, ξ_i = 1 for all i. We show an example of the self-tuned coherent chaotic state for a non-symmetric transfer function in S1 Fig.

Symmetry breaking in the self-tuned chaotic regime and transition to fixed point.

The example of Fig 4A illustrates dynamics which reside in the positive and negative coherent states with equal frequency, maintaining the h_i → −h_i symmetry of the underlying dynamic equations. We observe that in many realizations this symmetry is violated at the single trial level for sufficiently strong J₁, as demonstrated in Fig 5A.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Realization-dependent symmetry breaking in the self-tuned chaotic regime.

(A) Sample traces of coherent mode current, , (top) and histogram of values of (bottom) from a connectivity realization with real eigenvalue for J₁ = 40,60,80 increasing from left to right. Dynamics exhibit pronounced asymmetry. (B) Absolute value of the time-average coherent mode current, , as a function of J₁. Each colored line represents a single connectivity realization, averaged over 10 initial conditions. For many individual realizations, is significantly non-zero over a large range of values of J₁, while still not arriving at fixed point value (displayed by dashed line). We display the 37 realizations with real leading eigenvalue out of 100 total realizations from this set of trials. Thick black line shows average over those realizations. Dashed line shows . N = 8000 for all panels.

https://doi.org/10.1371/journal.pcbi.1006309.g005

We measure the asymmetry in a single trial as the absolute value of the time-averaged coherent activity, , and we find that asymmetry grows gradually with J₁ throughout the chaotic regime, and at different rates for different realizations (Fig 5B).

As is evident in Fig 5B, for many realization the level of asymmetry increases with J₁ until it reaches a maximal value of (dashed line). At this point, the system spends all the time at one of the possible states, suggesting a transition to fixed point.

Indeed, as has been previously reported by Garcia Del Molino et al [20], realizations of that have a real leading eigenvalue (see also Methods and S1 Appendix) yield a fixed point equation for the above linearized dynamics (Eq 14): (17) The fixed point requires so that all but the leading eigenmode decay to zero, and the resulting fixed point for is marginally stable.

Indeed, we find that realizations of with a real leading eigenvalue undergo a transition to a fixed point upon sufficient increase of J₁. Similar to the preceding chaotic state, the fixed point is highly coherent, with very small residuals. Furthermore, it exhibits the hallmarks of the self-tuned coherent state: the value of the coherent mode is close to , the residuals are aligned with the leading eigenvector of , and the the leading eigenvalue of the Jacobian matrix at the fixed point is close to zero independently of J1 (S2 Fig).

Oscillatory fluctuations in self-tuned chaos and transition to limit cycle.

The above symmetry breaking and transition to fixed point is observed only for some of the realizations of J. In most of the other cases, rather than symmetry breaking, we observe an increased oscillatory component in the chaotic dynamics. This is reflected in the presence of a large second peak in the normalized autocorrelation function, (Fig 6A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Realization-dependent oscillatory imprint on the self-tuned chaotic regime.

(A) Sample traces (top), and normalized autocorrelation (bottom) of coherent mode current for a connectivity realization with complex eigenvalue for J₁ = 25,30,35 increasing from left to right. Dynamics exhibit pronounced oscillatory power and the autocorrelation exhibits a pronounced peak near the same frequency that will dominate the limit cycle for larger J₁. (B) Height of second peak of the autocorrelation of the coherent mode input as a function of J₁. Each colored line represents a single connectivity realization, averaged over 10 initial conditions. For many realizations, there is a significant second peak in the autocorrelation over a long range of values of J₁ well before a limit cycle is reached. We display the 63 realizations which had complex leading eigenvalue out of 100 in this set of trials. Thick black line shows average over those realizations. (C) Observed period of oscillatory chaos vs phase of leading eigenvalue for 181 realizations from which we were able to measure an oscillatory period with chaotic fluctuations (out of 196 realizations with complex leading eigenvalue in this set of trials. In order to confine to realizations and values of J₁ that yielded chaos, we restrict to those with second peak of autocorrelation less than 0.8. These had average height of second peak over all realizations: 0.5). Dotted line shows prediction from theory: . The bulk of realizations are very well predicted although a notable fraction are not. The median error of prediction was 7.75 (average period over these realizations: 231, std: 212). N = 8000 for all panels.

https://doi.org/10.1371/journal.pcbi.1006309.g006

The origin of these oscillations can be traced to the nature of the leading eigenvalues of . As also reported in [20], if the leading eigenvalue of , λ₁, is complex then there is no fixed point solution. Instead, as we show in S1 Appendix, if we assume that undergoes a limit cycle with period T, we find that the period must satisfy (18) and additionally the period-average of must be equal to the critical value (see S1 Appendix).

Indeed, we observe that in the case of a with complex leading eigenvalue, most realizations exhibit oscillatory fluctuations. The height of the second peak in the autocorrelation, , grows gradually with J₁ with a realization-dependent rate (Fig 6B). Furthermore, for most individual realizations the period of the dominant oscillatory peak in the autocorrelation, even within the chaotic regime, is well-predicted by the above theoretical prediction for T (Eq 18, Fig 6C).

For individual realizations, there is a sufficiently large J₁ beyond which the second peak of reaches 1 and the dynamics transition to a pure limit cycle (S3 Fig). We note that some realizations with a real leading eigenvalue also exhibit oscillatory components in their chaotic dynamics for certain values of J₁, which we presume relate to complex subleading eigenvalues, but these do not exhibit a transition to pure limit cycle.

Realization-dependence and system-size scaling of the transition out of chaos.

We now consider the critical value, , of the strength of structured connectivity that yields a transition out of chaos. In contrast to the case without row balance, we find that in the row-balanced network the transition out of chaos occurs at values of J₁ scaling at least as . However, the particular value of varies considerably across realizations (S2 and S3 Figs).

The case of a real leading eigenvalue λ₁ ≈ g and the associated transition to fixed point provides a starting point for analyzing the transition out of chaos (a similar argument is made in Methods and S1 Appendix for the case of complex leading eigenvalue). In the limit of small residual currents, δh_i ≪ 1, the fixed point equation for the coherent mode current (Eq 4) is . Applying the fixed point requirements derived above that and , where u⁽¹⁾ is the leading eigenvector of , we find (as also reported in [20]) that (19) Because the above fixed point assumes where u⁽¹⁾ has norm 1, we must have be no larger than O(1). Therefore we expect that the critical value, , for a given realization will require yielding a sufficiently small value of in Eq 19. This suggests that scales roughly as inverse of the overlap between the leading eigenvector, u⁽¹⁾, and the output mode, ν. Indeed, we show numerically that the critical value for fixed g is negatively correlated with |ν^T u⁽¹⁾| (Fig 7A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Sufficiently strong structure yields transition out of chaos despite row balance.

(A) Scatterplot of the logarithm of transitional value, , vs the absolute value of the projection of the output mode, ν, on the leading eigenvector, u⁽¹⁾ for 300 connectivity realizations with N = 8000. r² = 0.29. (B) Fraction of realizations displaying chaotic activity as a function of the rescaled structured connectivity: . With this scaling the curve appears to be independent of N.

https://doi.org/10.1371/journal.pcbi.1006309.g007

We next ask how scales with system size. Since the typical value of |ν^T u⁽¹⁾| is O(1) we would naively expect that the typical transition might occur for . However, we observe numerically that the fraction of realizations exhibiting chaotic dynamics for a given value of J₁ appears to scale as and not as as expected (Fig 7B). For finite N chaos appears to be lost for some with , and the particular is highly realization-dependent. An analytical derivation of the actual value of requires a more comprehensive study of the network’s stability.

From our numerical work it appears that in the limit of large system size chaos persists for all values of J₁, for almost all realizations. Indeed for a network with N = 16000, for example, we increase J₁ up to values around 100 and observe that chaos persists for most realizations and coexists with a very high degree of spatial coherence. The coherence measure, χ, reaches values higher than 0.96 even as fluctuations persist (S4 Fig). Thus we conjecture that in the limit N → ∞, for almost all realizations, increasing J₁ indefinitely will yield self-tuned coherent chaos with χ → 1.

Multiple modes of coherence

We generalize our model in order to construct a network with multiple modes of coherent activity. In this extension we take the structural component, M, to be a low-rank matrix comprised of the sum of outer products between input modes, ξ^(k), and output modes, ν^(k). We require that the subspace spanned by the input modes be orthogonal to the subspace spanned by the output modes. Using singular value decomposition we can find orthogonal bases for each of these subspaces, such that without loss of generality we can additionally assume that each pair of input modes and each pair of output modes are orthogonal. In sum we assume ξ^(k) ⊥ ξ^(j) and ν^(k) ⊥ ν^(j) for all j ≠ k, and ξ^(k) ⊥ ν^(j) for all j, k.

We write the structured component of connectivity as (20) where we have introduced a parameter J₀ that controls the overall strength of the structured connectivity, and a set of parameters w_k satisfying that determine the relative weight of the different modes.

We can extend our schematic representation and think of each row vector as a separate output mode connected in a feed-forward-like manner to the input mode, ξ^(k), (Fig 8A) and then decompose the dynamics into the residual dynamics identical to the above (Eq 5) and the dynamics of the coherent activity along each separate input mode (by approximation analogously to Eq 4): (21) where each separate coherent current is given by .

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Coherent chaos along multiple modes.

(A) Schematic of network with three coherent modes displaying effective output modes, ν^(k), and input modes, ξ^(k), each of which are orthogonal to all others. (B) Matrix of Pearson correlation coefficients between firing rates, ϕ_i, of pairs of neurons in a network with three coherent modes and J₁ = 1. (C) Sample activity trace displaying sample single neuron activities and in thicker lines, three coherent mode activities, . (D) Generalized coherence, χ^(k), as a function of J₁ for 2,3,4 modes in the regime of passive coherence. Dashed line displays theory. (E) Sample activity traces display extreme coherence for two coherent modes. (F) Generalized coherence for two coherent modes with row balance random connectivity (in blue) as a function of J₁ extends to near complete coherence, while chaos persists. Compare network with one coherent mode (in black dashed line). Bar shows standard deviation over 100 realizations. For panels (B), (C), (E) N = 2000. For panels (D) and (F) N = 4096.

https://doi.org/10.1371/journal.pcbi.1006309.g008

The analytical results found above for the regime of passive coherence can be directly extended to the case of multiple modes. In particular, in the limit where J₀ ≪ g, the separate coherent modes are independent of each other and driven passively by the residual fluctuations (Fig 8C) such that (22) where is the autocorrelation function of the kth coherent mode.

The resulting covariance matrix, C_ij ≡ 〈ϕ_i ϕ_j〉, has a low-rank structure which is shaped by the input modes (Fig 8B). In particular, by Taylor expanding ϕ_i around δh_i: (23) Since the residuals are on average uncorrelated and the distinct input modes are orthogonal, we find using Eq 22 that to leading order: (24)

We generalize our measure of coherence to measure the d-dimensional coherence, or the fraction of total power which is shared along the d input mode directions: (25) and we find that for J₀ ≪ g and finite d (26)

Numerically, this prediction holds well for up to for at least up to d = 4 as we show in Fig 8D, and we expect it to hold for larger d as well.

We note that in the regime of passive coherence, just as in the case of a single coherent mode, we can relax the restrictions on ξ^(k) and ϕ: Our results hold for ξ^(k) any norm vector orthogonal to ν^(j) ∀j and ξ^(j) ∀j ≠ k, and also for ϕ any sigmoidal function.

In addition, we can generalize row balance by subtracting the weighted row-average for each input mode such that every ξ^(k) will reside in the null space of the new connectivity matrix, : (27)

For d > 2 this generalized row balance does not appear to preserve chaotic fluctuations and instead fixed points or limit cycles appear for J₀ ∼ O(1).

Intriguingly, for d = 2 we observe that with generalized row balance the chaotic regime persists as the structured connectivity is strengthened and the dynamics become increasingly coherent. The dynamics display switching-like activity in which at any time one of the two coherent modes appears to be near the critical value h^c while the other mode is near zero (Fig 8E). It appears that the generalized coherence approaches 1 while chaos persists (Fig 8F) such that we conjecture that just as in the case of d = 1, here too in the large N limit χ⁽²⁾ → 1 as J₀ → ∞.

Discussion

Coherent fluctuations are prevalent in cortical activity ranging in spatial scale from shared variability in membrane potential and spiking in local circuits to global signals measured across the scalp via EEG or across voxels via fMRI [3, 4, 9, 11]. Constructing a model that intrinsically generates coherent fluctuations has been a challenge to theorists.

We have studied the intrinsic generation of coherent chaotic dynamics in recurrent neural networks. Our model consists of rate-based neurons whose recurrent connections include a structured component in addition to random connections. The structured component is a low-rank connectivity matrix consisting of outer products between orthogonal pairs of vectors which allow local fluctuations to be summed along an output mode, amplified and projected to an input mode resulting in coherent fluctuations. The orthogonality of input and output mode effectively embeds a purely feedforward structure within the recurrent connectivity, thus avoiding feedback of the coherent fluctuations along the input mode.

In the regime where the structured component is weak, the local synaptic currents are effectively uncoupled from the coherent mode activity and their dynamics are similar to that of a random network with no structured component at all. The local fluctuations are summed by the structured component of connectivity to drive the coherent mode, which follows those fluctuations in a passive manner. In this regime of passive coherent chaos we derive a perturbative dynamical mean-field theory following [15, 19] which shows that the coherence grows linearly with the ratio of the strength of structured connectivity to the random connectivity. We show that this analysis extends to multiple modes of coherent activity yielding a finite-rank covariance pattern for the coherent fluctuations.

For moderate strength of structured connectivity the network exhibits significant realization-dependence and most realizations transition to either a fixed point or a limit cycle. A realization-dependent theory of these transitions is beyond the scope of this work. We add a row-balance constraint, placing the input mode in the null-space of the random connectivity matrix, and we observe that this constraint preserves chaos, reduces the variability between realizations, and increases the level of coherence.

With row balance, increased strength of structured connectivity yields a crossover to a distinct regime of self-tuned coherent chaos. In this regime the network undergoes Up-Down-like switching between two states each of which are characterized by slow, highly coherent dynamics. We show how row balance facilitates this regime by enabling the coherent mode to act as a dynamic gain to the dynamics of the residual currents. Consequently, intermittent marginal dynamics emerge as the coherent mode self-adjusts to one of two critical values. Interestingly the crossover to this self-tuned coherent regime happens for moderate strength of structured connectivity, independently of network size.

In the regime of self-tuned coherent chaos, realization-dependent qualitative differences begin to emerge with increasing strength of structured component, J₁. For realizations of the row-balanced random connectivity with real leading eigenvalue, symmetry-breaking emerges such that individual initial conditions yield trajectories that spend more time near one of the critical values of the coherent mode than the other. For realizations with complex leading eigenvalue, oscillatory fluctuations begin to emerge. The frequency of these oscillations is well predicted by the phase of the leading eigenvalue. Note that we have not addressed the question of the necessary scaling of J₁ for the emergence of realization-dependence in the chaotic regime for the limit of large system size.

As structured connectivity is further strengthened chaos persists even as coherence continues to increase until the dynamics are dominated almost entirely by the one dimensional fluctuations of the coherent mode. For a finite network, above some critical strength of the structured component the system converges to either a fixed point or a limit cycle, depending on the leading eigenvalue of the row-balanced random connectivity. Our numerical work suggests that the critical strength of structured connectivity grows with the system size, most likely scaling as (at least as ). Hence we conjecture that for the scaling of the strength of structured connectivity presented here, as the network size diverges coherent chaos persists independent of J₁ for most realizations, and the level of coherence can be brought arbitrarily close to 1.

Importantly, in the regime of weak structured connectivity and passive coherence some of the assumptions of our model can be loosened. First, in this regime row balance on the random connectivity is not necessary. Additionally, we need not require the input mode be binary but rather any general pair of orthogonal vectors can serve as input and output mode. Moreover we can loosen the restriction on orthogonality: a random pair of vectors can be used without qualitative impact on the dynamics presented here because the contribution of the realization-dependent overlap between the two vectors in this regime will be negligible relative to the typical contribution from the full-rank random connectivity, J.

On the other hand, achieving self-tuned and highly coherent chaos requires the network be finely-tuned to a high degree. The orthogonality of the input and output modes is not enough in order to achieve highly coherent chaos because of interactions between the random and structured components of connectivity. We therefore constrain the random component to satisfy row balance by ensuring that the input mode of the structured connectivity be in the null-space of the random component of connectivity. In addition to row balance, the self-tuned coherence regime depends on the choice of the single-neuron transfer function and input-mode vector. In the case where the transfer function is an odd function, such as tanh used throughout the main text, the input mode can be any binary vector. Otherwise, high coherence is achieved only for a uniform input mode, ξ_i = 1∀i. In the latter case, the theory developed here predicts coherent fluctuations that switch between two non-symmetric values of the coherent mode, corresponding to the two points where where the slope of the transfer function equals , and we have verified this numerically (S1 Fig).

An interesting question is whether the particular structure of the connectivity matrix in our model can be achieved by a biologically plausible synaptic learning rule. Prior studies of sequence generation have constructed learning rules that yield connectivity which is comprised of outer-products of random vectors [24, 25] and these could form the basis for learning the necessary orthogonal rank-one structure. Plausible learning rules for yielding balanced excitation-inhibition dynamics [26–28] could potentially provide a foundation for learning row balance. It is thus plausible that the constraints of our model can be achieved by an appropriate synaptic learning, especially for the more robust regime of mild coherence. On the other hand, it is unclear to us whether the high degree of fine tuning required for the self-tuned coherence regime can be achieved by a biologically plausible learning rule. Investigating candidates of appropriate learning rules for generating coherent chaos, is beyond the scope of this work.

In the case of a uniform input mode the model can be constructed as an excitation-inhibition network, for example with half the neurons defined as excitatory by setting ν_i = 1 and the other half defined as inhibitory via ν_i = −1 (or a larger inhibitory value to compensate for a smaller fraction of inhibitory neurons). From this perspective the coherent fluctuations, in particular in the regime of passive coherence, can be understood within the framework of dynamic excitation-inhibition balance [13]. In this case the pair of balance equations are degenerate and constrain only the mean excitatory population rate to be nearly equal the inhibitory rate, but otherwise leave the overall mean rate unconstrained. Local residual fluctuations yield only small differences in mean population rates, thus leaving the balance satisfied, but these small differences drive significant coherent fluctuations because of the strong balanced connectivity. In the general setting of excitation-inhibition balance the pair of balance equations fully determine the mean rates to leading order and no coherent fluctuations are possible without introducing shared fluctuations in the external drive. We note that excitation-inhibition networks in the literature have sometimes been constructed yielding degenerate balance equations. As we have shown here, such choices have dramatic impact on the dynamics and the results should not be assumed to be generalizable.

In parallel to our study a pre-print has been published which explores a very similar model [29]. The authors observe a similar phenomenon as the self-tuned coherence studied here, and attempt to explain it by an iterative numerical solution of locally Gaussian dynamic mean-field equations. They do not make the role of row balance clear in their analysis. In contrast, we have focused on analytical solutions in the limit of both weak and strong structured connectivity, deriving a perturbative dynamic mean-field solution for the regime of weak structure and passive coherence. As we have shown here row balance is critical for moderate structure and self-tuned coherent chaos. Additionally we have shown that a full understanding of the properties of the highly coherent regime requires a realization-dependent mean-field analysis. In particular, we have explained that the leading eigenvalue of the row-balanced random connectivity matrix impacts qualitative features of the chaotic dynamics, yielding either broken symmetry or oscillatory fluctuations. Furthermore the critical strength of structured connectivity that leads to a transition to either fixed point or limit cycle is correlated with the extent of overlap between the leading eigenvector and the output mode.

As mentioned in the main text and introduction, a previous study also explored the case of a single orthogonal E-I structured component [20]. They derived the fixed point and limit cycle solutions which we reviewed here, but did not focus on the chaotic regime and they did not discuss the role of network size in the transition out of chaos. Our focus here was on the chaotic regime, both the emergence of coherence for small structured connectivity and the imprint of the non-chaotic regime on the chaotic dynamics for moderate structured connectivity.

A separate study has claimed to observe coherent activity in excitation-inhibition networks of spiking neurons [30]. A study of the dynamics of spiking neurons is beyond the scope of our work, although we would conjecture that coherent activity would arise with orthogonal, rank-one E-I structure in that setting as well.

Previous work has shown how shared inputs from external drive can drive correlated fluctuations in excitation-inhibition networks [17, 18]. In our current work, in the context of rate neurons, we show that such correlated fluctuations can be generated internally by a recurrent network without external drive. In order to avoid either saturation or pure oscillations the coherent activity mode must not drive itself through a feedback loop. In order to achieve this it is necessary that the structured component embed an effectively-feedforward projection between a pair of orthogonal modes. In parallel, Darshan et al [31] have developed a theory for internally generated correlations in excitation-inhibition networks of binary units. The underlying principle is similar: the recurrent connectivity embeds a purely feedforward structure.

We note that the structured component of connectivity in our network is non-normal. The dynamics of non-normal matrices have drawn a fair amount of interest with suggested functional impact on working memory [32, 33] and amplification [34]. Non-normal matrices embed feedforward structure within recurrent connectivity, and the resulting dynamics even in a linear system are not fully determined by the eigenspectrum but depend on the structure of the corresponding eigenvectors [35]. It has been shown that E-I networks are generally non-normal, and that rank one E-I structure amplifies small differences between excitatory and inhibitory rates driving a large common response [34, 36]. This amplification is related to the way in our network, small fluctuations of the residuals are summed along the output mode and drive coherent fluctuations along the input mode in the regime of passive coherence, but these fluctuations are internally driven by the non-linear dynamics whereas the dynamics in those previous studies were linear. As pointed out in [37], a structured component such as in our model is purely feedforward and can be considered an extreme case of non-normality as it has only zero eigenvalues and therefore all the power in its Schur decomposition is in the off-diagonal. The results here depend on this property and cannot be extended to connectivity with only partial feedforward structure.

Rate model dynamics with a rank-one structured component have been studied in depth recently [38, 39]. Since these works focused on time-averaged activity and not fluctuations they did not observe coherent activity in the case of an outer product of a pair of orthogonal vectors as studied here. These works also differed in that the strength of the structured connectivity was scaled as . This scaling is similar to our limit of weak structured connectivity, J₁ ≪ 1, and guarantees that dynamic mean-field theory holds in the limit of large system size, but in that scaling coherent fluctuations will appear only as a finite-size correction.

It has been previously observed [21] and then proven [22] that adding an orthogonal outer-product to a random matrix generates realization-dependent outliers in the eigenspectrum, and furthermore that these outliers are be removed by row balance. It has been previously observed that performing such a subtraction has significant impact on the resulting dynamics [20, 40]. Yet the relationship between the change in eigenspectrum and the dynamics has not been made clear beyond the basic observations regarding the stability of a fixed-point at zero. Here we suggest that the impact of row balance on the chaotic dynamics is not directly related to the eigenspectrum but that this adjustment should be thought of as effectively subtracting the coherent-mode activity from each individual neuron, thus -preventing feedback loops to the coherent mode. We show that row balance enables the emergence of slow residual dynamics with the coherent mode playing the role of dynamic gain, and that it is crucial for the emergence of self-tuned chaos and highly coherent dynamics.

In conclusion we have presented a simple model which generates coherent chaos in which macroscopic fluctuations emerge through the interplay of random connectivity and a structured component that embeds a feedforward connection from an output mode to an orthogonal input mode.

Methods

Supporting information

Acknowledgments

We thank Rainer Engelken for advice on numerically calculating Lyapunov exponents, and Yu Hu for very helpful comments on the manuscript.