Adaptive scales of integration and response latencies in a critically-balanced model of the primary visual cortex.

The primary visual cortex (V1) integrates information over scales in visual space, which have been shown to vary, in an input-dependent manner, as a function of contrast and other visual parameters. Which algorithms the brain uses to achieve this feat are largely unknown and an open problem in visual neuroscience. We demonstrate that a simple dynamical mechanism can account for this contrast-dependent scale of integration in visuotopic space as well as connect this property to two other stimulus-dependent features of V1: extents of lateral integration on the cortical surface and response latencies.

Introduction

Stimuli in the natural world have quantitative characteristics that vary over staggering ranges. Our nervous system evolved to parse such widely-ranging stimuli, and research into how the nervous system can cope with such ranges has led to considerable advances in our understanding of neural circuitry. For example at the sensory transduction level, the physical magnitudes encoded into primary sensors, such as light intensity, sound pressure level and olfactant concentration, vary over exponentially-large ranges, leading to the Weber-Fechner law [1]. As neuronal firing rates cannot vary over such large ranges, the encoding process must compress physical stimuli into the far more limited ranges of neural activity that represent them. These observations have stimulated a large amount of research into the mechanisms underlying nonlinearly compression of physical stimuli in the nervous system. Of relevance to our later discussion is the nonlinear compression of sound intensity in the early auditory pathways [2-4], where it has been shown that poising the active cochlear elements on a Hopf bifurcation leads to cubic-root compression.

But other characteristics besides the raw physical magnitude still vary hugely. The wide range of spatial extents and correlated linear structures present in visual scenery [5-7] leads to a more subtle problem, if we think of the visual areas as fundamentally limited by corresponding anatomical connectivity. Research into this problem has been focused on elucidating the nature of receptive fields of neurons in the primary visual cortex (V1) [8-13]. Studies have found that as the contrast of a stimulus is decreased, the receptive field [14, 15] size or area of spatial summation in visual space increases (Fig 1) [12, 13, 16, 17]. As an example of contextual modulation of neuronal responses, this problem has naturally received theoretical attention [18-20]. However, current literature does not describe this phenomenon as structurally integral to the neural architecture but rather either highlight a different set of features or the contextual modulations are explicitly written in an ad hoc fashion. Our aim is to develop a model which displays this phenomenon structurally, as a direct consequence of the neural architecture. In our proposed models, multiple length scales emerge naturally without any fine tuning of the system’s parameters. This leads to length-tuning curves similar to the ones measured in Kapadia et al. over the entire range (Fig 1) [12].

Fig 1

Reprinted from [12] under a CC BY license, with permission from National Academy of Sciences, U.S.A., original copyright (1999). Measurements of single-neuron responses in the V1 area of monkeys to optimally oriented bars of light of different lengths and contrasts.

Panels a and b are measurements from two distinct neurons. The units of length along the horizontal axis are in minutes of arc. The solid, dotted, and dashed curves represent bars of light of 50% contrast, 15% contrast, and 50% contrast embedded in a textured background, respectively. The dashed curves are irrelevant to the focus of this paper.

The findings of Kapadia et al. demonstrate that receptive fields in V1 are not constant but instead grow and shrink, seemingly beyond naive anatomical parameters, according to stimulus contrast. The “computation” being carried out is not fixed but is itself a function of the input. Let us examine this distinction carefully. There are numerous operations in image processing, such as Gaussian blurs or other convolutional kernels, whose spatial range is fixed. It is very natural to imagine neural circuitry having actual physical connections corresponding to the nonzero elements of a convolutional kernel, and in fact a fair amount of effort has been expended trying to identify actual synapses corresponding to such elements [21, 22]. There are, however, other image-processing operations, such as floodfill (the “paint bucket”) whose spatial extent is entirely dependent on the input; the problem of “binding” of perceptual elements is usually thought about in this way, and mechanisms posited to underlie such propagation dynamics include synchronization of oscillations acting in a vaguely paint-bucket-like way [23-25]. This dichotomy is artificial because these are only the two extremes of a potentially continuous range. While the responses of neurons in V1 superficially appear to be convolutional kernels, their strong dependence on input characteristics, particularly the size of the receptive field, demonstrates a more complex logic in which spatial extent is determined by specific characteristics of the input. What is the circuitry underlying this logic?

Neurons in the primary visual cortex are laterally connected to other neurons on the cortical surface and derive input from them. Experiments have shown that the spatial extent on the cortical surface from which neurons derive input from other neurons through such lateral interactions varies with the contrast of the stimulus [26]. In the absence of stimulus contrast, spike-triggered traveling waves of activity propagate over large areas of cortex. As contrast is increased, the waves become weaker in amplitude and travel over increasingly small distances. These experiments suggest that the change in spatial summation area with increasing stimulus contrast may be consistent with the change in the decay constants of the traveling wave activity. However, no extant experiment directly links changes in summation in visual space to changes in integration on the cortical surface, and no explicit model of neural architecture has been shown to simultaneously account for, and thus connect, the input-dependence of spatial summation and lateral integration in V1. The latter one is our aim, and a crucial clue will come from the input-dependence of latencies.

Recently, a critically-balanced network model of cortex was proposed to explain the contrast dependence of functional connectivity [27]. It was shown that in the absence of input, the model exhibits wave-like activity with an infinitely-long ranged susceptibility, while in the presence of input, perturbed network activity decays exponentially with an attenuation constant that increases with the strength of the input. These results are in direct agreement with Nauhaus et al. [26].

We will now demonstrate that a similar model also leads to adaptive scales of spatial integration in visual space. Our model makes two key assumptions. The first is a local, not just global, balance of excitation and inhibition across the entire network; all eigenmodes of the network are associated with purely imaginary eigenvalues. It has been shown that such a critically-balanced configuration can be achieved by simulating a network of neurons with connections evolving under an anti-Hebbian rule [28]. The second key assumption is that all interactions in the network are described by the connectivity matrix; nonlinearities do not couple distinct neurons in the network.

In dynamical systems theory, the existence of purely imaginary eigenvalues implies the existence of an invariant subspace of the activity known as a center manifold [29]. In contrast to hyperbolic fixed points [29], where the linearization of the system fully describes the topological structure of the local solution, dynamics on center manifolds are not dominated by the linearization. This leads to complex nonlinear behavior where nonlinear terms and input parameters play a crucial role in determining the properties of the system such as relaxation timescales and correlation lengths [27]. In the model presented in this paper, the center manifold is full dimensional, and thus the rich, complex behavior we will be discussing is not surprising. We postulate that, in general, neural systems utilize center manifolds in order to flexibly integrate sensory input. Our aim in this paper is not to provide a detailed neuroanatomical and physiological model of V1, but rather to construct a toy model which provides an existence proof that center manifold dynamics can account for and connect three input-dependent computational properties of V1. This approach of constructing a minimalist toy model to explain how a given mechanism can lead to a particular set of properties is common practice in theoretical physics and is the underlying philosophical approach of several well known theoretical neuroscience models, e.g. Wilson-Cowan equations, Hopfield networks, and Kuramoto models [30, 31].

There are a number of examples of dynamical criticality in neuroscience, including experimental studies in motor cortex [32], theoretical [33] and experimental studies [34] of line attractors in oculomotor control, line attractors in decision making [35], Hopf bifurcation in the auditory periphery [2–4, 36, 37], olfactory system [38], and theoretical work on regulated criticality [39]. More recently, Solovey et al. [40] performed stability analysis of high-density electrocorticography recordings covering an entire cerebral hemisphere in monkeys during reversible loss of consciousness. Performing a moving vector autoregressive analysis of the activity, they observed that the eigenvalues crowd near the critical line. During loss of consciousness, the numbers of eigenmodes at the edge of instability decrease smoothly but drift back to the critical line during recovery of consciousness.

Dynamical criticality is distinct from statistical criticality [41] which is related to the statistical mechanics of second-order phase transitions. It has been proposed that neural systems [42], and more generally biological systems [43], are statistically critical in the sense that they are poised near the critical point of a phase transitions [44, 45]. Statistical criticality is characterized by power law behavior such as avalanches [46-48] and long-range spatiotemporal correlations [49, 50]. While both dynamical criticality and statistical criticality have had success in neuroscience, their relation is still far from clear [28, 43, 51].

We also examine the dynamics of the system and show that its activity exponentially decays to a limit cycle over multiple timescales, which depend on the strength of the input. Specifically, we find that the temporal exponential decay constants increase with increasing input strength. This result agrees with single-neuron studies which have found that response latencies in V1 decrease with increasing stimulus contrast [12, 52–54]. We now turn to describing our model.

Methods

Let be the activity vector for a network of neurons which evolve in time according to the normal form equation: In this model, originally proposed by Yan and Magnasco [27], neurons interact with one another through a skew-symmetric connectivity matrix A. The cubic-nonlinear term in the model is purely local and does not couple the activity states of distinct neurons, while the external input to the system may depend on time and have a complex spatial pattern.

The original model considered a 2-D checkerboard topology of excitatory and inhibitory neurons. For theoretical simplicity and computational ease, we will instead consider a 1-D checkerboard layout of excitatory and inhibitory neurons which interact through equal strength, nearest neighbor connections (Fig 2). In this case, A = (−1) s(δ + δ), where i, j = 0, 1, …, N − 1 and s is the synaptic strength. Boundary conditions are such that the activity terminates to 0 outside of the finite network.

Fig 2

Simplest connectivity matrix A.

A finite line of excitatory and inhibitory neurons. White nodes represent excitatory neurons. Black nodes are inhibitory. All connections have strength of the same magnitude.

We are specifically interested in the time-asymptotic response of the system, but explicitly integrating the stiff, high-dimensional ODE in (1) is difficult. Fortunately, we can bypass numerical integration methods by assuming periodic input of the form I(t) = Fe, where and look for solutions X(t) = Ze, where . Substituting these into (1), we find that:

And define g(Z) to be equal to the right hand side of (2).

The solution of (2) can be found numerically by using the multivariable Newton-Raphson method in : where and are the concatenations of the real and imaginary parts of Z and g, respectively. J is the Jacobian of with respect to

Results

Following the lead of previous work [18] and experimental studies [55], we assume the input strength from lateral geniculate nucleus to V1 to be a linear function of the stimulus contrast. To then test how the response of a single neuron in our network varies with both the contrast and length of the stimulus, we select a center neuron at index c and then calculate, for a range of input strengths, the response of the neuron as a function of input length around it. Formally, for each input strength level , we solve (3) for: where k = 0, …, N − 1, describes the spatial shape of the input, and 2l+ 1 is the length of the input in number of neurons. The response of the center neuron is taken as the modulus of Z, and we focus on the case where ω is an eigenfrequency of A and v the corresponding eigenvector.

The results for a 1-D checkerboard network of 64 neurons is shown in Fig 3. Here we fix a center neuron and sweep across a small range of eigenfrequencies ω of A. The curves from bottom to top correspond to an ascending order of base-2 exponentially distributed input strengths C = 2. For all eigenfrequencies, the peak of the response curves shift towards larger input lengths as the input strength decreases. In fact, for very weak input, the response curves rise monotonically over the entire range of input lengths without ever reaching a maximum in this finite network. This is in contrast to the response curves corresponding to strong input, which always reach a maximum but, depending on the eigenfrequency, exhibit varying degrees of response suppression beyond the maximum. This is consistent with variability of response suppression in primary visual cortex studies [12, 13]. In Fig 3, eigenfrequencies ω = 1.92, 1.96, 1.99 show the greatest amount of suppression while the others display little to none.

Fig 3

Length-response curves for different eigenfrequencies.

In each panel, which corresponds to a different eigenfrequency of A, we plot the response of a neuron (with index N/2) as a function of input length for a group of exponentially distributed input strengths. The blue arrow in the first plot indicates the direction of increasing input strength. The length of the input is recorded in number of neurons and the response is taken as the modulus of the amplitude of the time-asymptotic stable limit cycle.

To understand why certain eigenfrequencies lead to suppression, we fix the eigenfrequency to be ω = 1.92 and examine the response curves of different center neurons. The response of four center neurons (labeled by network position) and the modulus of the eigenfrequency’s corresponding eigenvector are plotted in Fig 4. The center neurons closest to the zeros of the eigenvector experience the strongest suppression for long line lengths. Neuron 38 closer to the peak of the eigenvector’s modulus experiences almost zero suppression. This generally holds for all eigenvectors and neurons in the network as all eigenvectors are periodic in their components with an eigenvalue-dependent spatial frequency. The periodicity of the eigenvectors arise from the fact that A2, which shares the same eigenvectors as A, is a circulant matrix.

Fig 4

Length-response curves of different neurons.

The top plot depicts the modulus of the eigenvector corresponding to eigenfrequency ω = 1.92. The 4 panels below are plots of the length-response curves for 4 different neurons in the network. The position of the neurons relative to the shape of the eigenvector are noted with the gray bars and arrows.

To strengthen the connection between model and neurophysiology, one can consider a critically-balanced network with an odd number of neurons so that 0 is now an eigenfrequency of the system. In our model, input associated with the 0-eigenmode represents direct current input to the system which is what neurophysiologists utilize in experiments; the visual input is not flashed [12, 13]. Contrary to the even case, long range connections must be added on top of the nearest neighbor connectivity in order to recover periodic eigenvectors and hence suppression past the response curves maximums.

Next, we show that the network not only selectively integrates input as a function of input strength but also operates on multiple time scales which flexibly adapt to the input. This behavior is not surprising given that in the case of a single critical Hopf oscillator, the half width of the resonance, the frequency range for which the oscillator’s response falls by a half, is proportional to the forcing strength of the input, where Γ is the half-width F the input strength [2]. Thus, decay constants in the case of a single critical oscillator should grow with the input forcing strength as .

Assuming input Fe, as described above, the network activity x(t), given by (1), decays exponentially in time to a stable limit cycle, X(t) = Ze. This implies that for any neuron i in the network, |x(t)| = e− f(t)+ |Z| during the approach to the limit cycle. We therefore plot log(||x(t)| − |Z||) over the transient decay period and estimate the slope of the linear regimes. We do this for a nearly network size input length (input length = 29, N = 32) and a range of exponentially distributed input strengths. In Fig 5, we plot representative transient periods of a single neuron corresponding to 3 input strengths: 2−10, 2−4, and 22. For weak input there is a fast single exponential decay regime (red) that determines the system’s approach to the stable limit cycle. As we increase the input, however, the transient period displays two exponential decay regimes: the fast decay regime (red) which was observed in the presence of weak input and a new slow decay regime (blue) immediately preceding the stable limit cycle. For very large input strength, the slow decay regime becomes dominant. The multiple decay regimes is a surprising result which doesn’t appear in the case of a single critical Hopf oscillator.

Fig 5

Input-strength-dependent timescales.

For neuron i in the network, we plot log(||x(t)| − |Z||) as a function of time for 3 different input strengths. Linear regions correspond to exponential decay. In the presence of weak input, a fast decay regime (red) guides the dynamics towards a stable limit cycle. For the intermediate input strength, a new, distinct slow decay regime appears (blue), which becomes dominant for strong input.

We estimate the exponential decay constants as a function of input strength and plot them on a log-log scale in Fig 6. The red circles correspond to the fast decay regime, while the blue circles correspond to the slow decay regime, which becomes prominent for large forcings. We separately fit both the slow and fast decay regimes with a best fit line. Unsurprisingly, the slopes of the lines are equal and approximately . Thus, the decay constants grow with the input as ∝ , where F is the input strength. This implies that the system operates on multiple timescales dynamically switching from one to another depending on the magnitude of the forcing. Larger forcings lead to faster network responses.

Fig 6

Input-strength dependence of exponential decay constants.

The temporal exponential decay constants for a range of input strengths is depicted above. A fast decay regime (red circles) is accompanied by a slow decay regime (blue circles) at large input strengths. Each decay regime is separately approximated by a least squares line (dashed lines) in log-log space.

In this paper, we consider a line of excitatory and inhibitory neurons, but our results hold equally well for a ring of neurons with periodic boundary conditions and appropriately chosen long range connections. Ring networks have extensively been studied as a model of orientation selectivity in V1 [56-62]. In agreement with recent findings [63], the critically-balanced ring network exhibits surround suppression in orientation space when long range connections are added on top of nearest neighbor connectivity.

Conclusion

We have shown that a simple dynamical system poised at the onset of instability exhibits an input-strength-dependent scale of integration of the system’s input and input-strength-dependent response latencies. This finding strongly complements our previous results showing that a similar nonlinear process with fixed, nearest neighbor network connectivity leads to input-dependent functional connectivity. This system is thus the first proposed mechanism that can account for contrast dependence of spatial summation, functional connectivity, and response latencies. In this framework, these three characteristic properties of signal processing in V1 are intrinsically linked to one another. As our model is just a toy model of center manifold dynamics, we do not suggest that a ring or 1-D line topology is necessarily present in V1 anatomy; although, if the brain does indeed utilize center manifolds in the processing of sensory input, it might be the case that the full high dimensional phase space of cortical dynamics could be reduced to simple, low dimensional structures on the center manifold.

The theory of V1 dynamics presented in this paper makes two testable predictions. The first is the specific form of the relationship between spatial and temporal frequencies of neural activity in V1. In physics, the relationship between the spatial and temporal frequencies in a system is known as the dispersion relation. Yan and Magnasco [27] have shown that the dispersion relation of the system considered in this paper and described by Eq (1) is elliptical, c2 k2+ ω2 = 1, where k and ω are the spatial and temporal frequencies, respectively, and c is a constant. If our theory is correct, multielectrode array recordings in V1 should reveal elliptic dispersion relations. Unfortunately, we are unaware of any studies that have examined dispersion relations in V1. Our theory also makes testable predictions regarding temporal response latencies in V1. In particular, our theory implies that temporal decay constants in V1 should increase as a power law with the contrast level. The predicted exponent for the power law is 2/3, which could be tested for in single or multielectrode recordings. Experiments could also test for the presence of multiple relaxation timescales, which our model predicts.