{'fraction_non_alphanumeric': 0.03585911685190748, 'fraction_numerical': 0.011827876239110843, 'mean_word_length': 4.598696657557284, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural wave interference and intrinsic tuning in distributed excitatory-inhibitory networks

We developed a model of cortical computation that implements key features of cortical circuitry and is capable of describing propagation of neural signals between cortical locations in response to spatially distributed stimuli. The model is based on the canonical neural circuit that consists of excitatory and inhibitory cells interacting through reciprocal connections, with recurrent feedback. The canonical circuit is used as a node in a distributed network with nearest neighbor coupling between the nodes. We find that this system is characterized by intrinsic preference for spatial frequency. The value of preferred frequency depends on the relative weights of excitatory and inhibitory connections between cells. This balance of excitation and inhibition changes as stimulus contrast increases, which is why intrinsic spatial frequency is predicted to change with contrast in a manner determined by stimulus temporal frequency. The dynamics of network preference is consistent with properties of the cortical area MT in alert macaque monkeys.


Introduction
We sought to develop a framework for modeling visual cortical processes satisfying a number of requirements: Biological realism: Include key features of cortical circuity that consists of excitatory and inhibitory neurons connected reciprocally and affected by recurrent feedback.
Distributed dynamics: Allow for propagation of signals between cortical locations, in order to predict responses to stimuli distributed spatially and temporally.
Mathematical tractability: Essential properties of dynamics should be tractable by methods of analysis and not numerical simulation alone.
Accordingly, we developed a spatially distributed model based on the canonical excitatory-inhibitory circuit introduced by Cowan (1972, 1973). The canonical circuit is illustrated in Figure 1A.
It consists of two units ("cells"): excitatory and inhibitory (E and I), connected reciprocally and recurrently. Previous studies of this circuit suggested that it can serve as a universal motif for dynamic models of cortical neural networks (e.g., Tsodyks et al., 1997;Ozeki et al., 2009;Ahmadian et al., 2013;Jadi and Sejnowski, 2014;Rubin et al., 2015). 1 1 Of particular note is the discovery by Tsodyks et al. (1997) that inhibitoryinhibitory connections are essential for providing stability in the otherwise unstable neural networks with recurrent excitation. This work showed that the architecture with recurrent excitation and inhibitory-inhibitory stabilization allows for intriguing variety of dynamic regimes in the network.
In our spatially distributed model, the canonical circuit serves as a network motif: a repetitive "node" in a larger network ( Figure 1B). This network architecture is a conservative generalization of the canonical circuit: from the spatially local system of Wilson and Cowan to a spatially distributed system suitable for modeling neural response to stimuli distributed both spatially and temporally (Savel'ev and Gepshtein, 2014). That is, the excitatory and inhibitory units of each node are each connected to both the excitatory and inhibitory units of the neighboring nodes, which is a form of connectivity least removed from the canonical model.
We used two analytical approaches. The first approach was inspired by the method of solving linear differential equations using Green's functions (Riley et al., 2006), commonly used in the analysis of mechanical and electrical systems (Challis and Sheard, 2003). A Green's function is a response of the system to impulse activation: a signal localized spatially or temporally. Having found the Green's function of the system, responses to distributed activation are represented by convolving the Green's function with the distributed activation. This approach revealed that the spatially distributed canonical circuit has a number of interesting emergent properties. One of these properties is A. The canonical Wilson-Cowan circuit contains one excitatory unit (E) and one inhibitory unit (I), each with recurrent feedback, and connected reciprocally (i.e., with the excitatory cell connected to the inhibitory cell, and vice versa). In the figure, the lines ending with small filled circles and with Tjunctions represent, respectively, excitatory and inhibitory connections between the cells.
B. The canonical circuits are arranged as "nodes" in a chain with nearest-neighbor coupling. Indices l of the nodes (notated as E l and I l ) indicate node locations in the chain. Currents i E (l) and i I (l) are the inputs into each node generated by the stimulus. the intrinsic preference of the system to the spatial frequency of stimulation. The system responds to spatially periodic stimuli of a certain spatial frequency more vigorously than other frequencies.
In the second approach, we investigated how this intrinsic stimulus preference changes at high stimulus contrasts, where the system becomes nonlinear. We found that the intrinsic spatial frequency of the system strongly depends on the relative weights of excitation and inhibition in cell connections. We also found that this balance of excitatory and inhibitory weights changes as luminance contrast of stimuli increases, which is why system preferences are predicted to change with contrast, in a manner determined by stimulus temporal frequency. The dynamics of system preference revealed in this study departs from the traditional view of neuronal selectivity, but our predictions happen to be strikingly similar to the dynamics of neuronal preferences discovered recently in a study of the Middle Temporal (MT) cortical area of alert macaque monkeys (Pawar et al. 2017(Pawar et al. , 2018Pawar et al., in preparation).

Approach
Consider a system of Wilson-Cowan equations for a chain of "nodes" each containing an excitatory cell ("E") and an inhibitory cell ("I"), as shown in Figure 1: The variables r E and r I represent the firing rates of the excitatory and inhibitory cells, τ E represents the relaxation time of excitation (in units of the relaxation time of inhibition), l is the node index in the chain, and g E , g I are sigmoid functions. W E and W I are positive weights: W E = w EE r E (l) +w EE r E (l + 1) +w EE r E (l − 1) − w EI r I (l) +w EI r I (l + 1) +w EI r I (l − 1) + i E (l, t), W I = w IE r E (l) +w IE r E (l + 1) +w IE r E (l − 1) − w I I r I (l) +w I I r I (l + 1) +w I I r I (l − 1) where w andw represent the weights of connections respectively within and between the nodes. Assuming that stimulus inputs i E , i I and cell responses r E , r I vary slowly on the scale of inter-node distance, we replace the discrete node index l by a continuous spatial variable x, to obtain: ∂ 2 r I ∂x 2 . By expanding the inverse sigmoid functions g −1 E and g −1 I in the Taylor series up to the third order, we obtain from the above the following pair of coupled partial differential equations: where K are the "interaction constants" reflecting the strengths of connections between the excitatory and inhibitory parts of the network (K EE = w EE + 2w EE , K EI = w EI + 2w EI , etc), and D are the "diffusion constants" reflecting the strengths of connections between the nodes (D EE =w EE , D EI =w EI , etc) and thus are responsible for spatial propagation of excitatory and inhibitory influences through the network. The Taylor expansion coefficients are β (second order) and γ (third order). Among these coefficients, γ and β can be different for excitatory (β E , γ E ) and inhibitory (β I , γ I ) cells; they define the degree of nonlinearity in the system. Parameter α describes how the input current i is divided between the excitatory and inhibitory cells. That is, in a system activated by the spatiotemporal stimulus j(x, t), the input currents

Point-source waves
First, consider how the system responds to a "point" stimulus j(x, t) = δ(x) that activates a singe node of the chain. This question requires that we solve system (2) with zero j(x, t), separately for x < 0 and For the analysis of system response to localized stimuli, we notate responses as r E (x) = G E (x) and r I = G I (x). We find the solutions separately for the left and right sides of the stimulus and then match the solutions at x = 0 so that: 2 We then substitute the solutions for x < 0 and x > 0 into (3) and obtain: where sign(x) = x/|x| for x = 0 and sign(0) = 0, and where 3 3 Notice that k 2 n = b for d = 0 which will serve a useful reference for Equation 13.
Equation 4 determines a spatial oscillation (illustrated in Figure 2) with the spatial frequency of k n . The period of spatial oscillation is The term λ determines the rate of spatial decay of system response. For very small λ (λ 1 and thus Spatial frequency k n can be called the natural or intrinsic spatial frequency of the system.
For stimuli more complex than a delta function, the general static solution in linear approximation can be written as: Figure 2: Response of the distributed circuit to a small spatial stimulus ("point stimulus") that activates a single node of the circuit.
The point activation propagates through the chain and forms a steady-state spatial pattern distributed across the chain. This pattern is an instance of neural wave and it has the form of spatially damped oscillation. suggesting how our model is related to the standard description of visual neural mechanisms in terms of linear filters (Campbell and Robson, 1968;Koenderink and van Doorn, 1987;DeValois and DeValois, 1988;Graham, 1989;Watson and Ahumada, 2016). To account for nonlinear properties of neural responses, the standard approach assumes that the stage of linear filters is followed by a separate nonlinear stage. In contrast to the standard picture, a single mechanism is used in our approach to account for both linear and nonlinear regimes, without postulating a separate nonlinear process. Indeed, by analogy of our system with a system of coupled oscillators (e.g., §28 in Landau and Lifshitz, 1969), it is expected that the intrinsic frequency of our system should change with contrast, as we show just below.

Distributed periodic stimulation
Next we studied system responses to higher contrasts, where nonlinear effects become significant, and so we should consider terms γ E , γ I , β E , β I . In writing the solutions of (2) we have found, however, that our main results (10-13) about the interaction of stimulus contrast and intrinsic spatial frequency of the circuit do not change qualitatively when terms β are ignored. We therefore present results of this analysis in a simplified form, omitting the quadratic terms in Taylor expansions. And to obtain steady-state solutions, in this section we also ignore time derivatives and the time course of stimulus (which we do consider in the next section). We derive system response to periodic stimuli of the form: where k is spatial frequency. We find the solutions in the form where constants E and I represent the amplitudes of excitatory and inhibitory activations. While keeping only the main harmonics in the solution, we obtain: The system (9) can be solved iteratively: where C = j 0 /µ is the effective stimulus contrast. 4 The new sys-4 As we show below, some of the coefficients used in the system (10-11) have an immediate interpretation in terms of intrinsic preference of the neural chain.
In particular, b determines the location of the peak of selectivity at low stimulus contrast and zero temporal frequency (12), and d determines the degree of selectivity (i.e., the "width" of tuning to spatial frequency). The curves in each panel are iterative solutions of (10-11). Each curve represents activation of one excitatory cell (E 0 in Figure 1B) at one stimulus contrast, plotted as a function of stimulus spatial frequency. Stimulus spatial frequency is normalized to the resonant spatial frequency of the network measured at the lowest tested contrast (the red curves).
In general, increasing stimulus contrast reveals different magnitudes of the resonant spatial frequency: falling resonant frequency in a system dominated by inhibition (A); rising resonant frequency in a system dominated by excitation (B).
tem (10-11) allows us to derive results of the (n + 1)th iterations for E n+1 and I n+1 using results of the nth iterations for E n and I n . The iterative procedure converges for the values of coefficients described below; the results are plotted in Figure 3. Figure 3 is a plot of iterative solutions of system (10-11) with the coefficients (ξ E , ξ I ) set to (1, 0.2) in panel B and (0.2, 1) in panel C. Effective stimulus contrast is varied by setting C = j 0 /µ to five values (0.005, 0.01, 0.02, 0.03, 0.04) for which the responses of 100 iterations are plotted against the spatial frequency k, normalized to the resonance frequency, k res at the lowest contrast (shown in red). The plots reveal that cell tuning shits with contrast: toward lower k in the inhibition-dominated network (ξ I > ξ E , panel B) and toward higher k in the excitation-dominated network (ξ E > ξ I , panel C). Other parameters are b=1, d=0.1, Notice that in the system (10-11), σ affects the bandwidth of response, and ξ shifts the location of the maximum of E on the dimension of spatial frequency k. The magnitude of k at which E reaches its maximum is the characteristic (or resonant) spatial frequency of the network, k = k res . Notice also that the maximum of E in this system is found near the minimum of the denominator of (10-11). Conditions of this minimum can be written explicitly: Remarkably, the terms for I and E in (12) have opposite signs. To appreciate this result, suppose that the activation of excitatory and inhibitory cells within a node have similar magnitudes A = E = I. We can therefore rewrite (12) as Recall that b, ξ I , and ξ E are constants and that A is an increasing function of contrast. It follows that the second term of (13), which is A 2 (ξ I − ξ E ), will increase with contrast when ξ I > ξ E and the resonant frequency k res will decrease ( Figure 3A). And when ξ I < ξ E , the second term of (13) will decrease with contrast and k res will increase ( Figure 3B). In summary, the network's intrinsic preference for stimulus spatial frequency (k res ) will generally change with stimulus contrast: k res will increase with contrast in a system dominated by excitation, and k res will decrease with contrast in a system dominated by inhibition.

Distributed spatiotemporal stimulation
Now we consider dynamic stimuli: drifting luminance gratings characterized in terms of both spatial and temporal frequencies, which we represent by j(x, t) = j 0 cos(ωt − kx), where ω is temporal frequency. Here we concentrate on low luminance contrasts, where the nonlinear terms can be omitted. As we show below, results of this analysis have a form that allows us to make predictions for the system behavior also at high luminance contrasts. We therefore seek solutions of (2) in the following form: In contrast to the in-phase response of the system to purely spatial stimuli, the response to dynamic stimuli is characterized by temporal delays relative to the responses to purely spatial stimuli. This means that responses to spatiotemporal stimuli have both sine and cosine components, even when the stimulus only contains cosine components. By substituting (14) in (2) with β = 0 and γ = 0, we find the following solutions for the amplitudes of activation of the excitatory (E c and E s ) and inhibitory (I c and I s ) cells: where Neural wave interference and intrinsic tuning 9 and where the denominator of every solution is the same: with ∆D = D EE − D I I and ∆K = K I I − K EE + 2. In spite of its complexity, equation (15) can be plotted numerically, as we did in panels A and C of Figure 4. The plots reveal a family of the familiar tuning functions.
As before, we find the intrinsic spatial frequency of the system by looking for the system's maximum response. The latter is expected where the denominator approaches its minimum. Since our goal is to understand how the maximum is controlled by temporal frequency ω, in the following we focus on the analysis of the denominator in different regions of the parameters space.
To begin, notice that the first term of the denominator H depends exclusively on spatial frequencies, the last term depends exclusively on temporal frequencies, while the middle term depends on the interaction of spatial and temporal frequencies. This structure suggests that the system can be studied using three distinct approaches: (1) finding the natural spatial frequency of the system (i.e., its maximum response) under variation of spatial frequency k, defined by equation ∂H/∂k = 0 for a fixed temporal frequency ω, (2) finding the natural temporal frequency of the system under variation of temporal frequency ω, defined by equation ∂H/∂ω = 0 for a fixed spatial frequency k, and (3) finding the natural tuning speed of the system for certain combination of k and ω, and under joint variation of k and ω, defined by equation ∂H/∂k = ∂H/∂ω = 0.
Here we pursue first two of these approaches, in the next two sections.

Variation of spatial frequency
Taking the derivative ∂H/∂k | k res (ω) = 0 results in the following equation: Notice that for ω = 0 the above equation reduces to reproducing our previous result for purely spatial stimuli (ω = 0) in equation (12) with the nonlinear terms removed. Figure 4: Intrinsic tuning to stimulus spatial frequency is predicted to depend on stimulus temporal frequency.
A. Response functions are plotted at different values of stimulus temporal frequency ω when the system is in Regime 1. The peak of each tuning function is the intrinsic spatial frequency of the system for the value of ω displayed in the corresponding box. Intrinsic spatial frequency k res increases with ω.
B. The plot at left, in red, describes Regime 1. Here intrinsic spatial frequency of the system increases with ω, approaching a vertical asymptote from the left. The plot at right, in black, describes Regime 2 (as in panel C), in which intrinsic spatial frequency decreases with ω, approaching the same vertical asymptote from the right.
C. Response functions are plotted as in panel A for the system in Regime 2.
Here intrinsic spatial frequency decreases with ω.
For the sake of understanding how stimulus temporal frequency (ω) affects resonant spatial frequency (k res ) of the neural chain, it is useful to write b as k 2 res (ω = 0), which we represent by k 2 res (0). For nonzero ω, we rewrite equation (16) by approximating k res (ω) ≈ k res (0) in the denominator of (16) so we can study the influence of ω of k res : Given that ω 2 is positive, the sign of the denominator of (17) will determine how ω affects the system behavior. In other words, varying ω can change k res two ways, implying two spatial regimes: [S1] When the denominator of (17) is negative, the numerator must be positive, implying that k res (ω) 2 − k 2 res (0) > 0. Here k res increases with ω, as shown is Figure 4A, approaching an asymptote where the denominator is zero. This regime is represented in Figure 4B by the red curve that approaches the asymptote (here represented by the vertical dashed line) from the left side, which is the side of low k res .
[S2] When the denominator of (17) is positive, the numerator must be negative, implying that k res (ω) 2 − k 2 res (0) < 0. In this regime, k res decreases with ω as shown in Figure 4C. This regime is represented in Figure 4B by the black curve, which approaches the same asymptote as above, yet from the right side, which is the side of high k res .
We know from our previous analysis (12-13), that the shift of tuning to higher spatial frequencies is observed at higher stimulus contrasts. It is therefore reasonable to expect that regime S2 should be found at higher stimulus contrasts.
This prediction has been recently corroborated in a physiological study of neuronal selectivity in the Middle Temporal (MT) cortical area of alert macaque monkeys. Pawar et al. (2017Pawar et al. ( , 2018 found that spatial frequency tuning of MT neurons depends on the luminance contrast and temporal frequency of stimulation. As contrast increased, preferences of cortical neurons shifted towards lower or higher spatial frequencies, more often the latter than the former. The changes in spatial frequency preference were largest at low temporal frequency, and they decreased as temporal frequency increased. Notably, preferred spatial frequencies (1) increased with temporal frequency at low stimulus contrast, and (2) decreased with temporal frequency at high stimulus contrast, in agreement with the pattern of predictions summarized in Figure 4B (Pawar et al., in preparation).

Variation of temporal frequency
Next consider the case when temporal frequency varies and spatial frequency does not. Here the maximum of r E is reach at a certain "natural" or "intrinsic" temporal frequency of the system. We estimate the natural temporal frequency using the condition ∂H/∂ω| ω=ω res = 0 and obtain: Each of the three terms on the right side of (18) can be positive or negative, as described in Table 1, producing qualitatively different patterns of behavior described in the table. Table 1: Regimes of intrinsic temporal frequency (ω res ) of the system discovered by varying temporal frequency ω for a fixed spatial frequency k of the stimulus (using the second approach listed on page 9).
The maximum of r E occurs at a non-zero temporal frequency ω even at zero spatial frequency k, i.e., ω res (0) = 0. In this regime, ω res monotonically increases as k increases. T2 positive negative positive Here ω res (0) = 0 as above, The initially positive ω res decreases, so the maximum of r E (ω) can disappear, within an interval of k. At larger k, ω res increases with k. T3 positive positive negative Again, ω res (0) = 0, but ω res increases and reaches a maximum, after which it drops to zero at a certain k. No maximum of r E (ω) occurs at larger k. No resonance occurs at k = 0, but it appears at some k = 0 when ω res = 0, and then increases with increasing k.  The intrinsic temporal frequency ω res of the system depends on stimulus spatial frequency k.
Here ω res is plotted as a function on k in green and blue, together with the two spatial resonance curves k res in red and black. The latter two curves are reproduced from Figure 4B. See text for detail.
The conditions listed in row T6 correspond to the conditions used to predict the interaction of k res and ω illustrated in Figure 4. A key property of this prediction is reproduced in Figure 5 (red and black curves) together with the prediction derived from row T6 (green and blue curves): The red and black curves were obtained using the first approach described on page 9, which is finding the natural spatial frequency of the system k res for a fixed temporal frequency ω.
The green and blue curves were obtained using the second approach described on page 9, which is finding the natural temporal frequency of the system ω res for a fixed spatial frequency k.
Of these results, the red and green curves correspond to low stimulus contrasts while the black and blue curves correspond to high stimulus contrasts. Accordingly, the intersection of the red and green curves, and also the intersection of the black and blue curves, correspond to the points of preferred stimulus speed of the system.

Discussion
We developed a model of spatially distributed cortical computation based on the spatially localized excitatory-inhibitory circuit originally proposed by Wilson & Cowan ( Figure 1A). In the distributed model, the localized circuits serve as "nodes" in a chain with the neighboring excitatory (E) and inhibitory (I) cells related through E-E, E-I, I-I, and I-E connections, i.e., forming a fully connected network with nearestneighbor coupling ( Figure 1B).

Intrinsic stimulus preference.
We investigated response properties of this distributed system by first modeling its response to a "point" stimulus that activates a single node in the network. The point activation propagates through the chain of nodes and forms a steady-state spatial pattern distributed across the chain. A typical steady-state pattern has the form of a spatially damped oscillation (Figure 2). The spatial frequency of this oscillation is the natural or intrinsic spatial frequency of the distributed network, determined by weights of excitatory and inhibitory connections between the cells within and between the nodes.
Neural interference. Network response to a spatially distributed stimulus can be understood in terms of spatial interference between the point-processes elicited by the stimulus on multiple nodes: At low stimulus contrasts, where system behavior is linear, the interference amounts to superposition of the point-generated waves. Spatially periodic stimuli elicits maximum response when the spatial frequency of the stimulus matches the intrinsic spatial frequency of the network. And at high stimulus contrasts, where system behavior is nonlinear, predicting system response requires more complex analysis. We developed an analytical solution to this problem (Equation 9) allowing one to solve the system iteratively and determine its stimulus preferences at any contrast.
Dynamics of stimulus preference. This analysis revealed that the preferred spatial frequency k res of the network should change as stimulus contrast increases, towards higher or lower frequency than at low contrasts. The direction of change depends on the balance of excitatory or inhibitory processes in the network (Equation 13). In particular, k res is predicted to increase with contrast in the network dominated by excitation, and to decrease with contrast in the network dominated by inhibition ( Figure 3).
The amount of change of k res is predicted to depend on stimulus temporal frequency ω. At low stimulus contrast, k res increases with ω, and at high stimulus contrast, k res decreases with ω ( Figure 4). The latter result can be described in terms of k res at low and high stimulus contrasts converging at a certain "attractor" value of k res as ω increases. The "attractor" value is represented in Figure 4 by the vertical asymptote.
Remarkably, the dynamics of neural preferences similar to that predicted by our model (Figure 4) was discovered recently in a study of neuronal preferences in the cortical area MT of alert macaque monkeys (Pawar et al. 2017(Pawar et al. , 2018Pawar et al., in preparation).