Global error signal guides local optimization in mismatch calculation

Domain: biology
Published: 2026-03-12
Editorial outcome: major_revision

Abstract
Corollary discharge denotes internal signals about the expected sensory consequences of one’s own actions, leading to reduced sensory responses to self-generated stimulation. To investigate the underlying neural circuit mechanism, here we introduce a biologically plausible three-factor learning rule, where a global signal guides the updating of local inhibitory synapses to enable the computation of mismatches between stimulus and their predictions. We show that our network model, endowed with both positive and negative prediction-error neurons, accounts for the salient physiological observations of motor-visual and motor-auditory mismatch responses in mice. Moreover, the model predicts that learning induces a bimodal distribution in activity correlation with stimulus and prediction, supported by our analysis of neural data from a recent experiment. Together, these results link global modulation to local learning for predictive error computation in the sensory areas, and predict how disrupting inhibition impairs mismatch computation in specific ways.

Introduction
Since Hermann von Helmholtz, an outstanding question in Psychology and Neuroscience has been to elucidate how self-generated movement is monitored in the brain through an internal signal called corollary discharge or efference copy, that conveys (top-down) expectation of sensory consequences of one’s own actions, to be compared with the actual stimulus induced (bottom-up) input in a sensory system [17] (Figure 1a). For instance, we make rapid eye movements several times per second. Yet, our vision is stable thanks to corollary discharge associated with saccades [64]. When we walk, we do not perceive that the world is moving, even though the visual input to our eyes suggests motion [41]. Similarly, when playing an instrument in a band, we may tune out the sound of our own actions relative to those of others [11]. What neural mechanisms underlie this selective suppression? Answering this question is of interest for basic neuroscience as well as psychiatry, since deficits in corollary discharges are believed to be a root cause of hallucination and delusion associated with mental disorders [24, 15].

In recent years, mismatch computation between a sensory response and corollary discharge of self-motion has been cast in the general framework of predictive coding [37]. Originally proposed for perception, predictive coding highlights the idea that a stimulus-induced response in a sensory circuit is compared with an internally generated expectation/prediction signal about that stimulus, and is suppressed when expectation is increased [56]. Applied to motor-sensory feedback, the motor cortex sends a copy of the movement signal to sensory areas, which serves as a local prediction of the incoming stimulus. If the prediction and actual input match, their signals cancel out, leading to reduced or no sensory perception (Figure 1b). For example, we may press a play button and hear the expected tune. But when the prediction and stimulus mismatch, the sensory area responds, signaling a prediction error (Figure 1c). This mismatch can occur in two ways: either a sound is heard without pressing the button, or we press the button while no sound follows. Both types of mismatch elicit stronger neural responses in mice than expected stimuli [2, 4]. Moreover, suppression in the expected case depends on experience, suggesting a critical role for neural plasticity.

Though the mathematical description of mismatch is as simple as computing the absolute difference between stimulus and prediction signals, it is not clear how the brain implements this operation (Figure 1d). One challenge is that neurons in sensory areas exhibit low baseline firing rates [18, 51, 57], making them highly responsive to excitatory input but relatively insensitive to inhibition. A straightforward solution is to introduce two separate neuronal populations: positive prediction error (pPE) neurons, which respond to stimulus minus prediction, and negative prediction error (nPE) neurons, which respond to prediction minus stimulus [37]. Another challenge is that long-range projections in the brain are predominantly excitatory. Therefore, to compute a difference between stimulus and prediction signals, local interneurons must be involved to convert long-range excitatory inputs into local inhibition.

However, accurate prediction error computation also requires well-tuned local connectivity. With naive connectivity (Figure 1e), stimulus and motor-related inputs can vary across individual neurons. Yet, these neurons receive dense inhibition from local interneuron populations [36], suggesting that each neuron receives a comparable level of inhibition. As a result, expected stimulus input cannot be effectively canceled by prediction-driven inhibition, and vice versa. This mismatch implies the need for a learning mechanism that tunes local inhibitory connections to accurately align inhibition with its corresponding excitatory drive, such that relayed stimulus inhibition targets prediction-driven excitation and relayed prediction inhibition targets stimulus-driven excitation. (Figure 1f). Notably, excitatory plasticity may not play a major role, as overall excitatory responses to motion or visual input are not significantly altered after learning [2]. Taken together, successful prediction error computation likely depends on inhibitory plasticity. This form of plasticity has been widely observed in cortex [72, 39, 13] and is modulated by neuromodulators [12]. In particular, noradrenaline signals broadcast by the locus coeruleus may facilitate local circuit optimization [35, 33], though the underlying mechanism remains unknown.

In this study, we introduce a biologically plausible three-factor learning rule, inspired by [26], that acts on inhibitory synapses from interneurons to pyramidal cells and enables local circuits to compute prediction errors. We show that this local rule converges to the same optimal connectivity as that learned via gradient descent. We test the rule in models composed of either simplified ReLU or biophysically realistic neuron models. The resulting inhibitory connectivity reproduces key experimental observations, including the selective motion-modulated auditory response [2] and the motion-visual mismatch response [4]. Moreover, we demonstrate that this learned connectivity is functionally optimal: perturbing interneuron activity, either by excitation or inhibition, consistently disrupts the population-level mismatch response. Finally, the model predicts a bimodal distribution of stimulus- and prediction-related correlation across neurons, a feature supported by reanalysis of experimental data.

Results
Three-factor learning rule optimizes prediction-error computation

We implemented a biologically feasible three-factor learning rule to align inhibition relayed from the prediction (or stimulus) with excitation driven by the stimulus (or prediction), such that population responses are suppressed when inputs are predicted, but enhanced during mismatches:

\Delta W \sim c(t)(R-R_0)r (1)

This learning rule depends on the firing rates of the presynaptic neuron r and postsynaptic neurons R, along with a global third factor c(t) representing the mismatch between prediction and stimulus. Inspired by recent studies [35, 33], this third factor may be mediated by noradrenaline, released from the locus coeruleus. Remarkably, the proposed rule converges to the same optimal solution as the traditional gradient descent algorithm, despite relying only on biologically plausible global error signals. The mathematical proof of this equivalence, valid for any number of ReLU neurons, is presented in the Supplementary Notes. We further validate the theory through simulations in a three-dimensional setting (i.e., 3D movement in the home cage leads to a predictable 3D visual flow change) and by replacing ReLU neurons with a realistic conductance-based model. Here, we provide an intuitive explanation in a simplified case that only includes n=2 RELU neurons to illustrate the necessary conditions required for our model.

This simplified case is designed to capture the minimal requirements for computing a prediction error. First, at least two excitatory neurons (n= 2) are needed, since a single neuron cannot encode bidirectional deviations from a low baseline firing rate. Likewise, two interneurons are included to represent distinct inhibitory populations. Second, one excitatory neuron is assumed to receive stronger sensory (stimulus-driven) input, whereas the other receives stronger top-down (prediction) input; an analogous arrangement applies to the two interneurons. Third, we retain only the dominant excitatory synapse onto each pyramidal neuron and each interneuron, neglecting weaker projections for simplicity. Fourth, all neurons are modeled as rectified linear units (ReLUs). These assumptions together define a minimal circuit consisting of two excitatory neurons (n = 2) and their corresponding interneurons (Figure 2a). Last, we assume both neurons project equally to an output neuron whose activity encodes the difference between stimulus and prediction. The objective is that, through learning, these two neurons evolve into one pPE and one nPE neuron, as shown in (Figure 2a).

We then trained the inhibitory connections from two interneurons IS, IP to both the pPE and nPE neurons to minimize the difference between the output neuron’s activity and the true error, defined as the absolute difference between the stimulus and prediction. The optimal solution in this setting requires the connectivity configured as in Figure 1d. In this case, the pPE neuron only receives relayed prediction inhibition but none of the relayed stimulus inhibition, and its firing rate reflects the positive prediction error RpPE = [S − P]+. Conversely, the nPE neuron only receives relayed stimulus inhibition but not prediction inhibition, and its firing rate reflects the negative prediction error RnPE =[P − S]+.

Learning how to compute prediction error accurately requires a model or an animal to use the prediction signal to cancel the incoming stimulus in the expected cases when both the prediction and stimulus signals are present (Expected set FE : |S−P| \approx 0). However, this result may be from training in both the expected case (FE and the mismatch set FMM : |S−P| > 0). Recent experiments have suggested that the neuromodulation signal may be different in the FE and FMM [35, 33]. To understand the updating difference in these two sets, we first examine the behavior of the gradient descent algorithm within the expected training set FE. Within the expected training set FE, the total inhibitory input converges to match the total excitation (Figure 2b, black arrow). Yet, in this condition, it does not matter whether the inhibition originates from the prediction or the stimulus, resulting in a family of equally valid solutions rather than a unique one. This family of solutions is referred to as the slow learning manifold in the following. Indeed, the optimal response in the expected case is zero, indicating that excessive inhibition can suppress neural activity trivially, thereby satisfying the objective without yielding a meaningful solution. To resolve this ambiguity, training on the mismatch set FMM, where stimulus and prediction inputs are not matched, is necessary. In this case, the source of inhibition becomes critical for accurate mismatch computation. As a result, the synaptic weights move along the slow learning manifold to the unique optimal solution (Figure 2b, blue arrow), where the circuit correctly associates inhibition with its corresponding excitatory drive.

The same optimal solution emerges from our three-factor learning rule, assuming expected training samples are more frequent than mismatched ones, as is typical in experimental paradigms. The third factor, c(t), presumably mediated by a neuromodulator such as noradrenaline [35], controls the sign of learning. Indeed, [62] demonstrated that different concentrations of neuromodulators could control whether paired pre- and post-synaptic activity led to long-term potentiation or depression. During training on expected samples, the rule behaves like Hebbian learning, with weight updates proportional to the product of pre- and post-synaptic activity, causing total inhibition to match excitation and thereby defining a slow-learning manifold. For mismatched samples, the rule becomes anti-Hebbian, with weight updates negatively correlated with the product of activity, driving weights toward the optimal configuration. Although inhibition temporarily deviates from the slow-learning manifold, the dominance of expected samples pulls total inhibition toward the manifold, keeping it close throughout learning. This alternation guides convergence along a zigzag trajectory (Figure 2c).

Next, we test our algorithm with a network of n = 40 ReLU neurons. At each learning step, the true prediction error is sampled from a truncated Gaussian distribution (Figure 2d), and the third factor c(t) reflects the true error. Here, the third factor is only required to differentiate expected cases from mismatched ones by its sign, while the exact form may vary. Here, we choose a linear function for simplicity, while other choice of the functions works (see below). For analytical convenience, we first assume that the sum of stimulus and prediction inputs is constant across neurons. Under this assumption, different neurons share the same slow learning manifold, which is defined by the total excitatory input, and they can be ordered along a one-dimensional axis according to their input affinity to the stimulus (Figure 2e). We further assume that the utility of mismatch signals saturates beyond a certain threshold (see Methods), which causes synaptic weights to different postsynaptic neurons to converge at distinct positions along the slow learning manifold. However, when considering two or more stimuli modulated by the same prediction signal, the sum of stimulus and prediction inputs will naturally differ across neurons, since enforcing identical sums would collapse the representation of distinct stimuli. Indeed, our learning rule still works, where the slow learning manifold differs for each neuron (Figure 2f). The corresponding activity and synaptic weights across neurons are shown in Supplementary Figure 1.

So far, we have shown that the proposed learning rule minimizes neuronal responses in the one-dimensional (1D) expected condition, where the prediction strength matches the stimulus strength (e.g., movement on a treadmill produces corresponding visual flow in one direction). However, animals operate in a three-dimensional (3D) environment, where free movement (e.g., a mouse exploring a cage) generates visual-flow changes along multiple axes. To test the generality of our learning rule, we extended the model to include three stimulus–prediction pairs, such that changes in each stimulus are canceled by changes in the corresponding prediction (Supplementary Figure 2a; see Methods). All stimuli and predictions project to the same excitatory neuron pool, but their synaptic strengths are shuffled to remove correlations between any pair. We further assume that distinct stimuli and prediction signals are relayed through separate interneuron (IN) pools. This assumption is consistent with the reported IN selectivity [50, 78] and/or the plasticity of excitatory synapses onto INs [43, 28]. In this scenario, the expected condition is defined as the L2-norm of the error vector being below a threshold, while the mismatch condition includes all remaining cases (Supplementary Figure 2b). Similarly, the sign of the third factor reverses between expected and mismatch conditions. As before, we sampled the prediction error at each learning step from a truncated 3D Gaussian distribution. After learning, we observed qualitatively similar outcomes: neural responses in the expected condition were suppressed, whereas those in the stimulus-only or prediction-only conditions were amplified (Supplementary Figure 2c, d). This arises because the relayed inhibition aligns with the excitation of the same pair (e.g., Supplementary Figure 2e, f). Similar cancellation patterns were observed for the other two stimulus–prediction pairs (Supplementary Figure 2g, h). In summary, although limited experimental evidence supports this approach in high-dimensional cases as far as we know, our model offers a biologically plausible method for calculating high-dimensional prediction error.

Prediction response depends on context after learning

After establishing our learning rule, we examine neural responses to stimulus and prediction signals before and after learning. From this point on, we use more realistic models for both pyramidal neurons and interneurons (see Methods), though ReLU neurons yield similar results (Supplementary Figure 1). Our pyramidal neuron model consists of two compartments and includes a dendritic voltage variable. With weak coupling between the somatic and dendritic compartments, the model captures both subtractive and divisive dendritic inhibition, depending on the level of dendritic depolarization (Supplementary Figure 3). Here, we use the model with strong dendro-somatic coupling, reflecting a model of Layer 2/3 pyramidal neurons. We further select the single-neuron parameters based on data recorded from Layer 2/3 in the mouse primary visual cortex [27]. Neurons receive a baseline input that maintains low but nonzero firing rates. Each neuron receives input from both the stimulus and prediction signals, with varying strengths. They also receive inhibition that is relayed through local interneuron populations. Given the high density of local inhibitory connectivity [36], we assume that the strength of inhibition is similar across neurons. For simplicity, there is no recurrent connectivity among pyramidal neurons or from pyramidal neurons to interneurons in this model.

Before learning, the majority of neurons exhibit stronger responses when both the stimulus and prediction are present, compared to when only one input is provided (Figure 3a). Here, we assume each neuron receives the same total input in the expected case, and sort neurons based on their affinity to the stimulus. At the population level, the response is highest when both inputs are present. Under our learning rule, inhibitory connectivity evolves from uniform inhibition into selective inhibition, such that excitatory input from the stimulus is cancelled by inhibition relayed from the prediction (Figure 3b, top), and excitatory input from the prediction is cancelled by inhibition relayed from the stimulus (Figure 3b, bottom).

After learning, population neuronal activity reflects prediction error computation (Figure 3c). The left-most neuron receives the strongest stimulus excitation and the least inhibition from the stimulus, resulting in stronger responses when the stimulus is present alone. The right-most neuron shows a similar pattern when only the prediction input is provided. However, when both stimulus and prediction inputs are present, all neurons exhibit low firing rates, slightly above baseline. At the population level, the average response in the expected condition, when both stimulus and prediction inputs are present, is lower than the response to either input alone. Furthermore, comparing activity distributions before and after learning, the average activity decreases only in the expected condition (Figure 3d). The model has similar qualitative results when each neuron receives a different total input in the expected case (Supplementary Figure 4).

From another perspective, the net effect of the prediction signal on population responses becomes context-dependent after learning. Before learning, adding prediction increases activity, both with the stimulus (Figure 3a, blue vs. grey) and without it (Figure 3a, red vs. green). After learning, however, prediction reduces activity when comparing the expected case to the stimulus-only case (Figure 3c, red vs. green), indicating that prediction-modulated inhibition is learned and context-specific. We will examine this further in the next section.

Here, the exact form of the third factor c(t) is not important. To directly test this, we use a piecewise-estimated c˜(t) instead to guide the learning and reach the same qualitative results (Supplementary Figure 5). In this case, the c˜(t) requires only an approximation of the global error by the animal.

Since the neurons on the left respond selectively to positive prediction error, [S − P]+, after learning, we refer to them as pPE neurons. The same applies to the neurons on the right, which respond to negative prediction error, and are referred to as nPE neurons. Along the x-axis of input strength, the neuron identity gradually transitions from pPE (neuron 1 to neuron 20) to nPE (neuron 21 to neuron 40). In the following sections, the representative pPE and nPE neurons correspond to neuron 1 and neuron 40, respectively.

Prediction input selectively modulates the expected stimulus

In many situations, we selectively suppress expected stimuli in complex environments, while leaving the representation of other stimuli unchanged. Recent experiments have shown that this suppression is highly specific: only the response to the expected sound associated with self-generated motion is reduced, whereas responses to other, non-predicted sounds remain unaffected [4, 3]. Does the prediction error computation learned by our model show the same specificity to expected stimuli?

We assess this selectivity by introducing a probed stimulus in addition to the expected stimulus. The excitation of the probed stimulus is generated by shuffling that of the expected stimulus. Importantly, the prediction excitation to each neuron is also shuffled, such that it is uncorrelated with both the excitation of the expected and the probed stimuli (Figure 4a).

During training, as shown previously, the synaptic weights gradually align the inhibition relayed from the stimulus with the prediction input (Supplementary Figure 6a). In contrast, inhibition relayed from the prediction aligns specifically with the expected stimulus, but not with the probed stimulus (Figure 4b). This is because the error signal driving the third-factor learning rule reflects only the mismatch between the expected stimulus and the prediction, and does not incorporate the probed stimulus.

After learning, in the absence of prediction input, the response distributions to the expected and probed stimuli are similar (Figure 4c). In contrast, when both stimulus and prediction inputs are present, strong inhibition emerges selectively for the expected stimulus, but not the probed stimulus (Figure 4d, Supplementary Figure 6b).

To quantify the prediction-modulated response, we used the same modulation index (MI; see Methods) as in [4, 3] to measure activity changes within the same neuron population (Figure 4e). The MI for most neurons for the expected stimulus is close to minus one, indicating strong selective inhibition. In contrast, the MI for the probed stimulus shows a broader distribution centered above zero, consistent with weaker or absent modulation. This difference in MI between expected and probed stimuli was also observed experimentally [4, 3]. One notable difference is that, in the experiments, the MI for the probed stimulus is closer to zero or slightly negative. This discrepancy may be explained by a global, non-specific inhibitory effect from motor to auditory cortex, potentially mediated by Parvalbumin-expressing (PV) interneurons [4].

It is also important for the sensory cortex to maintain accurate stimulus processing in the absence of prediction input. We find that discriminability, measured using a form of the Fisher information metric (see Methods), is comparable between the expected and probed stimuli (Supplementary Figure 6c).

Model reproduces motor-visual mismatch response after learning

In the previous section, we examined how neural responses depended on whether a stimulus was predicted. Here, we investigate a complementary scenario in which the stimulus is absent when a motor signal is present (Figure 5a). Recent studies [2, 35] used a motor-visual mismatch protocol to explore neural responses when expected sensory input was omitted. In this paradigm, dark-reared mice were passively trained in a virtual environment where their movement on a treadmill controlled the visual flow of a grating pattern (coupled training, CT). In parallel, a control group of mice received the same visual input generated from the CT animals, but the visual flow was decoupled from their own movement (non-coupled training, NT). As a result, the CT group may learn that the visual flow was a consequence of their movement, but not the NT group. After training, Mismatch events were introduced by abruptly pausing the visual flow during locomotion for both groups. In these experiments, a larger fraction of neurons in the CT group responded to the mismatch, and the average mismatch response was significantly stronger in the CT group compared to the NT group.

We simulate the same protocol in our model, both before and after learning. During the expected period, both prediction and stimulus inputs are active, whereas during the mismatch period, the stimulus input is paused while the prediction input remains. In our model, both nPE and pPE neurons exhibit reduced activity during the expected condition after learning (Figure 5b), consistent with previous results (Figure 3d). However, these two neuron types respond differently to the mismatch signal. pPE neurons show a negative mismatch response, indicating that they are normally excited by the stimulus. In contrast, nPE neurons show a positive mismatch response, suggesting that the stimulus inhibits their activity. At the population level, more neurons are excited than inhibited during mismatch (Figure 5c). Before learning (Figure 5d), the strongest excitatory and inhibitory mismatch responses at the single-neuron level are similar in magnitude. On average (Figure 5e), responses in the expected condition (stimulus with prediction) are comparable to those in the mismatch condition (prediction only), indicating that mismatch signals are not readily distinguishable. After learning (Figure 5f), however, the strongest excitatory mismatch responses (e.g., neuron 40) are approximately ten times larger than the strongest inhibitory responses (e.g., neuron 1). This asymmetry arises from two factors: first, responses during the expected condition are strongly suppressed; second, nPE neurons become selectively inhibited by stimulus-driven inhibition, resulting in stronger disinhibition when stimulus input is absent.

Mismatch response under perturbation

If the learned connectivity is optimal in computing prediction error, then any perturbations of it are suboptimal and should do worse. To test this, we probed the model by perturbing the activity of interneuron populations, mimicking a photostimulation experiment (Figure 6a).

By inhibiting the IS population, the pyramidal neuron population becomes disinhibited (Figure 6b). However, because inhibition from IS is aligned with prediction-driven excitation, nPE neurons experience stronger disinhibition than pPE neurons. This asymmetry is reflected in the increasing activity across neurons along the neuron index, which is sorted by the stimulus affinity. As a result, nPE neurons no longer compute prediction error accurately: they exhibit elevated responses in the expected condition. Consequently, the change of population response stops reflecting the prediction error. In the opposite scenario, excitation of the IS population leads to widespread inhibition across the pyramidal neuron population (Figure 6c). Under this condition, nPE neurons are strongly inhibited in both the expected and mismatch cases. However, accurate mismatch computation requires nPE neurons to show increased activity in the mismatch condition, and this differential response is lost under excessive inhibition. As a result, nPE neurons fail to compute the prediction error accurately in this case as well. When plotting the mismatch response as a function of perturbation strength (Figure 6d; additional examples shown in Supplementary Figure 7a, b), we observe that performance is optimal at zero perturbation, indicating that any deviation impairs the network’s ability to compute prediction error. This degradation in performance is primarily driven by nPE neurons, rather than pPE neurons (Figure 6e), consistent with the fact that nPE neurons receive stronger inhibition from the IS population.

In addition to perturbing the IS population, perturbations of the IP population (Figure 6f) yield qualitatively similar outcomes. However, the specific mechanisms underlying performance degradation differ across conditions. When the IP population is inhibited, pPE neurons become disinhibited in the expected condition (Figure 6g, Supplementary Figure 7c), resulting in spurious responses that distort prediction error signaling. When the IP population is excited, the example pPE and nPE neurons appear less affected (Figure 6g), but nPE neurons that receive both prediction and stimulus inhibition (e.g., neurons 20 to 30) become excessively inhibited and cease contributing to the population-level prediction error signal (Supplementary Figure 7d).

The impairment of predictive coding by NMDA receptor blockade is widely observed across species [38, 21, 71]. To test this in our model, we block NMDA channels on all excitatory connections. As expected, this impairs mismatch computation (Figure 6h), due to reduced excitatory drive across neurons, leading to a general decrease in response amplitude (Figure 6i, j; Supplementary Figure 7e, f). Despite the differing mechanisms, the same principle holds: any perturbation disrupts prediction error computation.

This concept of optimization can be further tested in a conjugated task, where a mismatch signal is introduced by pausing the prediction input in the model (Figure 7a). Such a mismatch could arise experimentally through optogenetic suppression of axons projecting from motor areas to the visual cortex, thereby disrupting the prediction signal.

In contrast to the previous scenario, pPE neurons exhibit a strong positive response to the mismatch, while nPE neurons show a modest negative response (Figure 7c). Across the neuronal population, the mismatch response gradually decreases from positive to negative values, although the average response remains positive (Figure 7d). When examining mismatch responses as a function of perturbation strength, the population response is again optimized in the unperturbed condition: when neither the IS nor the IP populations, nor the NMDA channels, are manipulated (Figure 7e, g, i). The mechanisms underlying performance degradation under perturbation are similar to those described earlier (Figure 7f, h, j, Supplementary Figure 8). For example, inhibiting the IP population leads to strong disinhibition of pPE neurons in the expected condition, reducing the difference between the expected and stimulus-only responses (Figure 7f, Supplementary Figure 8a). As a result, the population-level mismatch response is diminished (Supplementary Figure 8b). This mechanism is analogous to the effect of inhibiting the IS population in the original mismatch task.

Bimodal correlation distribution emerges after learning

Our model predicts that, after learning, prediction-driven excitation should be balanced by stimulus-driven inhibition, and stimulus-driven excitation should be balanced by prediction-driven inhibition. These relationships between excitation and inhibition can be revealed by examining the correlation between neuronal responses and varying levels of stimulus and prediction input. Here, we focus on how these correlations change as a result of learning.

We first fix the prediction input strength while randomly varying the stimulus input each second (Figure 8a). In our model, before learning, a subset of neurons exhibit stronger responses when both stimulus and prediction inputs are present, compared to either input alone (Figure 3a; approximately from neuron 15 to 25). Indeed, two example neurons (Figure 8b; neuron 15 and neuron 25) show positive correlations with stimulus input. After learning, however, almost no neurons maintain elevated responses when both inputs co-occur (Figure 3c). After learning, neuron 25 switches from a positive to a negative correlation with stimulus strength (Figure 8b). We then reverse the setup by fixing the stimulus input while randomly varying the prediction input (Figure 8c). Before learning, both example neurons again show positive correlations with the prediction strength. After learning, neuron 15 switches to a negative correlation with prediction strength (Figure 8d).

We next vary the stimulus and prediction input strengths simultaneously and compute each neuron’s response correlation with stimulus and prediction input separately. To improve statistical reliability, we use a larger model with N = 200 neurons and a continuous simulation duration of 200 seconds. In Figure 8e, each dot represents a single neuron, with its correlation with stimulus input on the x-axis and correlation with prediction input on the y-axis. The color indicates the mismatch response after learning. The dots lie approximately on a unit circle, suggesting that the variance in neural responses can be fully explained by the combination of stimulus and prediction input.

We define the angular deviation θ for each neuron as the difference between the angle of its correlation vector in polar coordinates and π/4. A value of θ = 0 indicates equal correlation with stimulus and prediction input. We further weight each θ by the radius of the corresponding point in polar coordinates and plot the weighted distribution of θ. Before learning, θ follows a unimodal distribution centered near zero (Figure 8f). After learning, however, the correlation points cluster in the second and fourth quadrants, consistent with a functional split between nPE and pPE populations (Figure 8g). This separation is further reflected in the emergence of a bimodal distribution in the weighted distribution of θ (Figure 8h).

Finally, we test whether a similar bimodal correlation structure emerges in experimental data. Previously in Figure 5, we showed that our model reproduces neuronal responses observed in both non-coupled training (NT) and coupled training (CT) mice during a motor–visual mismatch protocol [2]. Here, we analyze data from the same study, focusing on neuronal responses in NT and CT mice under the open-loop condition, where movement and visual flow are decoupled.

As in the model, we compute the correlation between neuronal responses and either visual flow (stimulus) or locomotion speed (prediction), and represent these in polar coordinates (Figure 8i, k). We then plot the weighted distribution of the angular difference θ as defined earlier (Figure 8j, l). In the NT group, the distribution of θ is unimodal and centered near zero, indicating similar correlation with stimulus and prediction inputs. In contrast, the CT group shows a bimodal distribution, suggesting that the functional separation into pPE and nPE neurons emerges only with coupled sensorimotor experience, consistent with our model predictions. This bimodality is further confirmed by fitting the angular data to Gaussian mixture models (GMM) with one or two components (see Methods; Figure 8j, l, Supplementary Figure 9).

Methods
Simplified Neural Circuit

We begin with a simplified neural circuit with ReLU neurons to illustrate the mechanism of our three-factor learning rule. Here, all the variables are dimensionless and are updated instantaneously without temporal dynamics.

For the i-th ReLU unit, activity is given by:

R_i = [\lambda_S u_S + \lambda_P u_P - W_S u_S - W_P u_P]_+ (2)

For simplicity, we omit the subindex i in the following. Within the above equation, λS and λP denote the excitatory input strengths from stimulus and prediction, respectively. The stimulus and prediction input are denoted by uS and uP, or simply S and P when unambiguous. WS and WP represent the inhibitory synaptic weights from interneuron populations IS and IP, whose activity is the same as the stimulus and prediction input uS and uP. The operator [x]+ denotes a rectifier (ReLU) function, where [x]+ = max(0, x). The synaptic weights WS and WP are plastic and updated according to our three-factor learning rule.

In the case of homeostasis, λS + λP = 1. Thus, for i-th neuron, we can write λS = λi and λP = 1 − λi. In the case of two neurons Nneu = 2, we set λ1 = 1, λ2 = 0. In the case of forty neurons Nneu = 40, the values of λi are linearly spaced from 1 to 0 across the population.

In the case without homeostasis, we shuffled the λPi across neurons such that the total input λSi + λPi is no longer constrained to equal 1.

Three-factor learning rule

Our learning mechanism can optimize the mismatch calculation, as proven in the Supplementary Notes.

The inhibitory connections are initiated with Wj(t0) for both interneuron populations j = S, P with some jitter, mimicking the average synaptic strength from a population that should be similar across post-synaptic neurons. Again, we omit the subscript i in the following derivations.

These inhibitory synapses are updated with a three-factor learning rule as follows:

\Delta W_j = \alpha c(t)(R - R_0)r_j, (3)

where j = S, P. Here, R denotes the activity of the postsynaptic pyramidal neuron, and rj represents the steady-state activity of the presynaptic interneuron. R0 is a small target firing rate shared by all neurons, and α is the learning rate.

The third factor, c(t), is defined as:

c(t) = c_0 - |u_S - u_P|, (4)

such that when uS = uP, c(t) = 1, and when |uS − uP| = 1, c(t) = −1. The threshold c0 of 0.5 is chosen for simplicity. Also, our choice of c0 implies that a rough estimation of the prediction error is sufficient to serve as the third factor, since any input with c(t) > 0 is considered as an expected case. We further show that an estimated piece-wise third factor is sufficient to guide the local learning. In this case,

c˜(t) = 1, if |u_S - u_P| < 0.2; (5)

c˜(t) = -1, if |u_S - u_P| > 0.8; (6)

c˜(t) = 0, others. (7)

During training, ϵ = uS − uP is drawn from a truncated normal distribution (Figure 2d) with mean zero and standard deviation σϵ = 0.5, reflecting the assumption that training samples are dominated by expected cases where uS ≈ uP.

The mismatch modulation saturates when the mismatch signal becomes sufficiently large, and is defined as:

c(t)eff = c(t) [1 - R¯ / (\kappa R_0)]+, if c(t) < 0, (9)

where κ is a scaling factor that sets the saturation level of the mismatch utility, and R¯ is the average activity of the neuron population.

In the case of realistic neuron models, synaptic weight updates are further capped at high conductance values:

\Delta W_j^eff = \Delta W_j \cdot (\eta W_0 - W_j) / (\eta - 1), (10)

where η > 1 is the saturation factor for synaptic weights, and W0 is the initial averaged synaptic strength.

Generalization in the 3-dimension

We next tested whether the proposed learning rule generalizes to a scenario in which the model need to use a three-dimensional (3D) prediction signal to cancel a 3D stimulus input. The simulation results are shown in Supplementary Figure 2.

In this case, the i-th ReLU unit, activity is given by:

R_i = [\sum_{j=1}^3 (\lambda_i^{Sj} u^{Sj} + \lambda_i^{Pj} u^{Pj} - W_i^{Sj} u^{Sj} - W_i^{Pj} u^{Pj})]_+ (11)

For each timestep, the maximum of the pair νj = max(uSj, uPj) is randomly drawn from a linear distribution from 0 to 1. Then, the error for the pair ϵj = uSj − uPj is drawn from a truncated Gaussian distribution with mean zero and standard deviation σϵ = 0.4 and further scaled by the maximum input νj. The value of uSj and uPj is calculated based on the νj and the error ϵj. Each stimulus-prediction pair is generated independently. Please note that ϵj is only used in generating the training samples, while the model has no access to the error signal in each dimension.

The third factor remains as a scaler, which is calculated based on the L2-norm of the difference between the 3D-stimulus and the 3D-prediction at each time step:

ϵ = ||(u^{S1}, u^{S2}, u^{S3}) - (u^{P1}, u^{P2}, u^{P3})||_2 (12)

c(t) = c_0 - ϵ, (13)

with c0 = 0.3.

Pyramidal neuron and interneuron models

In this work, we use a simplified hybrid dendritic model. For each pyramidal neuron, the dendritic membrane potential Vd and somatic firing rate Rs are updated at each time step as follows:

Cd dVd/dt = Id + Id,leak - gc(Vd - Vs) + IbAP, (14)

τs dRs/dt = -Rs + f(Is - gc(Vs - Vd)), (15)

where Is and Id are the total synaptic currents into the soma and dendrite compartments, respectively. The dendritic leak current is given by Id,leak = −gd,leak(Vd − Vd,rest), where Vd,rest is the resting dendritic potential and gd,leak is the leak conductance. The soma and dendrite are coupled via conductance gc, resulting in a coupling current Isd = gc(Vd − Vs).

The dendritic compartment also receives input from backpropagating action potentials (bAPs), which we model as a rate-averaged current following [42]: IbAP = −gc(Vd − VbAP)tbAPRs, where VbAP is the averaged membrane potential during a spike, tbAP is the duration of the bAP, and Rs is the firing rate at the soma.

As observed in [42], the membrane potential fluctuations at the soma are considerably smaller (10 mV) than those at the dendrite (50 mV). To simplify the computation of the coupling current between somatic and dendritic compartments, we approximate the real somatic potential in simulations by the time-averaged somatic potential between action potentials, defined as Vs,eff ≡ Vs0 = ⟨Vs⟩t.

The activation function f that converts somatic input current to firing rate is taken from [1]:

f(I) = \Delta V / (\tau(Vth - Vreset) (1 - exp(-a \Delta V / \sigma_V))), (16)

where ΔV = I/gl + Vl − Vth.

The model parameters are selected based on prior experimental work [53, 42]. The coupling conductance gc is set to a small value for neurons with extensive dendritic trees, such as Layer 5 pyramidal cells, where our single-neuron model reproduces shunting dendritic inhibition (Supplementary Figure 3a). In contrast, gc is set higher for neurons with more compact dendritic arbors, such as Layer 2/3 pyramidal neurons, where dendrites are closer to the soma and dendritic inhibition exhibits primarily subtractive effects (Supplementary Figure 3b).

We use a point-neuron model to describe the average activity of each interneuron population:

τ drj/dt = -rj + f(I), (17)

where j = S, P denotes the stimulus- or prediction-driven interneuron population, and f(I) is the same activation function defined previously.

Synapse model

We use a conductance-based synaptic model in our simulations, following [75]. The AMPA and GABA synapses are described by:

Isyn = -gsyn s_syn (VP - Vsyn), (18)

d s_syn / dt = -s_syn / τsyn + Rpre, (19)

where syn = A, G denotes AMPA and GABA synapses, respectively. Here, gsyn is the maximal synaptic conductance, s_syn is the gating variable representing the fraction of open channels, Vsyn is the synaptic reversal potential (determined by whether the synapse is excitatory or inhibitory), τsyn is the synaptic time constant, and Rpre is the presynaptic firing rate.

The NMDA synapse model includes both a voltage-dependent magnesium block, fMg(V), and a saturating gating variable, sN:

IN = -gN sN (V - VE) fMg(V), (20)

fMg(V) = 1 / (1 + exp(-0.062V) / 3.57) (21)

d sN / dt = -sN / τN + (1 - sN) γ Rpre. (22)

Here, gN is the maximal NMDA conductance, VE is the excitatory reversal potential, τN is the NMDA time constant, γ is a saturation scaling factor, and Rpre is the presynaptic firing rate. The function fMg(V) captures the voltage-dependent magnesium block characteristic of NMDA receptors.

The stimulus input to both pyramidal neurons and the first interneuron population IS is modeled as an excitatory conductance to the soma, denoted by gsoma. When illustrating the behavior of individual pyramidal neurons in Supplementary Figure 3, the dendritic compartment additionally receives excitatory and inhibitory input conductances, denoted by gdend. Because our model employs conductance-based synapses with distinct reversal potentials for excitation and inhibition, the total input current is not simply a linear difference between excitatory and inhibitory conductances. This formulation allows us to capture more complex and dynamic behaviors in the circuit.

Input and network connectivity

We include only the necessary connections required for the model. A schematic of the network architecture is shown in Figure 1e. The maximum synaptic conductance from interneuron population j to pyramidal neuron i is denoted by Wji, where j = S, P corresponds to the stimulus- and prediction-driven interneuron populations, respectively.

In the simulations, synapses originating from the same presynaptic interneuron population share a common gating variable, denoted by sj. The total inhibitory conductance received by a postsynaptic pyramidal neuron i from both interneuron populations is given by:

g_i = \sum_{j \in S,P} W_i^j s_j, (23)

where sj is the gating variable for population j.

Our model captures pyramidal neuron activity at both the single-cell and population levels. The network consists of Nneu pyramidal neurons, with two distinct interneuron populations relaying stimulus- and prediction-driven inhibition, respectively.

In the case where the sum of stimulus and prediction inputs is the same across neurons, the excitatory synaptic inputs from the stimulus and prediction are anti-correlated. The stimulus input strength is linearly spaced from a maximum value gSmax to zero. Thus, for the i-th neuron, the stimulus input strength is defined as:

gSi = gSmax (Nneu - i) / Nneu (24)

With input strength u ∈ [0, 1], the total excitatory current from the stimulus is:

ISi = -u gSi (Vs0 - VE). (25)

Similarly, the stimulus input to the IS interneuron population is defined by the synaptic strength parameter gIS. In the perturbation experiments, the perturbation level is normalized by this baseline input strength. A perturbation value of ISperturb = -1 indicates that the IS population receives no excitation when the stimulus is presented, whereas ISperturb = 1 indicates that the IS population receives double the baseline excitation.

The top-down input is modeled as a presynaptic predictor neuron firing at a rate of 5(1 + uP) Hz, where uP ∈ [0, 1] represents the prediction strength. This yields a predictor firing rate in the range [5, 10] Hz. The synaptic weights from the predictor neuron to the dendritic compartments of pyramidal neurons are linearly spaced from zero to a maximum value gEdmax. For the i-th neuron, the top-down connectivity strength is given by gEdi = gEdmax i / Nneu. Furthermore, the excitatory synaptic input is composed of a mixture of NMDA and AMPA receptors, with a fraction κ mediated by NMDA synapses and the remaining (1 − κ) by AMPA synapses.

The input conductance is computed by substituting the predictor neuron’s firing rate into Equations 19 and 22. Specifically, the dynamics of the AMPA and NMDA gating variables become:

d sA / dt = -sA / τA + 5(1 + uP), (26)

and

d sN / dt = -sN / τN + (1 - sN) γ · 5(1 + uP). (27)

The input current to pyramidal neuron i depends on both the gating variables sj and the dendritic membrane potential Vdi:

IEdi = -gAi sA (Vdi - VE) - gNi sN (Vdi - VE) fMg(Vdi), (28)

where gAi = κ gEdi and gNi = (1 − κ) gEdi.

Similarly, the prediction input to the IP population is defined by the input strength gIP. In the perturbation experiments, the degree of perturbation is normalized by the input strength to the IS population, gIS, so that equivalent perturbation magnitudes are comparable across the IS and IP populations.

We further include background stimulus conductance to all neurons in the model, ensuring that each cell maintains a similar low, but non-zero, baseline firing rate.

In the realistic model, we shuffle the synaptic weights from the predictor neuron to pyramidal cells, such that the input strengths from the stimulus and prediction are uncorrelated. Naturally, the sum of stimulus and prediction inputs is different across neurons.

We do not include lateral connectivity between pyramidal neurons in this model. While lateral inhibition may aid in differentiating responses between nPE and pPE neurons, it cannot contribute to canceling stimulus or prediction inputs in the expected condition, where activity is low across the population. For simplicity, we omit these connections from the model.

All parameters used in this study are listed in the Supplementary Table. All code will be available on GitHub upon publication acceptance.

Data analysis

In the motion-modulated auditory response analysis (Figure 3), we quantify the selective modulation effect of the prediction signal using the same modulation index (MI) as defined in [4]. If the prediction increases the response to the stimulus, the MI is positive; if it suppresses the response, the MI is negative. Specifically, for each cell, let rs denote the response to the unmodulated stimulus and rm the response to the modulated stimulus. We compute the angular deviation from the diagonal (indifference) line as: θ = arctan(rm / rs) − π/4. The resulting θ ∈ [−π, π] is linearly rescaled to the range [−2, 2] for visualization clarity. An MI of −1 indicates that the response to the modulated stimulus is zero. Only cells with responses greater than 2 × (1 + 1.28), corresponding to a statistically significant response (p = 0.1) compared to the baseline rate r = 2, are included in the analysis.

The discriminability between stimulus pairs is computed to assess whether the sensory area can continue to process stimuli accurately. Discriminability is quantified using a form of the Fisher information metric [16], specifically the squared d-prime, summed across neurons:

F = \sum_{i=1}^{Nneu} (rAi - rBi)^2 / (rAi + rBi) (29)

where A and B denote two different stimulus conditions. With this formulation, the Fisher information can be interpreted as the discriminability of an optimal linear decoder [46].

In the motion-vision mismatch test (Figure 4), we examine the emergence of a bimodal distribution in the correlation diagram following learning. To generate this distribution, we randomize the strengths of the stimulus and prediction inputs and run 200 trials continuously in the model. For each neuron, we calculate the correlation between its response and the inputs strengths uS and uP. We then compute an angular difference, θ, in polar coordinates, where θ = 0 indicates equal correlation with stimulus and prediction. This angular metric is further weighted by the radial distance in the polar plane to emphasize neurons with larger response variance. A bimodal distribution of θ after learning reflects the emergence of both pPE and nPE neuronal populations. Because neuronal responses in the model are solely determined by the stimulus and prediction input strengths, the response variance is almost entirely explained by these two variables. Consequently, the data points lie approximately on a ring manifold in the correlation space.

We further apply the same analysis to data collected from C57BL/6J mice of either sex (n = 6 per group) under both coupled training (CT) and non-coupled training (UT) conditions, originally published in [2]. The total number of cells under examination is N = 939 for the CT group, and N = 690 for the NT group. We demonstrate that the bimodal distribution emerges only in the coupled training condition by fitting the angular data to Gaussian Mixture Models (GMMs) with one or two components (see Methods; Figure 8j, l, Supplementary Figure 9).

To perform this analysis on weighted data, we first construct a pseudo-dataset. For each cell in polar coordinates with angle θi and radius ri, we replicate θi by 100 × ri times to reflect the weighting. Only data with θi ∈ [−π/2, π/2] are included, in order to avoid spurious peaks near θ = −π. All replicated values are collected to form the pseudo-dataset. The resulting pseudo-dataset is then fitted to Gaussian Mixture Models using the MATLAB function fitgmdist, with either one or two components and default parameter settings.

To evaluate whether one-component or two-component fitting is superior, we use the Bayesian Information Criterion (BIC) [61]. Typically, a difference of ΔBIC > 10 is considered strong evidence in favor of the model with the lower BIC. In our analysis, fitting two-component GMMs to both pseudo CT and NT datasets yields ΔBIC > 10, favoring the two-component model in both cases (Supplementary Figure 9). However, when we overlay the fitted GMMs on the pseudo-datasets, only the CT datasets show a true bimodal distribution. In contrast, the uncoupled condition shows a single dominant peak, suggesting that the means of the two Gaussian components are nearly identical (Supplementary Figure 9a, blue dashed line).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors study a neuronal encoding model of prediction error (PE), which results and includes both positive and negative encoding PE neurons. An objective function is proposed to derive the distribution of positive and negative PE neurons. A 3-factor learning rule is convincingly shown to produce a solution comparable to true gradient descent over the objective/cost function. This is shown both using analytical derivations and simulations. Several predictions of the model are tested in neural data.\nOverall, the idea is interesting, but I think that the model lacks realism, as explained next. I think that the paper is better suited for a more theoretically oriented audience.\nThe authors start to motivate their problem by defining the cost function in Eq. 32. This states that the average activity of two neurons (populations) should be the PE signal |S-P|. When one neuron is assumed to have S inputs and the other P inputs, plus inhibition, the optimal result is one neuron encoding positive PE and the other negative PE, respectively. However, this derivation assumes that the two neurons have the same tuning curve. The case of heterogeneous slopes in their tuning curves has not been considered in this basic theoretical setting (e.g., just two neurons with different slopes), and this strongly limits the interest of the theoretical formulation of the problem. An equivalent formulation to Eq. 32 but with heterogeneous population is not provided as far as I can see.\nMy second major concern is that the learning rate c(t) after Eq. 49 requires exactly computing |S-P|. However, the simple average of activities of a heterogenous population will not equal |S-P|. A decoder would need to be trained to approximate that quantity, so that it can be used in the 3-factor learning rule. If I understand well, the authors have not studied how |S-P| could computed/approximated by heterogeneous populations, nor how it could be used through learning by the 3-factor learning rule, limiting the generality of the results.\nThe derivation that follows Eq. 32 is not very interesting or novel in my opinion. Here the authors take as result that to minimize the cost in Eq. 32 one needs a dataset with both expected and mismatched examples. This is hardly an interesting result because to fit any functon(e.g., the signal |S-P|) based on two base functions (e.g., the activity of two neuronal types) one needs to sufficiently span the range of possible values of the functions’ arguments. Therefore, I am not sure about the novelty of this result beyond saying that to fit a function one needs to have a broad range of values to do so.\nFinally, what would be the objective/cost function in the case of temporal mismatch between S and P? The mathematical analysis seems to be restricted to the case of simultaneously presented S and P.\nThe results on Sec. 0.0.1 are interesting and non-trivial. The analysis on the 3-factor learning rule is elegant and seems to be complete.\nMinor\n-I was missing a section in the Appendix where all the details of the models could be found at once. For instance, in the Methods, the simplified neuron model is explained, but there is no explanation there on how the “u” are defined, chosen and used. In the same place, “two neuron” should be “two neurons”.\n-“One caveat” paragraph in Appendix is repeated twice.", "reviewer_2": "The paper is a well-formulated model of how inhibitory plasticity onto excitatory neurons can lead to the excitatory neurons responding only when a stimulus mismatches a prediction signal that could be an efference copy from motor areas. The three-factor plasticity rule (based on presynaptic activity, postsynaptic activity, and mismatch signal), which has experimental motivation, is local (so feasible) and can produce the optimal solutions in a manner similar to gradient descent. While there are some quirks in the model (two compartments but a soma with a membrane potential that never changes) I don’t think they influence the overall results and claims which are quite general given the learning rules. There is good connection with experimental data.\nMy comments are to help clarity, particularly around the figures which I think would be more helpful if more self-contained.\nSpecific Comments:\nIt is worth a mention or discussion of if and how this study is related to the mismatch negativity observed in EEG studies of humans, or more generally prediction errors arising from sequences of stimuli rather than self-generated motor plans.\nThe mention of homeostasis is rather strange and needs more explanation about what is assumed. Homeostasis is for long-term firing rates of neurons in response to all the stimuli they receive over a significant period of time. Also homeostatic set points are different for each cell, ranging over many orders of magnitude. The assumption that the sum of two specific input strengths is equal to 1 (p.17 after Eq 1) seems very different from any obvious assumption of homeostasis. To be honest, the point of the comparison of the situation with/without homeostasis (Fig 2E/F) is unclear here. Is there a prediction? Or is it just a demonstration that the weight changes together follow the same manifold when using the 3-factor rule as that followed in an optimization procedure?\nFig 2A is confusing as it shows inhibitory connections going to both types of prediction error neuron from both types of interneuron, while excitatory input goes to just one prediction error neuron from either stimulus/prediction, in conflict with Figure 1D/F, where the inhibition ends up matching the opposite type of excitation only. Should the labels W^P and W^S be identical to the pPE and nPE neurons in this figure? Maybe this is before learning and there are two distinct W^P and two W^S values which end up being (1,0) and (0,1) in this toy example, but as shown it is confusing. Perhaps, as in Fig 1, show a before and after so it is clear what is achieved? And then perhaps show what happens to the connections to the nPE unit so the symmetry is clear. I think I eventually understood it, but the idea that we are just selecting two extreme cells from the range implicit in Fig 1E/F would help, and then focusing only on the inputs to one of those two cells in Fig 2.\nFig 4A, lower panel: why are the prediction inputs and expected stimulus inputs not matched by the strength of opposing inhibitory inputs? The different variations in line widths indicating different strengths of inputs that is not present for the Probe is confusing.\np.14: the focus on NT “and CT mice … in which movement and visual flow are not coupled”. A bit more explanation is needed here as in the CT condition, what exactly is coupled? i.e., what is the difference between CT and NT if not a change in the coupling between movement and visual flow? I assume it is something like “not naturally coupled and rendered uncorrelated in the NT condition but artificially coupled in the CT condition?”\np.14 “an unsigned prediction error” could be explained (signal proportional to absolute value of the error) if only in parentheses, but perhaps compared/contrasted in a paragraph with the RPE literature used in ML.\np.15 “these cumbersome assumptions” – I think the word “cumbersome” sounds rather too negative when addressing a model that has for example perhaps/arguably some more biologically realistic features (inhibitory synapses on the soma, rather than only on the dendrites as in the model presented here – or inclusion of somatic conductance that does not change the somatic membrane potential). And of course, biological cells are pretty cumbersome in their number of features that may warrant modeling.\np.16 “perturbing I^S”: I presume here it is “increasing I^S” not just any perturbation? Or perhaps you mean the spatial shuffle here, not a shift in magnitude?", "reviewer_3": "Meng and Wang study an important question in neuroscience: what learning rules acting with locally-available information at individual synapses can improve system-wide performance on sensory and cognitive tasks. They focus on prediction: the system aims to predict what will happen next, and learns from mismatches between those predictions and incoming stimuli. They derive a simple set of learning rules that rely on global error signals (presumably communicated to all the synapses by the release of neuropeptides) as a \"third factor\" that interacts with pre- and post-synaptic activity to drive changes in synaptic strength. They demonstrate the efficacy of this rule in simple settings with scalar stimuli and predictions and show that their learning rules succeed in causing the networks to learn. They also show that perturbations to the learned networks reduce performance. Finally, they show that there is a bimodal distribution in correlations between neurons in the learned model (reflected the positive and negative prediction error populations), which is recapitulated by analysis of data from the Keller lab.\nOverall, there is lots to like about this paper. I am, however, disappointed by the fact that the model only considers scalar stimuli and predictions. (Please correct me if I am mistaken about this fact.) Thatt simplification makes good senses as a starting point for an initial analysis. However, if the model is to explain biological learning, it must work in the context of right high-dimensional stimuli (e.g., natural scenes) and models that work in the scalar stimulus case often fail in this relevant case of high-dimensional stimuli. I have other, smaller suggestions for the authors. But my biggest point is that they should show that their learning rules work (and reproduce experimental data) in the relevant naturalistic conditions.\nMinor points:\n1\\ Please include your learning rules near the start of the Results section. Those rules are a core result of your work and they are needed in order to understand the rest of the results you show. I understand the temptation to keep equations out of the main text (to make it more accessible to non-mathematical folks), but in this case the paper is much less readable as a result of that choice.\n2\\ The authors should do a more thorough search through the theory literature. For example, Machens, Deneve and others published a somewhat-similar (though not identical) predictive coding idea about a decade ago (cited below) that also succeeded in explaining the random-looking spike trains emitted by cortical neurons. And more recently, many papers from the Bogacz lab have implemented predictive coding with biologically plausible learning rules. The authors should clarify how their work advances the field beyond the previous state of the art.\nBoerlin, M., Machens, C.K. and Denève, S., 2013. Predictive coding of dynamical variables in balanced spiking networks. PLoS computational biology, 9(11), p.e1003258.", "reviewer_4": "This manuscript demonstrates that a three-factor learning rule first introduced by Gerstner and colleagues naturally leads to differentiated patterns of prediction error in neural networks. Specifically, this learning rule when acting on synapses from inhibitory to excitatory neurons leads to emergence of two kinds of neurons: positive prediction neurons that receive stimulus input minus prediction, and negative prediction neurons that receive prediction minus stimulus. These bimodal distribution of responses are important in neural circuits given the rectified properties of neural responses that should be able to encode both the cases where stimulus arrives without the expected reward and where reward arrives unexpectedly. The results are technically valid, but the writing of the paper can be made more exciting to emphasize the impact of the findings. Perhaps, the authors can make predictions for the relative abundances of positive and negative prediction error neurons and compare these to neural data and whether the relative number of these neurons is expected to change with learning. For example, one imagine that situations where negative prediction errors occur should be less common and this may be reflected in the relative numbers of these neurons after learning."}, "author_response": "The authors' response to the reviewers is provided in the text following the reviewer reports in the original document, which includes detailed point-by-point rebuttals to each reviewer's comments."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for addressing my questions.\nI still think that the theory exactly applies only for N=2 and in the case of identical tuning curves.\nAs a result, the theory is not complete, even though the authors use very simple relu models and no temporal dynamics in the theory part.\nAs the authors acknowledge, there is not much realism in the learning rule requiring a specific form for c(t). Although approximations are shown, the theory do not really tell what approximations can be made and when they would fail.\nI still think that this paper is interesting, but I also think that it is appropriate for a more technical audience after providing a full theory covering the missing parts described above.", "reviewer_2": "I think that with all the revisions in place, this is now a very nice paper.", "reviewer_3": "The authors did a thorough job of responding to my critiques and those of the other reviewers (as far as I can tell). The inclusion of 3-D stimuli is a very big improvement to the paper. I would have liked to see an extension to arbitrary number of stimulus dimensions, but that could be left for a future study.", "reviewer_4": "The manuscript was marginally improved but not in a substantive way. The message remains rather incremental and feels better suited for a theory journal. The introduction about sensory-motor mismatch is placed. The problem considered here -- that of predictive coding -- is more general than the sensory-motor mismatch, and there are no specific predictions for the sensory-motor mismatch. The advance here is to implement learning with two kinds of neurons, as in on and off neurons in the sensory coding."}, "author_response": "The authors provided a final set of responses to the reviewers' second-round comments, addressing the remaining concerns regarding the theoretical completeness and the scope of the study."}]

Decision letter