Coherent chaos in a recurrent neural network with structured connectivity

We present a simple model for coherent, spatially correlated chaos in a recurrent neural network. Networks of randomly connected neurons exhibit chaotic fluctuations and have been studied as a model for capturing the temporal variability of cortical activity. The dynamics generated by such networks, however, are spatially uncorrelated and do not generate coherent fluctuations, which are commonly observed across spatial scales of the neocortex. In our model we introduce a structured component of connectivity, in addition to random connections, which effectively embeds a feedforward structure via unidirectional coupling between a pair of orthogonal modes. Local fluctuations driven by the random connectivity are summed by an output mode and drive coherent activity along an input mode. The orthogonality between input and output mode preserves chaotic fluctuations by preventing feedback loops. In the regime of weak structured connectivity we apply a perturbative approach to solve the dynamic mean-field equations, showing that in this regime coherent fluctuations are driven passively by the chaos of local residual fluctuations. When we introduce a row balance constraint on the random connectivity, stronger structured connectivity puts the network in a distinct dynamical regime of self-tuned coherent chaos. In this regime the coherent component of the dynamics self-adjusts intermittently to yield periods of slow, highly coherent chaos. The dynamics display longer time-scales and switching-like activity. We show how in this regime the dynamics depend qualitatively on the particular realization of the connectivity matrix: a complex leading eigenvalue can yield coherent oscillatory chaos while a real leading eigenvalue can yield chaos with broken symmetry. The level of coherence grows with increasing strength of structured connectivity until the dynamics are almost entirely constrained to a single spatial mode. We examine the effects of network-size scaling and show that these results are not finite-size effects. Finally, we show that in the regime of weak structured connectivity, coherent chaos emerges also for a generalized structured connectivity with multiple input-output modes.

We write the full dynamics without row balance: and we defineh Applying these definitions to the full dynamics (and noting that ν T φ = ν T δφ), the exact coherent mode dynamics and by subtracting these from the full dynamics of h, the decomposed dynamics are: whereĴ ≡ I − ξξ T N J as introduced in the main text in Eqns (5) and (6).
We observe that the constraint ξ T δh = 0 must be satisfied automatically by the residual dynamics (Eqn 5), and this can be confirmed by verifying that In the regime where J 1 1 the φ j are nearly uncorrelated and therefore ξ T N Jφ ∼ O 1 √ N and can be ignored.
This yields the approximate coherent dynamics presented in Eqn (4) of the main text.
In the regime with strong structured connectivity we must consider this term in Eqn 4. To that end we write φ =φξ + δφ and also write the transformation of the input mode via the random matrix as where a is a realization-dependent scalar and ξ ⊥ is a realization-dependent vector orthogonal to ξ. That yields coherent mode dynamics and residual dynamics dδh dt = −δh +Ĵδφ +φξ ⊥ (9) As J 1 increases and the fluctuations in the coherent activity,φ (t), drive feedback in two ways. First of all, J maps the coherent activity back along the input mode driving direct feedback to the coherent currenth via the term a φ .
Secondly, J maps the coherent activity in a realization-dependent direction, ξ ⊥ , orthogonal to the input mode.
This drives the residual activity fluctuations δh via the termφξ ⊥ , and this biasing of the residual fluctuations may in turn generate feedback to the coherent current through the output mode via ν T δφ.
Both of these feedback terms are realization dependent, and both of them are canceled via the row balance sub- which yields exact coherent mode dynamics and residual dynamics dδh dt = −δh +Ĵδφ And in this case the residual dynamics are again uncorrelated so that ξ T N Jδφ ∼ O 1 √ N and can again be ignored in the coherent mode dynamics.
Note thatφ no longer drives feedback to either the residual or the coherent dynamics. Nevertheless the dynamics are still coupled in both directions as δφ depends onh.
2 Perturbative Dynamic Mean-Field Theory in the Limit of Weak Structured Connectivity We derive the dynamic mean-field equations in the limit of small J 1 using a perturbative approach. We write the mean-field dynamics of the residuals as and the coherent component as where η i and m are assumed to be uncorrelated, mean-zero Gaussians. For general J 1 the assumption of Gaussianity fails, therefore we assume J 1 1.
The autocorrelation of η i is given by where we have introduced [] as the notation for averaging over realizations. We assume that φ j is independent of J ij and so the terms j = k have average zero over realizations and we get The autocorrelation of m is given by And again the j = k terms fall in the realization average so that Next we define the autocorrelation of the residuals and the autocorrelation of the coherent current and we can follow previous work [3,2] and write the dynamic mean-field equations for ∆ δ (τ ) as and for∆ (τ ) as Next we note that for J 1 g we assume that h 1 so we have Therefore we have that to leading order and then following previous results [3,2] we can write this to leading order as an integral over Gaussians: Note that it is possible to compute the sub-leading correction term as well, but for our purposes this is unnecessary.
The dynamic equation for∆ (τ ) is identical to that for ∆ δ (τ ) except with J 1 in place of g, so we conclude (as presented in Eqn (9) of the main text) that the resulting leading order autocorrelation of the coherent mode is Thus for J 1 g fluctuations in the coherent input are driven passively by the random source which is generated self-consistently by the residual fluctuations, and the resulting autocorrelation of the coherent mode is simply a scaled version of the autocorrelation of the residuals.
It is worth noting that for J 1 ∼ g the assumption of Gaussianity is broken due to the cross-correlations between the φ j .
3 Analysis of the Limit of Strong Structured Connectivity with Row

Balance
In the limit of large J 1 we assume δh i 1, and approximate φ j ≈ φ ξ jh + φ h δh j , where we have made use of the symmetry of the transfer function and the binary restriction on ξ j . Note that this linearization clearly holds without symmetric transfer function for the case of uniform ξ j = 1 as well.
Using the random connectivity with row balance constraint,J, and following the exact decomposition above (Eqns 5 and 4) this yields dynamical equations: In this regimeh acts as a dynamic gain on the local synaptic currents through φ h . Givenh the equation for the residual currents is linear and therefore their dynamics can be decomposed in the eigenbasis of the matrix where P ξ = I − ξξ T N is the projection matrix onto the the subspace orthogonal to ξ.
We observe a fine-point not noted in [1]: It may seem intuitive that the eigenvalues ofJ determine the dynamics. In fact, as we show these eigenvalues are identical to those ofĴ. However, had we ignored the constraint ξ T δh = 0 then the residual dynamics would have been determined by J and its eigenvalues, and these are not the same as those ofĴ.
We claim thatĴ = P ξ J andJ = JP ξ have the same eigenvalues. Suppose λ is an eigenvalue ofĴ with associated eigenvector u, then u must be orthogonal to ξ. If Ju is orthogonal to ξ as well, thenJu =Ĵu = λu, and we are done. Otherwise we can write Ju = λu + aξ, and thusJ (λu + aξ) = Jλu = λ (λu + aξ) so that λ is also an eigenvalue ofJ. Suppose now that λ is an eigenvalue ofJ. Again if the associated eigenvector is orthogonal to ξ then it is also an eigenvector ofĴ with eigenvalue λ and we are done. Otherwise we write the eigenvector ofJ as u + aξ and then we haveJ (u + aξ) = Ju = λ (u + aξ). ThereforeĴu = λu.
We write the eigenvectors as u (i) withĴu (i) = λ i u (i) . We write the vector of residual current as δh = i c i u (i) and note that as mentioned above u (i) ⊥ ξ, so that the constraint ξ T δh = 0 is satisfied. This yields dynamics The only (marginally) stable, non-zero fixed point is achieved with c 1 = 0 and c i = 0 for all i > 1. And the fixed-point equation is This fixed point only exists if λ 1 is real, and yields a fixed-point requirement forh * : In order to close the loop we turn to the fixed point equation for the coherent dynamics. Ignoring the term φ h ξ T N Jδh, which yields an O 1 √ N correction we find: which in turn, using δh * = c 1 u (1) , yields a solution to leading order for c 1 : J1ν T u (1) , as reported in [1].
If λ 1 is complex there is no fixed point but rather a limit-cycle solution to the dynamics of the complex-valued c 1 exists with δh (t) = Re c 1 (t) u (1) , and c i = 0 for all other eigenmodes. Assumingh (t) is periodic with period T , we can separate variables and integrate Eqn 30 in order to find c 1 (t) is given by for t ≤ T , where Φ (t) =´t 0 ds φ h (s) . Writing c 1 (0) = c 0 1 exp (iθ 0 ) and using Reλ 1 ≈ g, this gives c 1 (t) = c 0 1 exp (−t + gΦ (t)) exp (i (θ 0 + Imλ 1 Φ (t))) A limit cycle in phase withh (t) means c 1 (T ) = c 1 (0) and this requires that both gΦ (T ) = T and also Im [λ 1 ] Φ (T ) = 2π. From the first requirement we find that the average value of φ over a period must be the critical value: And combining the second requirement yields an expression for the period (Eqn ??, as reported in [1] as well): We can further write a self-consistency expression forh (t) by taking c 1 (t) as given by Eqn 35 and integrate over the coherent mode dynamics: which yieldsh Without loss of generality we assume thath (0) =h c = φ −1 1 g , then This is analogous to the fixed point equation forh * and c * 1 . In both cases the requirement that δh i = c 1 u (1) i 1 requires that c 1 be maximally O (1) and motivates our conjectures about the realization-dependence and system-size scaling of the transition out of chaos. In particular, we expect and confirm numerically that the critical value of J 1 for transition to either fixed point or limit cycle is inversely proportional to ν T u (1) and grows with network size (see main text and Fig 7).
In the case of complex leading eigenvalue, simulations confirm that a projection of the full synaptic current dynamics into the coherent mode and the leading eigenvector plane (consisting of real and imaginary parts of u (1) ) accounts for well over 0.99 of the total variance. Even restricting ourselves to the variance of the residuals, δh i , we find that 0.98 of the variance is restricted to the leading eigenvector plane (Fig S3).
For N = 4000 we simulate 219 realizations of random connectivity with complex leading eigenvalue and find that for sufficiently large J 1 all but one of these realizations yield highly oscillatory dynamics with period predicted nearly perfectly by theory (Fig S3).
We note that in the limit of large N we expect that the typical size of the imaginary component of the leading eigenvalue, λ 1 , shrinks such that the typical period grows. These longer period oscillations are characterized by squarewave-like shape in which the dynamics of the coherent component slows around the critical valueh c = φ −1 1 Reλ1 , which is identical to the fixed-point value ofh when λ 1 is real. (Fig S3) The fraction of realizations with real leading eigenvalue in the large N limit has not been calculated analytically to our knowledge. We find numerically that this fraction appears to saturate roughly around 3 10 for N 8000.