Integrator readout

The first readout scheme, the integrator readout (IR), first filters the summed activity of the readout neurons, which are a subset of the BCN. This set of readout neurons can be divided into three subsets ($$ S^{\mathrm{\text{RS}}}$$, $$ S^{\mathrm{\text{FS}}}$$, and $$ S^{\mathrm{\text{SOM}}}$$), which are a random selection from the RS, FS, and SOM populations, respectively (see Fig 2A). Unless otherwise indicated, the size of the three readout sets is $$ N_{\mathrm{\text{read}}}^{\mathrm{\text{RS}}}=1000$$, $$ N_{\mathrm{\text{read}}}^{\mathrm{\text{FS}}}=100$$, and $$ N_{\mathrm{\text{read}}}^{\mathrm{\text{SOM}}}=100$$, respectively. The spike trains emitted by all neurons within the readout sets are linearly filtered by using the following dynamical equation:

To compute false positive and correct detection rates, a single lower decision boundary θ₋ is used by the IR, as depicted in Fig 3A. A detection event is registered if A_ir(t) falls at least once below the boundary θ₋ within the detection time window (0, T_w = 600ms). In catch trials, no stimulus is present. These catch trials are used to determine the false positive rate:

Fig 3

A detection event is defined by a lower threshold crossing for the integrator readout (IR) and by an upper threshold crossing for the differentiator readout (DR) and the differentator network readout (DNR).

We show here an example for the IR (A) and the DR (B). A: The black trace represents one realization of the IR readout activity during one catch trial (i.e. in the absence of stimulus). The red trace represents one realization of the readout activity in the presence of a stimulus. B: Same for the DR. The readout activity of the IR and DR are defined in Eqs (1) and (4), respectively. The IR reacts whenever the readout activity falls at least once below the detection boundary θ₋ (dashed-dotted line) within the detection time window (delimited by vertical dotted lines). The DR and the DNR react whenever the readout activity exceeds at least once the detection boundary θ₊ (dashed-dotted line) within the detection time window (delimited by vertical dotted lines).

In Fig 3A, the red trace is a stimulus trial which falls below the chosen threshold. Hence, this trial is considered a correct detection. The black trace represents a catch trial. Because it does not cross the threshold within the detection window, no false positive event is registered in this trial. Note that the noise at the single-trial level is large.

Differentiator readout

The differentiator readout (DR) first reads in the input from the network in the same way as the IR and it takes the difference between A_ir evaluated at two times separated by a lag ΔT. The result is then convolved with an exponential filter $$ {ℱ}_{{\tau }_f}(t)=\text{exp}(-t/{\tau }_f)/{\tau }_f$$ to reduce the noise

Trajectories computed from Eq (4) are used in combination with a simple threshold detector, as in the previous case. Note that the DR, however, uses an upper detection threshold θ₊, as shown in Fig 3B. Black and red trials plotted in Fig 3B represent again a catch and stimulus trial, respectively. They both exceed the chosen threshold at least once within the detection window. Therefore, both trials contribute to false positive and correct detection rates, respectively. These rates are defined as in the previous case, except that here an upper threshold is used:

There is nothing profound in the choice to endow the IR with a lower detection boundary and the DR with an upper detection boundary. The essential reason is that, in this way, we obtain a positive effect size (defined below), as it was observed in the experiments. In other words, this combination works for the readout schemes we consider here.

Differentiator network readout

The operation performed by the DR, i.e. the subtraction of a delayed copy of the readout activity, can be approximated by the simple network shown in Fig 2C. The differentiator network readout (DNR) consists of two populations: one readout population of 10000 RS neurons ($$ S^B$$) and one population of 2000 FS inhibitory neurons ($$ I$$). Each neuron within both populations receives the same number of feed-forward input connections from the three populations of the BCN as the size of the readout set of the respective cell type. More precisely, each neuron in the DNR receives input from $$ N_{\mathrm{\text{read}}}^{\mathrm{\text{RS}}}=1000$$ randomly chosen RS neurons, $$ N_{\mathrm{\text{read}}}^{\mathrm{\text{FS}}}=100$$ randomly chosen FS neurons, and $$ N_{\mathrm{\text{read}}}^{\mathrm{\text{SOM}}}=100$$ randomly chosen SOM neurons. Neurons in the readout population $$ S^B$$ evolve according to the same dynamical equation as RS neurons of the BCN, while neurons in $$ I$$ follow the same dynamical equations as the FS neurons of the BCN (see Methods).

One way to obtain a network that approximates a differentiator is that the mean input to the readout population $$ S^B$$ via the feed-forward inhibitory pathway from the BCN via $$ I$$ to $$ S^B$$ should cancel the direct feed-forward excitatory pathway input at a later time. To this end, the average weight of connections from $$ I$$ to $$ S^B$$, $$ J_{ei}^R$$, should be chosen such that a static change in the input from the direct pathway Δμ_e is compensated by the static change in the input from the inhibitory pathway $$ \Delta {\mu }_I$$ (see Fig 4). The value of $$ J_{ei}^R$$ that approximately satisfies the condition $$ \Delta {\mu }_e+\Delta {\mu }_I\approx 0$$ is computed from a linear-response calculation, (see p. 28). Importantly, the inhibitory pathway is given an additional transmission delay ΔT = 10 ms, so that changes in the input from the BCN to the DNR are counterbalanced at a later time. As a matter of fact, choosing $$ J_{ei}^R$$ such that the cancellation is not perfect grants an even better agreement with the experimental data, as shown below. More details on the DNR can be found in the Methods along with all parameter values.

Fig 4

Tuning of the differentiator readout network to implement the operation of the differentiator readout scheme.

A perturbation in the firing rate of the BCN (Δr) causes a perturbation in the mean input to the RS readout neurons, Δμ_e, and a perturbation in the firing rate of the inhibitory readout population $$ I$$. This change in firing rate causes a shift in the input from $$ I$$ to $$ S^B$$ ($$ \Delta {\mu }_I$$). The strength of the connection from $$ I$$ to $$ S^B$$ ($$ J_{ei}^R$$) is adjusted such that $$ \Delta {\mu }_e+\Delta {\mu }_I\approx 0$$. This cancellation reaches $$ S^B$$ with a time lag ΔT so that the readout network roughly implements the operation of the DR Eq (4).

The DNR activity is obtained by filtering the average firing rate of the readout neurons $$ S^B$$ with the same exponential filter used for the DR:

The DNR uses an upper detection boundary as the DR. Hence, false positive and hit rates are obtained in exactly the same way as done for the DR and depicted in Fig 3B):

Effect of stimulus duration

The effect of changing the stimulus duration will be considered first. To this end, stimuli of length 100, 200, and 400 ms are used (Fig 7A). The stimulus intensity is kept constant at 25% of the maximum current. In the experiment, when the stimulated cell was a RS neuron, the three stimuli evoked 6 ± 3, 11 ± 5, and 23 ± 10 spikes in the targeted neuron, respectively (Fig 7B, squares with different colors as in panel A). In the model, the number of evoked spikes was 7 ± 1, 12 ± 2, and 20 ± 5 spikes, respectively. The average number of evoked spikes generated by the model is within one standard deviation of the experimental data. However, the spread of the spike count distribution is smaller in the model, which is not surprising, considering the multiple possible noise sources that are not modeled, and that only some of the cellular parameters are randomly distributed in the model.

Fig 7

Experimental data for stimuli of equal intensity and different duration are compatible with simulation results from DR and DNR, but not from IR.

A: stimulation currents used in this figure. The color coding (red: 100 ms, green: 200 ms, blue: 400 ms) applies to all panels. B: number of spikes evoked in the targeted neuron vs. stimulus duration in the model (black circles and dashed line) and in the experiments (squares colored as in A). First row (C,F,I): deviation from the spontaneous value of the time-dependent mean readout activity normalized by the spontaneous standard deviation, defined in Eq (13). Second row (D,G,J): deviation from the spontaneous value of the time-dependent standard deviation of the readout activity normalized by the spontaneous standard deviation, defined in Eq (16). Third row (E,H,K): average effect size for the three stimuli in A. Open circles with error bars: experimental data (average over 119 RS cells and 2407 total trials); filled circles of corresponding color: model simulations. First column (C,D,E): IR. Second column (F,G,H): DR. Third column (I,J,K): DNR. Black line: catch trial condition.

The remaining panels, Fig 7C–7K give an overview of $$ {\hat{\mu }}_X(t)$$, $$ {\hat{\sigma }}_X(t)$$, and the effect size measured by all three detectors with stimuli of different duration. In all plots, the color coding is as in Fig 7A and the black line represents catch trials. The various panels are organized as follows: the first, second, and third column represent results obtained from the IR, DR, and DNR, respectively; the first row shows the time-dependent trial average (expressed as standard score) $$ {\hat{\mu }}_X(t)$$, the second row shows $$ {\hat{\sigma }}_X(t)$$, and the third row displays the effect size defined as in Eq (12) (filled dots) superimposed with the experimental results (open circles with error bars).

The average response of the IR activity, $$ {\hat{\mu }}_{\text{ir}}(t)$$, is shown in Fig 7C for the three stimuli. It is a reduction that shows its deepest dip around t = 40 ms and then relaxes back owing to the spike-frequency adaptation of RS and SOM-LTS populations. Since a large part of the IR activity A_ir(t) is driven by the RS population and the further input from the SOM population has the same average postsynaptic effect (see Fig 5A and 5B and recall that the effect of the SOM population is inhibitory), it is not surprising that the time course of $$ {\hat{\mu }}_{\text{ir}}(t)$$ upon 400 ms stimulation (Fig 7C, blue line) closely resembles that of the RS population, i.e. the mirrored SOM response, shown in Fig 5. The disinhibitory input from the FS readout population does have an antagonist effect, but it is not sufficient to remove the average decrease in the IR activity. The different roles of the three readout populations will be discussed in more detail below (see p. 19).

In the case of the 100 ms stimulus (Fig 7C, red line), the mean deviation goes back to zero shortly after the stimulus is turned off, while $$ {\hat{\mu }}_{\text{ir}}(t)$$ for the other stimulus durations (green and blue line) settles around -0.3 for the remainder of the respective stimulation time window. The time-dependent standard deviation $$ {\hat{\sigma }}_{\text{ir}}(t)$$, shown in Fig 7D, displays a mild increase, especially in the later part of the stimuli. The time courses of $$ {\hat{\mu }}_{\text{ir}}(t)$$ and $$ {\hat{\sigma }}_{\text{ir}}(t)$$ suggest that the readout activity has more chances of reaching a lower detection threshold for longer lasting stimuli. Indeed, the effect size measured by the IR strongly depends on the stimulus duration, as can be seen in Fig 7E (filled circles), which is in contrast with the experimental data (open circles with error bars).

When the DR is used, the picture changes rather drastically. The time-dependent mean $$ {\hat{\mu }}_{\text{dr}}(t)$$ displays two peaks and one trough in response to all three signals (Fig 7F), even though the two peaks partly overlap in the case of the 100 ms stimulus (red line). Because the DR considers differences in the IR readout activity, each peak corresponds to an upswing of $$ {\hat{\mu }}_{\text{ir}}(t)$$ and the trough to the initial sharp drop. The most prominent feature in $$ {\hat{\sigma }}_{\text{dr}}(t)$$ is again a mild increase in the later part of the simulation time window (Fig 7G). The main difference in the response to the three stimuli is the position of the last peak. Hence, it stands to reason that the DR activity A_dr(t) has similar chances to reach the upper barrier regardless of the signal length. Indeed, the effect size measured by the DR displays only a weak dependence on the signal duration, due to the mild increase in $$ {\hat{\sigma }}_{\text{dr}}(t)$$ (Fig 7H).

The time course of $$ {\hat{\mu }}_{\text{dnr}}(t)$$ qualitatively resembles that of $$ {\hat{\mu }}_{\text{dr}}(t)$$ (see Fig 7I), thus suggesting that the DNR does approximately operate as a differentiator. The most evident difference between $$ {\hat{\mu }}_{\text{dnr}}(t)$$ and $$ {\hat{\mu }}_{\text{dr}}(t)$$ is the value of the plateau between the two peaks, which is slightly below the zero level. This negative response is due to the chosen tuning of the recurrent inhibition within the DNR and partly compensates the increase in the time-dependent standard deviation $$ {\hat{\sigma }}_{\text{dnr}}(t)$$, which behaves similarly to the mean, except that it is above the zero level in the central part of all stimuli (Fig 7J), and thus, it slightly increases the detection chance of longer stimuli. The combination of the two effects leads to an effect size that is barely dependent on the stimulus duration and that agrees rather well with the experimental results (Fig 7K).

To summarize the results of Fig 7, the effect size measured by the IR is larger but strongly depends on the duration of the stimulus. The DR and the DNR detect the stimulus with a reliability that is essentially independent of the signal duration and that is of the same magnitude as the experimental data.

Effect of stimulus intensity

In the second experiment, we vary the firing rate of the stimulated cell by changing the current intensity while keeping the total area under the step stimulus, i.e. the injected charge, constant. As depicted in Fig 8A, the stimuli lasting 100 ms (red), 200 ms (green), and 400 ms have an intensity corresponding to 100%, 50%, and 25% of the maximum current, respectively. In this way, the total number of elicited spikes is approximately the same for each stimulus. In the experiment, the three stimuli evoked a firing rate of (109 ± 52) Hz, (54 ± 23) Hz, and (30 ± 10) Hz, respectively. In the model, the average evoked rates are (150 ± 25) Hz, (103 ± 20) Hz (50 ± 12) Hz. Note that the maximum current in the model is chosen such that the number of elicited spikes roughly matches the data of the previous experiment (dependence on stimulus duration), which were based on a different set of cells. Still, the evoked firing rates in the model are in a similar range as the experimental values.

Fig 8

Only results from DR and DNR are consistent with the experimental data for the stimulation with stimuli of intensity inversely proportional to the duration (signals are shown in A).

Panels B-J are organized as panels C-K of Fig 7 with the same color coding. Open circles with error bars in D, G, and J are the experimental average effect size measured from 55 RS cells (1469 total trials).

The remainder of Fig 8 reports all detection statistics arranged in the same way as before. The shape of $$ {\hat{\mu }}_{\text{ir}}(t)$$, the mean response of the IR to the three stimuli, is similar to the previous case, although the initial drop is much stronger here when the stronger stimuli are used (Fig 8B). A small peak in the time-dependent standard deviation $$ {\hat{\sigma }}_{\text{ir}}(t)$$ is observed right after the signal onset for the two stronger stimuli (Fig 8C). The pronounced initial response to the stronger signals partially compensates the shorter duration of the signal in terms of chances of reaching the detection barrier. As a result, the effect size measured by the IR for the two shorter, but stronger, stimuli is larger than in the previous case, and a clear dependence on the stimulus duration can be seen also in this case, as shown in Fig 8D.

The time-dependent mean of the DR activity, $$ {\hat{\mu }}_{\text{dr}}(t)$$, shown in Fig 8E, displays an initial drop followed by two peaks for each of the three stimuli. Here, however, both the first trough and the two subsequent peaks are more pronounced for signals of larger intensity. The behavior of $$ {\hat{\sigma }}_{\text{dr}}(t)$$ is similar to $$ {\hat{\sigma }}_{\text{ir}}(t)$$, and it marginally favors the detection of the longer signals (Fig 8F). As the most prominent features of $$ {\hat{\mu }}_{\text{dr}}(t)$$ have opposing effects and become stronger simultaneously upon growing stimulus intensity, the net effect on the effect size is barely noticeable (Fig 8G). Furthermore, the effect size is of magnitude comparable with the experimental observations.

Results obtained from the DNR are qualitatively similar to those from the DR. The principal differences are that both positive and negative deflections in $$ {\hat{\mu }}_{\text{dnr}}(t)$$ are more pronounced, and that the plateau between the two peaks is below the zero level (Fig 8H). The increase in $$ {\hat{\sigma }}_{\text{dnr}}(t)$$ is similar to that of $$ {\hat{\sigma }}_{\text{dr}}(t)$$ (Fig 8I). Finally, the effect size is essentially independent of the stimulus, as in the experiments (Fig 8J).

Hence, the effect measured by the DR and the DNR shows barely any dependence on the intensity of the stimulus, as it is observed in the experimental data. However, the effect size measured by the IR is larger and dependent on which of the three stimuli is used, which is not consistent with the data.

Effect of stimulus regularity

In the third and last in silico experiment, random stimuli will be used to evoke irregular spike trains. These stimuli, in accordance with the experimental procedure [6, 48], are a random shuffling of six current steps of length 10, 20, 40, 80, 160, and 90 ms, and with current intensity 100%, 50%, 25%, 12.5%, 6.25%, and -50% of the maximum current, respectively. In other words, each sequence consists of a random permutation of five positive (depolarizing) current steps with intensity inversely proportional to the duration and of one hyperpolarizing step, which inhibits the cell from firing. Two example signals are shown in Fig 9A. Note that stimuli are varied in each trial and not frozen. The irregular stimuli are constructed such that their total duration is 400 ms. The response probability to these stimuli will be compared to that of regular steps of 400 ms at 25% of the maximum current, which was used in both previous cases (plotted in blue). In the experiments, irregular current injections generated spike trains with an average firing rate of (24 ± 11) Hz and average CV of (1.1 ± 0.3). In the model, the average rate is (27 ± 5) Hz and the average CV (1.3 ± 0.3).

Fig 9

In the comparison irregular vs. regular stimulation, simulation results from DNR agree well with the experimental data, whereas results from the IR are most inconsistent with the data.

Blue lines and symbols refer to the 400 ms regular stimulus, as in the previous cases; Red and orange lines and symbols refer to two different random samples of 10000 irregular stimuli (two specific realizations are shown in A, the average stimulus is shown in B; for details see text). Black line is catch trial condition (no stimulus). First row (C,F,I): deviation from the spontaneous value of the time-dependent mean readout activity normalized by the spontaneous standard deviation, defined in Eq (13). Second row (D,G,J): deviation from the spontaneous value of the time-dependent standard deviation of the readout activity normalized by the spontaneous standard deviation, defined in Eq (16). Third row (E,H,K): average effect size. Open circles with error bars: experimental data (average over 62 RS cells and 1780 total trials); filled circles of corresponding color: model simulations. First column (C,D,E): IR. Second column (F,G,H): DR. Third column (I,J,K): DNR.

In Fig 9 we compare the simulation results obtained when irregular stimuli are used to those obtained from the 400 ms regular current injection. The response to the latter is shown once more in Fig 9 because it serves here as reference case but it will not be discussed in depth (see the discussion above). Results for irregular stimulation are based on two sets of irregular stimuli, constructed by choosing from all possible permutations with equal probability. Despite the large total number of trials used here (N_trials = 10000), it is advisable to compare results for two independent sets of stimuli because the large number of possible permutations (6!=720) implies that the number of trials per signal is limited and that finite-size fluctuations due to the particular choice of signals may be non-negligible. Results for the two independent sets of irregular stimuli are plotted in red and orange in Fig 9. We recall that $$ {\hat{\mu }}_X(t)$$ and $$ {\hat{\sigma }}_X(t)$$ are obtained by averaging over different realizations of the irregular stimuli and are not related to the two particular signals shown in Fig 9A.

The average signal is plotted in Fig 9B. By recalling that each irregular signal corresponds to one permutation of the step sequence, and recognizing that the reversed sequence is another valid and equally probable permutation, one can conclude that the average signal must possesses time-reversal symmetry. Indeed, it displays a peak just after the stimulus onset, followed by a mild trough and then by a plateau barely above the zero level. Just before the end of the stimulation time window, the symmetric trough followed by a peak can be seen. Accordingly, $$ {\hat{\mu }}_{\text{ir}}$$ shows two dips at the same time where the two peaks in the average signal are seen (Fig 9C). Note that although the first and second half of the average signal are perfectly symmetrical, the second dip in $$ {\hat{\mu }}_{\text{ir}}(t)$$ is more pronounced than the first.

Considering how different the average signal is from each particular realization of the irregular sequence, it may seem questionable to average over signals that provoke rather heterogeneous responses. However, the detector as well as the animals in the actual experiments do not know which realization of the irregular sequence is used in each trial. Therefore, this averaging ensemble, in which the stimulus is drawn in each trial, correctly represents the experimental situation and it makes sense to consider its time-dependent mean, as done above. The variability due to the particular realization of the input signal, which mostly averages out in the mean, reveals itself in an increased time-dependent standard deviation, $$ {\hat{\sigma }}_{\text{ir}}(t)$$ (Fig 9D), which is above the zero level in the entire stimulation time window and grows further towards the end of the stimulus. This increase of $$ {\hat{\sigma }}_{\text{ir}}(t)$$ above the zero level enhances the chances of reaching the detection threshold. As a result, the effect size upon irregular stimulation is as large as in the experiments, but not as large as that observed for the regular stimulus (Fig 9E filled dots), which is not consistent with the experimental observations, in which it is the other way around (Fig 9E, open circles with error bars).

The average DR activity in response to the irregular stimuli (Fig 9F, red and orange line) and to the 400 ms regular stimulus (blue line) are rather similar to each other. The main difference is that the initial trough and peak are somewhat smaller for irregular stimulation. Furthermore, a small dip is observed for irregular stimuli just before the last peak. Also the standard deviation $$ {\hat{\sigma }}_{\text{dr}}(t)$$ (Fig 9G) is similar in the two cases. Consequently, the average effect size upon irregular stimulation measured by the DR is similar for regular and irregular stimulation (Fig 9H, red and orange vs. blue full circles), which is in better agreement, but not completely consistent with the experimental observations.

The time-dependent mean DNR activity $$ {\hat{\mu }}_{\text{dnr}}(t)$$ (shown in Fig 9I) is qualitatively similar to $$ {\hat{\mu }}_{\text{dr}}(t)$$ although the peaks are more pronounced. Importantly, $$ {\hat{\sigma }}_{\text{dnr}}(t)$$ is generally larger than $$ {\hat{\sigma }}_{\text{dr}}(t)$$ in the case of irregular stimulation (Fig 9J) and displays a strong peak at the end. The larger $$ {\hat{\sigma }}_{\text{dnr}}(t)$$ sented (Fig 9K), which is as large as in the data.

In summary, Fig 9 shows that the DNR has the largest degree of consistency with the experimental observation that irregular stimuli are easier to detect than a regular current step, as opposed to the IR, which yields a smaller effect upon irregular stimulation than upon regular stimulation.

Roles of cell types and delay difference on readout performance

After having shown that the DNR is most consistent with the experimental results, we will now explore how the readout performance depends on the size of the three readout populations, i.e. the effect of each neuron class on the detectability of the single-neuron stimulation. The simulation results are summarized in Fig 10. For brevity, we will focus on the IR (first row in Fig 10) and DNR (second row in Fig 10), and use only the 400 ms stimulus (represented by blue circles) and irregular stimulation (red squares).

Fig 10

Increasing the size of RS and SOM readout sets enhances the detectability of single-cell stimulation, whereas the size of the FS readout set has an opposing effect.

Effect size obtained upon 400 ms regular (blue circles) and irregular (red squares) stimulation as a function of the size of each readout population. The upper row (A, C, E) shows results for the IR readout, the lower row (B, D, F) displays the effect size measured from the DNR readout. The reference parameters used in the previous figures are indicated by the vertical dashed lines; in each panel the size of the readout population indicated on the x-axis is varied while all other parameters are kept fixed at their default value (as in Figs 7–9).

First, we can systematically vary the size of the RS readout population. Intuitively, the effect size should increase, the more RS neurons are fed into the readout. Simulation results demonstrate that this is indeed the case both for the IR readout (Fig 10A) and the DNR readout (Fig 10B): the effect size grows monotonically with the size of the RS readout population.

Next, we can study the effect of a systematic variation of the FS readout population’s size. In Fig 10C and 10D), it can be seen how the effect size measured by the IR (DNR) decreases when the FS readout population is enlarged. The strong self-inhibition of FS cells tends to stabilize their firing rate. Hence, despite the strong SOM response within the BCN, the average firing rate of FS cells within the BCN is, on average, only slightly reduced during the single-cell stimulation. However, due to the strong output weights of FS cells, they have a considerable effect on the readout activity. Because they disinhibit the readout during the stimulation, the effect of FS cells on the readout is antagonist to that of RS cells.

The readout performance is improved by feeding more SOM neurons into the IR (Fig 10E). The readout performance of the DNR is generally enhanced by enlarging the number of SOM readout cells (Fig 10F). However, the effect size tends to saturate (for irregular stimulation) or even slightly decrease when the size the readout population approaches the entire SOM population (200 cells). One possible reason may be stronger cross-correlations between SOM cells.

As a final comment to Fig 10, we note that the DNR detects the irregular signal with a higher accuracy over the entire range of parameters, while the IR readout is almost always better at detecting the regular signal. In this respect, these findings are robust, and the DNR is a stronger candidate than the IR as a possible readout mechanism which is in line with the experimental results.

One further crucial parameter for the operation of the DNR is the delay difference of the two input pathways (ΔT). It is therefore important to check how robust our results are with respect to this parameter. As shown in Fig 11, the effect size measured upon irregular stimulation is larger than that observed for regular stimulation also for when ΔT is decreased below 10 ms, although it decreases in magnitude in both stimulus conditions. When ΔT is further decreased below 5 ms, the effect size reaches zero and becomes negative for the regular stimulus. Hence, the qualitative picture is robust for ΔT ⪆ 5 ms, as the weaker effect could be enhanced by leveraging, for instance, the dependencies seen in Fig 10.

Fig 11

Effect size measured by the DNR is robust for ΔT ⪆ 5 ms.

Effect size measured by the DNR as a function of the delay difference ΔT measured upon regular (blue circles) and irregular (red squares) stimulation.

We note here that it is difficult to try out all possible parameter combinations due to the high computational cost of simulating the network for multiple conditions and a large number of trials, which is needed to measure a possibly small effect size reliably. The simulation results shown in Figs 10 and 11 required about two weeks on a computing cluster with more than 200 cores.

Robustness of readout mechanisms against slow input non-stationarity

The above results show that the DNR is more consistent with the experimental results. However, the effect size measured by the IR is larger in magnitude in most cases, except for the case of irregular stimulation, in which the DNR is only slightly better. This observation raises the question of why the animal should opt for a detection strategy which tends to detect the signal less often, considering that experimental subjects were rewarded upon successful detection. One possible answer, which we will briefly explore in the following, is that the DNR might be more robust with respect to slow fluctuation of the background input. It is possible to get a feeling for these slow non-stationarities by looking at Fig 12A, which shows the spontaneous firing rate of some cells over long juxtacellular recording sessions. The strong changes in firing rate are mostly unrelated to any stimulation event. A readout that integrates spiking in a sliding time window might experience more difficulties in distinguishing the possibly small changes induced by the single-cell stimulation from the strong background variations, whereas a differentiator readout might still separate the timescales of stimulus and background noise.

Fig 12

The differentiator network readout (DNR) is more robust than the integrator readout (IR) against static readout noise, a proxy for slow modulations of the background network activity.

A: Firing rate of four different cells (filtered with a 1 sliding window) during juxtacellular recording sessions in the barrel cortex. The slow changes in firing rate are unrelated to stimulation events and give an example of slow non-stationarities present in vivo (note that the time window shown here is larger by orders of magnitude from the typical duration of a single trial in the detection experiments). B: Effect size measured in the model for the case of 400 ms stimulation on increasing σ_r,ext, the standard deviation of the global background input (input readout noise). Light blue diamonds represent results obtained from the integrator readout (IR), whereas green diamonds are results from the differentiator network readout (DNR). C: same as B, but for the case that irregular stimuli are used.

Although a thorough investigation of this hypothesis goes beyond the scope of this study, in this last subsection we will explore this conjecture by adding a “static modulation” to the network background input. The basic idea is that a slow modulation of the network input will look nearly constant, if observed for the total duration of one of our trials (≈ 2 s). More precisely, in each trial we randomly draw r_ext,th, one of the parameters that sets the global intensity of the external background noise (see Methods). Each sample is drawn from a lognormal distribution with constant mean (equal to 10 Hz, as in the simulations above) and standard deviation σ_r,ext (where σ_r,ext = 0 corresponds to the situation considered so far, i.e. no background noise modulation). Fig 12B shows the effect size measured by the IR (light blue diamonds) and DNR (green triangles) upon increasing σ_r,ext for the case of 400 ms regular stimulus. Both detectors are affected by the additional noise in the background input. However, the effect size measured by the IR drops faster and becomes eventually smaller than the effect measured by the DNR. For the case that an irregular stimulus is applied, Fig 12C shows that both detectors suffer a performance loss when σ_r,ext is increased. However, the performance loss is stronger for the IR than for the DNR, which measures a larger effect size over the entire range of σ_r,ext considered here.

Single-neuron properties and total input to neurons

We modeled all neurons as leaky integrate-and-fire point neurons [90]. The kth neuron follows the differential equation

If the kth neuron belongs to the FS population, its total input current I_total,k is just the sum of the external input and of the recurrent input, i.e. it reads:

The adaptation current in the last equation obeys [37, 38]:

Both Δa_k and τ_a,k are randomly drawn from a lognormal distribution with standard deviation equal to 20% of the mean. For RS neurons, the mean of the two distributions are τ_a,e = 100 ms and Δa_e = 0.3 nA, respectively; for SOM-LTS neurons they are τ_a,s = 50 ms and Δa_s = 0.2 nA, respectively.

With this choice of parameters, the strength of the spike-frequency adaptation roughly agrees with in vitro measurements from the layer IV of the rat barrel cortex [29, 36].

External input to the network

The external input encompasses one constant term and two excitatory Poissonian shot-noise processes:

The constant term is set to R_m,k I₀ = 10 mV for all neurons. The second term in Eq (28) represents the input from the thalamus, and the third mimics incoming spikes from surrounding cortical areas. Because the thalamus has a higher average firing rate, “thalamic” input spikes at times $$ t_{k,j,l}^{th}$$ occur at an average rate of r_ext,th = 10 Hz, while “cortical” input spikes at times $$ t_{k,p,q}^{bc}$$ have a lower rate of r_ext,bc = 2 Hz. The number of external input spike trains depends on the cell type. Experimental studies suggest that SOM cells, in contrast to RS and FS cells, receive only weak input from the thalamus and from distant brain areas [29, 40]. Therefore, if the kth neuron belongs to the SOM-LTS population, then the number of external inputs is set to zero C_ext,th,k = C_ext,bc,k = 0, whereas when k is a RS or a FS neuron, then C_ext,th,k = 500. Furthermore, dendrites of FS neurons tend to be more localized than those of pyramidal cells, i.e. to receive more input from local RS neurons and less from distant ones. Hence, the number of inputs mimicking the cortical surroundings is C_ext,bc,e = 2000 when k is a RS neuron, and C_ext,bc,i = 1000 when k is a FS neuron. Each input spike causes a PSP drawn independently from an exponential distribution with mean J_ext,e = 0.1 mV when k is a RS neuron, and from an exponential distribution with mean J_ext,i = 0.2 mV when k is a FS neuron, because both thalamic and cortical excitatory postsynaptic potential (EPSP) amplitudes are larger in FS cells than in RS cells [29].

Recurrent input to RS neurons

The recurrent input term I_rec,k(t) depends on the neuron type. If k is a RS cell, it is

Connections to neuron k originate from three sets of neurons: $$ P_e(k)$$, formed by C_ee = 300 randomly selected RS neurons, $$ P_i(k)$$, consisting of C_ei = 200 randomly selected FS neurons, and $$ P_s(k)$$, composed of C_es = 100 randomly selected SOM-LTS neurons. Hence, the connection probability of RS-to-RS synapses is C_ee/N_e = 15%, of FS-to-RS and of SOM-LTS-to-RS is C_ei/N_i = C_es/N_s = 50%, consistent with the experimental observations that the connections between RS cells are sparse whereas those between RS and inhibitory cells are dense [27, 29, 41, 91–93]. Transmission delays are drawn uniformly in the range 0.5 ms to 1.0 ms [93]. All synaptic weights in Eq (29) undergo short-term depression (STD)

Recurrent input to FS neurons

The recurrent input to a FS neuron reads:

Effective model for gap junctions

Both FS and SOM neurons in the rat somatosensory cortex are coupled by gap junctions [32, 34, 40, 94]. In a simplified picture, gap junctions act as a passive conductance coupling between the membrane voltage of two neurons. The standard way of mimicking the effect of a gap junction between neuron k and neuron ℓ would be the following additional current for neuron k [95, 96]:

Recurrent input to SOM-LTS neurons

Finally, the recurrent input to a SOM-LTS neuron is

The three terms in Eq (33) represent the input from excitatory RS neurons, from inhibitory FS neurons, and from gap junctions, respectively. Gap-junctions are modeled in the same way as for FS neurons: their amplitudes and delays are drawn from the same distributions. The first term in Eq (33) is the input from C_se = 1000 randomly chosen RS neurons (connection probability C_se/N_e = 50%). These are the only connections that undergo short-term facilitation instead of depression, and for which random transmission failures were modeled (details on the model below). The static baseline amplitudes J_ki of each synapse are drawn independently from an exponential distribution and have mean J_se = 0.1 mV. The second term in Eq (33) represents the input from C_si = 100 randomly selected FS neurons (connection probability C_si/N_i = 25%). These connections have the average maximum strength J_si = 0.25 mV, undergo short-term depression and obey Eq (30). Chemical synapses between SOM neurons are infrequent and weak [29, 40] and were omitted for simplicity.

Model of short-term depression

Except for those connecting RS to SOM neurons, all chemical synapses in the model undergo short-term depression. Each weight J_kj(t) has a time dependence described by

The parameter values chosen to model strong depression (all connections depicted in blue in Fig 1) are τ_D = τ_D,s = 150 ms and U_se = U_se,s = 0.2. With this choice, the eighth PSP of a 40 pre-synaptic regular spike train is about one half of the maximal amplitude [29]. Most chemical synapses in the barrel cortex are depressing ([29, 43, 69]). However, inhibitory synapses originating from SOM-LTS neurons and terminating onto RS neurons show only weak depression or slight facilitation. Here, these connections are modeled as mildly depressing (Fig 1, light blue) by setting τ_D = τ_D,w = 50 ms and U_se,w = 0.05. For simplicity, also SOM-to-FS connections were given the same STD parameters.

Short-term facilitation and transmission failures

Excitatory synapses from RS neurons to SOM-LTS neurons (depicted in red in Fig 1) are strongly facilitating ([29, 41, 42]). If the parameter U_se in Eq (36) is turned into a dynamical variable, u(t), facilitating synapses can be described [45, 99]. The amplitude of the PSPs is proportional to the product R(t)u(t). Considering the connection from the RS neuron i to the SOM-LTS neuron k, the time evolution of the synaptic amplitude is described by (note that the conventions have been slightly changed with respect to ref. [45]):

RS-to-SOM synapses stand out from all other synapses considered here because of a much higher occurrence of transmission failures at low presynaptic firing rates (the average failure rate is ≈10% for RS-to-RS synapses, ≈ 5% for synapses to and from FS neurons, and ≳ 50 for RS to SOM-LTS synapses [29]). However, the failure rate of RS-to-SOM-LTS synapses decreases to ≈ 10% upon repeated stimulation at 40 (failure rates for other synapses weakly depend on the presynaptic firing rate [29]). Here, transmission failures are modeled only for RS-to-SOM synapses via a stochastic binary variable S(p_f):

In the last equation, p_f,rest = 0.5 is the baseline failure rate. Upon each presynaptic spike, the failure rate decreases by G(p_f, Δp_f, p_min) and relaxes back to the baseline value with the time constant τ_f = 250 ms. The size of each downward jump is Δp_f = 0.1 but is constrained to values above Δp_f = p_min = 0.1, a condition which is imposed by the piecewise linear function

In the end, the synaptic weight from the RS neuron i to the SOM-LTS neuron k obeys the following equation:

If the average effect of synaptic failures is taken into account, a 40 Hz presynaptic stimulation causes the eighth PSP to be about eight times larger than the first, which is in a reasonable qualitative agreement with the strong amplification measured in vitro [29, 42].