By systematically varying the degree of modular specialization in the feedforward projections (modularity parameter, $m$, see Materials and methods and Figure 1), we can control the segregation of stimulus-specific pathways across the network and investigate how it influences the characteristics of neural representations as the signal propagates. If the feedforward projections are unstructured or moderately structured ($m\lesssim 0.8$), information about the input fails to permeate the network, resulting in a chance-level reconstruction accuracy in the last sub-network, SSN₅, even in the absence of noise (see Figure 1b, c). However, as $m$ approaches a switching value $m_{\mathrm{switch}}\approx 0.83$, there is a qualitative transition in the system’s behavior, leading to a consistently higher reconstruction accuracy across the sub-networks (Figure 1b–e), regardless of the amount of noise added to the signal (Figure 1f, g).

Beyond this transition point, reconstruction accuracy improves with depth, that is the signal is more accurately represented in SSN₅ than in the initial sub-network, SSN₀, with an effective accuracy gain of over 40% (Figure 1d, g). While the addition of noise does impair the absolute reconstruction accuracy in all cases (see Figure 1—figure supplement 1), the denoising effect persists even if the input is severely corrupted (${\sigma }_{\xi }=3$, see Figure 1f, g). This is a counter-intuitive result, suggesting that topographic modularity is not only necessary for reliable communication across multiple populations (see Zajzon et al., 2019), but also supports an effective denoising effect, whereby representational precision increases with depth, even if the signal is profoundly distorted by noise.

The sequential denoising effect observed beyond the transition point $m_{\mathrm{switch}}\approx 0.83$ results in an increasingly accurate input encoding through progressively more precise internal representations. In general, such a phenomenon could be achieved either through noise suppression, stimulus-specific response amplification or both. In this section, we examine these possibilities by analyzing and comparing the input-driven dynamics of the different sub-networks. The strict segregation of stimulus-specific sub-populations in SSN₀ is only fully preserved across the system if $m=1$, in which case signal encoding and transmission primarily rely on this spatial segregation. Spiking activity across the different SSNs (Figure 2a) demonstrates that the system gradually sharpens the segregation of stimulus-specific sub-populations; indeed, in systems with fully modular feedforward projections, activity in the last sub-network is concentrated predominantly in the stimulated sub-populations. This effect can be observed in both excitatory (E) and inhibitory (I) populations, as both are equally targeted by the feedforward excitatory projections. The sharpening effect consists of both noise suppression and response amplification (Figure 2b), measured as the relative firing rates of the non-stimulated ${\nu }_5^{\mathrm{NS}}/{\nu }_0^{\mathrm{NS}}$ and stimulated sub-populations ${\nu }_5^{\mathrm{S}}/{\nu }_0^{\mathrm{S}}$, respectively. For ,${\displaystyle m<m_{\mathrm{s}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{c}\mathrm{h}}}$. noise suppression is only marginal and responses within the stimulated pathways are not amplified (${\displaystyle {\nu }_5^{\mathrm{S}}/{\nu }_0^{\mathrm{S}}<1}$).

Figure 2 with 1 supplement see all

Download asset Open asset

Figure 2—source data 1

Code and data for Figure 2 and related figure supplements.

https://cdn.elifesciences.org/articles/77009/elife-77009-fig2-data1-v2.zip

Download elife-77009-fig2-data1-v2.zip

Mean-field analysis of the stationary network activity (see Materials and methods and Appendix B) predicts that the firing rates of the stimulus-specific sub-populations increase systematically with modularity, whereas the untuned neurons are gradually silenced (Figure 2c, left). At the transition point $m_{\mathrm{switch}}\approx 0.83$, mean firing rates across the different sub-networks converge, which translates into a globally uniform signal encoding capacity, corresponding to the zero-gain convergence point in Figure 1d, g. As the degree of modularity increases beyond this point, the self-consistent state is lost again as the functional dynamics across the network shifts toward a gradual response sharpening, whereby the activity of stimulus-tuned neurons become increasingly dominant (Figure 2a–c). The effect is more pronounced for the deeper sub-networks. Note that the analytical results match well with those obtained by numerical simulation (Figure 2c, right).

In the limit of very deep networks (up to 50 SSNs, Figure 2d) the system becomes bistable, with rates converging to either a high-activity state associated with signal amplification or a low-activity state driven by the background input. The transition point is observed at a modularity value of $m=0.83$, matching the results reported so far. Below this value, elevated activity in the stimulated sub-populations can be maintained across the initial sub-networks (<10), but eventually dies out; the rate of all neurons decays and information about the input cannot reach the deeper populations. Importantly, for $m=0.83$, the transition toward the high-activity state is slower. This allows the input signal to faithfully propagate across a large number of sub-networks ($\approx 15$), without being driven into implausible activity states.

The departure from the balanced activity in the initial sub-networks can be better understood by zooming in at the synaptic level and analyzing how topography influences the synaptic input currents. The segregation of feedforward projections into stimulus-specific pathways breaks the symmetry between excitation and inhibition (see Figure 3a) that characterizes the balanced state (Haider et al., 2006; Shadlen and Newsome, 1994), for which the first two sub-networks were tuned (see Materials and methods). E/I balance is thus systematically shifted toward excitation in the stimulated populations and inhibition in the non-stimulated ones. Neurons belonging to sub-populations associated with the active stimulus receive significantly more net overall excitation, whereas the other neurons become gradually more inhibited. This disparity grows not only with modularity but also with network depth. Overall, across the whole system, increasing modularity results in an increasingly inhibition-dominated dynamical regime (inset in Figure 3a), whereby stronger effective inhibition silences non-stimulated populations, thus sharpening stimulus/feature representations by concentrating activity in the stimulus-driven sub-populations.

Figure 3

Download asset Open asset

Figure 3—source data 1

Code and data for Figure 3.

https://cdn.elifesciences.org/articles/77009/elife-77009-fig3-data1-v2.zip

Download elife-77009-fig3-data1-v2.zip

To gain an intuitive understanding of these effects from a dynamical systems perspective, we linearize the network dynamics around the stationary working points of the individual populations (Tetzlaff et al., 2012) in order to obtain the effective connectivity $W$ of the system (see Materials and methods and Appendix B). The effective impact of a single spike from a presynaptic neuron $j$ on the firing rate of a postsynaptic neuron $i$ (the effective weight $w_{ij}\in W$) is determined not only by the synaptic efficacies $J_{ij}$, but also by the statistics of the synaptic input fluctuations to the target cell $i$ that determine its excitability (see Materials and methods, Equation 6). This analysis reveals that there is an increase in the effective synaptic input onto neurons in the stimulated sub-populations as a function of modularity (Figure 3b). Conversely, non-stimulated neurons effectively receive weaker excitatory (and stronger inhibitory) drive and become increasingly less responsive (see Figure 3a, b). The role of topographic modularity in denoising can thus be understood as a transient, stimulus-specific change in effective connectivity.

For low and moderate topographic precision ($m\lesssim 0.83$), denoising does not occur as the effective weights are sufficiently similar to maintain a stable E/I balance across all populations and sub-networks (Figure 3a, b), resulting in a relatively uniform global dynamical state (indicated in Figure 3c by a constant spectral radius for $m\lesssim 0.83$, see also Materials and methods) and stable linearized dynamics (${\displaystyle \rho (W)<1}$).

However, as the feedforward projections become more structured, the system undergoes qualitative changes: after a weak transient ($0.83\lesssim m\lesssim 0.85$) the spectral radius $\rho$ in the deep SSNs expands due to the increased effective coupling to the stimulated sub-population (Figure 3b); the spectral radius eventually ($m\gtrsim 0.85$) contracts with increasing modularity (Figure 3c, d). Given that $\rho$ is determined by the variance of $W$, that is heterogeneity across connections (Rajan and Abbott, 2006), this behavior is expected: most weights are in the non-stimulated pathways, which decrease with larger $m$ and network depth (Figure 3b). Strong inhibitory currents (Figure 3a) suppress the majority of neurons, thereby reducing noise, as demonstrated by the collapse of the bulk of the eigenvalues toward the center for larger $m$ (Figure 3d). Indicative of a more constrained state space, this contractive effect suggests that population activity becomes gradually entrained by the spatially encoded input along the stimulated pathway, whereas the responses of the non-stimulated neurons have a diminishing influence on the overall behavior.

By biasing the effective connectivity of the system, precise topography can thus modulate the balance of excitation and inhibition in the different sub-networks, concentrating the activity along specific pathways. This results in both a systematic amplification of stimulus-specific responses and a systematic suppression of noise (Figure 2b). The sharpness/precision of topographic specificity along these pathways thus acts as a critical control parameter that largely determines the qualitative behavior of the system and can dramatically alter its responsiveness to external inputs.

How can the system generate and maintain the elevated inhibition underlying such a noise-suppressing regime? On the one hand, feedforward excitatory input may increase the activity of certain excitatory neurons in ${\mathrm{E}}_{\mathrm{i}}$ of sub-network ${\mathrm{SSN}}_{\mathrm{i}}$, which, in turn, can lead to increased mean inhibition through local recurrent connections. On the other hand, denoising could depend strongly on the concerted topographic projections onto ${\mathrm{I}}_{\mathrm{i}}$. Such structured feedforward inhibition is known to play important functional roles in, for example, sharpening the spatial contrast of somatosensory stimuli (Mountcastle and Powell, 1959) or enhancing coding precision throughout the ascending auditory pathways (Roberts et al., 2013).

To investigate whether recurrent activity alone can generate sufficiently strong inhibition for signal transmission and denoising, we maintained the modular structure between the excitatory populations and randomized the feedforward projections onto the inhibitory ones ($m=0$ for ${\mathrm{E}}_{\mathrm{i}}\to {\mathrm{I}}_{\mathrm{i}+1}$, compare top panels of Figure 4a, b). This leads to unstable firing patterns in the downstream sub-networks, characterized by significant accumulation of synchrony and increased firing rates (see bottom panels of Figure 4a, b and Figure 4—figure supplement 1a, b). These effects, known to result from shared pre-synaptic excitatory inputs (see e.g. Shadlen and Newsome, 1998; Tetzlaff et al., 2003; Kumar et al., 2008a), are more pronounced for larger $m$ and network depth (see Figure 4—figure supplement 1). Compared with the baseline network, whose activity shows clear spatially encoded stimuli (sequential activation of stimulus-specific sub-populations [Figure 4a, bottom]), removing structure from the projections onto inhibitory neurons abolishes the effect and prevents accurate signal transmission.

Figure 4 with 1 supplement see all

Download asset Open asset

Figure 4—source data 1

Code and data for Figure 4 and related figure supplements.

https://cdn.elifesciences.org/articles/77009/elife-77009-fig4-data1-v2.zip

Download elife-77009-fig4-data1-v2.zip

These effects of unstructured inhibitory projections are so marked that they can be observed even if a single set of projections is modified: this can be seen in Figure 4c, where only the ${\mathrm{E}}_4\to {\mathrm{I}}_5$ connections are randomized. It is worth noting, however, that the excessive synchronization that results from unstructured inhibitory projections (Figure 4c, bottom left, no additional input condition) can be easily counteracted by driving ${\mathrm{I}}_5$ (the inhibitory population that receives only unstructured projections) with additional uncorrelated external input. If strong enough (${\nu }_{\mathrm{X}}^{+}\approx 10\mathrm{s}\mathrm{p}\mathrm{k}/\sec$), this additional external drive pushes the inhibitory population into an asynchronous regime that restores the sharp, stimulus-specific responses in the excitatory population of the corresponding sub-network (see Figure 4c, bottom right, and Figure 4—figure supplement 1c).

These results emphasize the control of inhibitory neurons’ responsiveness as the main causal mechanism behind the effects reported. Elevated local inhibition is strictly required, but whether this is achieved by tailored, stimulus-specific activation of inhibitory sub-populations, or by uncorrelated excitatory drive onto all inhibitory neurons appears to be irrelevant and both conditions result in sharp, stimulus-tuned responses in the excitatory populations.

We have demonstrated that, by controlling the different sub-networks’ operating point, the sharpness of feedforward projections allows the architecture to systematically improve the quality of internal representations and retrieve the input structure, even if profoundly corrupted by noise. In this section, we investigate the robustness of the phenomenon in order to determine whether it can be entirely ascribed to the topographic projections (a structural/architectural feature) or if the particular choices of models and model parameters for neuronal and synaptic dynamics contribute to the effect.

To do so, we study two alternative model systems on the signal denoising task. These are structured similar to the baseline system explored so far, comprising separate sequential sub-networks with modular feedforward projections among them (see Figure 1 and Materials and methods), but vary in total size, neuronal and synaptic dynamics. In the first test case, only the models of synaptic transmission and corresponding parameters are altered. To increase biological verisimilitude and following Zajzon et al., 2019, synaptic transmission is modeled as a conductance-based process, with different kinetics for excitatory and inhibitory transmission, corresponding to the responses of $\mathrm{AMPA}$ and ${\mathrm{GABA}}_{\mathrm{a}}$ receptors, respectively, see Materials and methods and Supplementary file 3 for details. The results, illustrated in Figure 5a, demonstrate that task performance and population activity across the network follow a similar trend to the baseline model (Figures 1 and 2a, b). Despite severe noise corruption, the system is able to generate a clear, discernible representation of the input as early as SSN₂ and can accurately reconstruct the signal. Importantly, the relative improvement with increasing modularity and network depth is retained. In comparison to the baseline model, the transition occurs for a slightly different topographic configuration, $m_{\mathrm{switch}}\approx 0.85$, at which point the network dynamics converges toward a low-rate, stable asynchronous irregular regime across all populations, facilitating a linear firing rate propagation along the topographic maps (Figure 5—figure supplement 1).

Figure 5 with 1 supplement see all

Download asset Open asset

Figure 5—source data 1

Code and data for Figure 5 and related figure supplements.

https://cdn.elifesciences.org/articles/77009/elife-77009-fig5-data1-v2.zip

Download elife-77009-fig5-data1-v2.zip

The second test case is a smaller and simpler network of nonlinear rate neuron models (see Figure 5b and Materials and methods) which interact via continuous signals (rates) rather than discontinuities (spikes). Despite these profound differences in the neuronal and synaptic dynamics, the same behavior is observed, demonstrating that sequential denoising is a structural effect, dependent on the population firing rates and thus less sensitive to fluctuations in the precise spike times. Moreover, the robustness with respect to the network size suggests that denoising could also be performed in smaller, localized circuits, possibly operating in parallel on different features of the input stimuli.

Despite their ubiquity throughout the neocortex, the characteristics of structured projection pathways is far from uniform (Bednar and Wilson, 2016), exhibiting marked differences in spatial precision and specificity, aligned with macroscopic gradients of cortical organization. This non-uniformity may play an important functional role supporting feature aggregation (Hagler and Sereno, 2006) and the development of mixed representations (Patel et al., 2014) in higher (more anterior) cortical areas. Here, we consider two scenarios in the baseline (current-based) model to examine the robustness of our findings to more complex topographic configurations.

First, we varied the size of stimulus-tuned sub-populations (parametrized by $d_{\mathrm{i}}$, see Materials and methods) but kept them fixed across the network. For small sub-populations and intermediate degrees of topographic modularity, the activity along the stimulated pathway decays with network depth, suggesting that input information does not reach the deeper SSNs (see Figure 6a and Figure 6—figure supplement 1). These results place a lower bound on the size of stimulus-tuned sub-populations below which no signal propagation can occur, as reflected by the negative gain in performance for $d=0.01$ (Figure 6b). Whereas denoising is robust to variation around the baseline value of $d=0.1$ that yielded perfect partitioning of the feedforward projections (see Supplementary Materials), an upper bound may emerge due to increasing overlap between the maps ($d=0.2$ in Figure 6b). In this case, the activity may ‘spill over’ to other pathways than the stimulated one, corrupting the input representations and hindering accurate transmission and decoding. This can be alleviated by reduced or no overlap (as in Figure 6a), in which case signal propagation and denoising is successful for larger map sizes (${\displaystyle {\nu }_5^{\mathrm{S}}/{\nu }_0^{\mathrm{S}}>1}$ also for ${\displaystyle d>0.1}$). We thus observe a trade-off between map size, overlap and the degree of topographic precision that is required to accurately propagate stimulus representations (see Discussion).

Figure 6 with 1 supplement see all

Download asset Open asset

Figure 6—source data 1

Code and data for Figure 6 and related figure supplements.

https://cdn.elifesciences.org/articles/77009/elife-77009-fig6-data1-v2.zip

Download elife-77009-fig6-data1-v2.zip

Second, we took into account the fact that these structural features are known to vary with hierarchical depth resulting in increasingly larger sub-populations and, consequently, increasingly overlapping stimulus selectivity (Smith et al., 2001; Patel et al., 2014; Bednar and Wilson, 2016). To capture this effect, we introduce a linear scaling of map size with depth ($d_{i+1}=\delta +d_i$ for $i\ge 1$, see Materials and methods). The ability of the circuit to gradually clean the signal’s representation is fully preserved, as illustrated in Figure 6c. In fact, for intermediate modularity (${\displaystyle m<0.9}$) broadening the projections can further sharpen the reconstruction precision (compare curves for $\delta =0.02$ and $\delta =0$).

Taken together, these observations demonstrate that a gradual denoising of stimulus inputs can occur entirely as a consequence of the modular wiring between the subsequent processing circuits. Importantly, this effect generalizes well across diverse neuron and synapse models, as well as key system properties, making modular topography a potentially universal circuit feature for handling noisy data streams.

The results so far indicate that the modular topographic projections, more so than the individual characteristics of neurons and synapses, lead to a sequential denoising effect through a joint process of signal amplification and noise suppression. To better understand how the system transitions to such an operating regime, it is helpful to examine its macroscopic dynamics in the limit of many sub-networks (Toyoizumi, 2012; Cayco-Gajic and Shea-Brown, 2013; Kadmon and Sompolinsky, 2016). We apply standard mean-field techniques (Fourcaud and Brunel, 2002; Helias et al., 2013; Schuecker et al., 2015) to find the asymptotic firing rates (fixed points across sub-networks) of the stimulated and non-stimulated sub-populations as a function of topography (Figure 2d). For this, we can approximate the input μ to a group of neurons as a linear function of its firing rate $\nu$ with a slope $\kappa$ that is determined by the coupling within the group and an offset given by inputs from other groups of neurons (orange line in Figure 7a). With an approximately sigmoidal rate transfer function, the self-consistent solutions are at the intersections marked in Figure 7a.

Figure 7 with 1 supplement see all

Download asset Open asset

Figure 7—source data 1

Code and data for Figure 7 and related figure supplements.

https://cdn.elifesciences.org/articles/77009/elife-77009-fig7-data1-v2.zip

Download elife-77009-fig7-data1-v2.zip

Formally, all neurons in the deep sub-networks of one topographic map form such a group as they share the same firing rate (asymptotic value). The coupling $\kappa$ within this group comprises not only recurrent connections of one sub-network but also modular feedforward projections across sub-networks. For small modularity, the group is in an inhibition-dominated regime (${\displaystyle \kappa <0}$) and we obtain only one fixed point at low activity (Figure 7a, left). Importantly, the firing rate of this fixed point is the same for stimulated and non-stimulated topographic maps. Any influence of input signals applied to SSN₀ therefore vanishes in the deeper sub-networks and the signal cannot be reconstructed (fading regime). As topographic projections become more concentrated (larger $m$), $\kappa$ changes sign and gradually leads to two additional fixed points (as conceptually illustrated in Figure 7a and quantified in Figure 7b by numerically solving the self-consistent mean-field equations, see also Appendix B): an unstable one (red) that eventually vanishes with increasing $m$ and a stable high-activity fixed point (black). The bistability opens the possibility to distinguish between stimulated and non-stimulated topographic maps and thereby reconstruct the signal in deep sub-networks: in the active regime beyond the critical modularity threshold (here $m\ge m_{\mathrm{crit}}=0.76$), a sufficiently strong input signal can drive the activity along the stimulated map to the high-activity fixed point, such that it can permeate the system, while the non-stimulated sub-populations still converge to the low-activity fixed point. Note that this critical modularity represents the minimum modularity value for which bistability emerges. It typically differs from the actual switching point ${\displaystyle m_{switch}}$, which additionally depends on the input intensity.

In the potential energy landscape $U$ (see Materials and methods), where stable fixed points correspond to minima, the bistability that emerges for more structured topography $m\ge m_{\mathrm{crit}}=0.76$ can be understood as a transition from a single minimum at low rates (Figure 7c, inset) to a second minimum associated with the high-activity state (Figure 7c). Even though the full dynamics of the spiking network away from the fixed point cannot be entirely understood in this simplified potential picture (see Appendix B), qualitatively, more strongly modular networks cause deeper potential wells, corresponding to more attractive dynamical states and higher firing rates (see Figure 9—figure supplement 2).

Because the intensity of the input signal dictates the rate of different populations in the initial sub-network SSN₀ (Figure 7d), it also determines, for any given modularity, whether the rate of the stimulated sub-population is in the basin of attraction of the high-activity (see Figure 7e, solid markers and arrows) or low-activity (dashed, blue marker and arrow) fixed point. Denoising, and therefore increasing signal reconstruction, is thus achieved by successively (across sub-networks) pushing the population states toward the self-consistent firing rates.

As reported above, for the baseline network and (standard) input ($\lambda =0.05$) used in Figures 1 and 2, the switching point between low and high activity is at $m=0.83$ (blue markers in Figure 7d, f). Stronger input signals move the switching point toward the minimal modularity $m=0.76$ of the active regime (black markers in Figure 7d, f), while weaker inputs only induce a switch at larger modularities (gray markers in Figure 7d, f).

Noise in the input simply shifts the transition point to the high-activity state in a similar manner, with more modular connectivity required to compensate for stronger jitter (Figure 7g). However, as long as the mean firing rate of the stimulated sub-population in SSN₀ is slightly higher than that of the non-stimulated ones (up to 0.5 spks/sec), it is sufficient to position the system in the attracting basin of the high-rate fixed point and the system is able to clean the signal representation. This indicates a remarkably robust denoising mechanism.

In addition to properties of the input, the critical modularity marking the onset of the active regime is also influenced by neuronal and connectivity features. To build some intuition, it is helpful to consider the sigmoidal activation function of spiking neurons (Figure 8a). The nonlinearity of this function prohibits us from obtaining quantitative, closed-form analytical expressions for the critical modularity and requires a numerical solution of the self-consistency equations (Figure 7b). However, since the continuous rate model shows a qualitatively similar behavior to the spiking baseline model (see Section ‘A generalizable structural effect’), we can study a fully analytically tractable model with piecewise linear activation function (Figure 8a, b) to expose the dependence of the critical modularity on both neuron and network properties (see detailed derivations in Appendix B).

Figure 8 with 2 supplements see all

Download asset Open asset

Figure 8—source data 1

Code and data for Figure 8 and related figure supplements.

https://cdn.elifesciences.org/articles/77009/elife-77009-fig8-data1-v2.zip

Download elife-77009-fig8-data1-v2.zip

In this simple model, the output is zero for inputs below ${\mu }_{\min }=15$ and at maximum rate ${\nu }_{\max }=150$ for inputs above ${\mu }_{\max }=400$. In between these two bounds, the output is linearly interpolated $\nu (\mu )={\nu }_{\max }(\mu -{\mu }_{\min })/({\mu }_{\max }-{\mu }_{\min })$. As discussed before, successful denoising is achieved if the non-stimulated sub-populations are silent, ${\nu }^{\mathrm{NS}}=0$, and the stimulated sub-populations are active, ${\displaystyle {\nu }^{\mathrm{S}}>0}$. Note that in the following we focus on this ideal scenario representing perfect denoising, but, in principle, intermediate solutions with ${\displaystyle {\nu }^{\mathrm{S}}\gg {\nu }^{\mathrm{N}\mathrm{S}}>0}$ may also occur and could still be considered as successful denoising. Analyzing for which neuron, network and input properties this scenario is achieved, we obtain multiple conditions for the modularity that need to be fulfilled.

The first condition illustrates the dependence of the critical modularity on the neuron model (Figure 8c, purple horizontal line)

where $N_{\mathrm{C}}$ is the number of stimulus-specific sub-populations and $\alpha \le 1$ (typically with a value of 0.25) represents the (reduced) noise ratio in the deeper sub-networks, with $\alpha$ scaling the noise and $1-\alpha$ scaling the feedforward connections (see Materials and methods). This is necessary to ensure that the total excitatory input to each neuron is consistent across the network. In particular, the critical modularity depends on the dynamic range of input ${\mu }_{\max }-{\mu }_{\min }$ and output ${\nu }_{\max }$. The condition represents a lower bound on the modularity required for denoising. Importantly, while it depends on the effective coupling strength ${\displaystyle \mathcal{J}}$, the noise ratio $\alpha$ and the number of maps $N_{\mathrm{C}}$ (see Materials and methods), it does not depend on the nature of the recurrent interactions (E/I ratio) and the strength of the external background input. In addition, we find two additional critical values of the modularity (cyan and green curves in Figure 8c–e), both of which do depend on the strength of the external background input ${\nu }_{\mathrm{X}}$ and the recurrent connectivity (E/I ratio $\gamma g$):

Depending on the external input strength ${\nu }_{\mathrm{X}}$, these are either upper or lower bounds. In the denominator of these expressions, the total input (recurrent and external) is compared to the limits of the dynamic range of the neuron model. The cancellation between recurrent and external inputs in the inhibition-dominated baseline model typically yields a total input within the dynamic range of the neuron, such that modularity in feedforward connections can decrease the input of the non-stimulated sub-populations to silence them, and increase the input of the stimulated sub-populations to support their activity. The competition between the excitatory and inhibitory contributions ensures that the total input does not lead to a saturating output activity. Thus, for inhibitory recurrence, denoising can be achieved at a moderate level of modularity over a large range of external background inputs (shaded black and hatched regions in Figure 8c), which demonstrates a robust denoising mechanism even in the presence of changes in the input environment.

In contrast, if recurrent connections are absent, strong inhibitory external background input is required to counteract the excitatory feedforward input and achieve a denoising scenario (Figure 8d). Fixed points at non-saturated activity ${\displaystyle {\nu }^{\mathrm{S}}>0}$ are also present for low excitatory external input, but unstable due to the positive recurrent feedback. This is because in networks without recurrence, there is no competition between the recurrent input and the external and feedforward inputs. As a result, the input to both the stimulated and non-stimulated sub-populations is typically high, such that modulation of the feedforward input via topography cannot lead to a strong distinction between the pathways as required for denoising. In these networks, one typically observes high activity in all populations. A similar behavior can be observed in excitation-dominated networks (Figure 8e), where the inhibitory external background input must be even stronger to compensate the excitatory feedforward and recurrent connectivity and reach a stable denoising regime.

Note that inhibitory external input is not in line with the excitatory nature of external inputs to local circuits in the brain and is therefore biologically implausible. One way to achieve denoising in excitation-dominated networks for excitatory background inputs would be to shift the dynamic range of the activation function (see Figure 8—figure supplement 1), which is, however, not consistent with the biophysical properties of real neurons (distance between threshold and rest as compared to typical strengths of postsynaptic potentials). In summary, we find that recurrent inhibition is crucial to achieve denoising in biologically plausible settings.

These results on the role of recurrence and external input can be transferred to the behavior of the spiking model. While details of the fixed point behavior depend on the specific choice of the activation function, Figure 8f, h shows that there is also no denoising regime for the spiking model in case of no or excitation-dominated recurrence and a biologically plausible level of external input. Instead, one finds high activity in both stimulated and non-stimulated sub-populations, as confirmed by network simulations (Figure 8g, i). Figure 8—figure supplement 2 further confirms that even reducing the external input to zero does not avoid this high-activity state in both stimulated and non-stimulated sub-populations for ${\displaystyle m<1}$.

The analysis considered in the sections above is restricted to a system driven with a single external stimulus. However, to adequately understand the system’s dynamics, we need to account for the fact that it can be concurrently driven by multiple input streams. If two simultaneously active stimuli drive the system (see illustration in Figure 9a), the qualitative behavior where the responses along the stimulated (non-stimulated) maps are enhanced (silenced) is retained if the strength of the two input channels is sufficiently different (Figure 9b, top panel). In this case, the weaker stimulus is not strong enough to drive the sub-population it stimulates toward the basin of attraction of the high-activity fixed point. Consequently, the sub-population driven by this second stimulus behaves as a non-stimulated sub-population and the system remains responsive to only one of the two inputs, acting as a WTA circuit. If, however, the ratio of stimulus intensities varies, two active sub-populations may co-exist (Figure 9b, center) and/or compete (bottom panel), depending also on the degree of topographic modularity.

Figure 9 with 2 supplements see all

Download asset Open asset

Figure 9—source data 1

Code and data for Figure 9 and related figure supplements.

https://cdn.elifesciences.org/articles/77009/elife-77009-fig9-data1-v2.zip

Download elife-77009-fig9-data1-v2.zip

To quantify these variations in macroscopic behavior, we focus on the dynamics of SSN₅ and measure the similarity (correlation coefficient) between the firing rates of the two stimulus-specific sub-populations as a function of modularity and ratio of input intensities ${\lambda }_2/{\lambda }_1$ (see Materials and methods and Figure 9c). In the case that both inputs have similar intensities but the feedforward projections are not sufficiently modular, both sub-populations are activated simultaneously (Co-Ex, red area in Figure 9c). This is the dynamical regime that dominates the earlier sub-networks. However, this is a transient state, and the Co-Ex region gradually shrinks with network depth until it vanishes completely after approximately 9–10 SSNs (see Figure 9d).

For low modularity, the system settles in the single stable state associated with near-zero firing rates, as illustrated schematically in the energy landscape in Figure 9e, (1) (see Materials and methods, Appendix B, and Supplementary Materials for derivations and numerical simulations). Above the critical modularity value, the system enters one of two different regimes. For ${\displaystyle m>0.84}$ and an input ratio below 0.7 (Figure 9c, gray area), one stimulus dominates (WTA) and the responses in the two populations are uncorrelated (Figure 9b, top panel). Although the potential landscape contains two minima corresponding to either population being active, the system always settles in the high-activity attractor state corresponding to the dominating input (Figure 9e, (2)).

If, however, the two inputs have comparable intensities and the topographic projections are sharp enough (${\displaystyle m>0.84}$), the system transitions into a different dynamical state where neither stimulus-specific sub-population can maintain an elevated firing rate for extended periods of time. In the extreme case of nearly identical intensities (${\lambda }_2/{\lambda }_1\ge 0.9)$ and high modularity, the responses become anti-correlated (Figure 9b, bottom panel), that is the activation of the two stimulus-specific sub-populations switches, as they engage in a dynamic behavior reminiscent of WLC between multiple neuronal groups (Lagzi and Rotter, 2015; Rost et al., 2018). The switching between the two states is driven by stochastic fluctuations (Figure 9e, (3)). The depth of the wells and width of barrier (distance between fixed points) increase with modularity (see Figure 9e, (4) and Figure 9—figure supplement 2), suggesting a greater difficulty in moving between the two attractors and consequently fewer state changes. Numerical simulations confirm this slowdown in switching (Figure 9f).

We wish to emphasize that the different dynamical states arise primarily from the feedforward connectivity profile. Nevertheless, even though the synaptic weights are not directly modified, varying the topographic modularity does translate to a modification of the effective connectivity weights (Figure 3b). The ratio of stimulus intensities also plays a role in determining the dynamics, but there is a (narrow) range (approximately between 0.75 and 0.8) for which all 3 regions can be reached through sole modification of the modularity. Together, these results demonstrate that topography can not only lead to spatial denoising but also enable various, functionally important network operating points.

Until now, we have considered continuous but piecewise constant, step signals, with each step lasting for a relatively long and fixed period of $200\mathrm{ms}$. This may give the impression that the denoising effects we report only works for static or slowly changing inputs, whereas naturalistic stimuli are continuously varying. Nevertheless, sensory perception across modalities relies on varying degrees of temporal and spatial discretization (VanRullen and Koch, 2003), with individual (sub-)features of the input encoded by specific (sub-)populations of neurons in the early stages of the sensory hierarchy. In this section, we will demonstrate that denoising is robust to the temporal properties of the input and its encoding, as we relax many of the assumptions made in previous sections.

We consider a sinusoidal input signal, which we discretize and map onto the network according to the depiction in Figure 10a. This approach is similar to previous works, for instance it can mimic the movement of a light spot across the retina (Klos et al., 2018). By varying the sampling interval $dt$ and number of channels $k$, we can change the coarseness of the discretization from step-like signals to more continuous approximations of the input. If we choose a high sampling rate ($dt=1\mathrm{ms}$) and sufficient channels ($k=40$), we can accurately encode even fast changing signals (Figure 10b). Given that each input-driven SSN is inhibition-dominated and therefore close to the balanced state, the network exhibits a fast tracking property (van Vreeswijk and Sompolinsky, 1996) and can accurately represent and denoise the underlying continuous signal in the spiking activity (Figure 10c, top). This is also captured by the readout, with the tracking precision increasing with network depth (Figure 10c, bottom). In this condition, there is a performance gain of up to 50% in the noiseless case (Figure 10d, top) and similar values for varying levels of noise (Figure 10d, bottom).

Figure 10 with 1 supplement see all

Download asset Open asset

Figure 10—source data 1

Code and data for Figure 10 and related figure supplements.

https://cdn.elifesciences.org/articles/77009/elife-77009-fig10-data1-v2.zip

Download elife-77009-fig10-data1-v2.zip

Note that due to the increased number of input channels (40 compared to 10) projecting to the same number of neurons in SSN₀ as before $(800)$, for the same ${\sigma }_{\xi }$ the effective amount of noise each neuron receives is, on average, four times larger than in the baseline network. Moreover, the task was made more difficult by the significant overlap between the maps ($N_{\mathrm{C}}=20$) as well as the resulting decrease in neuronal input selectivity. Nevertheless, similar results were obtained for slower and more coarsely sampled signals (Figure 10e–g).

We found comparable denoising dynamics for a large range of parameter combinations involving the map size, number of maps, number of channels, and signal complexity. Although there are limits with respect to the frequencies (and noise intensity) the network can track (see Figure 10—figure supplement 1), these findings indicate a very robust and flexible phenomenon for denoising spatially encoded sensory stimuli.

We consider a feedforward network architecture where each sub-network (SSN) is a balanced random network (Brunel, 2000) composed of $N=10000$ homogeneous LIF neurons, grouped into a population of $N^{\mathrm{E}}=0.8N$ excitatory and $N^{\mathrm{I}}=0.2N$ inhibitory units. Within each sub-network, neurons are connected randomly and sparsely, with a fixed number of $K_{\mathrm{E}}=ϵN^{\mathrm{E}}$ local excitatory and $K_{\mathrm{I}}=ϵN^{\mathrm{I}}$ local inhibitory inputs per neuron. The sub-networks are arranged sequentially, that is the excitatory neurons ${\mathrm{E}}_{\mathrm{i}}$ in ${\mathrm{SSN}}_{\mathrm{i}}$ project to both ${\mathrm{E}}_{\mathrm{i}+1}$ and ${\mathrm{I}}_{\mathrm{i}+1}$ populations in the subsequent sub-network ${\mathrm{SSN}}_{\mathrm{i}+1}$ (for an illustrative example, see Figure 1a). There are no inhibitory feedforward projections. Although projections between sub-networks have a specific, non-uniform structure (see next section), each neuron in ${\mathrm{SSN}}_{\mathrm{i}+1}$ receives the same total number of synapses from the previous SSN, $K_{\mathrm{FF}}$.

In addition, all neurons receive $K_{\mathrm{X}}$ inputs from an external source representing stochastic background noise. For the first sub-network, we set $K_{\mathrm{X}}=K_{\mathrm{E}}$, as it is commonly assumed that the number of background input synapses modeling local and distant cortical input is in the same range as the number of recurrent excitatory connections (see e.g. Brunel, 2000; Kumar et al., 2008b; Duarte and Morrison, 2014). To ensure that the total excitatory input to each neuron is consistent across the network, we scale $K_{\mathrm{X}}$ by a factor of $\alpha =0.25$ for the deeper SSNs and set $K_{\mathrm{FF}}=\left(1-\alpha \right)K_{\mathrm{E}}$, resulting in a ratio of 3:1 between the number of feedforward and background synapses.

Within each SSN, each neuron is assigned to one or more of $N_{\mathrm{C}}$ sub-populations SP associated with a specific stimulus ($N_{\mathrm{C}}=10$ unless otherwise stated). This is illustrated in Figure 1a for $N_{\mathrm{C}}=2$. We choose these sub-populations so as to minimize their overlap within each ${\mathrm{SSN}}_{\mathrm{i}}$, and control their effective size $C_{\mathrm{i}}^{\beta }=d_{\mathrm{i}}{N}^{\beta },\beta \in [\mathrm{E},\mathrm{I}]$, through the scaling parameter $d_{\mathrm{i}}\in [0,1]$. Depending on the size and number of sub-populations, it is possible that some neurons are not part of any or that some neurons belong to multiple such sub-populations (overlap).

We study networks composed of LIF neurons with fixed voltage threshold and static synapses with exponentially decaying postsynaptic currents or conductances. The sub-threshold membrane potential dynamics of such a neuron evolves according to:

where ${\tau }_{\mathrm{m}}$ is the membrane time constant, and $R{I}^{\beta }$ is the total synaptic input from population $\beta \in [E,I]$. The background input $I^{\mathrm{X}}$ is assumed to be excitatory and stochastic, modeled as a homogeneous Poisson process with constant rate ${\nu }_{\mathrm{X}}$. Synaptic weights $J_{\mathrm{ij}}$, representing the efficacy of interaction from presynaptic neuron $j$ to postsynaptic neuron $i$, are equal for all realized connections of a given type, that is, $J_{\mathrm{EE}}=J_{\mathrm{IE}}=J$ for excitatory and $J_{\mathrm{EI}}=J_{\mathrm{II}}=gJ$ for inhibitory synapses. All synaptic delays and time constants are equal in this setup. For a complete, tabular description of the models and model parameters used throughout this study, see Supplementary files 1–5.

Following previous works (Zajzon et al., 2019; Duarte and Morrison, 2014), we choose the intensity of the stochastic input ${\nu }_{\mathrm{X}}$ and the E–I ratio $g$ such that the first two sub-networks operate in a balanced, asynchronous irregular regime when driven solely by background input. This is achieved with ${\nu }_{\mathrm{X}}=12\mathrm{spikes}/\mathrm{s}$ and $g=-12$, resulting in average firing rates of $\sim 3\mathrm{spikes}/\mathrm{s}$, coefficient of variation ($C{V}_{\mathrm{ISI}}$) in the interval $[1.0,1.5]$ and Pearson cross-correlation (CC) ≤0.01 in ${\mathrm{SSN}}_0$ and ${\mathrm{SSN}}_1$.

In Section ‘A generalizable structural effect’ we consider two additional systems, a network of LIF neurons with conductance-based synapses and a continuous firing rate model. The LIF network is described in detail in Zajzon et al., 2019. Spike-triggered synaptic conductances are modeled as exponential functions, with fixed and equal conduction delays for all synapses. Key differences to the current-based model include, in addition to the biologically more plausible synapse model, longer synaptic time constants and stronger input (see also Zajzon et al., 2019 and Supplementary file 3 for the numerical values of all parameters).

The continuous rate model contains $N=3000$ nonlinear units, the dynamics of which are governed by:

where $\mathit{x}$ represents the activation and $\mathit{r}$ the output of all units, commonly interpreted as the synaptic current variable and the firing rate estimate, respectively. The rates $r_{\mathrm{i}}$ are obtained by applying the nonlinear transfer function $\tanh \left(x_{\mathrm{i}}\right)$, modified here to constrain the rates to the interval ${\displaystyle [0,1]}$ is the neuronal time constant, ${\mathit{b}}^{\mathrm{rec}}$ is a vector of individual neuronal bias terms (i.e., a baseline activation), and $J$ and $J^{\mathrm{in}}$ are the recurrent (including feedforward) and input weight matrices, respectively. These are constructed in the same manner as for the spiking networks, such that the overall connectivity, including the input mapping onto ${\mathrm{SSN}}_0$, is identical for all three models. Input weights are drawn from a uniform distribution, while the rest follow a normal distribution. Finally, $\mathit{\xi }$ is a vector of $N$ independent realizations of Gaussian white noise with zero mean and variance scaled by ${\sigma }_{\mathrm{X}}$. The differential equations are integrated numerically, using the Euler–Maruyama method with step $\delta t=1\mathrm{ms}$, with specific parameter values given in Supplementary file 5.

We evaluate the system’s ability to recover a simple, continuous step signal from a noisy variation using linear combinations of the population responses in the different SSNs (Maass et al., 2002). This is equivalent to probing the network’s ability to function as a denoising autoencoder (Bengio et al., 2013).

To generate the $N_{\mathrm{C}}$-dimensional input signal $u\left(t\right)$, we randomly draw stimuli from a predefined set $S=\{S_1,S_2,\mathrm{\dots },S_{N_{\mathrm{C}}}\}$ and set the corresponding channel to active for a fixed duration of 200 ms (Figure 1a, left). This binary step signal $u\left(t\right)$ is also the target signal to be reconstructed. The effective input is obtained by adding a Gaussian white noise process with zero mean and variance ${\sigma }_{\xi }^2$ to $u\left(t\right)$, and scaling the sum with the input rate ${\nu }_{\mathrm{in}}$. Rectifying the resulting signal leads to the final form of the continuous input signal $z\left(t\right)={\left[{\nu }_{\mathrm{in}}\left(u\left(t\right)+\xi \left(t\right)\right)\right]}_{+}$. This allows us to control the amount of noise in the input, and thus the task difficulty, through a single parameter ${\sigma }_{\xi }$.

To deliver the input to the circuit, the analog signal $z\left(t\right)$ is converted into spike trains, with its amplitude serving as the rate of an inhomogeneous Poisson process generating independent spike trains. We set the scaling amplitude to ${\nu }_{\mathrm{in}}=K_E\lambda {\nu }_{\mathrm{X}}$, modeling stochastic input with fixed rate $\lambda {\nu }_{\mathrm{X}}$ from $K_E=800$ neurons. If not otherwise specified, $\lambda =0.05$ holds, resulting in a mean firing rate below 8 spks/sec in ${\mathrm{SSN}}_0$ (see Figure 2c).

Each input channel $k$ is mapped onto one of the $N_{\mathrm{C}}$ stimulus-specific sub-populations of excitatory and inhibitory neurons in the first (input) sub-network ${\mathrm{SSN}}_0$, chosen according to the procedure described above (see also Figure 1a). This way, each stimulus $S_{\mathrm{k}}$ is mapped onto a specific set of sub-populations in the different sub-networks, that is, the topographic map associated with $S_{\mathrm{k}}$.

For each stimulus in the sequence, we sample the responses of the excitatory population in each ${\mathrm{SSN}}_{\mathrm{i}}$ at fixed time points (once every ms) relative to stimulus onset. We record from the membrane potentials $V_m$ as they represent a parameter-free and direct measure of the population state (Duarte et al., 2018; Uhlmann et al., 2017). The activity vectors are then gathered in a state matrix ${\displaystyle X_{{\mathrm{S}\mathrm{S}\mathrm{N}}_{\mathrm{i}}}\in {\mathbb{R}}^{N^{\mathrm{E}}\times T}}$, which is then used to train a linear readout to approximate the target output of the task (Lukoševičius and Jaeger, 2009). We divide the input data, containing a total of 100 stimulus presentations (yielding ${\displaystyle T=20,000}$ samples), into a training and a testing set (80/20%), and perform the training using ridge regression (L2 regularization), with the regularization parameter chosen by leave-one-out cross-validation on the training dataset.

Reconstruction performance is measured using the normalized root mean squared error (NRMSE). For this particular task, the effective delay in the build-up of optimal stimulus representations varies greatly across the sub-networks. In order to close in on the optimal delay for each ${\mathrm{SSN}}_{\mathrm{i}}$, we train the state matrix $X_{{\mathrm{SSN}}_{\mathrm{i}}}$ on a larger interval of delays and choose the one that minimizes the error, averaged across multiple trials.

In Section ‘Reconstruction and denoising of dynamical inputs’, we generalize the input to a sinusoidal signal $x\left(t\right)=\sin \left(a\cdot t\right)+\cos \left(b\cdot t\right)$, with parameters $a$ and $b$. From this, we obtain $u\left(t\right)$ through the sampling and discretization process described in the respective section, and compute the final input $z\left(t\right)={\left[{\nu }_{\mathrm{in}}\left(u\left(t\right)+\xi \left(t\right)\right)\right]}_{+}$ as above.

To better understand the role of structural variations on the network’s dynamics, we determine the network’s effective connectivity matrix $W$ analytically by linear stability analysis around the system’s stationary working points (see Appendix B for the complete derivations). The elements $w_{ij}\in W$ represent the integrated linear response of a target neuron $i$, with stationary rate ${\displaystyle {\nu }_i}$, to a small perturbation in the input rate ${\displaystyle {\nu }_j}$ caused by a spike from presynaptic neuron $j$. In other words, $w_{ij}$ measures the average number of additional spikes emitted by a target neuron $i$ in response to a spike from the presynaptic neuron $j$, and its relation to the synaptic weights is defined by Tetzlaff et al., 2012; Helias et al., 2013:

Note that in Figure 3 we ignore the contribution $\tilde{\beta }$ resulting from the modulation in the input variance ${\sigma }_{\mathrm{j}}^2$ which is significantly smaller due to the additional factor ${\displaystyle 1/{\sigma }_i\sim \mathcal{O}(1/\sqrt{N})}$. Importantly, the effective connectivity matrix $W$ allows us to gain insights into the stability of the system by eigenvalue decomposition. For large random coupling matrices, the effective weight matrix has a spectral radius $\rho ={\max }_k\left(\mathrm{Re}\left\{{\lambda }_{\mathrm{k}}\right\}\right)$ which is determined by the variances of $W$ (Rajan and Abbott, 2006). For inhibition-dominated systems, such as those we consider, there is a single negative outlier representing the mean effective weight, given the eigenvalue ${\lambda }_{\mathrm{k}}^{*}$ associated with the unit vector. The stability of the system is thus uniquely determined by the spectral radius $\rho$: values smaller than unity indicate stable dynamics, whereas ${\displaystyle \rho >1}$ lead to unstable linearized dynamics.

For the mean-field analysis, the $N_{\mathrm{C}}$ sub-populations in each sub-network can be reduced to only two groups of neurons, the first one comprising all neurons of the stimulated SPs and the second one comprising all neurons in all non-stimulated SPs. This is possible because (1) the firing rates of the excitatory and inhibitory neurons within one SP are identical, owing to homogeneous neuron parameters and matching incoming connection statistics, and (2) all neurons in non-stimulated SPs have the same rate ${\nu }^{\mathrm{NS}}$ that is in general different from the rate of the stimulated SP ${\nu }^{\mathrm{S}}$. Here we only sketch the main steps, with a detailed derivation given in Appendix B.

The mean inputs to the first sub-network can be obtained via

where $\gamma =K_{\mathrm{I}}/K_{\mathrm{E}}$ and ${\displaystyle \mathcal{J}={\tau }_{\mathrm{m}}{K}_{\mathrm{E}}J}$. Both equations are of the form

where $\kappa$ is the effective self-coupling of a group of neurons with rate $\nu$ and input μ, and $I$ denotes the external inputs from other groups. Equation 8 describes a linear relationship between the rate $\nu$ and the input μ. To find a self-consistent solution for the rates ${\nu }^{\mathrm{S}}$ and ${\nu }^{\mathrm{NS}}$, the above equations need to be solved numerically, taking into account in addition the f–I curve $\nu \left(\mu \right)$ of the neurons that in the case of LIF model neurons also depends on the variance ${\sigma }^2$ of inputs. The latter can be obtained analogous to the mean input μ (see Appendix B). Note that for general nonlinearity $\nu \left(\mu \right)$ there is no analytical closed-form solution for the fixed points.

Starting from ${\mathrm{SSN}}_1$, networks are connected in a fixed pattern such that the rate ${\nu }_i$ in ${\mathrm{SSN}}_i$ also depends on the excitatory input from the previous sub-network ${\mathrm{SSN}}_{i-1}$ with rate ${\nu }_{i-1}$. For a fixed point, we have ${\nu }_i={\nu }_{i-1}$ (Toyoizumi, 2012). In this case, we can effectively group together stimulated/non-stimulated neurons in successive sub-networks and re-group equations for the mean input in the limit of many sub-networks, obtaining the simplified description (details see Appendix B)

The scaling terms of the firing rates incorporate the recurrent and feedforward contributions from the stimulated and non-stimulated groups of neurons. They depend solely on some fixed parameters of the system, including modularity $m$ (see Appendix B). Importantly, Equations 9 and 10 and have the same linear form as (Equation 8) Equation 8 and can be solved numerically as described above. Again, for general nonlinear $\nu \left(\mu \right)$ there is no closed-form analytical solution, but see below for a piecewise linear activation function $\nu \left(\mu \right)$. The numerical solutions for fixed points are obtained using the root finding algorithm root of the scipy.optimize package (Virtanen et al., 2020). The stability of the fixed points is obtained by inserting the corresponding firing rates into the effective connectivity Equation 6. On the level of stimulated and non-stimulated sub-populations, the effective connectivity matrix reads

from which we obtain the maximum eigenvalue $\rho$, which for stable fixed points must be smaller than 1.

The structure of fixed points for the stimulated sub-population (see discussion in ‘Modularity as a bifurcation parameter’) can furthermore be intuitively understood by studying the potential landscape of the system. The potential $U$ is thereby defined via the conservative force $F=-\frac{dU}{d{\nu }^{\mathrm{S}}}=-{\nu }^{\mathrm{S}}+\nu (\mu ,{\sigma }^2)$ that drives the system toward its fixed points via the equation of motion $\frac{d{\nu }^{\mathrm{S}}}{dt}=F$ (Wong and Wang, 2006; Litwin-Kumar and Doiron, 2012; Schuecker et al., 2017). Note that μ and ${\sigma }^2$ are again functions of ${\nu }^{\mathrm{S}}$ and ${\nu }^{\mathrm{NS}}$, where the latter is the self-consistent rate of the non-stimulated sub-populations for given rate ${\nu }^S$ of the stimulated sub-population, ${\nu }^{\mathrm{NS}}={\nu }^{\mathrm{NS}}\left({\nu }^{\mathrm{S}}\right)$ (details see Appendix B).

In Figure 9, we consider two stimuli $S_1$ and $S_2$ to be active simultaneously for 10 s. Let $S{P}_1$ and $S{P}_2$ be the two corresponding SPs in each sub-network. The firing rate of each SP is estimated from spike counts in time bins of 10 ms and smoothed with a Savitzky-Golay filter (length 21 and polynomial order 4). We compute a similarity score based on the correlation between these rates, scaled by the ratio of the input intensities ${\lambda }_2/{\lambda }_1$ (with ${\lambda }_1$ fixed). This scaling is meant to introduce a gradient in the similarity score based on the firing rate differences, ensuring that high (absolute) scores require comparable activity levels in addition to strong correlations. To ensure that both stimuli are decodable where appropriate, we set the score to 0 when the difference between the rate of $S{P}_2$ and the non-stimulated SPs was <1 spks/sec ($S{P}_1$ had significantly higher rates). The curves in Figure 9c mark the regime boundaries: coexisting (Co-Ex) where score is >0.1 (red curve); WLC where score is <−0.1 (blue); WTA (gray) and where the score is in the interval (−0.1, 0.1), and either ${\displaystyle {\lambda }_2/{\lambda }_1<0.5}$ holds or the score is 0. While the Co-Ex region is a dynamical regime that only occurs in the initial sub-networks (Figure 9d), the WTA and WLC regimes persist and can be understood again with the help of a potential $U$, which is in this case a function of the rates of the two SPs (details see Appendix B).

All numerical simulations were conducted using the Neural Microcircuit Simulation and Analysis Toolkit (NMSAT) v0.2 (Duarte et al., 2017), a high-level Python framework for creating, simulating and evaluating complex, spiking neural microcircuits in a modular fashion. It builds on the PyNEST interface for NEST (Gewaltig and Diesmann, 2007), which provides the core simulation engine. To ensure the reproduction of all the numerical experiments and figures presented in this study, and abide by the recommendations proposed in Pauli et al., 2018, we provide a complete code package that implements project-specific functionality within NMSAT (see Data availability) using NEST 2.18.0 (Jordan et al., 2019).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

This section expands on the limitations arising from the definitions of topographic modularity and map sizes used in this study. By imposing a fixed connection density on the feedforward connection matrices, the projection probabilities between neurons tuned to the same $\left(p_{\mathrm{c}}\right)$ and different $\left(p_0\right)$ stimuli are uniquely determined by the modularity $m$ and the parameter $d_0$ and $\delta$, which control the size of stimulus-specific sub-populations (see Materials and methods). For notational simplicity, here we consider the merged excitatory and inhibitory sub-populations tuned to a particular stimulus in a given sub-network ${\displaystyle {\mathrm{S}\mathrm{S}\mathrm{N}}_{\mathrm{i}}}$, with a total size $C_{\mathrm{i}}=C_{\mathrm{i}}^{\mathrm{E}}+C_{\mathrm{i}}^{\mathrm{I}}$.

Under the constraints applied in this work, the total density of a feedforward adjacency matrix between ${\mathrm{SSN}}_{\mathrm{i}}$ and ${\mathrm{SSN}}_{\mathrm{i}+1}$ can be computed as:

where $U_0^{\mathrm{i}}$ and $U_{\mathrm{c}}^{\mathrm{i}}$ are the number of realizable connections between similarly and differently tuned sub-populations, respectively. Since $U_{\mathrm{c}}^{\mathrm{i}}=N^2-U_0^{\mathrm{i}}$, we can simplify the notation and focus only on $U_0^{\mathrm{i}}$. We distinguish between the cases of non-overlapping and overlapping stimulus-specific sub-populations:

where each potential synapse is counted only once, regardless of whether the involved neurons belong to any or multiple overlapping sub-populations. This ensures consistency with the definitions of the probabilities $p_{\mathrm{c}}$ and $p_0$. Alternatively, we can express $U_0^{\mathrm{i}}$ as:

For the case with no overlap, we can derive an additional constraint on the minimum sub-populations size $C_{\mathrm{i}}$ for the required density ${\sigma }_{\mathrm{i}}$ to be satisfied, which we define in relation to the total number of sub-populations ${\displaystyle N_C}$:

The equality holds in the case of $m=1$ and all-to-all feedforward connectivity between similarly tuned sub-populations, that is, $p_{\mathrm{c}}=1$.

For an analytical investigation of the role of topographic modularity on the network dynamics, we used mean-field theory (Fourcaud and Brunel, 2002; Helias et al., 2013; Schuecker et al., 2015). Under the assumptions that each neuron receives a large number of small amplitude inputs at every time step, the synaptic time constants ${\tau }_{\mathrm{s}}$ are small compared to the membrane time constant ${\tau }_{\mathrm{m}}$, and that the network activity is sufficiently asynchronous and irregular, we can make use of theoretical results obtained from the diffusion approximation of the LIF neuron model to determine the stationary population dynamics. The equations in this section were partially solved using a modified version of the LIF Meanfield Tools library (Layer et al., 2020).