We use stochastic neural field theory to analyze the stimulus-dependent tuning of neural variability in ring attractor networks. We apply perturbation methods to show how the neural field equations can be reduced to a pair of stochastic nonlinear phase equations describing the stochastic wandering of spontaneously formed tuning curves or bump solutions. These equations are analyzed using a modified version of the bivariate von Mises distribution, which is well-known in the theory of circular statistics. We first consider a single ring network and derive a simple mathematical expression that accounts for the experimentally observed bimodal (or M-shaped) tuning of neural variability. We then explore the effects of inter-network coupling on stimulus-dependent variability in a pair of ring networks. These could represent populations of cells in two different layers of a cortical hypercolumn linked via vertical synaptic connections, or two different cortical hypercolumns linked by horizontal patchy connections within the same layer. We find that neural variability can be suppressed or facilitated, depending on whether the inter-network coupling is excitatory or inhibitory, and on the relative strengths and biases of the external stimuli to the two networks. These results are consistent with the general observation that increasing the mean firing rate via external stimuli or modulating drives tends to reduce neural variability.
Consider a pair of mutually coupled ring networks labeled j = 1,2. Let uj (θ, t) denote the activity at time t of a local population of cells with stimulus preference θ ∈ [−π, π) in network j. Here θ could represent the direction preference of neurons in area-middle temporal cortex (MT) [10], the orientation preference of V1 neurons, after rescaling θ → θ/2 [26,27], or a coordinate in spatial working memory [22,28,29]. For concreteness, we will refer to θ as a direction preference. The variables uj evolve according to the neural field equations [21,22,30,31] τ d u 1 (θ, t) = [ - u 1 (θ, t) + ∫ - π π J 1 (θ - θ ′) f (u 1 (θ ′, t) ) d θ ′ + ϵ ∫ - π π K 1 (θ - θ ′) f (u 2 (θ ′, t) ) d θ ′ + ϵ h 1 (θ) ] d t + 2 ϵ d W 1 (θ, t) (1a) τ d u 2 (θ, t) = [ - u 2 (θ, t) + ∫ - π π J 2 (θ - θ ′) f (u 2 (θ ′, t) ) d θ ′ + ϵ ∫ - π π K 2 (θ - θ ′) f (u 1 (θ ′, t) ) d θ ′ + ϵ h 2 (θ) ] d t + 2 ϵ d W 2 (θ, t) (1b) where ϵ is a constant scale factor (see below), Jj (θ − θ′) is the distribution of intra-network connections between cells with stimulus preferences θ′ and θ in network j, Kj (θ − θ′) is the corresponding distribution of inter-network connections to network j, and hj (θ) is an external stimulus. The firing rate function is assumed to be a sigmoid f (u) = f 0 1 + e - γ (u - η), (2) with maximal firing rate f0, gain γ and threshold η. The final term on the right-hand side of each equation represents external additive noise, with Wj (θ, t) a θ-dependent Wiener process. In particular, E [ d W j (θ, t) ] = 0, E [ d W i (θ, t) d W j (θ ′, s) ] = δ i, j C j (θ - θ ′) δ (t - t ′) d t d t ′, (3) where δ (t) is the Dirac delta function and δi, j is the Kronecker delta function. For concreteness, we will take C (θ) = aδ (θ) + b cos (θ) for constants a, b. For b ≠ 0, the noise is colored in θ (which is necessary for the solution to be spatially continuous) and white in time. (One could also take the noise to be colored in time by introducing an additional Ornstein-Uhlenbeck process. For simplicity, we assume that the noise processes in the two networks are uncorrelated, which would be the case if the noise were predominantly intrinsic. Correlations would arise if some of the noise arose from shared fluctuating inputs. For a discussion of the effects of correlated noise in coupled ring networks see [22].) The external stimuli are taken to be weakly biased inputs of the form ϵ h j with h j = h ¯ j cos (θ - θ ¯ j), (4) where θ ¯ j is the location of the peak of the input (stimulus bias) and h ¯ j is the contrast. Finally, the time-scale is fixed by setting the time constant τ = 10 msec. The maximal firing rate f0 is taken to be 100 spikes/sec. The weight distributions are 2π-periodic and even functions of θ and thus have cosine series expansions. Following [21], we take the intra-network recurrent connections to be J j (θ - θ ′) = J ¯ j cos (θ - θ ′), (5) which means that cells with similar stimulus preferences excite each other, whereas those with sufficiently different stimulus preferences inhibit each other. It remains to specify the nature of the inter-network connections. As we have already mentioned, we consider two different network configurations (see Fig 1): (A) a vertically connected two layer or laminar model and (B) a horizontally connected single layer model. In model A, the inter-network weight distribution is taken to have the form K j (θ - θ ′) = E j + K ¯ j cos (θ - θ ′), (6) which represents vertical coupling between the layers. We also assume that both layers receive the same stimulus bias, that is, θ ¯ 1 = θ ¯ 2 = θ ¯ in Eq (4). In model B, the inter-network weight distribution represents patchy horizontal connections, which tend to link cells with similar stimulus preferences [32–35]. This is implemented by taking K j (θ - θ ′) = K ¯ j δ (θ - θ ′). (7) Now the two networks can be driven by stimuli with different biases so that θ ¯ 1 ≠ θ ¯ 2. Note that in order to develop the analytical methods of this paper, we scale the internetwork coupling, the noise terms and the external stimuli in Eq (1) by the constant factor ϵ. Taking 0 < ϵ ≪ 1 (weak noise, weak inputs and weak inter-network coupling) will allow us to use perturbation methods to derive explicit parameter-dependent expressions for neural variability. We do not claim that cortical networks necessarily operate in these regimes, but use the weakness assumption to obtain analytical insights and make predictions about the qualitative behavior of neural variability. In the case of weak inter-network connections, the validity of the assumption is likely to depend on the source of these connections. For example, in model B, they arise from patchy horizontal connections within superficial or deep layers of cortex, which are known to play a modulatory role [36]. On the other hand, vertical connections between layers are likely to be stronger than assumed in our modeling analysis, at least in the feedforward direction [37]. Finally, the weak stimulus assumption depends on a particular view of how cortical neurons are tuned to stimuli, which is based on the theory of ring attractor networks, see the Discussion.
Let us begin by considering stimulus-dependent neural variability in a single ring network evolving according to the stochastic neural field equation d u (θ, t) = [ - u (θ, t) + ∫ - π π J (θ - θ ′) f (u (θ ′, t) ) d θ ′ + ϵ h (θ) ] d t + 2 ϵ d W (θ, t), (12) where E [ d W (θ, t) ] = 0, E [ d W (θ, t) d W (θ ′, t ′) ] = C (θ - θ ′) δ (t - t ′) d t d t ′, (13) with J (θ) = J ¯ cos θ, h (θ) = h ¯ cos θ, C (θ) = a δ (θ) + b cos θ. A clear demonstration of the suppressive effects of an external stimulus can be seen from direct numerical simulations of Eq (12), see Fig 2. In the absence of an external stimulus, the center-of-mass (phase) of the bump diffuses on the ring, whereas it exhibits localized fluctuations when a weakly-biased stimulus is present. Clearly, the main source of neural variation is due to the wandering of the bump, which motivates the amplitude phase decomposition given by Eq (9). Applying the perturbation analysis of Materials and Methods yields a one-network version of the phase Eq (10), which takes the form d β (t) = - ϵ Λ sin β (t) d t + 2 ϵ D d w (t), (14) with Λ = h ¯ / A and D = C ¯ / 2 A 2, where A is the amplitude of the bump for ϵ = 0. Eq (14) is known as a von Mises process, which can be regarded as a circular analog of the Ornstein-Uhlenbeck process on a line, and generates distributions that frequently arise in circular or directional statistics [45]. The von Mises process has been used to model the trajectories of swimming organisms [46,47], oscillators in physics [48], bioinformatics [49], and the data fitting of neural population tuning curves [50]. (Nonlinear stochastic phase equations analogous to (14) also arise in models of ring attractor networks with synaptic heterogeneities, which have applications to spatial working memory [23,51,52]). Introduce the probability density p (β, t | β 0,0 ] d β = P [ β < β (t) < β + d β | β (0) = β 0 ]. This satisfies the forward Fokker-Planck equation (dropping the explicit dependence on initial conditions) ∂ p (β, t) ∂ t = ∂ ∂ β [ ϵ Λ sin (β) p (β, t) ] + ϵ D ∂ 2 p (β, t) ∂ β 2 (15) for β ∈ [−π, π] with periodic boundary conditions p (−π, t) = p (π, t). It is straightforward to show that the steady-state solution of Eq (15) is the von Mises distribution p (β) = M (β; 0, κ), κ = h ¯ ϵ A D, (16) with M (β; β *, κ) ≔ 1 2 π I 0 (κ) exp (κ cos (β - β *) ). (17) Here I0 (κ) is the modified Bessel function of the first kind and zeroth order (n = 0), where I n (κ) = 1 2 π ∫ - π π exp (κ cos θ) cos (n θ) d θ. Sample plots of the von Mises distribution are shown in Fig 3. One finds that M (β; β*, κ) → 1/2π as κ → 0; since κ ∼ h ¯ this implies that in the absence of an external stimulus one recovers the uniform distribution of pure Brownian motion on the circle. On the other hand, the von Mises distribution becomes sharply peaked as κ → ∞. More specifically, for large positive κ, M (β; β *, κ) ≈ 1 2 π σ 2 e - (β - β *) 2 / 2 σ 2, σ 2 = κ - 1. (18) We thus have an explicit example of the noise suppression of fluctuations by an external stimulus, since σ 2 ∝ 1 / h ¯. (We are assuming that the time for the distribution of the stochastic phase variable to reach steady-state is much shorter than the time for the amplitude-phase decomposition (9) to break down. This can be proven rigorously using variational methods for sufficiently small ϵ, since the time for a large transverse fluctuation becomes exponentially large [44]). Moments of the von Mises distribution are usually calculated in terms of the circular moments of the complex exponential x = eiβ = cos β + i sin β. The nth circular moment is defined according to μ n = ⟨ z n ⟩ κ, β * = ∫ - π π z n M (β; β *, κ) d β = I n (κ) I 0 (κ) e i n β *. (19) In particular, ⟨ cos β ⟩ κ, β * = I 1 (κ) I 0 (κ) cos β *, ⟨ sin β ⟩ κ, β * = I 1 (κ) I 0 (κ) sin β *. (20) We can use these moments to explore stimulus-dependent variability in terms of the stochastic wandering of the bump or tuning curve. That is, consider the leading order approximation u (θ, t) ≈ A cos (θ + β (t) ), with β (t) evolving according to the von Mises SDE (14). Trial-to-trial variability can be captured by averaging the solution with respect to the stationary von Mises density (16). First, ⟨ U ⟩ (θ) = A ∫ - π π cos (θ + β) M (β, 0, κ) d β = A [ ⟨ cos β ⟩ κ, 0 cos θ - ⟨ sin β ⟩ κ, 0 sin θ ] ≔ A (κ) cos θ, A (κ) = A I 1 (κ) I 0 (κ). (21) from Eq (20). Hence, the mean amplitude A (κ) is given by the first circular moment of the von Mises distribution, see inset of Fig 3. When κ = 0 (zero external stimulus), the amplitude vanishes due to the fact that the random position of the bump is uniformly distributed around the ring. As the stimulus contrast h ¯ increases the wandering of the bump is more restricted and A (κ) monotonically increases. Second, ⟨ U 2 ⟩ (θ) = A 2 ∫ - π π cos 2 (θ + β) M (β, 0, κ) d β = A 2 2 ∫ - π π (1 + cos (2 [ θ + β ]) ) M (β, 0, κ) d β = A 2 2 [ 1 + ⟨ cos 2 β ⟩ κ, 0 cos 2 θ - ⟨ sin 2 β ⟩ κ, 0 sin 2 θ ] = A 2 2 [ 1 + I 2 (κ) I 0 (κ) cos 2 θ ]. It follows that the variance is var (U) = A 2 2 [ 1 + I 2 (κ) I 0 (κ) cos 2 θ - 2 (I 1 (κ) I 0 (κ) cos θ) 2 ] = A 2 2 { 1 - (I 1 (κ) I 0 (κ) ) 2 - [ (I 1 (κ) I 0 (κ) ) 2 - I 2 (κ) I 0 (κ) ] cos 2 θ } (22) In Fig 4 (a), we show example plots of the normalized variance var (U) /A2 as a function of the parameter κ, which is a proxy for the input amplitude h ¯, since κ ∝ h ¯. It can be seen that our theoretical analysis reproduces the various trends observed in [10]: (i) a global suppression of neural variability that increases with the stimulus contrast; (ii) a directional tuning of the variability that is bimodal; (iii) a peak in the suppression of cells at the preferred directional selectivity. One difference between our theoretical results and those of [10] is that, in the latter case, the directional tuning of the variance is not purely sinusoidal. Part of this can be accounted for by noting that we consider the variance of the activity variable u rather than the firing rate f (u). Moreover, for analytical convenience, we take the synaptic weight functions etc. to be first-order harmonics. In Fig 4 (b) we show numerical plots of the variance in the firing rate, which exhibits the type of bimodal behavior found in [10] when the ring network operates in the marginal regime. We now turn to a pair of coupled ring networks that represent vertically connected layers as shown in Fig 1 (a) (model A), with inter-network weight distribution (6). For analytical tractability, we impose the symmetry conditions A1 = A2 = A and K ¯ 1 = K ¯ 2 = K ¯. However, we allow the contrasts of the external stimuli to differ, h ¯ 1 ≠ h ¯ 2. Also, without loss of generality, we set θ ¯ 1 = θ ¯ 2 = 0. Eq (10) then reduce to the form, see Materials and methods d β 1 = - ϵ Λ 1 sin (β 1) d t - ϵ K ¯ sin (β 1 - β 2) d t + 2 ϵ d w 1 (t), (23a) d β 2 = - ϵ Λ 2 sin (β 2) d t - ϵ K ¯ sin (β 2 - β 1) d t + 2 ϵ d w 2 (t), (23b) with E [ d w j (t) ] = 0, E [ d w j (t) d w k (t ′) ] = δ j k D j δ (t - t ′) d t ′ d t. (24) Given our various simplifications, we can rewrite Eq (23) in the more compact form d β j = - ϵ ∂ Φ (β 1, β 2) ∂ β j d t + 2 ϵ d w j (t), j = 1,2 (25) where Φ is the potential function Φ (β 1, β 2) = - Λ 1 cos (β 1) - Λ 2 cos (β 2) - K ¯ cos (β 1 - β 2). (26) Introduce the joint probability density p (β 1, β 2, t | β 1,0, β 2,0, 0 ] d β 1 d β 2 = P [ β j < β j (t) < β j + d β j, j = 1,2 | β j (0) = β j, 0, j = 1,2 ]. This satisfies the two-dimensional forward Fokker-Planck equation (dropping the explicit dependence on initial conditions) ∂ p (β 1, β 2, t) ∂ t = ϵ ∑ j = 1,2 ∂ ∂ β j [ ∂ Φ (β 1, β 2) ∂ β j p (β 1, β 2, t) ] + ϵ ∑ j = 1,2 D j ∂ 2 p (β 1, β 2, t) ∂ β j 2 (27) for βj ∈ [−π, π] and periodic boundary conditions p (−π, β2, t) = p (π, β2, t), p (β1, −π, t) = p (β1, π, t). The existence of a potential function means that we can solve the time-independent FP equation. Setting time derivatives to zero, we have ∑ j = 1,2 ∂ J j ∂ β j = 0, J j = ϵ ∂ Φ ∂ β j p + ϵ D j ∂ p ∂ β j, where Jj is a probability current. In the stationary state the probability currents are constant, but generally non-zero. However, in the special case D1 = D2 = D, then there exists a steady-state solution in which the currents vanish. This can be seen by rewriting the vanishing current conditions as J j = ϵ p ∂ ∂ β j [ Φ + ϵ D ln p ] = 0. This yields the steady-state probability density, which is a generalization of the von Mises distribution, p (β 1, β 2) = N - 1 e - Φ (β 1, β 2) / ϵ D = N - 1 exp (κ 1 cos (β 1) + κ 2 cos (β 2) + χ cos (β 1 - β 2) ) ≔ M 2 (β 1, β 2; κ 1, κ 2, χ), (28) where κ j = h ¯ j ϵ A j D ≥ 0, χ = K ¯ ϵ D, and N is the normalization factor N (κ 1, κ 2, χ) = ∫ - π π ∫ - π π exp (κ 1 cos (β 1) + κ 2 cos (β 2) + χ cos (β 1 - β 2) ) d β 1 d β 2. (29) The distribution M2 (β1, β2; κ1, κ2, χ) is an example of a bivariate von Mises distribution known as the cosine model [49]. The normalization factor can be calculated explicitly to give N (κ 1, κ 2, χ) = (2 π) 2 [ I 0 (κ 1) I 0 (κ 2) I 0 (χ) + 2 ∑ s = 1 ∞ (- 1) s I s (κ 1) I s (κ 2) I s (χ) ]. (30) The corresponding marginal distribution for β1 is p (β 1) = ∫ - π π p (β 1, β 2) d β 2 = N (κ 1, κ 2, χ) - 1 2 π I 0 (κ 13 (β 1) ) exp (κ 2 cos (β 1) ), (31) where κ 13 (β) 2 = κ 1 2 + χ 2 + 2 κ 1 χ cos β. An analogous result holds for the marginal density p (β2). We now summarize a few important properties of the cosine bivariate von Mises distribution [49]: We will assume that the vertical connections are maximal between neurons with the same stimulus preference so that K ¯ ≥ 0 and χ ≥ 0. It then follows that p (β1, β2) is unimodal. Moreover, from Eq (32) we have Σ = 1 κ 1 κ 2 + χ (κ 1 + κ 2) (κ 2 + χ χ χ κ 1 + χ). (33) For zero inter-network coupling (χ = 0), we obtain the diagonal matrix Σ = diag (κ 1 - 1, κ 2 - 1) and we recover the variance of the single ring networks, that is, var (β j) = κ j - 1; there are no interlaminar correlations. On the other hand, for χ > 0 we find two major effects of the interlaminar connections. First, the vertical coupling reduces fluctuations in the phase variables within a layer. This is most easily seen by considering the symmetric case κ1 = κ2 = κ for which Σ = 1 κ (κ + 2 χ) (κ + χ χ χ κ + χ). (34) Clearly, var (β j) = 1 κ κ + χ κ + 2 χ < κ - 1. (35) (This result is consistent with a previous study of the effects of inter-network connections on neural variability, which focused on the case of zero stimuli and treated the bump positions as effectively evolving on the real line rather than a circle [22]. In this case, inter-network connections can reduce the variance in bump position, which evolves linearly with respect to the time t.) The second consequence of interlaminar connections is that they induce correlations between the phase β1 (t) and β2 (t). Having characterized the fluctuations in the phases β1 (t) and β2 (t), analogous statistical trends will apply to the trial-to-trial variability in the tuning curves. This follows from making the leading-order approximation uj (x, t) ∼ A cos (θ + βj (t) ), and then averaging the βj with respect to the bivariate von Mises density M2 (β1, β2; κ1, κ2, χ). In the large κj regime, this could be further simplified by averaging with respect to the bivariate normal distribution under the approximations cos (β) ≈ 1 − β2/2 and sin β ∼ β. Both the mean and variance of the tuning curves are similar to the single ring network, see Eqs (21) and (22): ⟨ U 1 ⟩ (θ) = A ⟨ cos β 1 ⟩ cos θ, (36) and var (U 1) = A 2 2 [ 1 - ⟨ cos β 1 ⟩ 2 - [ ⟨ cos β 1 ⟩ 2 - ⟨ cos 2 β 1 ⟩ ] cos 2 θ ] (37) Their dependence on the coupling strength χ and input parameter κ1 = κ2 = κ is illustrated in Fig 5. Finally, ⟨ U 1 (θ) U 2 (θ ′) ⟩ = A 2 ∫ - π π ∫ - π π cos (θ + β 1) cos (θ ′ + β 2) M 2 (β 1, β 2; κ, κ, χ) d β 1 d β 2 = A 2 ∫ - π π ∫ - π π (cos θ cos β 1 - sin θ sin β 1) × (cos θ ′ cos β 2 - sin θ ′ sin β 2) M 2 (β 1, β 2; κ, κ, χ) d β 1 d β 2 = A 2 (cos θ cos θ ′ ⟨ cos β 1 cos β 2 ⟩ + sin θ sin θ ′ ⟨ sin β 1 sin β 2 ⟩) = A 2 (cos θ cos θ ′ ⟨ cos β 1 cos β 2 ⟩ + sin θ sin θ ′ ⟨ sin β 1 sin β 2 ⟩). so that inter-network covariance take the form ⟨ U 1 (θ) U 2 (θ ′) ⟩ - ⟨ U 1 (θ) ⟩ ⟨ U 2 (θ ′) ⟩ = A 2 cos θ cos θ ′ [ ⟨ cos β 1 cos β 2 ⟩ - ⟨ cos β 1 ⟩ ⟨ cos β 2 ⟩ ] + A 2 sin θ sin θ ′ ⟨ sin β 1 sin β 2 ⟩. In particular, for θ = θ′ we have ⟨ U 1 (θ) U 2 (θ) ⟩ - ⟨ U 1 (θ) ⟩ ⟨ U 2 (θ) ⟩ = A 2 2 [ ⟨ sin β 1 sin β 2 ⟩ + ⟨ cos β 1 cos β 2 ⟩ - ⟨ cos β 1 ⟩ ⟨ cos β 2 ⟩ ] - A 2 2 [ ⟨ sin β 1 sin β 2 ⟩ - ⟨ (cos β 1 cos β 2 ⟩ - ⟨ cos β 1 ⟩ ⟨ cos β 2 ⟩) ] cos (2 θ). (38) The resulting correlation tuning curve behaves in a similar fashion to the variance, see Fig 5 (c), where corr (U 1, U 2) = ⟨ U 1 (θ) U 2 (θ) ⟩ - ⟨ U 1 (θ) ⟩ ⟨ U 2 (θ) ⟩ var (U 1) var (U 2). (39) (Note that our definition of the cross-correlation function differs from that used, for example, by Churchland et al [9]. These authors consider the covariance matrix of simultaneous recordings of spike counts obtained using a 96-electrode array. The matrix is then decomposed into a network covariance matrix and a diagonal matrix of private single neuron noise. Our definition involves pairwise correlations between the activity of two distinct populations. Nevertheless, consistent with the findings of Churchland et al. [9], we find that the cross-correlations decrease in the presence of a stimulus). The above qualitative analysis can be confirmed by numerical simulations of the full neural field Eq (1), as illustrated in Fig 6 (a) –6 (d) for a pair of identical ring networks. In Fig 6 (e) –6 (h), we show corresponding results for the case where network 2 receives a weaker stimulus than network 1 (h ¯ 1 = 2 and h ¯ 2 = 0. 5). In the absence of interlaminar connections, the phase of network 2 fluctuates much more than the phase of network 1. When interlaminar connections are included, fluctuations are reduced, but network 2 still exhibits greater variability than network 1. This latter result is consistent with an experimental study of neural variability in V1 [24], which found that neural correlations were more prominent in superficial and deep layers of cortex, but close to zero in input layer 4. One suggested explanation for these differences is that layer 4 receives direct feedforward input from the LGN. Thus we could interpret network 1 in model A as being located in layer 4, whereas network 2 is located in a superficial layer, say. Our final example concerns a pair of coupled ring networks that represent horizontally connected hypercolumns within the same superficial layer, say, as shown in Fig 1 (b) (model B), with inter-network weight distribution (7). Again, for analytical tractability, we impose the symmetry conditions A1 = A2 = A and K ¯ 1 = K ¯ 2 = K ¯. However, unlike model A, we take the contrasts to be the same, h ¯ 1 = h 2 ¯ = h ¯, but allow the biases of the two inputs to differ, θ ¯ 1 ≠ θ ¯ 2. Eq (10) become, see Materials and methods d β 1 = - ϵ Λ sin (β 1 + θ ¯ 1) d t + ϵ K (β 1 - β 2) d t + 2 ϵ d w 1 (t), (40a) d β 2 = - ϵ Λ sin (β 2 + θ ¯ 2) d t + ϵ K (β 2 - β 1) d t + 2 ϵ d w 2 (t), (40b) with wj (t) given by Eq (24) and K (β) = 2 K ¯ Γ ∫ - π π f ′ (A cos (θ - β) ) sin (θ - β) f (A cos (θ) ) d θ. (41) We can rewrite K (β) in the form K (β) = - 2 K ¯ A | Γ | ∂ ϕ ∂ β, ϕ (β) = f (A cos (θ - β) ) f (A cos (θ) ). (42) Note that ϕ (−β) = ϕ (β) and thus ϕ′ (−β) = −ϕ′ (β). A sample plot of the potential ϕ (β) is shown in Fig 7 (a), together with an approximate curve fitting based on a von Mises distribution. For the given firing rate parameters η = 0. 5 and γ = 4, the unperturbed bump amplitude is A ≈ 1. 85. As in the case of model A, we can rewrite Eq (40) in the more compact form d β j = - ϵ ∂ Ψ (β 1, β 2) ∂ β j d t + 2 ϵ d w j (t), j = 1,2 (43) where Ψ is the potential function Ψ (β 1, β 2) = - Λ cos (β 1 + θ ¯ 1) - Λ cos (β 2 + θ ¯ 2) - K ¯ ϕ (β 1 - β 2), (44) and we have absorbed the factor 2/ (A|Γ|) into the constant K ¯. The corresponding two-dimensional forward Fokker-Planck equation is ∂ p (β 1, β 2, t) ∂ t = ϵ ∑ j = 1,2 ∂ ∂ β j [ ∂ Ψ (β 1, β 2) ∂ β j p (β 1, β 2, t) ] + ϵ ∑ j = 1,2 D j ∂ 2 p (β 1, β 2, t) ∂ β j 2 (45) for βj ∈ [−π, π] and periodic boundary conditions p (−π, β2, t) = p (π, β2, t), p (β1, −π, t) = p (β1, π, t). Following the analysis of model A, if D1 = D2 = D then the stationary density takes the form p (β 1, β 2) = M - 1 e - Ψ (β 1, β 2) / ϵ D = M - 1 exp (κ cos (β 1 + θ ¯ 1) + κ cos (β 2 + θ ¯ 2) + χ ϕ (β 1 - β 2) ), (46) where κ = h ¯ ϵ A D ≥ 0, χ = K ¯ ϵ D, and M is a normalization factor. Long-range horizontal connections within superficial layers of cortex are mediated by the axons of excitatory pyramidal neurons. However, they innervate both pyramidal neurons and feedforward interneurons so that they can have a net excitatory or inhibitory effect, depending on stimulus conditions [36,53,54], More specifically, they tend to be excitatory at low contrasts and inhibitory at high contrasts. Suppose that ring network 1 represents a hypercolumn driven by a stimulus h ¯ cos θ and network 2 represents a hypercolumn driven by a stimulus h ¯ cos (θ - θ ¯), see Fig 1 (b). In Fig 7 (b) and 7 (c) we plot how the normalized maximal mean and variance of network 1 (at θ = ±π/2) varies with the directional bias θ ¯ of the input to network 2. We also show the baseline mean and variance in the absence of horizontal connections (χ = 0). It can be seen that the mean and variance covary in opposite directions. In particular, for inhibitory horizontal connections (χ < 0) the variance is facilitated relative to baseline when the two stimuli have similar biases (θ ¯ ≈ 0) and is suppressed when they are sufficiently different (θ ¯ ≈ ± π). The converse holds for excitatory horizontal connections (χ > 0). In the Discussion, these results will be explored within the context of surround modulation.
In order to utilize perturbation methods, we assumed that the ring networks were driven by weakly biased stimuli. This assumption depends on a particular view of how cortical neurons are tuned to stimuli. Consider the most studied example, which involves orientation tuning of cells in V1. The degree to which recurrent processes contribute to the receptive field properties of V1 neurons has been quite controversial over the years [55–58]. The classical model of Hubel and Wiesel [59] proposed that the orientation preference and selectivity of a cortical neuron in input layer 4 arises primarily from the geometric alignment of the receptive fields of thalamic neurons in the lateral geniculate nucleus (LGN) projecting to it. (Orientation selectivity is then carried to other cortical layers through vertical projections). This has been confirmed by a number of experiments [60–64]. However, there is also significant experimental evidence suggesting the importance of recurrent cortical interactions in orientation tuning [65–71]. One issue that is not disputed is that some form of inhibition is required to explain features such as contrast-invariant tuning curves and cross-orientation suppression [58]. The uncertainty in the degree to which intracortical connections contribute to orientation tuning of V1 neurons is also reflected in the variety of models. In ring attractor models [26,27,72,73], the width of orientation tuning of V1 cells is determined by the lateral extent of intracortical connections. Recurrent excitatory connections amplify weakly biased feedforward inputs in a way that is sculpted by lateral inhibitory connections. Hence, the tuning width and other aspects of cortical responses are primarily determined by intracortical rather than thalamocortical interconnections. On the other hand, in push-pull models, cross-orientation inhibition arises from feedforward inhibition from interneurons [62,74]. Finally, in normalization models, a large pool of orientation-selective cortical interneurons generates shunting inhibition proportional in strength to the stimulus contrast at all orientations [75]. In the end, it is quite possible that are multiple circuit mechanisms for generating tuned cortical responses to stimuli, which could depend on the particular stimulus feature, location within a feature preference map, and cortical layer [58]. Surround modulation (SM) refers to the phenomenon in which stimuli in the surround of a neuron’s receptive field (RF) modulate the neuron’s response to stimuli simultaneously presented inside the RF. SM is a fundamental property of sensory neurons in many species and sensory modalities, and is thought to play an important role in contextual image processing. As with mechanisms of orientation tuning, there is considerable debate over whether feedforward or intracortical circuits generate SM, and whether this results from increased inhibition or reduced excitation [19,36,53,54,76–82]. SM has been characterized in many species, commonly using circular grating patches of increasing radius or grating patches confined to the RF surrounded by annular gratings, and varying systematically the grating parameters. Modulatory effects are typically quantified in terms of changes in the mean firing rates of single neurons recorded from the center. Some of the main features of SM in V1 are as follows (see [36] and references therein): (i) SM is spatially extensive. For example, in primates, modulatory effects from the surround (both facilitatory and suppressive) can be evoked at least 12. 5 degrees away from a neuron’s RF center. (ii) SM is tuned to specific stimulus parameters. The strongest suppression is induced by stimuli in the RF and surround of the same orientation, spatial frequency, drift direction, and speed, and weaker suppression or facilitation is induced by stimuli of orthogonal parameters (e. g. , orthogonally oriented stimuli or stimuli drifting in opposite directions). (iii) SM is contrast dependent. Surround stimulation evokes suppression when the center and surround stimuli are of high contrast, but can be facilitatory when they are of low contrast. One way to interpret the results of model B is to treat networks 1 and 2 as hypercolumns driven by center and surround stimuli, respectively. SM is then mediated by the horizontal connections that can have a net excitatory or inhibitory effect, depending on stimulus conditions. Here, for simplicity, we impose the sign of the horizontal connections by hand. However, one could develop a more detailed model that implements the switch between excitation and inhibition using, for example, high threshold interneurons [54]. The major prediction of our analysis is that whenever the surround modulation suppresses (facilitates) the center firing rate, the corresponding variance is facilitated (suppressed). One of the main simplifications of our neural field model is that we do not explicitly distinguish between excitatory and inhibitory populations. This is a common approach to the analysis of neural fields, in which the combined effects of excitation and inhibition are incorporated using, for example, Mexican hat functions [83–85]. In the case of the ring network, the spontaneous formation of population orientation tuning curves or bumps is implemented using a cosine function, which represents short-range excitation and longer-range inhibition around the ring. We note, however, that the methods and results presented in this paper could be extended to the case of separate excitatory and inhibitory populations, as well as different classes of interneuron, as has been demonstrated elsewhere for deterministic neural fields [27,54]. One major difference between scalar and E-I neural fields is that the latter can also exhibit time-periodic solutions, which would add an additional phase variable associated with shifts around the resulting limit cycle. The effects of noise on limit cycle oscillators can be analyzed in an analogous fashion to wandering bumps [86,87]. We also note that neural variability in a two-population (E-I) stabilized supralinear network has been analyzed extensively using linear algebra [20]. Another possible extension of our work would be to consider higher-dimensional neural fields. For example, one could replace the ring attractor on S1 by a spherical attractor on S2. In the latter case, marginally stable modes would correspond to rotations of the sphere. (Mathematically speaking, this corresponds to the action of the Lie group SO (3) rather than SO (2) for the circle.) One could generalize the Fourier analysis of the ring network by using spherical harmonics, as previously shown for deterministic neural field models of orientation and spatial frequency tuning in V1 [88,89]. One could also consider a planar neural field with Euclidean-symmetric weights, for which marginally stable modes would be generated by the Euclidean group of rigid body transformations of the plane (translations, rotations and reflections.) However, this example is more difficult since the marginally stable manifold is non-compact, and one cannot carry out a low-dimensional harmonic reduction. In order to obtain analytical results, one has to use Heaviside rate functions [30,90]. A third possible extension would be to develop a more detailed model of the laminar structure of cortex. Roughly speaking, cortical layers can be grouped into input layer 4, superficial layers 2/3 and deep layers 5/6 [37,91–93]. They can be distinguished by the source of afferents into the layer and the targets of efferents leaving the layer, the nature and extent of intralaminar connections, the identity of interneurons within and between layers, and the degree of stimulus specificity of pyramidal cells. In previous work, we explored the role of cortical layers in the propagation of waves of orientation selectivity across V1 [94], under the assumption that deep layers are less tuned to orientation. This suggests considering coupled ring networks that differ in their tuning properties. Another modification would be to consider asymmetric coupling between layers, both in terms of the range of coupling and its strength. Interestingly, the properties of SM also differ across cortical layers, suggesting different circuits and mechanisms generating SM in different layers. More specifically, surround fields in input layer 4 are smaller than in other layers, and SM is weaker and untuned for orientation. Moreover, SM is stronger and more sharply orientation-tuned in superficial layers compared to deep layers [36]. Therefore, it would be interesting to consider coupled ring networks that combine models A and B. One final comment is in order. Neural variability in experiments is typically specified in terms of the statistics of spike counts over some fixed time interval, and compared to an underlying inhomogeneous Poisson process. Often Fano factors greater than one are observed. In this paper, we worked with stochastic firing rate models rather than spiking models, so that there is some implicit population averaging involved. In particular, we focused on the statistics of the variables uj (x, t), which represent the activity of local populations of cells rather than of individual neurons, with f (uj) the corresponding population firing rate [30]. This allowed us to develop an analytically tractable framework for investigating how neural variability depends on stimulus conditions within the attractor model paradigm. In order to fit a neural field model to single-neuron data, one could generate spike statistics by taking f (uj) to be the rate of an inhomogeneous Poisson process. Since f (uj) is itself stochastic, this would result in a doubly stochastic Poisson process, which is known to produce Fano factors greater than unity [95]. Moreover, the various phenomena identified in this paper regarding stimulus-dependent variability would carry over to a spiking model, at least qualitatively. However, one should not expect a mean-field reduction to capture everything in a spiking model. For example, multivariate doubly stochastic Poisson processes can have correlations between their spike times in addition to the correlations induced by shared rate fluctuations. Spiking network models typically do produce these spike timing correlations that are not captured by most mean-field reductions, even those that account for correlated firing rate fluctuations [13,15,96–98]. These correlations could, in turn, affect auto-correlation and firing rates in the network.
First, suppose that there are no external inputs, no inter-network coupling (J12 = J21 = 0), and no noise (ϵ = 0). Each network can then be described by a homogeneous ring model of the form ∂ u (θ, t) ∂ t = - u (θ, t) + ∫ - π π J (θ - θ ′) f (u (θ ′, t) ) d θ ′. (47) Let J (θ) = J ¯ cos θ and consider the trial solution u (θ, t) = U (θ) with U (θ) an even, unimodal function of θ centered about θ = 0. This could represent a direction tuning curve in MT ( (in the marginal regime) or a stationary bump encoding a spatial working memory. It follows that U (θ) satisfies the integral equation U (θ) = J ¯ ∫ - π π cos (θ - θ ′) f (U (θ ′) ) d θ ′. (48) Substituting the cosine series expansion cos (θ - θ ′) = cos (θ) cos (θ ′) + sin (θ) sin (θ ′) (49) into the integral equation yields the even solution U (θ) = A cos θ with the amplitude A satisfying the self-consistency condition A = J ¯ ∫ - π π cos (θ) f (U (θ) ) d θ = J ¯ g (A). (50) The amplitude Eq (50) can be solved explicitly in the large gain limit γ → ∞, for which f (u) → H (u − κ), where H is the Heaviside function [21]. That is, A = 1 + κ ± 1 - κ, corresponding to a marginally stable large amplitude wide bump and an unstable small amplitude narrow bump, consistent with the original analysis of Amari [90]. On the other hand, at intermediate gains, there exists a single stable bump rather than an unstable/stable pair of bumps, see Fig 8. Linear stability of the stationary solution can be determined by considering weakly perturbed solutions of the form u (θ, t) = U (θ) + ψ (θ) eλt for |ψ (θ) | ≪ 1. Substituting this expression into Eq (47), Taylor expanding to first order in ψ, and imposing the stationary condition (48) yields the infinite-dimensional eigenvalue problem [27] (λ + 1) ψ (θ) = ∫ - π π J (θ - θ ′) f ′ (U (θ ′) ) ψ (θ ′) d θ ′. (51) This can be reduced to a finite-dimensional eigenvalue problem by applying the expansion (49): (λ + 1) ψ (θ) = A cos (θ) + B sin (θ), (52) where A = J ¯ ∫ - π π cos (θ) f ′ (U (θ) ) ψ (θ) d θ, B = J ¯ ∫ - π π sin (θ) f ′ (U (θ) ) ψ (θ) d θ. (53) Substituting Eqs (52) into (53) then gives the matrix equation [21] (λ + 1) (A B) = J ¯ (I [ cos 2 θ ] I [ cos θ sin θ ] I [ cos θ sin θ ] I [ sin 2 θ ]) (A B), (54) where for any periodic function v (θ) I [ v (θ) ] = ∫ - π π v (θ) f ′ (U (θ) ) d θ. (55) Integrating Eq (50) by parts shows that for A ≠ 0 I [ sin 2 θ ] = ∫ - π π sin 2 θ f ′ (U (θ) ) d θ = 1 / J ¯. Hence, exploiting the fact that I is a linear functional of v, I [ cos 2 θ ] = I [ 1 - sin 2 θ ] = I [ 1 ] - I [ sin 2 θ ] = I [ 1 ] - 1 / J ¯. Finally, integration by parts establishes that I [ cos θ sin θ ] = ∫ - π π cos θ sin θ f ′ (U (θ) ) d θ = - ∫ - π π sin θ f (U (θ) ) d θ = 0, since U (θ) is even. Eq (54) now reduces to (λ + 1) (A B) = J ¯ (I [ 1 ] - 1 / J ¯ 0 0 1) (A B), (56) which yields the pair of solutions λ 0 = 0, λ e = 2 [ J ¯ ∫ 0 π f ′ (U (θ) ) d θ - 1 ]. (57) The zero eigenvalue is a consequence of the fact that the bump solution is marginally stable with respect to uniform shifts around the ring; the generator of such shifts is the odd function sinθ. The other eigenvalue λe is associated with the generator, cosθ, of expanding or contracting perturbations of the bump. Thus linear stability of the bump reduces to the condition λe < 0. This can be used to determine the stability of the pair of bump solutions in the high-gain limit [21]. (Note that there also exist infinitely many eigenvalues that are equal to −1, which form the essential spectrum. However, since they lie in the left-half complex λ-plane, they do not affect stability). A variety of previous studies have shown how breaking the underlying translation invariance of a homogeneous neural field by introducing a nonzero external input stabilizes wave and bump solutions to translating perturbations [21,99–102]. For the sake of illustration, suppose that h (θ) = h ¯ cos (θ) in the deterministic version of Eq (1). This represents a weak θ-dependent input with a peak at θ = 0. Extending the previous analysis, one finds a stationary bump solution U (θ) = A cos θ + ϵ h ¯ cos θ, with A satisfying the implicit equation A = J ¯ ∫ - π π cos θ f (A cos θ + ϵ h ¯ cos θ) d θ. Again, this can be used to determine both the width and amplitude of the bump in the high-gain limit. Furthermore, the above analysis can be extended to establish that, for weak inputs, the bump is stable (rather than marginally stable) with respect to translational shifts [21]. The amplitude phase decompositions (βj, vj) defined by Eq (9) are not unique, so additional mathematical constraints are needed, and this requires specifying the allowed class of functions of vj (the appropriate Hilbert space). We will take take vj ∈ L2 (S1), that is, vj (θ) is a periodic function with ∥ v j∥ 2 = 〈 v j, v j 〉 = ∫ - π π v j (θ) 2 d θ < ∞. Substituting the decomposition into the stochastic neural field Eq (1) and using Ito’s lemma gives [103] U 1 ′ (θ + β 1) d β 1 + 1 2 U 1 ′ ′ (θ + β 1) d β 1 2 + ϵ d v 1 (θ, t) = [ - U 1 (θ + β 1) - ϵ v 1 (θ, t) + ∫ - π π J 1 (θ - θ ′) f (U 1 (θ ′ + β 1) + ϵ v 1 (θ ′, t) ) d θ ′ ] d t + [ ∫ - π π ϵ K 1 (θ - θ ′) f (U 2 (θ ′ + β 2) + ϵ v 2 (θ ′, t) ) d θ ′ + ϵ h 1 (θ) ] d t + 2 ϵ d W 1 (θ, t). and U 2 ′ (θ + β 2) d β 2 (t) + 1 2 U 2 ′ ′ (θ + β 2) d β 2 2 + ϵ d v 2 (θ, t) = [ - U 2 (θ + β 2) - ϵ v 2 (θ, t) + ∫ - π π J 2 (θ - θ ′) f (U 2 (θ ′ + β 2) + ϵ v 2 (θ ′, t) ) d θ ′ ] d t + [ ∫ - π π ϵ K 2 (θ - θ ′) f (U 1 (θ ′ + β 1) + ϵ v 1 (θ ′, t) ) d θ ′ + ϵ h 2 (θ) ] d t + 2 ϵ d W 2 (θ, t). Introduce the series expansions v j = v j, 0 + ϵ v j, 1 + O (ϵ), Taylor expanding the nonlinear function F, imposing the stationary solution (48), and dropping all O (ϵ) terms. This gives [21,30], after dropping the zero index on vj, 0, ϵ d v 1 (θ, t) = ϵ L β 1 1 v 1 (θ, t) d t + ϵ K ^ 1 (θ + β 2) d t + ϵ h 1 (θ) d t + 2 ϵ d W 1 (θ, t), - U 1 ′ (θ + β 1) d β 1 (58a) ϵ d v 2 (θ, t) = ϵ L β 2 2 v 2 (θ, t) d t + ϵ K ^ 2 (θ + β 1) + ϵ h 2 (θ) d t + 2 ϵ d W 2 (θ, t) - U 2 ′ (θ + β 2) d β 2, (58b) where L β j are the following linear operators L β j v (θ, t) = - v (θ, t) + ∫ - π π J j (θ - θ ′) f ′ (U j (θ ′ + β) ) v (θ ′, t) d θ ′, (59) and K ^ 1 (θ + β) = ∫ - π π K 1 (θ - θ ′) f (U 2 (θ ′ + β) ) d θ ′, K ^ 2 (θ + β) = ∫ - π π K 2 (θ - θ ′) f (U 1 (θ ′ + β) ) d θ ′. (60) It can be shown that the operator L 0 j has a 1D null space spanned by U j ′ (θ). The fact that U j ′ (θ) belongs to the null space follows immediately from differentiating Eq (48) with respect to θ. Moreover, U j ′ (θ) is the generator of uniform translations around the ring, so that the 1D null space reflects the marginal stability of the bump solution. (Marginal stability of the bump means that the linear operator L 0j has a simple zero eigenvalue while the remainder of the discrete spectrum lies in the left-half complex plane. The spectrum is discrete since S1 is a compact domain.) This then implies a pair of solvability conditions for the existence of bounded solutions of Eq (58a), namely, that dvj is orthogonal to all elements of the null space of the adjoint operator L β j j †. The corresponding adjoint operator is L β j † v (θ, t) = - v (θ, t) + f ′ (U j (θ + β) ) ∫ - π π J j (θ - θ ′) v (θ ′, t) d θ ′. (61) Let V j (θ) span the 1D adjoint null space of L 0 †. Now taking the inner product of both sides of Eq (58a) with respect to V j (θ + β j) and using translational invariance then yields the following SDEs to leading order: d β 1 = ϵ H 1 (β 1) d t - ϵ K 1 (β 1 - β 2) d t + 2 ϵ d w 1 (t), (62a) d β 2 = ϵ H 2 (β 1) d t - ϵ K 2 (β 2 - β 1) d t + 2 ϵ d w 2 (t), (62b) where H j (β) = Γ j - 1 ∫ - π π V j (θ) h j (θ - β) d θ, (63) for Hj (β + 2π) = Hj (β), K j (β) = Γ j - 1 ∫ - π π V j (θ) K ^ j (θ + β) d θ, (64) and Γ j = ∫ - π π V j (θ) U j ′ (θ) d θ, (65) Here wj (t) are scalar independent Wiener processes, E [ d w j (t) ] = 0, E [ d w j (t) d w k (t ′) ] = δ j, k D j δ (t - t ′) d t ′ d t, with D j = 1 Γ j 2 ∫ - π π ∫ - π π V j (θ) V j (θ ′) C j (θ - θ ′) d θ ′ d θ. (66) Note that stochastic phase equations similar to (62) were previously derived in [21,22], except that the functions Hj (β) and K j (β) were linearized, resulting in a system of coupled Ornstein-Uhlenbeck (OU) processes: d β 1 = ϵ ν 1 β 1 d t - ϵ r 1 (β 1 - β 2) d t + 2 ϵ d w 1 (t), (67a) d β 2 = ϵ ν 2 β 2 d t - ϵ r 2 (β 2 - β 1) d t + 2 ϵ d w 2 (t), (67b) for constant coefficients ν1, ν2, r1, r2. Properties of one-dimensional OU processes were then used to explore how the variance in the position of bump solutions depended on inter-network connections and statistical noise correlations. However, it should be noted that the variables βj (t) are phases on a circle (rather than positions on the real line), so that the right-hand side of Eq (67) should involve 2π-period functions. Therefore, the linear approximation only remains accurate on sufficiently short times scales for which the probability of either of the phases winding around the circle is negligible. In order to illustrate this point, consider an uncoupled OU process evolving according to d β j = ϵ ν j β j d t + 2 ϵ d w j (t). A standard analysis shows that [103] ⟨ β j (t) ⟩ = β 0 e - ν j t, ⟨ β j (t) 2 ⟩ - ⟨ β j (t) ⟩ 2 = ϵ D j ν j [ 1 - e - 2 ν j t ]. In particular, the variance approaches a constant ϵD/2νj in the large t limit. The corresponding density is given by the Gaussian ρ (β, t | β 0,0) = 1 2 π ϵ D j [ 1 - e - 2 ν j t ] / ν j exp (- (β - β 0 e - k t) 2 2 ϵ D j [ 1 - e - 2 ν j t ] / ν j). Although the linear approximation is sufficient if one is interested in estimating the diffusivity Dj, which was the focus of [21,22], it does not yield the correct steady-state distribution on the ring in the limit t → ∞. Indeed, for vj → 0, the density of the OU process converges point-wise to zero, whereas ρ (β, t) → 1/2π on the ring. In our paper, we are interested in the full steady-state densities rather than just the diffusivities Dj. In order to determine the functions Hj and K j we need to obtain explicit expressions for the null vectors V j. We will take h j (θ) = h ¯ j cos (θ - θ ¯ j). Applying the expansion (49) to the adjoint equation L 0 j † V j = 0 with L 0 j † defined by Eq (61), we can write [21] V j (θ) = f ′ (U j (θ) ) [ C j cos θ + S j sin θ ], with C j = J ¯ j ∫ - π π cos θ V j (θ) d θ, S j = J ¯ j ∫ - π π sin θ V j (θ) d θ. Substituting the expression for V j (θ) into the expressions for Cj and Sj then leads to a matrix equation of the form (56) with λ = 0. Since I [ 1 ] ≠ 1, it follows that Cj = 0 so that, up to scalar multiplications, V j (θ) = f ′ (U j (θ) ) sin θ, U (θ) = A j cos θ. (68) Now substituting V (θ) into Eq (63), we have H j (β) = Γ j - 1 ∫ - π π V j (θ) h j (θ - β) d θ = h ¯ j Γ ∫ - π π f ′ (U j (θ) ) sin θ cos (θ - θ ¯ j - β) d θ = h ¯ j Γ ∫ - π π f ′ (U j (θ) ) sin θ [ cos θ cos (β + θ ¯ j) + sin θ sin (β + θ ¯ j) ] d θ = - Λ j sin (β + θ ¯ j), (69) with Λ j = - h ¯ j Γ j ∫ - π π f ′ (U j (θ) ) sin 2 θ d θ. (70) We have used the fact that f″ (Uj (θ) ) is an even function of θ, so that the coefficient for cos (β + θ ¯ j) is zero. The constant Γj can be calculated from Eq (65): Γ j = ∫ - π π V j (θ) U j ′ (θ) d θ = - A j ∫ - π π f ′ (U j (θ) ) sin 2 θ d θ < 0. (71) It follows that Λ j = h ¯ j A j > 0. (72) The calculation of K j (β) depends on whether we consider model A or model B, see Fig 1. From Eqs (6), (60) and (64), we have for model A K 1 (β) = Γ 1 - 1 ∫ - π π V 1 (θ) [ ∫ - π π K 1 (θ - θ ′) f (U 2 (θ ′ + β) ) d θ ′ ] d θ, = Γ 1 - 1 ∫ - π π f ′ (U 1 (θ) ) sin θ [ ∫ - π π [ E 1 + K ¯ 1 cos (θ - θ ′) ] f (U 2 (θ ′ + β) ) d θ ′ ] d θ = Γ 1 - 1 [ ∫ - π π f ′ (U 1 (θ) ) sin 2 θ d θ ] [ ∫ - π π K ¯ 1 sin θ ′ f (U 2 (θ ′ + β) ) d θ ′ ] = - 1 A 1 [ ∫ - π π K ¯ 1 sin (θ ′ - β) f (U 2 (θ ′) ) d θ ′ ] = K ¯ 1 A 1 sin β [ ∫ - π π cos θ ′ f (A 2 cos (θ ′) ) d θ ′ ] ≡ K ¯ 1 A 2 A 1 sin β, (73a) where we have used the stationary condition (8), and K 2 (β) = K ¯ 2 A 2 sin β [ ∫ - π π cos θ ′ f (A 1 cos (θ ′) ) d θ ′ ] ≡ K ¯ 2 A 1 A 2 sin β. (73b) Similarly, from Eqs (7), (60) and (64), we have for model B K 1 (β) = Γ 1 - 1 ∫ - π π f ′ (U 1 (θ) ) sin θ [ ∫ - π π K ¯ 1 δ (θ - θ ′) f (U 2 (θ ′ + β) ) d θ ′ ] d θ = 2 K ¯ 1 Γ 1 ∫ - π π f ′ (U 1 (θ) ) sin θ f (U 2 (θ + β) ) d θ = 2 K ¯ 1 Γ 1 ∫ - π π f ′ (U 1 (θ - β) ) sin (θ - β) f (U 2 (θ) ) d θ. (74a) Similarly, K 2 (β) = 2 K ¯ 2 Γ 2 ∫ - π π f ′ (U 2 (θ - β) ) sin (θ - β) f (U 1 (θ) ) d θ. (74b) Finally, from Eq (66), the diffusion coefficients Dj become D j = 1 Γ j 2 ∫ - π π ∫ - π π C j (θ - θ ′) f ′ (U j (θ) ) f ′ (U j (θ ′) ) sin θ sin θ ′ d θ ′ d θ. (75) One finds that the diffusivities decreases as the spatial correlation lengths increase. For example, in the case of spatially homogeneous noise (C j (θ - θ ′) = C ¯ j), Dj = 0 since f′ (Uj (θ) ) is even. On the other hand, for spatially uncorrelated noise (C j (θ - θ ′) = C ¯ j δ (θ - θ ′) ), we have D j = C ¯ j Γ j 2 ∫ - π π sin 2 θ [ f ′ (U j (θ) ] 2 d θ > 0. (76) In Results we take C j (θ - θ ′) = C ¯ j cos (θ - θ ′) so that D j = 1 Γ j 2 ∫ - π π ∫ - π π C ¯ j cos (θ - θ ′) f ′ (A j cos (θ) ) f ′ (A j cos (θ ′) ) sin θ sin θ ′ d θ ′ d θ = C ¯ j Γ j 2 [ ∫ - π π f ′ (A j cos (θ) ) sin 2 θ d θ ] 2 = C ¯ j 2 A j 2. (77) All numerical simulations were performed in Matlab. One dimensional numerical simulations were performed using a forward Euler method scheme in time and a trapezoidal rule for integration in θ. Time steps were taken to be Δt = 0. 001, and orientation steps Δθ = 0. 01π.
A topic of considerable current interest concerns the neural mechanisms underlying the suppression of cortical variability following the onset of a stimulus. Since trial-by-trial variability and noise correlations are known to affect the information capacity of neurons, such suppression could improve the accuracy of population codes. One of the main candidate mechanisms is the suppression of noise-induced transitions between multiple attractors, as exemplified by ring attractor networks. The latter have been used to model experimentally measured stochastic tuning curves of directionally selective middle temporal (MT) neurons. In this paper we show how the stimulus-dependent tuning of neural variability in ring attractor networks can be analyzed in terms of the stochastic wandering of spontaneously formed tuning curves or bumps in a continuum neural field model. The advantage of neural fields is that one can derive explicit mathematical expressions for the second-order statistics of neural activity, and explore how this depends on important model parameters, such as the level of noise, the strength of recurrent connections, and the input contrast.
perturbation theory
neural networks
neuroscience
learning and memory
cognitive neuroscience
mathematics
cognition
network analysis
algebra
memory
interneurons
quantum mechanics
neuronal tuning
computer and information sciences
animal cells
physics
working memory
cellular neuroscience
cell biology
linear algebra
neurons
biology and life sciences
cellular types
physical sciences
cognitive science
eigenvalues
Stochastic neural field model of stimulus-dependent variability in cortical neurons