Introduction

Information about the sensory world is represented in the brain through the dynamics of neural population activity^1,2. One prominent theory about the principles that may guide the design of neural computations for sensory function is efficient coding^3,4,5. This theory posits that neural computations are optimized to maximize the information that neural systems encode about features of sensory stimuli while at the same time limiting the metabolic cost. Efficient coding has been highly influential, especially in visual neuroscience and computational vision^6,7,8,9, and has been developed to become a normative theory of how networks are organized and designed to optimally process natural sensory stimuli in visual ^10,11, auditory¹² and olfactory sensory pathways¹³.

The first normative neural network models^4,11 designed with efficient coding principles had at least two major levels of abstractions. First, information was assumed to be processed in a purely feedforward manner, whereas information processing in real neural circuits often involves recurrent or feedback computations. Second, neural dynamics was greatly simplified, ignoring the spiking nature of neural activity. Instead, in biological networks considerable amount of information are encoded or transmitted only through the millisecond-precise timing of spikes^{14,15,16,17,18,19,20}. Also, these earlier works mostly considered encoding of static sensory stimuli, whereas the sensory environment changes continuously at multiple timescales and the dynamics of neural networks encodes these temporal variations of the environment ^21,22,23,24.

Recent years have witnessed a considerable effort and success in laying down the mathematical tools and methodology to understand how to formulate efficient coding theories of neural networks with much more biological realism²⁵. This work has established the incorporation of recurrent connectivity^26,27, of spiking neurons, and of time-varying stimulus inputs^{28,29,30,31,32,33,34,35}. In these models, the efficient coding principle has been implemented by designing networks whose activity minimizes the encoding accuracy (the error between a desired representation and a linear readout of network’s activity) subject to a constraint on the metabolic cost of processing (proportional to the total number of spikes fired by a population of neurons). This double objective is captured by a loss function that trades-off encoding accuracy and metabolic cost. The minimization of the loss function is performed through a greedy approach, by assuming that a neuron will emit a spike only if this will decrease the loss. This, in turn, yields a set of leaky integrate-and-fire (LIF) neural equations which govern the network dynamics²⁸, which can also include biologically plausible non-instantaneous synaptic delays^36,35,34. These previous implementations, however, had neurons that did not respect Dale’s law. In recent work³⁷, we further extended the biological plausibility of these models by analytically deriving how to implement efficient coding in networks of spiking neurons that respect Dale’s law. These networks take the form of generalized leaky integrate-and-fire (gLIF) models of excitatory (E) and inhibitory (I) neurons endowed with spike-triggered adaptation^38,39,40, which can provide highly accurate predictions of spike times in biological networks⁴¹. Efficient spiking models have thus the potential to provide a unifying theory of neural coding through spiking dynamics of E-I circuits ^42,37 with elements that are fully biologically plausible and potentially interpretable as biophysical variables.

However, despite the major progress described above, as well the progress provided by other studies of efficient coding with spikes ^{31,43,44,33,29}, we still lack a thorough characterization of which structural, coding, biophysical and dynamical properties of excitatory-inhibitory recurrent networks of spiking neurons are directly related to efficient coding principles. Previous studies only rarely made predictions that could be quantitatively compared against experimentally measurable properties of biological neural networks. As a consequence, we still do not know which, if any, fundamental properties of cortical networks emerge directly from imposing efficient coding.

To address the above questions, we analyze systematically our biologically plausible efficient coding model of E and I neurons that respect Dale’s law³⁷ to make concrete predictions about experimentally measurable structural, coding and dynamical features of neural activity that arise from efficient coding. We systematically investigated how experimentally measurable emergent dynamical properties, such as firing rates, trial-to-trial spiking variability of single neurons, E-I balance⁴⁵ and noise correlations, relate to optimally-efficient coding. We further analyze how the organization of the connectivity arising by imposing efficient coding relates to the anatomical and effective connectivity recently reported in visual cortex, which suggests competition between excitatory neurons with similar stimulus tuning. We found that several key and robustly found empirical properties of cortical circuits match the predictions of our efficient coding network, lending support to the notion that efficient coding may be a design principle that has shaped the evolution of cortical circuits and that may be used to conceptually understand and interpret them.

Results

Assumptions and emergent properties of the efficient E-I network derived from first principles

We study the properties of a spiking neural network in which the dynamics and structure of the network are analytically derived starting from first principles of efficient coding of sensory stimuli. The model relies on a number of assumptions, described next.

The network responds to M time-varying features of a sensory stimulus, s_k(t) (e.g., for a visual stimulus, contrast, orientation, etc) received as inputs from an earlier sensory area (e.g., retina). We model features as independent Ornstein–Uhlenbeck (OU) processes (see Methods). The network’s objective is to compute a leaky integration of sensory features, a relevant computation of cortical sensory areas⁴⁶. The target representations of the network, x_k(t), are defined as

with τ a characteristic integration time-scale (Fig. 1A).

Figure 1.

The network is composed of two neural populations of excitatory (E) and inhibitory (I) neurons, defined by their postsynaptic action which respects Dale’s law. For each population, y ∈ {E, I}, we define a population readout of each feature, , as a filtered weighted sum of spiking activity of neurons in the population,

Table 1.

where is the spike train of neuron i of type y and are its decoding weights for features k = 1, …, M. As a result of the optimization, the decoding weights of the neurons are equivalent to neuron’s tuning parameters to the stimulus features (see Methods;⁴²). We draw tuning parameters from a normal distribution with zero mean and SD for E neurons and for I neurons. We assume that every neuron encodes multiple (M > 1) stimulus features⁴⁷ and define the vector as the tuning vector of neuron i (Fig. 1B).

Unlike previous approaches^28,48, we hypothesize that E and I neurons have distinct normative objectives and define cell-type specific loss functions relative to the activity of the E and I neuron types. To implement at the same time, as requested by efficient coding, the constraints of faithful stimulus representation with limited computational resources⁴⁹, we define the loss functions of E and I population as a weighted sum of a time-dependent encoding error and time-dependent metabolic cost:

We refer to β, the parameter controlling the relative importance of the metabolic cost over the encoding error, as the metabolic constant of the network. We hypothesize that population readouts of E neurons, , track the target representations, x_k(t), and the population readouts of I neurons, , track the population readouts of E neurons , by minimizing the squared error between these quantities³⁷ (see also^28,48 for related approaches). Furthermore, we hypothesize the time-resolved metabolic cost to be proportional to the estimate of a momentary firing rate of the neural population. We thus define the variables of loss functions in Eq. 3 as

where is the low-pass filtered spike train of neuron i (single neuron readout) with time constant . We then impose the following condition for spiking: a neuron shall emit a spike only if this decreases the loss function of its population in the immediate future. The condition for spiking also includes a noise term accounting for sources of stochasticity in spike generation⁵⁰, including the effect of unspecific inputs from the rest of the brain.

By assuming that each neuron emits a spike at time t only if this decreases the loss function of its population (Eq. 3), we derived the dynamics and network structure of a spiking network that instantiates efficient coding (Fig. 1C, see Methods). The derived dynamics of the subthreshold membrane potential and obey the equations of the generalized leaky integrate and fire (gLIF) neuron

where , and are synaptic current, spike-triggered adaptation current and non-specific external current, respectively, R_m is the membrane resistance and is the resting potential. This dynamics is complemented with a fire-and-reset rule: when the membrane potential reaches the firing threshold ϑ^y, a spike is fired and is set to the reset potential V^reset,y. The analytical solution in Eq. (5) holds for any number of neurons (with at least 1 neuron in each population) and predicts an optimal spike pattern to encode the presented external stimulus.

The synaptic currents in E neurons, , consist of feedforward currents, obtained as stimulus features s_k(t) weighted by the tuning weights of the neuron, and of recurrent inhibitory currents. Synaptic currents in I neurons, , consist of recurrent excitatory and inhibitory currents (Fig. 1C).

The optimization of the loss function also yields structured recurrent connectivity (Fig. 1D). The synaptic strength between two neurons is proportional to their tuning similarity if the tuning similarity is positive; otherwise the synaptic weight is set to zero (Fig. 1E,F) to ensure that Dale’s law is respected. This also sets the overall connection probability to 0.5. (For a study of how efficient coding would be implemented if the above Dale’s law constraint was removed and each neuron is free to have either an inhibitory or excitatory effect depending on the postsynaptic target, see Supplementary Text 1). Neurons with opposite tuning have low connection probability, consistent with experimental results^51,52,53 (Fig. 1D). Note that the structured recurrent connectivity leads to both E and I cells being stimulus-tuned, even though I cells do not receive feedforward inputs (Fig. 1C). The spike-triggered adaptation current of neuron i in population y, , has the dynamics of its low-pass filtered spike train . This current realizes spike-frequency adaptation or facilitation depending on the difference between the time constants of population and single neuron readout (see Results section “Weak spike-triggered adaptation optimizes network efficiency”). Finally, external currents have a constant mean, that depends on the parameter β, plus fluctuations that depend on the noise in the condition for spiking with intensity σ. Importantly, the relative weight of the metabolic cost over the encoding error controls the operating regime of the network biophysically plausibly, by modulating the mean of the external current. (Previous studies interpreted changes of the metabolic constant β as changes to the firing thresholds, which has less biophysical plausibility^36,33) (see section “State-dependent coding and dynamics are controlled by the metabolic cost on spiking”).

To summarize, the analytical derivation of an optimally efficient network includes gLIF neurons^{54,41,40,55,56}, a distributed code with mixed selectivity to the input stimuli, spike-triggered adaptation current, structured synaptic connectivity, and an operating regime controlled by the metabolic constant β.

The equations for the E-I network of gLIF neurons in Eq. (5) optimize the loss functions at any given time and for any set of parameters. In particular, the network equations have the same analytical form for any positive value of the metabolic constant β. To find a set of parameters that optimizes the overall performance, we defined a performance measure as the average over time and trials of the loss function. We then optimized the parameters by setting the metabolic constant β such that the encoding error weights 70 % and the metabolic error weights 30 % of the total performance, and by choosing all other parameters such as to minimize numerically our network performance measure (see Methods). The numerical optimization was performed by simulating a model of 400 E and 100 I units, a network size relevant for computations within one layer of a cortical microcolumn⁵⁷. The set of model parameters that optimized network efficiency is detailed in Table 1. Unless otherwise stated, in all simulations we will use the optimal parameters of Table 1 and only vary those parameters detailed in the figure axes.

With optimally efficient parameters, population readouts closely tracked the target signals (Fig. 1G, M=3, R² = [0.95, 0.97] for E and I neurons, respectively). When stimulated by our 3-dimensional time-varying feedforward input, the optimal E-I network provided a precise and unbiased estimator of the multi-dimensional and time-dependent target signal (Fig. 1H).

Next, we examined the emergent dynamical properties of an optimally efficient E-I network. The distribution of firing rates was well described by a log-normal distribution (Fig. 1I, left). Neurons fired irregularly, with mean coefficient of variation (CV) slightly smaller than 1 (Fig. 1I, right; CV= [0.97, 0.95] for E and I neurons, respectively). We assessed E-I balance in single neurons through two complementary measures. First, we calculated the average (global) balance of E-I currents by taking the time-average of the net sum of currents⁵⁸. Second, we evaluated the instantaneous⁵⁹ (also termed detailed⁴⁵) E-I balance using the Pearson correlation (ρ) of E and I currents received by a single neuron over time (see Methods).

We observed a strong average E-I balance (indicated by a net sum of synaptic inputs close to zero, with only a weak residual of inhibition in both E and I cells (Fig. 1J). Furthermore, we found a moderate instantaneous balance, stronger in I compared to E cell type (Fig. 1K-L, ρ = [0.44, 0.25], for I and E neurons, respectively). The presence of instantaneous balance between E and I synaptic currents within single neurons has been reported in cortical data^59,60.

Competition across neurons with similar stimulus tuning emerging in efficient spiking networks

We next explored coding properties emerging from recurrent synaptic interactions between E and I populations in the optimally efficient networks.

An approach that has recently provided empirical insight into local recurrent interactions between neurons is measuring the effective connectivity with cellular resolution, by photostimulating individual neurons and measuring the effect of such perturbation on other neurons in the network. Recent effective connectivity experiments photostimulated single E neurons in primary visual cortex and measured its effect on neighbouring neurons, finding that the photostimulation of an E neuron led to a decrease in firing rate of similarly tuned close-by neurons⁶¹. This effective lateral inhibition²⁶ between E neurons with similar tuning to the stimulus implements competition between neurons for the representation of stimulus features (termed feature-specific competition⁶¹).

To assess how E-I interactions shape coding in efficient networks, we simulated photostimulation experiments in these networks. We performed such experiments in the absence of the feedforward input to insure all effects are only due to the recurrent processing and not to feedforward processing. We stimulated a randomly selected single “target” E neuron and measured the change in the instantaneous firing rate from the baseline firing rate, Δz_i(t), in all the other I and E neurons (Fig. 2A, left). The photo-stimulation was modeled as an application of a constant depolarising current with a strength parameter, a_p, proportional to the distance between the resting potential and the firing threshold (a_p = 0 means no stimulation, while a_p = 1 indicates photostimulation at the firing threshold). We quantified the effect of the simulated photostimulation of a target E neuron on other E and I neurons, distinguishing neurons with either similar or different tuning with respect to the target neuron (Fig. 2A, right; Supplementary Fig. S2).

Figure 2.

The photostimulation of the target E neuron increased the instantaneous firing rate of similarly-tuned I neurons and reduced that of other similarly-tuned E neurons (Fig. 2B, Supplementary Fig. S2). We quantified the effective connectivity as the difference between the time-averaged firing rate of the recorded cell in presence or absence of the photostimulation of the targeted cell, measured during perturbation and up to 50 ms after. We found positive effective connectivity on I and negative effective connectivity on E neurons with similar tuning to the stimulated neuron, with a positive correlation between tuning similarity and effective connectivity on I neurons and a negative correlation on E neurons (Fig. 2C). As we varied the strength of the photostimulation, the firing rate of the target neuron increased proportionally to the photostimulation strength, as did the effect of the perturbation on I and E neurons with similar tuning to the target neuron (Fig. 2D, Supplementary Fig. S2). As we varied the time window of photostimulation, we found that the effective connectivity converges within a time window of about 300 ms (Fig. 2E). We confirmed these effects of photostimulation in presence of a weak feedforward input (Supplementary Fig. S2), similar to the experiments of Ref⁶¹ in which photostimulation was applied during the presentation of visual stimuli with weak contrast.

In summary, lateral excitation of I neurons and lateral inhibition of E neurons with similar tuning is an emerging coding property of the efficient E-I network. Lateral excitation and inhibition leads to competition between neurons with similar tuning to stimulus features, comparable to that found in the visual cortex^61,62. An intuitive summary of how this mechanism is implemented is that the E neuron that fires first activates I neurons with similar tuning. In turn, these I neurons inhibit all similarly tuned E neurons (Fig. 2A, right), preventing them to generate redundant spikes and encoding the sensory information that has already been encoded by the first spike. Suppression of redundant spiking allows efficient coding because it reduces the metabolic cost without compromising on encoded information³⁶.

To explore further the consequences of E-I interactions for stimulus encoding, we next investigated the dynamics of lateral inhibition in a network driven by the feed-forward sensory input but without perturbing neurons. In this case, shared feedforward inputs s_k(t) create a particular pattern of voltage correlations in E-E neuronal pairs, where voltage correlations linearly depend on the tuning similarity (Fig. 2F, left). The feedforward inputs are shared across neurons and weighted by the tuning parameters of E neurons. For this reason, they cause strong positive voltage correlations between E-E neuronal pairs with very similar tuning and strong negative correlations between pairs with very different (opposite) tuning (Fig. 2F, top-left). Voltage correlations between E-E pairs vanished regardless of tuning similarity when we made the inputs independent across neurons (Fig. 2F, top-middle), showing the relation between tuning similarity and voltage correlation occurs because of shared feedforward inputs. In contrast to E neurons, I neurons do not receive feedforward inputs and are driven only by similarly tuned E neurons (Fig. 2A, right). This causes positive voltage correlations in I-I neuronal pairs with similar tuning and vanishing correlations in neurons with different tuning (Fig. 2F, bottom-left). Such dependence of voltage correlations on tuning similarity disappears when removing the structure from the E-I synaptic connectivity (Fig. 2F, bottom-right).

Although membrane potentials could be strongly correlated or anti-correlated depending on tuning similarity (Fig. 2F, left), the coordination of spike timing of pairs of E neurons (measured with cross-correlograms or CCGs) was very weak (Fig. 2G-H). For I-I neuronal pairs, the peaks of CCGs were stronger than those observed in E-E pairs, but they were present only at very short lags (lags < 1 ms), and the same was true for E-I pairs. Additionally, noise correlations measured as Pearson correlation on spike counts in trials with the same stimulus (r_SC) had values distributed around zero (Fig. 2H). These findings lead to two conclusions. First, recurrent interactions of the efficient E-I network wipe away the effect of membrane potential correlations to produce largely uncorrelated spiking output, consistently with the efficient coding hypothesis of reducing redundancy in cases of low noise^3,6. Second, such precise cancelling of correlations between voltages and the spiking output reflects the millisecond precision of information processing in efficient E-I networks.

The effect of structured connectivity on coding efficiency and neural dynamics

The analytical solution of the optimally efficient E-I network predicts that recurrent synaptic weights are proportional to the tuning similarity between neurons. We here investigated the role of such efficient connectivity structure by comparing the behavior of an efficiently structured network with a similar but randomly structured E-I network of the type studied in previous works^63,64,23. We removed the connectivity structure by randomly permuting synaptic weights across neuronal pairs. We either randomized connections within a single connectivity type (E-I, I-I or I-E) or within all these three connectivity types at once (“all”). Such procedure destroys the relationship between tuning similarity and synaptic strength as in Fig. 1F while it preserves Dale’s law and the overall distribution of connectivity weights. We found that randomizing the connectivity structure significantly altered neural dynamics and coding (Fig. 3A-H). The structure in E-I and in I-E connectivity has a major effect on efficient coding. Randomizing E-I and I-E connectivity led to several-fold increases in the encoding error as well as to significant increases in the metabolic cost (Fig. 3A-B). In particular, with unstructured E-I connectivity the network failed completely to encode the target with I population (Fig. 3C).

Figure 3.

Unstructured E-I and I-E connectivity also yielded an increase of the variance in the membrane potentials (Fig. 3D) and firing rate in E neurons (Fig. 3E), while pulling the average net synaptic inputs towards inhibition (Fig. 3F) and removing the instantaneous balance (Fig. 3G). Together, these findings suggest a shift from mean-driven to fluctuation-driven spiking activity as the connectivity structure is removed. The structure of E-I connectivity was also found to be crucial for the linear relation between voltage correlations and tuning similarity in pairs of I neurons (Fig. 3H, magenta). Interestingly, we found no effect of connectivity structure on the variability of spiking of single neurons, with both structured and unstructured networks showing strong variability (Supplementary Fig. S3), suggesting that the variability of spiking is independent of the connectivity structure.

Randomizing I-I connectivity was less detrimental to the coding efficiency as it led to a slightly higher encoding error, but to a lower metabolic cost, and still allowed for a relatively good tracking of target signals in both cell types (Fig. 3C, “permuted I to I”). Contrary to randomization of the E-I and I-E connectivity, shuffling I-I connectivity decreased the variance of the membrane potential, decreased the firing rate in E neurons and increased instantaneous balance in E neurons. Thus it had opposite effects compared to shuffling of E-I and I-E connectivity. To understand if there was a minimal connectivity structure necessary for efficient coding, we also removed the connectivity structure only partially, keeping like-to-like connectivity structure and removing all structure beyond like-to-like. This manipulation only had very modest effects on network’s coding and almost no effect on neural dynamics (Supplementary Fig. S3), thus showing that like-to-like structure of connectivity is largely sufficient to achieve efficient coding.

Finally, we analyzed how the structure in recurrent connectivity influences lateral inhibition that we observed in efficient (structured) networks (see Fig. 2A-E). We found that the dependence of lateral inhibition on tuning similarity vanish when the connectivity structure is fully removed (Fig. 3I, right), thus showing that connectivity structure is necessary for lateral inhibition. While networks with unstructured E-I and I-E connectivity still show inhibition in E neurons upon single neuron optostimulation (because of the net inhibitory effect of recurrent connectivity; Supplementary Fig. S4), this inhibition was largely unspecific to tuning similarity. Unstructured connectivity decreased the correlation between tuning similarity and effective connectivity from r = [0.31, − 0.54] in E and I neurons in a structured network to r = [0.02, − 0.13] and r = [0.57, 0.11] in networks with unstructured E-I and I-E connectivity, respectively (Fig. 3I, first and third from the left). Removing the structure in I-I connectivity, in contrast, increased the correlation between effective connectivity and tuning similarity in E neurons (r = [0.30, − 0.65], Fig. 3I, second from the left), showing that lateral inhibition takes place irrespective of the I-I connectivity structure. Furthermore, a partial removal of connectivity structure where we only removed the connectivity structure beyond like-to-like had smaller effects on lateral inhibition (Supplementary Fig. S4), thus confirming that like-to-like connectivity pattern is sufficient for lateral excitation/inhibition in I and E neurons.

While optimally structured connectivity predicted by efficient coding is biologically plausible, it may be difficult to realise it exactly on a synapse-by-synapse basis in biological networks. We verified the robustness of the model to small deviations from the optimal synaptic weights by adding a random jitter, proportional to the synaptic strength, to all synaptic connections (see Methods). The encoding performance and neural dynamics were barely affected by such perturbation, demonstrating that the network is robust against random perturbations of the optimal synaptic weights (Supplementary Fig. S3).

In summary, we found that some aspects of recurrent connectivity structure, such as the like-to-like organization of E-I and I-E connectivity, are crucial to achieve efficient coding. Instead, for other aspects there is considerable flexibility; the organization of I-I connectivity is less crucial, as is the connectivity structure beyond like-to-like, and adding small perturbations to optimal weights has only minor effects. Structured E-I and I-E, but not I-I connectivity, is necessary for a robust dependence of lateral inhibition on tuning similarity.

Weak spike-triggered adaptation optimizes network efficiency

We next investigated the role of spike-triggered adaptation current,, that emerges from the optimally efficient solution (Eq. 5). This current provides a within-neuron feedback triggered by each spike, with time constant equal to that of the single neuron readout (E neurons) and (I neurons). The strength of the current is proportional to the difference in inverse time constants of single neuron and population readouts, , and it is thus absent in previous studies assuming that these time constants are equal^{29,28,33,31,44,42}.

Depending on the sign of the difference of time constants, this spike-triggered current is negative, giving spike-triggered adaptation³⁹, if the single-neuron readout has longer time constant than the population readout , or positive, giving spike-triggered facilitation, if the opposite is true (Table 2). We expected that network efficiency would benefit from spike-triggered adaptation (in short, adaptation), because accurate encoding requires fast temporal dynamics of the population readouts, to capture fast fluctuations in the target signal, while we expect a slower dynamics in the readout of single neuron’s firing frequency,, a process that could be related to homeostatic regulation of single neuron’s firing rate^65,66. Measuring performance of a simulated E-I network, we indeed found that optimal coding efficiency is achieved with weak adaptation in both cell types, and in particular in regimes where the adaptation is stronger in E compared to I neurons (Fig. 4A). We note that adaptation in E neurons promotes efficient coding because it enforces every spike to be error-correcting, while a spike-triggered facilitation in E neurons would lead to additional spikes that might be redundant and reduce network efficiency. Contrary to previously proposed model of adaptation in LIF neurons³⁸, here strength and the time constant of adaptation are not independent, but they both depend on , with larger yielding both longer and stronger adaptation.

Table 2.

Figure 4.

To gain further insights on how adaptation influences network performance, we set the adaptation in one cell type to 0 and vary the strength of adaptation in the other cell type by varying the time constant of the single neuron readout. In the absence of adaptation in I neurons , adaptation in E neurons resulted in an increase of the encoding error in E neurons and a decrease in I neurons (Fig. 4B, top). Conversely, adaptation in I neurons (with no adaptation in E neurons) was harmful for the efficiency of the model, as it led to an increase in the encoding error in both cell types (Fig. 4B, bottom).

Firing rates and variability of spiking were sensitive to the strength of adaptation. As expected, adaptation in E neurons caused a decrease in the firing levels in both cell types (Fig. 4C). In contrast, adaptation in I neurons decreased the firing rate in I neurons, but increased the firing rate in E neurons, due to a decrease in the level of inhibition. Furthermore, adaptation decreased the variability of spiking, in particular in the cell type with strong adaptation (Fig. 4D), a well-known effect of spike-triggered adaptation in single neurons⁶⁷.

Instantaneous balance of synaptic currents predicts network efficiency better than the average E-I balance

Next, we tested the capability of instantaneous and average E-I balance to predict the efficiency of the network. Measuring average balance and instantaneous balance of synaptic inputs from electrophysiology recordings is possible^59,60,58, while measuring efficiency from empirical data is challenging. The estimation of network efficiency requires the comparison between typically unknown network’s target representations and the population readouts. The estimation of the population readout, in turn, requires an estimation of decoding weights and the knowledge of spiking dynamics from a complete neural network.

We focused the analysis on regimes with adaptation, because these regimes gave better performance. In regimes with adaptation, time constants of single neuron readout influenced the average imbalance (Fig. 4E) as well as the instantaneous balance (Fig. 4F) in E and I cell type. The average balance was precise (with the net synaptic current close to 0) with strong adaptation in E neurons, and it got weaker when increasing the adaptation in I neurons (Fig. 4E). However, regimes with precise average balance in both cell types coincided with suboptimal efficiency (compare Fig. 4A, right and E).

To test how well the average imbalance and the instantaneous balance of synaptic inputs predict network efficiency, we concatenated the column-vectors of the measured average loss and of the average imbalance in each cell type and computed the Pearson correlation between these quantities. The correlation between the average imbalance and the average loss was weak in the E cell type (R = 0.16) and close to zero in the I cell type (R = 0.02), suggesting almost no relation between efficiency and average imbalance in the E cell type. In contrast, the average loss was negatively correlated with the instantaneous balance in both E (R = − 0.35) and in I cell type (R = − 0.45), showing that instantaneous balance of synaptic inputs is positively correlated with network efficiency.

When measured for varying levels of spike-triggered adaptation, unlike the average balance of synaptic inputs, the instantaneous balance is therefore a reliable predictor of network efficiency.

State-dependent coding and dynamics are controlled by the metabolic cost on spiking

In our derivation of efficiency objectives, we obtained non-specific external current (in the following, non-specific current), described by the term and comprising mean and fluctuations (see Methods). Non-specific current captures the ensemble of all synaptic currents that are unrelated and un-specific with respect to the stimulus features. This non-specific term collates effects of synaptic currents from neurons untuned to the stimulus^68,69, as well as synaptic currents from other brain areas. This term can also be conceptualized as the “background” synaptic activity that is thought to provide a large fraction of all synaptic inputs to both E and I neurons in cortical networks⁷⁰, and which may modulate feedforward-driven responses by controlling how far is typically the membrane potential from the firing threshold⁷¹. Likewise, in our model, the external current does not directly convey information about the feedforward input features, but influences the operating regime of the network. The mean of the non-specific external currents is proportional to the metabolic constant β and its fluctuations reflect the noise that we assumed in the condition for spiking. Since β governs the trade-off between encoding error and metabolic cost (Eq. 3), higher values of β imply that more importance is assigned to the metabolic efficiency than to coding accuracy, yielding a reduction in firing rates. In the expression for the non-specific synaptic current, we found that the mean of the current is negatively proportional to the metabolic constant β (see Methods). The non-specific current is typically depolarizing, meaning that increasing β yields a weaker non-specific current and increases the distance between mean membrane potential and the firing threshold. Thus, an increase of the metabolic constant is expected to create a network state that is less responsive to the feedforward signal.

We found the metabolic constant β to significantly influence the spiking dynamics (Fig. 5A). The optimal efficiency was achieved for non-zero levels of the metabolic constant (Fig. 5B). The metabolic constant modulated the firing rate as expected, with the firing rate decreasing with the increasing of the metabolic constant (Fig. 5C, top). It also modulated the variability of spiking, as increasing the metabolic constant decreased the variability of spiking in single neurons (Fig. 5C, bottom). Furthermore, it modulated the average imbalance and the instantaneous balance in opposite ways: larger values of β led to regimes that had stronger average balance, but weaker instantaneous balance (Fig. 5D). We note that, even with suboptimal values of the metabolic constant, the neural dynamics remained within biologically relevant ranges.

Figure 5.

The fluctuation part of the non-specific current, modulated by the noise intensity σ, that we added in the definition of spiking rule for biological plausibility (see Methods), strongly affected the neural dynamics as well (Fig. 5E). The optimal performance was achieved with non-vanishing noise levels (Fig. 5F) and the beneficial effect of the noise in the non-specific current arose from its impact on the instantaneous E-I balance. While the average firing rate of both cell types, as well as the variability of spiking in E neurons, increased with noise variance (Fig. 5G), the average and instantaneous balance of synaptic currents exhibited a non-linear behavior as a function of noise variance (Fig. 5H). Due to decorrelation of membrane potentials by the noise, instantaneous balance decreased with increasing noise variance (Fig. 5H, bottom). Some level of noise in the non-specific inputs is therefore necessary to establish the optimal level of instantaneous E-I balance. Interestingly, single neurons manifest significant levels of spiking variability already in the absence of noise in the non-specific inputs (Fig. 5H, bottom), indicating that the recurrent network dynamics generates substantial variability even in absence of variability in the external current. Variability in absence of noise demonstrates the intrinsic chaotic behavior of the network⁷².

In summary, non-specific external currents derived in our optimal solution have a major effect on coding efficiency and on neural dynamics. The noise in the external current is particularly important to obtain optimal levels of the instantaneous E-I balance in I neurons.

Optimal ratio of E-I neuron numbers and of the mean I-I to E-I synaptic efficacy coincide with biophysical measurements

Next, we investigated how coding efficiency and neural dynamics depend on the ratio of the number of E and I neurons (N^E : N^I or E-I ratio) and on the relative synaptic strengths between E-I and I-I connections.

Efficiency objectives (Eq. 3) are based on population, rather than single-neuron activity. Our efficient E-I network thus realizes a computation of the target representation that is distributed across multiple neurons (Fig. 6A). We predict that, if number of neurons within the population decreases, neurons have to fire more spikes to achieve an optimal population readout because the task of tracking the target signal is distributed among fewer neurons. To test this prediction, we varied the number of I neurons while keeping the number of E neurons constant. As predicted, a decrease of the number of I neurons (and thus an increase in the ratio of the number of E to I neurons) caused a linear increase in the firing rate of I neurons, while the firing rate of E neurons stayed constant (Fig. 6B, top). However, the variability of spiking and the average synaptic inputs remained relatively constant in both cell types as we varied these ratios (Fig. 6B, bottom, C), indicating a compensation for the change in the ratio of E-I neuron numbers through adjustment in the firing rates. These results are consistent with the observation in neuronal cultures of a linear change in the rate of postsynaptic events but unchanged postsynaptic current in either E and I neurons for variations in the E-I neuron number ratio⁷³.

Figure 6.

The ratio of the number of E to I neurons had a significant influence on coding efficiency. We found a unique minimum of the encoding error of each cell type, while the metabolic cost increased linearly with the ratio of the number of E and I neurons (Fig. 6D). We found the optimal ratio of E to I neuron numbers to be in range observed experimentally in cortical circuits (Fig. 6D, bottom, black arrow, N^E : N^I = 3.75 : 1;⁷⁴). Due to the linear increase of the cost with the ratio of the number of E and I neurons (Fig. 6D, bottom, green), strong weighting of the error predicted higher ratios (Fig. 6E, bottom). Also the encoding error (RMSE) alone, without considering the metabolic cost, predicted optimal ratio of the number of E to I neurons within a plausible physiological range, N^E : N^I = [3.75 : 1, 5.25 : 1], with stronger weightings of the encoding error by I neurons predicting higher ratios (Fig. 6E, top).

Next, we investigated the impact of the strength of excitatory and inhibitory synaptic efficacy (EPSPs and IPSPs). In our model, the mean synaptic efficacy is fully determined by the distribution of tuning parameters (see Methods). As evident from the expression for the population readouts (Eq. 2), the amplitude of tuning parameters (which are also decoding weights) determines the amplitude of jumps of the population readout caused by spikes (Fig. 6F). The stronger the amplitude of these weights, the larger is the average impact of spikes on the population signals.

We parametrized the distribution of decoding weights as a normal distributions centered at zero, but allowed the standard deviation (SD) of distributions relative to E and I neurons ( and ) to vary across E and I cell type. With such parametrization, we were able to analytically evaluate the mean E-I, I-I and I-E synaptic efficacy (see Methods). We found that in the optimally efficient network, the mean E-I and I-E synaptic efficacy is exactly balanced.

We next searched for the optimal ratio of the mean I-I to E-I efficacy as the parameter setting that maximizes network efficiency. Network efficiency was maximized when such ratio was about 3 to 1 (Fig. 6G). Our results predict the maximum E-I and I-E synaptic efficacy, averaged across neuronal pairs, of 0.75 mV, and the maximal I-I efficacy of 2.25 mV, values that are consistent with empirical measurements in the primary sensory cortex^75,52,53.

Similarly to the ratio of E-I neuron numbers, a change in the ratio of mean E-I to I-E synaptic efficacy was compensated for by a change in firing rates, with stronger I-I synapses leading to a decrease in the firing rate of I neurons (Fig. 6H). Conversely, weakening the E-I and I-E synapses resulted in an increase in the firing rate in E neurons (Supplementary Fig. S5). This is easily understood by considering that weakening the E-I and I-E synapses activates less strongly the lateral inhibition in E neurons (Fig. 2) and thus leads to an increase in the firing rate of E neurons. We also found that single neuron variability remained almost unchanged when varying the ratio of mean I-I to E-I efficacy (Fig. 6H, bottom) and the optimal ratio corresponded with previously found optimal levels of average and instantaneous balance of synaptic inputs (Fig. 6I). The instantaneous E-I balance monotonically decreased with increasing ratio of I-I to E-I efficacy (Fig. 6I, bottom, Supplementary Fig. S5).

In summary, our analysis suggests that optimal coding efficiency is achieved with four times more E neurons than I neurons and with mean I-I synaptic efficacy about 3 times stronger than the E-I and I-E efficacy. The optimal network has less I than E neurons, but the impact of spikes of I neurons on the population readout is stronger, also suggesting that spikes of I neurons convey more information.

Dependence of efficient coding and neural dynamics on the timescales and dimensionality of the stimulus

We finally investigated how the network’s behavior depends on the timescales and dimensionality of the input stimulus features. We manipulated the stimulus timescales by changing the time constant of the Ornstein-Uhlenbeck (O-U) process. The network efficiently encoded stimulus features when their time constants varied between 1 and 200 ms, with stable encoding error, metabolic cost (Fig. 7A) and neural dynamics (Supplementary Fig. S6).

Figure 7.

Finally, we tested how the network’s behavior changed when we varied the number of stimulus features M processed by the network. The encoding error of E (RMSE^E) and I neurons (RMSE^I) had a minimum at 3 and 4 stimulus features, respectively (Fig.7B, top), while the metabolic cost increased monotonically with the number of features (Fig.7B, bottom). The number of features that optimized network efficiency (the average loss) ranged between M = [1, 4]. With strong weighting of the error (g_L ≥ 0.89), the optimal number of features was M = 4, and with strong weighting of the cost, (g_L < 0.27), the optimal number of features was M = 1. It is intriguing that the optimal encoding performance, when assuming the weighting for the error is stronger than for the cost, is achieved not for a single stimulus feature, but for 3 or 4 independent features. Increasing the number of features beyond the optimal number resulted in a monotonic increase in firing rates for both cell types and in a contrasting effect on average and instantaneous balance, as it increased the average E-I balance and weakened the instantaneous balance (Supplementary Fig. S6).

In sum, we found the optimal network efficiency in presence of several (3 or 4) stimulus features, and a surprising ability of the network to accurately encode stimuli on a wide range of timescales.

Advantages of E-I versus one cell type model architecture for coding efficiency and robustness to parameter variations

Neurons in the brain are either excitatory or inhibitory. To understand how differentiating E and I neurons benefits efficient coding, we compared the properties of our efficient E-I network with an efficient network with a single cell type (1CT). The 1CT model is a simplification of the E-I model (see Supplementary Text 1) and has been derived and analyzed in previous studies^{29,28,36,33,44,42}. We compared the average encoding error (RMSE), the average metabolic cost (MC), and the average loss (see Supplementary Text 2) of the E-I model against the one cell type (1CT) model. Compared to the 1CT model, the E-I model exhibited a higher encoding error and metabolic cost in the E population, but a lower encoding error and metabolic cost in the I population (Fig. 7C). The average loss of the E-I model was significantly smaller than that of the 1CT model when using the typical weighting of the error and the cost of g_L = 0.7 (Fig. 7D), as well as for the vast majority of other weightings (g_L ≤ 0.95; Supplementary Fig. S1).

We further compared the 1CT and E-I models in terms of the robustness of firing rates to changes in the metabolic constant. Consistently with previous studies^36,35, firing rates in the 1CT model were highly sensitive to variations in the metabolic constant (Fig. 7E, note the logarithmic scale on the y-axis), with a superexponential growth of the firing rate with the inverse of the metabolic constant in regimes with metabolic cost lower than optimal. This is in contrast to the E-I model, whose firing rates exhibited lower sensitivity to the metabolic constant, and never exceeded physiological limits (Fig. 5C). Because our E-I model does not incorporate a saturating input-output function as in³⁴ that would constrain the range of firing rates, the ability of the E-I model to maintain firing rates within biologically plausible limits emerges as a highly desirable dynamic property.

In summary, we found that the optimal E-I model is more efficient than the 1CT model. Beyond the performance of optimal models, the E-I model is advantageous with respect to the 1CT model also because it does not enter into states of physiologically unrealistic firing rates.

Discussion

We analyzed comprehensively the structural, dynamical and coding properties that emerge in networks of spiking neurons that implement optimally the principle of efficient coding. We demonstrated that efficient recurrent E-I networks form highly accurate and unbiased representations of stimulus features with biologically plausible parameters, biologically plausible neural dynamics, instantaneous E-I balance and like-to-like lateral inhibition. The network can implement efficient coding with stimulus features varying over a wide range of timescales and when encoding even multiple such features. Here we discussed the implications of these findings.

By a systematic study of the model, we determined the model parameters that optimize network efficiency. Strikingly, the optimal parameters (including the ratio between the number of E and I neurons, the ratio of I-I to E-I synaptic efficacy and parameters of non-specific currents) were consistent with parameters measured empirically in cortical circuits, and generated plausible spiking dynamics. This result lends credibility to the hypothesis that cortical networks might be designed for efficient coding and may operate close to optimal efficiency, as well as provides a solid intuition about what specific parameter ranges (e.g. higher numbers of E and than I neurons) may be good for. Efficient networks still exhibited realistic dynamics and reasonably efficient coding in the presence of moderate deviations from the optimal parameters, suggesting that the optimal operational point of such networks is relatively robust. We also found that optimally efficient analytical solution derives generalized LIF (gLIF) equations for neuron models³⁷. While gLIF^67,40 and LIF^63,64 models are reasonably biologically plausible and are widely used to model and study spiking neural network dynamics, it was unclear how their parameters affect network-level information coding. Our study provides a principled way to determine uniquely the parameter values of gLIF networks that are optimal for efficient information encoding. Studying the dynamics of gLIF networks with such optimal parameters thus provides a direct link between optimal coding and neural dynamics. Moreover, our formalism provides a framework for the optimization of neural parameters that can in principles be used not only for neural network models that study brain function but also for the design of artificial neuromorphic circuits that perform information coding computations^76,77.

Unlike in previous randomly-connected recurrent networks of LIF and gLIF spiking neurons,^63,64 in our efficient-coding solution, a highly structured E-I, I-I and I-E synaptic connectivity emerges as an optimal structural solution to support efficient coding. Our model generates a number of insights about the role of structured connectivity in efficient information processing. A first insight is that I neurons develop stimulus feature selectivity because of the structured recurrent connectivity. This is in line with recent reports of stimulus feature selectivity of inhibitory neurons, including in primary visual cortex^78,79,80. A second insight is that a network with structured connectivity shows stronger average and instantaneous E-I balance, as well as significantly lower variance in membrane potentials compared to an equivalent network with the same connections organized randomly. This implies that the connectivity structure determines the operating regime of the network. In particular, a network structured as in our efficient coding solution operates in a dynamical regime that is more stimulus-driven, compared to an unstructured network that is more fluctuation driven. A third insight is that the structured network exhibits a several-fold lower encoding error compared to unstructured networks and achieves this precision with lower firing rates. Network with structured recurrent connectivity creates more precise representations with less spikes and is therefore significantly more efficient compared to unstructured networks. Our analysis of the effective connectivity created by the efficient connectivity structure shows that this structure sharpens stimulus representations, reduces redundancy and increases metabolic efficiency by implementing feature-specific competition, that is a negative effective connectivity between E neurons with similar stimulus tuning, as proposed by recent theories³⁰ and experiments ^61,62 of computations in visual cortex.

Our perturbation experiments on single E neurons predict a negative like-to-like effective connectivity between E neurons with similar tuning, as found experimentally in the mouse primary visual cortex with 2-photon optogenetic perturbations of E neurons^61,62. This suggests that the effective connectivity found in mouse visual cortex could reflect efficient coding in visual cortex. Comparing effective connectivity in models and experiments is also useful for ruling in and out different theories of how efficient coding may be implemented in primary visual cortex. Earlier theories^4,11 found evidence for efficient coding in visual cortex and proposed that such efficient computations relied only on feedforward connectivity; thus they predicted null effective connectivity between visual neurons and were ruled out by the empirical effective connectivity measures^61,62. Our model, instead, implements efficient coding with recurrent interactions, suggesting a mechanism that is compatible with these empirical measures. Importantly, we made predictions for further optogenetics experiments that could better constraints models of visual cortical efficient coding. Previous studies⁶¹ optogenetically stimulated E neurons but did not determine whether the recorded neurons where excitatory or inhibitory. Our model predicts that stimulation of E neurons would increase firing in similarly tuned I neurons and decrease firing in similarly tuned E neurons. Our analysis confirms earlier model predictions⁸¹ that like-to-like connectivity between E and I neurons is necessary for lateral inhibition and competition between E neurons. Beyond like-to-like connectivity, our model predicts an optimally efficient connectivity where synaptic strength positively correlates with pair-wise tuning similarity, a connectivity pattern that was recently observed experimentally ⁸².

Our study determines how structured E-I connectivity affects the dynamics of E-I balancing and how this relates to information coding. Previous work³² proposed that the E-I balance in efficient spiking networks operates on a finer time scale than in classical balanced E-I networks with random connectivity⁶⁴. However, a theory to determine the exact levels of instantaneous E-I balance that is optimal for coding was lacking. Consistent with the general idea put forth in^32,31,48, we here showed that moderate levels of E-I balance are optimal for coding, and that too strong levels of instantaneous E-I balance are detrimental to coding efficiency. Our results predict that like-to-like structured E-I-E connectivity is necessary for optimal levels of temporal E-I balance. Finally, the E-I-E structured connectivity that we derived supports optimal levels of instantaneous E-I balance and causes desynchronization of the spiking output. Such intrinsically generated desynchronization is a desirable network property that in previously proposed models could only be achieved by the less plausible addition of strong noise to each neuron^31,35.

We found that our efficient network, optimizing the representation of a leaky integration of stimulus features, does not require recurrent E-E connections. Supporting this prediction, recurrent E-E connections were reported to be sparse in primary visual cortex⁸³), and the majority of E-E synapses in the visual cortex were suggested to be long-range⁸⁴. However, future studies could address the role of recurrent excitatory synapses, that were shown to emerge in efficient coding networks implementing computations beyond leaky integration such as linear mixing of features³⁷. Efficient networks with E-E connectivity show neural dynamics that goes well beyond the canonical case analyzed here and can potentially describe persistent network dynamics⁴⁴. Such networks would also allow to address whether biologically plausible efficient networks exhibit criticality, as suggested by⁸⁵. Finally, we note that efficient encoding might be the primary normative objective in sensory areas, while areas supporting high-level cognitive tasks such as decision-making might include other computational objectives such as efficient transmission of information downstream to generate reliable behavioral outputs^86,87,88,25.

Acknowledgements

V.K. and T.S. thank Tatiana Engel for her contribution to the discussion of results and for her comments on an earlier version of the manuscript. This project was supported by funding from Technische Universität Berlin (“Equal Opportunity Program” to VK), by Internal Research Funding of Technische Universität Berlin (to TS), by NIH Brain Initiative (grants U19 NS107464, R01 NS109961, R01 NS108410 to SP), and the Simons Foundation for Autism Research Initiative (SFARI; grant 982347 to SP).

Code availability

The complete computer code for reproducing the results is available as a Github repository [will be shared upon acceptance].

Methods

Overview of the current approach and of differences with previous approaches

In the following, we present a detailed derivation of the E-I spiking network implementing the efficient coding principle. The analytical derivation is based on previous works on efficient coding with spikes^28,36, and in particular on our recent work³⁷. While these previous works analytically derived feedforward and recurrent transmembrane currents in leaky integrate-and fire neuron models, these models did not contain any synaptic current unrelated to feedforward and recurrent processing. Non-specific synaptic current was suggested to be important for an accurate description of coding and dynamics in cortical networks⁷¹. In the model derivation that follows, we also derived non-specific external current from efficiency objectives.

As we mapped the efficient coding objective on biologically plausible neural implementations, we found that such implementations (with plausible biophysical parameters) requires a transmembrane current that is independent of feedforward and recurrent processing. We interpreted this current as non-specific external current (shortly, non-specific current), collating the ensemble of synaptic projections from other brain areas that are not directly involved in processing of feedforward stimulus features⁷⁰, as well as synaptic inputs from the local network from neurons that are not tuned to feedforward stimulus features⁶⁹. The mechanistic effect of the non-specific current is to regulate the distance to firing threshold, a role that is close to the notion of “background” synaptic activity in cortical neurons⁷¹.

Moreover, previous models on efficient coding did not thoroughly consider physical units of variables that were interpreted as biophysical quantities (such as membrane potentials, firing thresholds, etc.). As these biophysical variables were derived from computational variables (such as target signals and population readouts), it remained unclear how biophysical variables might acquire their physical units. Here, we assigned physical units to the computational variables and thus naturally endowed the model with physical units. The network developed here allows for a better compatibility of efficient spiking models with neurobiology compared to previous works on efficient coding with spikes. With this model, we aim to describe neural dynamics and computation in early sensory cortices such as the primary visual cortex in rodents, even though many principles of the model developed here could be relevant throughout the brain.

Introducing variables of the model

We consider two types of neurons, excitatory neurons E and inhibitory neurons I. We denote as N^E and N^I the number of E-cells and I-cells, respectively. The spike train of neuron i of type y ∈ {E, I}, i = 1, 2, …, N^y, is defined as a sum of Dirac delta functions,

where is the time of the α-th spike of that neuron, defined as a time point at which the membrane potential of neuron i crosses the firing threshold.

We define the readout of the spiking activity of neuron i of type y (in the following, “single neuron readout”) as a leaky integration of its spike train,

with λ_r denoting the inverse time constant. This way, the quantity represents an estimate of the instantaneous firing rate of neuron i.

We denote as s_k(t), k = 1, 2, …, M the set of M dynamical features of the external stimulus (in the following, stimulus features) which are transmitted to the network through a feedforward sensory pathway. The stimulus features have the unit of the square root of millivolt, . The k-th dimension of the target signal is then obtained through a leaky integration of the feedforward variable, s_k(t)²⁹, with inverse time constant λ, as

Furthermore, we define a linear population readout of the spiking activity of E and I neurons

with in units of . Here, each neuron i of type y is associated with a vector of M tuning parameters representing the readout weight of neuron i with respect to the M population readouts in Eq. 9. These readout weights can be combined in the M × N^y matrix . The rows of this matrix define the patterns of readout weights for each signal dimension k = 1, …, M.

Loss functions

We assume that the activity of a population y ∈ {E, I} is set so as to minimize a time-dependent encoding error and a time-dependent metabolic cost:

with β^y > 0 in units of mV the Lagrange multiplier which controls the weight of the metabolic cost relative to the encoding error. The time-dependent encoding error is defined as the squared distance between the targets and their estimates, and the role of estimates is assigned to the population readouts . In E neurons, the targets are defined as the target signals x_k(t), and their estimators are the population readouts of the spiking activity of E neurons,. In I neurons, the targets are defined as the population readouts of E neurons and their estimators are the population readouts of I neurons . Furthermore, the time-dependent metabolic cost is proportional to the squared estimate of the instantaneous firing rate, summed across neurons from the same population. Following these assumptions, we define the variables of loss functions in Eq. 10 as

We use a quadratic metabolic cost because it promotes the distribution of spiking across neurons²⁸. In particular, the loss function of I neurons, L^I (t) implies the relevance of the approximation: (see ϵ^I in the Eq. 11), which will be used in what follows.

When shall a neuron spike?

We minimize the loss function by positing that neuron i of type y ∈ {E, I} emits a spike as soon as its spike decreases the loss function of its population y in the immediate future³⁷. We also define t⁻ and t⁺ as the left- and right-sided limits of a spike time , respectively. Thus, at the spike time, the following jump condition must hold:

with in units of mV. Here, the arguments t⁻ and t⁺ denote the left- and right-sided limits of the respected functions at time t. Furthermore, we added a noise term on the right-hand side of the Eq. (12) in order to consider the stochastic nature of spike generation in biological networks⁵⁰. A convenient choice for the noise is the Ornstein-Uhlenbeck process obeying

where is a Gaussian white noise with auto-covariance function ⟨η_i(t)η_j(t′)⟩ = δ_ijδ(t − t′). The process has zero mean and auto-covariance function , with the variance of the noise.

By applying the condition for spiking in Eq. (12) using y = E and y = I, respectively, we get

According to the definitions in Eqs. (7) and (9), if neuron i fires a spike at time , it causes a jump of its own filtered spike train (but not of other neurons j ≠ i), as well as of the population readout of the population it belongs to. Therefore, when neuron i fires a spike, we have for a given neuron j and a given population readout k:

By inserting Eq. (15a)-(15b) in Eq. (12), we find that neuron i of type y should fire a spike if the following condition holds:

These equations tell us when the neuron i of type E and I, respectively, emits a spike, and are similar to the ones derived in previous works^37,28. In addition to what has been found in these previous works, we here also find that each term on the left- and right-hand side in the Eq 16a has the physical units of millivolts.

We note that the expression derived from the minimization of the loss function of E neurons in the top row of Eq. (16a) is independent of the activity of I neurons, and would thus lead to the E population being unconnected with the I population. In order to derive a recurrently connected E-I network, the activity of E neurons must depend on the activity of I neurons. We impose this property by using the approximation of estimates that holds under the assumption of efficient coding in I neurons (see ϵ^I in the Eq. 11),. This yields the following conditions:

We now define new variables and as proportional to the left- and the right-hand side of these expressions,

The variables and are interpreted as the membrane potential and the firing threshold of neuron i of cell type y ∈ {E, I}.

Dynamic equations for the membrane potentials

In this section we develop the exact dynamic equations of the membrane potentials for y ∈ {E, I} according to the efficient coding assumption. It is practical to use the vector notation and rewrite variables in Eq. (17) as

with the squared length of the tuning vector of neuron i of type y. We also rewrite Eq. (8)-(9) in vector notation as

with x(t) := [x₁(t), …, x_M (t)]^⊤ the vector of M target signals, the vector of estimates of cell type y, and the vector of spike trains for N^y cell type y ∈ {E, I}.

In the case of E neurons, the time-derivative of the membrane potential in Eq. (18), is obtained as

By inserting the dynamic equations of the target signal , its estimate (Eq. 19) and of the single neuron readout (Eq. 7 in the case y = E), we get

where in the last line we used the definition of from the Eq. (18).

In the case of I neurons, the time derivative of the membrane potential in Eq. (18) is

By inserting the dynamic equations of the population readouts of E neurons and of the I neurons (Eq. 19) and of the single neuron readout (Eq. 7 in the case y = I), we get

where in the last line we used the definition of from Eq. (18).

Leaky integrate-and-fire neurons

The terms on the right-hand-side in Eqs. (21) and (23) can be interpreted as transmembrane currents. The last term in these equations, , can be interpreted as a current instantaneously resetting the membrane potential upon reaching the firing threshold²⁸. Indeed, when the membrane potential reaches the threshold, it triggers a spike and causes a jump of the membrane potential by an amount −β^y; this realizes resetting of the membrane potential which is equivalent to the resetting rule of integrate-and-fire neurons^54,56. Thus, by taking into account the resetting mechanism and defining the time constants of population and single neuron readout τ := λ⁻¹ and , we can rewrite Eqs. (21) and (23) as a leaky integrate-and-fire neuron model,

In the Eq. 24 we wrote explicitly the terms , which correspond to the synaptic projections of N^x presynaptic neurons of type x to the postsynaptic neuron i of type y, with the quantity denoting the synaptic weight. We note that, in the case of I neurons, the element with j = i describes an autapse, i.e., a projection of a neuron with itself; this term is equal to , and thus contributes to the resetting of the neuron i.

Imposing Dale’s principle on synaptic connectivity

We now examine the synaptic terms in Eq. (24). As a first remark, we see that synaptic weights depend on tuning parameters . For the sake of generality we drew tuning parameters from a normal distribution with vanishing mean, which yielded both positive and negative values of . This has the desirable consequence that a spike of a neuron with a positive tuning parameter in signal dimension k, pulls the estimate, , up, while a spike of a neuron with pulls the estimate down, allowing population readouts to track both positive and negative fluctuations of the target signal on a fast time scale.

Another consequence of synaptic connectivity in the Eq. (24) is that the synaptic weight between a presynaptic neuron j of type x and a postsynaptic neuron i of type y is symmetric and depends on the similarity of tuning vectors of the presynaptic and the postsynaptic neuron: . The sign of this scalar product is positive between neurons with similar tuning and negative between neurons with different tuning (and zero when the two tuning vectors are orthogonal). Thus, for a presynaptic neuron j of type x, the synaptic weights of its outgoing connections can be both positive and negative, because some of its postsynaptic neurons have similar tuning to the neuron j while others have different tuning. This is inconsistent with Dale’s principle⁸⁹, which postulates that a particular neuron can only have one type of effect on postsynaptic neurons (excitatory or inhibitory), but never both. To impose this constraint in our model, we set synaptic weights between neurons with different tuning to zero. To this end, we define the rectified connectivity matrices,

with x, y ∈{E, I} and [a]₊ ≡ max(0, a) a rectified linear function. This manipulation is also plausible from a biological point of view, because in the cortex, the connection probability of neurons with very different (e.g. opposite) tuning is typically close to 0⁵¹. Since the elements of the matrix J^yx are all non-negative, it is the sign in front of the synaptic term in the Eq. (24) that determines the sign of the synaptic current between neurons i and j. The synaptic current is excitatory if the sign is positive, and inhibitory if the sign is negative.

It is also interesting to note that rectification affects the rank of connectivity matrices. Without rectification, the product in Eq. (25) yields a connectivity matrix with rank smaller or equal to the number of input features to the network, M, similarly as in previous works^29,43,44. Since typically the number of input features is much smaller than the number of neurons, i.e., M << N^y, this would give a low-rank connectivity matrix. However, rectification in Eq. (25), necessary to ensure Dale’s principle in presence of positive and negative tuning parameters, typically results in a substantial increase of the rank of the connectivity matrix.

Using the synaptic connectivity defined in Eq. (25), we rewrite the network dynamics from Eq. (24) as:

These equations express the neural dynamics which minimizes the loss functions (Eq. (10)) in terms of a generalized leaky integrate-and-fire model with E and I cell types, and are consistent with Dale’s principle.

In principle, it is possible to use the same strategy as for the E-I network to enforce Dale’s principle in model with one cell type (introduced by²⁸). To do so, we constrained the recurrent connectivity of the model with a single cell type from³⁶ by keeping only connections between neurons with similar tuning vectors and setting other connections to 0 (see Supplementary text). This led to a network of only inhibitory neurons, a type of network model which is less relevant for the description of biological networks.

Model with resting potential and an external current

In the model given by the Eq. (26) the resting potential is equal to zero. In order to account for biophysical values of the resting potential and to introduce an implementation of the metabolic constant that is consistent with neurobiology, we add a constant value to the dynamical equations of the membrane potentials , the firing thresholds and the reset potentials . This does not change the spiking dynamics of the model, as what matters to correctly infer the efficient spiking times of neurons is the distance between the membrane potential and the threshold.

Furthermore, in the same equations, the role of the metabolic constant β^y as a biophysical quantity is questionable. The metabolic constant β^y is an important parameter that weights the metabolic cost over the encoding error in the objective functions (Eq. 10). On the level of computational objectives, the metabolic constant naturally controls firing rates, as it allows the network to fire more or less spikes to correct for a certain encoding error. A flexible control of the firing rates is a desirable property, as gives the possibility to potentially capture different operating regimes of efficient spiking networks³⁶. In the spiking model we developed thus far (Eq. 26), similarly to previous efficient spiking models^36,33, the metabolic constant β^y controls the firing threshold. In neurobiology, however, strong changes to the firing threshold that would reflect metabolic constraints of the network are not plausible. We thus searched for an implementation of the metabolic constant β^y that is consistent with neurobiology.

The condition for threshold crossing of the neuron i can be written by Eq. (26) as

with c an arbitrary constant in units of millivolts. In Eq. (27) we added a constant c/2 and a resting potential on the left- and right-hand side of the firing rule. Moreover, we shifted the noise and the dependency on the parameter β from the firing threshold to the membrane potential. Thus, we assumed that the firing threshold is independent of the metabolic constant and the noise, and we instead assumed the dependence on the metabolic constant and noise in the membrane potentials.

We now define new variables for y ∈ {E, I}:

and rewrite the model in Eq. 26 in these new variables

where and are the independent Gaussian white noise processes defined in Eq. (13) above. We note that all terms on the right-hand side of Eq. (29) have the desired units of mV. The model in Eq. (29) is an efficient E-I spiking network with improved compatibility with neurobiology. We have expressed two new terms in the membrane potentials of E and I neurons, one dependent on the metabolic constant β^y and one on the noise that we assumed in the condition for spiking (see Eq. 12). We will group these two terms to define an external current, a current that is well known in spiking models of neural dynamics⁴⁰.

Efficient generalized leaky integrate-and-fire neuron model

Finally, we rewrite the model from Eq. (29) in a compact form in terms of transmembrane currents, and discuss their biological interpretation. The efficient coding with spikes is realized by the following model for the neuron i of type y ∈ {E, I}:

with R_m the current resistance. The leak current,

with τ = R_mC_m and C_m the capacitance of the neural membrane⁵⁴, arose by assuming the same time constant for the target signals x_k and estimates and (see Eq. 19). We see that the passive membrane time constant τ = λ⁻¹ can be traced back to the time constant of the population read-out in Eq. (9). The synaptic currents are defined as

where we note the presence of a feedforward current to E neurons,

which consist in a linear combination of the stimulus features s(t) weighted by the readout weights . The stimulus features can be traced back to the definition of the target signals in Eq. (8). This current emerges in E neurons, as a consequence of having the target signals x_k(t) in the loss function of the E population (see Eqs. 10-11). I neurons do not receive the feedforward current because their loss function does not contain the target signal.

The current providing within-neuron feedback triggered by each spike,

was recently recovered³⁷. This current has the kinetics of the single neuron readout (i.e., low-pass filtered spike train). Its sign depends on the relation between the time constant of the population readout τ = λ⁻¹ and single neuron readout , because the metabolic constant β^y is non-negative by definition (Eq. 10). If the single neuron readout is slower than the population readout, , within-neuron feedback is negative, and can thus be interpreted as spike-triggered adaptation. On the contrary, if the single neuron readout is faster than the population readout, , the within-neuron feedback is positive and can thus be interpreted as spike-triggered facilitation. In a special case where the time constant of the single neuron and population readout are assumed to be equal, within-neuron feedback vanishes.

Finally, we here derived the non-specific external current:

that captures the ensemble of non-specific synaptic currents received by each single neuron. The non-specific current has a homogeneous mean across all neurons of the same cell type, and a neuron-specific fluctuation. The mean of the non-specific current can be traced back to the weighting of the metabolic cost over the encoding error in model objectives (Eq. 10), while the fluctuation can be traced back to the noise intensity that we assumed in the condition for spiking (Eq. 12). The non-specific external current might arise because of synaptic inputs from other brain areas than the brain area that delivers feedforward projections to the E-I network we consider here, or it might result from synaptic activity of neurons that are part of the local network, but are not tuned to the feedforward input⁶⁹.

We also recall the fast and slower time scales of single neuron activity:

and the connectivity matrices

The structure of synaptic connectivity is fully determined by the similarity of tuning vectors of the presynaptic and the postsynaptic neurons ( and , respectively), while the distribution of synaptic connectivity weights is fully determined by the distribution of tuning parameters .

Stimulus features

We define stimulus features as a set of k = 1, …, M independent Ornstein-Uhlenbeck processes with vanishing mean, standard deviation σ_s and the correlation time τ_s,

If not mentioned otherwise, we use the following parameters: σ_s = 2 (mV)^1/2 and τ_s = 10 ms. Variables η_k(t) are independent Gaussian white noise processes with zero mean and covariance function ⟨η_k(t)η_l(t′)⟩ = δ_klδ(t − t′). These variables should not be confused with the Gaussian white noises in Eq. (29).

Parametrization of synaptic connectivity

In the efficient E-I model, synaptic weights are parametrized by tuning parameters through Eq. (25). The total number of synapses in the E-I, I-I and I-E connectivity matrices (including silent synapses with zero synaptic weight) is n_syn = 2N^EN^I + (N^I)², while the number of tuning parameters is n_w = M (N^E + N^I). Because the number of stimulus features M is expected to be much smaller than the number of E or I neurons, the number of tuning parameters n_w is much smaller than the number of synapses n_syn.

We can achieve a further substantial decrease in the number of free parameters by using a parametric distribution of tuning parameters . We have set the tuning parameters following a normal distribution and found that excellent performance can be achieved with random draws of tuning parameters from the normal distribution, thus without searching for a specific set of tuning parameters. This drastically decreased the number of free parameters relative to synaptic weights to only a handful of parameters that determine the distributions of tuning parameters.

Given M features, we sample tuning parameters,, with i = 1, …, N^y, y ∈ {E, I}, as random points uniformly distributed on a M -dimensional sphere of radius . We obtain this by sampling, for each neuron, a vector of M i.i.d. standard Gaussian random variables,, with , and normalizing the vector such as to have length equal to ,

This ensures that the length of tuning vectors in Eq. (32) is homogeneous across neurons of the same cell type, i.e., . Parameters and determine the heterogeneity (spread) of tuning parameters.

By combining Eq. (25) and Eq. (32), we obtain the synaptic weights, , as a function of the angle, , between the tuning vectors of presynaptic neurons, , and postsynaptic neurons,,

In the M = 3 dimensional case, we have that the distribution of the angle between two vectors is , with . Thus, the average strength of synaptic weights between the pre- and the postsynaptic population can be calculated as

Thus, the upper bound for the synaptic weight between cell types x and y is simply

From the Eq. (34), we have that the mean E-I connectivity is equal to the mean I-E connectivity . As we consider the ratio of the mean connectivity between E-I and I-I neurons, we find that it is given by the following:

Performance measures

Average encoding error and average metabolic cost

The definition of the time-dependent loss functions (Eq. 10) induces a natural choice for the performance measure: the mean squared error (MSE) between the targets and their estimators for each cell type. In the case of the E population, the time-dependent encoding error is captured by the variable ϵ^E(t) in the Eq. (11) and in case of I population it is captured by ϵ^I (t) defined in the same equation. We used the root MSE (RMSE), a standard measure for the performance of an estimator⁴⁰. For the cell type y ∈ {E, I} in trial q, the RMSE is measured as

with ⟨z_q(t)⟩_t,q denoting the time- and trial-average.

Following the definition of the time-dependent metabolic cost in the loss functions (Eq. 10), we measured the average metabolic cost in a trial q for the cell type y ∈ {E, I} as

with time-dependent metabolic cost κ^y(t) as in model’s objectives (Eq. 11) and ⟨z_q(t)⟩ _t,q the time- and trial-average. The square root was taken to have the same scale as for the RMSE (see Eq. 37).

The bias of the estimator

The MSE can be decomposed into the bias and the variance of the estimator. The time-dependent bias of estimates, y ∈ {E, I}, were evaluated for each time point over q = 1, …, Q trials. The time-dependent bias in input dimension k = 1, …, M is defined as

with ⟨z_q(t)⟩ _t,q the trial-averaged realization at time t. To have an average measure of the encoding bias, we averaged the bias of estimators over time and over input dimensions:

The averaging over time and input dimensions is justified because s_k(t) are independent realizations of the Ornstein-Uhlenbeck process (see Eq.31) with vanishing mean and with the same time constant, and variance across input dimensions.

Criterion for determining optimal model parameters

The equations of the E-I spiking network in Eqs. 30a-30h (Methods), derived from the instantaneous loss functions, give efficient coding solutions valid for any set of parameter values. However, to choose parameters values in simulated data in a principled way, we performed a numerical optimization of the performance function detailed below. Numerical optimization gave the set of optimal parameters listed in Table 1. When testing the efficient E-I model with simulations, we used the optimal parameters in Table 1 and changed only the parameters plotted in the figure axes on a figure-by-figure basis.

To estimate the optimal set of parameters θ = θ^*, we performed a grid search on each parameter θ_i while keeping all other parameters fixed as specified in Table 1. While varying the parameters, we measured a weighted sum of the time- and trial-averaged encoding error and metabolic cost. For each cell type y ∈ {E, I}, we computed

with ⟨z_q(t)⟩ _t,q the average over time and over trials and with ϵ^y(t) and κ^y(t) as in model’s objectives (Eq. 11).

To optimize the performance measure, we used a value of g_L = 0.7. The parameter g_L in the Eq. (40a) regulates the relative importance of the average encoding error over the average metabolic cost. Since the performance measure in Eq. (40a) is closely related to the average over time and trials of the instantaneous loss function (Eq. 10) where the parameter β regulates the relative weight of instantaneous encoding error over the metabolic cost, setting g_L is effectively achieved by setting β.

The optimal parameter set θ = θ^* reported in Table 1 is the parameter set that minimizes the sum of losses across E and I cell type

For visualization of the behavior of the average metabolic cost (Eq. 38) and average loss (Eq. 40a) across a range of a specific parameter θ_i, we summed these measures across the E and I cell type and normalized them across the range of tested parameters.

The exact dynamic and performance of our model depends on the realizations of random variables which describe the the tuning parameters , the Gaussian noise in the non-specific currents , and the initial conditions of the membrane potential , that were randomly drawn from a normal distribution in each simulation trial. To capture the performance of a “typical” network, we iterated the performance measures across trials with different realizations of these random variables, and averaged the performance measures across trials. We used 100 simulation trials for each parameter search.

Functional activity measures

Tuning similarity

The pair-wise tuning similarity was measured as the cosine similarity⁹¹, defined as:

with the length of the tuning vector in Euclidean space and α the angle between the tuning vectors and .

Cross-correlograms of spike timing

The time-dependent coordination of spike timing was measured with the cross-correlogram (CCG) of spike trains, corrected for stimulus-driven coincident spiking. The raw cross-correlogram (CCG) for neuron i of cell type y and neuron j of cell type x was measured as follows:

with q = 1, …, Q simulation trials with identical stimulus and T the duration of the trial. We subtracted from the raw CCG the CCG of trial-invariant activity. To evaluate the trial-invariant cross-correlogram, we first computed the peri-stimulus time histogram (PSTH) for each neuron as follows:

The trial-invariant CCG was then evaluated as the cross-correlation function of PSTHs between neurons i and j,

Finally, the temporal coordination of spike timing was computed by subtracting the correction term from the raw CCG:

Average imbalance of synaptic inputs

We considered time and trial-averaged synaptic inputs to each E and I neuron i in trial q, evaluated as:

with synaptic currents to E neurons and to I neurons as in Eq. (30c). Synaptic inputs were measured in units of mV. We reported trial-averages of the net synaptic inputs from the Eq. (43).

Instantaneous balance of synaptic inputs

We measured the instantaneous balance of synaptic inputs as the Pearson correlation of time-dependent synaptic inputs incoming to the neuron i. For those synaptic inputs that are defined as weighted delta-spikes (for which the Pearson correlation is not well defined; see Eq. 30c), we convolved spikes with a synaptic filter ,

where we used the expression for the feedforward synaptic current from the Eq. (30d). Note that the feedforward synaptic current is already already low-pass filtered (see Eq. 31). Using synaptic inputs from the Eq. 44, we computed the Pearson correlation of synaptic inputs incoming to single E neurons, for i = 1, …, N^E, and to single I neurons, for i = 1, …, N^I. The coefficients were then averaged across trials.

Perturbation of connectivity

To test the robustness of the model to random perturbations of synaptic weights, we applied a random jitter to optimally efficient recurrent synaptic connectivity weights. The random jitter was proportional to the synaptic weight, , where σ_J is the strength of the perturbation and are independent standard normal random variables. All three recurrent connectivity matrices (E-I, I-I and I-E) were randomly perturbed at once.

Computer simulations

We run computer simulations with Matlab R2023b (Mathworks). The membrane equation for each neuron was integrated with Euler integration scheme with the time step of dt = 0.02 ms.

The simulation of the E-I network with 400 E units and 100 I units for an equivalent of 1 second of neural activity lasted approximately 1.65 seconds on a laptop.

Supplementary material

Supplementary text 1: Derivation of the one cell type model

An efficient spiking model network with one cell type (1CT) has been developed previously²⁸, and properties of the 1CT model where the computation is assumed to be the leaky integration of inputs has been addressed in a number of previous studies^{29,43,36,33,42}. Compared to the efficient E-I model, the 1CT model can be seen as a simplification, and can be treated similarly to the E-I model, which is what we demonstrate in this section.

As the name of the model suggests, all neurons in the 1CT model are of the same cell type, and we have i = 1, …, N such neurons. We can then use the definitions in Eqs. (6) - (9) (now without the index y) and a loss function similar to the one in³⁶, but with only one (quadratic) regularizer

with β₁ > 0. The encoding error of the one cell type model minimizes the squared distance between the target signal x_k(t) and the estimate . As we apply the condition for spiking as for the E-I network (Eq. 12 without the index y) and follow the same steps as for the E-I network, we get

with ξ_i(t) the noise at the condition for spiking. Same as in the E-I model, we define the noise as an Ornstein-Uhlenbeck process with zero mean, obeying

where η_i is a Gaussian white noise and λ = τ⁻¹ is the inverse time constant of the process. We now define proxies of the membrane potential and the firing threshold as

Differentiating the proxy of the membrane potential u_i(t) and rewriting the model as an integrate-and-fire neuron, we get

We now proceed in the same way as with the E-I model and define new variables

In these new variables, we can rewrite the membrane equation of the 1CT model as follows:

Finally, we rewrite the model with a more compact notation of a leaky integrate-and-fire neuron model with transmembrane currents,

with currents

Note that the model with one cell type does not obey Dale’s law, since the same neuron sends to its postsynaptic targets excitatory and inhibitory currents, depending on the tuning similarity of the presynaptic and the postsynaptic neuron w_i and w_j (Eq. S.8b). In particular, if the pre- and postsynaptic neurons have similar selectivity , the recurrent interaction is inhibitory, and if the neurons have different selectivity , the interaction is excitatory. Simply put, neurons with similar selectivity inhibit each other while neurons with different selectivity excite each other³⁶.

Dale’s law can be imposed to the 1CT model the same way as in the E-I model. To do so, we removed synaptic interactions between neurons with different selectivity by rectifying the connectivity matrix,

However, this manipulation results in a network with only inhibitory recurrent synaptic interactions, and thus a network of only inhibitory neurons. Network with only inhibitory interactions is less relevant for the description of recurrently connected biological networks.

Supplementary text 2: Analysis of the one cell type model and comparison with the E-I model

We re-derived the 1CT model as a simplification of the E-I network (Supplementary Text 1, Supplementary Fig. S1A-B), with objective function of the same form as L^E and by allowing a single type of neurons sending both excitatory and inhibitory synaptic currents to their post-synaptic targets (Supplementary Fig. S1C). Similarly to the E-I model, also the 1CT model exhibits structured connectivity, with synaptic strength depending on the tuning similarity between the presynaptic and the postsynaptic neuron. Pairs of neurons with stronger tuning similarity (dissimilarity) have stronger mutual inhibition (excitation); see Supplementary Fig. S1D.

We compared the coding performance of the E-I model with that of a fully connected 1CT model. Both models received the same set of stimulus features and performed the same computation. In the 1CT model, tuning parameters were drawn from the same distribution as used for the E neurons in the E-I model. We used the same membrane time constant τ in both models, while the metabolic constants (β of the E-I model and β₁ of the 1CT model) and the noise intensity (σ of the E-I model and σ₁ of the 1CT model) were chosen such as to optimize the average loss for each model (Fig. 5B for E-I model, Supplementary Fig. S1F-G for 1CT model). Parameters of the 1CT model are listed in the Supplementary Table S1. A qualitative comparison of the E-I and the 1CT model showed that with optimal parameters, both models accurately tracked multiple target signals (Fig. 1G and Supplementary Fig. S1E).

To compare the performance of the E-I and the 1CT models also quantitatively, we measured the average encoding error (RMSE), metabolic cost (MC) and loss of each model. The RMSE and the MC in the 1CT model were measured as in Eq. 37 and 38, while the average loss of each model was evaluated as follows:

Unless mentioned otherwise, we weighted stronger the encoding error compared to the metabolic cost and used g_L = 0.7. Note that our comparison of the losses is conservative, because the metabolic cost is defined as a sum of activities across neurons (Eq. 38) and the total number of neurons in the E-I model (N^E + N^I) is larger than the number of neurons in the 1CT model (N^1CT = N^E).

Table S1.

Supplementary Figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.