Neuron modelPoint neuron

In the simulations with a postsynaptic neuron described by a single variable (point neuron), we implemented a LIF neuron with after-hyperpolarization (AHP) current and conductance-based synapses. The postsynaptic neuron’s membrane potential, u(t), evolved according to a first-order differential equation:

$$\begin_}}}}\frac}}u(t)}}}t}&=&-[u(t)-_}}}}]-_}}}}(t)[u(t)-_}}}}]+R_}}}}(t)\\ &&-_}}}}(t)[u(t)-_}}}}]-_}}}}_}}}}}(t)[u(t)-_}}}}_}}}}}]\\ &&-_}}}}(t)_}}}}(u(t))[u(t)-_}}}}],\end$$

(3)

where τm is the membrane time constant (τm = RC; leak resistance × membrane capacitance); urest is the resting membrane potential; gAHP(t) is the conductance of the AHP channel with reversal potential EAHP; Iext(t) is an external current used to mimic experimental protocols to induce excitatory plasticity; and gX(t) and EX are the conductance and the reversal potential of the synaptic channel X, respectively, with X = . Excitatory NMDA channels were implemented with a nonlinear function of the membrane potential, caused by a Mg2+ block, whose effect was simulated by the function:

$$_}}}}(u)=_}}}}\exp [_}}}}(u-_}}}})]\right)}^,$$

(4)

where aNMDA and bNMDA are parameters50. The AHP conductance was modeled as:

$$\frac}}_}}}}(t)}}}t}=-\frac_}}}}(t)}_}}}}}+_}}}}_}}}}(t),$$

(5)

where τAHP is the characteristic time of the AHP channel; AAHP is the amplitude of increase in conductance due to a single postsynaptic spike; and Spost(t) is the spike train of the postsynaptic neuron:

$$_}}}}(t)=\mathop\limits_\delta (t-_}}}}^),$$

(6)

where \(_}}}}^\) is the time of the kth spike of the postsynaptic neuron, and δ( ⋅ ) is the Dirac’s delta. The synaptic conductance was modeled as:

$$\frac}}_(t)}}}t}=-\frac_(t)}_}+\mathop\limits__(t)_(t),$$

(7)

where τX is the characteristic time of the neuroreceptor X. The sum on the right-hand side of equation (7) corresponds to presynaptic spike trains weighted by the synaptic strength wj(t). The presynaptic spike train of neuron j was modeled as:

$$_(t)=\mathop\limits_\delta \left(t-_^\right),$$

(8)

where \(_^\) is the time of the kth spike of neuron j. The postsynaptic neuron elicited an action potential whenever the membrane potential crossed a spiking threshold from below. We simulated two types of threshold: fixed or adaptive.

Fixed spiking threshold. A fixed spiking threshold was implemented as a parameter, uth. When the postsynaptic neuron’s membrane potential crossed uth from below, a spike was generated, and the postsynaptic neuron’s membrane potential was instantaneously reset to ureset and then clamped at this value for the duration of the refractory period, τref. All simulations with a single postsynaptic neuron were implemented with a fixed spiking threshold (Figs. 2–6, Extended Data Figs. 2, 3 and 5–7 and Supplementary Figs. 3 and 4), except the simulations in which the action potential was explicitly implemented (Extended Data Fig. 2c,g,k and Supplementary Figs. 2 and 3d; details in the Supplementary Modeling Note).

Adapting spiking threshold. For the simulations of the recurrent network, we used an adapting spiking threshold, uth(t). When the postsynaptic neuron’s membrane potential crossed uth(t) from below, a spike was generated, and the postsynaptic neuron’s membrane potential was instantaneously reset to ureset without any additional clamping of the membrane potential (the refractory period that results from the adapting threshold is calculated below). Upon spike, the adapting spiking threshold, uth(t), was instantaneously set to \(_}}}}^\), decaying back to its baseline according to:

$$_}}}}\frac}}_}}}}(t)}}}t}=-_}}}}(t)+_}}}}^,$$

(9)

where τth is the decaying time for the spiking threshold variable, and \(_}}}}^\) is the baseline for spike generation. The maximum depolarization of the membrane potential is linked to the reversal potential of NMDA, and, thus, the absolute refractory period can be calculated as:

$$_}}}}=_}}}}\ln \left(\frac_}}}}^-_}}}}^}_}}}}-_}}}}^}\right),$$

(10)

which is the time the adapting threshold takes to decay to the same value as the reversal potential of the NMDA channels.

Two-layer neuron

The two-layer neuron was simulated as a compartmental model with a spiking soma that receives input from NB dendritic branches. The soma was modeled as a LIF neuron and the dendrite as a leaky integrator (without generation of action potentials). Somatic membrane potential evolved according to:

$$\begin_}}}}\frac}}_}}}}(t)}}}t}&=&-[_}}}}(t)-_}}}}]-_}}}}(t)[_}}}}(t)-_}}}}]\\ &&-\mathop\limits_^_}}}}}_[_}}}}(t)-_(t)].\end$$

(11)

The soma of the two-layer neuron was similar to the point neuron (equation (3)); however, synaptic currents were injected on the dendritic tree, which interacted with the soma passively through the last term on the right-hand side of equation (11), Ji being the conductance that controls the current flow due to connection between the soma and the ith dendrite. In equation (11), ui(t) is the membrane potential of the dendritic branch i. When the somatic membrane potential, usoma(t), crossed the threshold, uth, from below, the postsynaptic neuron generated an action potential, being instantaneously reset to ureset and then clamped at this value for the duration of the refractory period, τref.

Dendritic compartments received presynaptic inputs as well as a sink current from the soma. The membrane potential of the ith branch, ui(t), evolved according to the following differential equation:

$$\begin_}}}}\frac}}_(t)}}}t}&=&-[_(t)-_}}}}]-_[_(t)-_}}}}(t)]\\ &&-_}}},i}(t)[_(t)-_}}}}]\\ &&-_}}}}_}}}},i}(t)[_(t)-_}}}}_}}}}}]\\ &&-_}}},i}(t)_}}}}(_(t))[_(t)-_}}}}].\end$$

(12)

Spikes were not elicited in dendritic compartments, but, due to the gating function HNMDA(u) and the absence of spiking threshold, voltage plateaus occurred naturally when multiple inputs arrived simultaneously on a compartment (Fig. 6a). We simulated two compartments (NB = 2) with the same coupling with the soma, Ji: one whose synapses changed according to the co-dependent synaptic plasticity model and one with fixed synapses that acted as a noise source.

Coupling strength as function of electrotonic distance. The crucial parameter introduced when including dendritic compartments was the coupling, Ji, between soma and the dendritic compartment i. Steady changes in membrane potential at the soma are attenuated at dendritic compartments, and this attenuation has been shown to decrease with distance. Without synaptic inputs and steady membrane potential at both soma and dendritic compartments, equations (11) and (12) are equal to zero, which results in:

where ai is the passive dendritic attenuation of the dendritic compartment i,

$$_=\frac}_-_}}}}}}_}}}}-_}}}}},$$

(14)

with \(}_}}}}\) being a constant steady state held at the soma and \(}_\) being the resulting steady state at the dendritic compartment i. The coupling between soma and the dendritic compartment i is a function of distance as follows:

$$_=_(d)=\frac_^}^},$$

(15)

where d* is a parameter that we fitted from experimental data from ref. 45 (Fig. 6c). We used this fitted parameter to approximate the distance to the soma in Fig. 6f and Extended Data Figs. 6 and 7 according to the soma–dendrite coupling strength used in our simulations.

Co-dependent synaptic plasticity model

The co-dependent plasticity model is a function on both spike times and input currents. We first describe how synaptic currents are accounted and then how excitatory and inhibitory plasticity models were implemented. We defined a variable Ej(t) to represent the process triggered by excitatory currents that influence plasticity at the synapse connecting a presynaptic neuron j to the postsynaptic neuron. We considered NMDA currents, which reflect influx of calcium into the postsynaptic cell, as the trigger for biochemical processes that are represented by the state of Ej(t). Its dynamics are described by the weighted sum (Gaussian envelope) of the synapse-specific filtered NMDA current, \(}_(t)\),

$$_(t)=\mathop\limits__^}}}}(\;j,k)}_(t),$$

(16)

where \(_^}}}}(j,k)\) is the function describing the effect of synapse k in the plasticity of synapse j (based on physical distance considering that both synapses are connected onto the same postsynaptic neuron; details below). The synapse-specific filtered NMDA current dynamics are given by:

$$_}}}}\frac}}}_(t)}}}t}=-}_(t)-_}}},j}(t)_}}}}(u(t))\left[u(t)-_}}}}\right],$$

(17)

where τE is the characteristic time of the excitatory trace; u(t) is the postsynaptic membrane potential (dendritic membrane potential for the two-layer neuron model); and gNMDA,j(t) is the conductance of the jth excitatory synapse connected onto the postsynaptic neuron, with dynamics given by:

$$\frac}}_}}},j}(t)}}}t}=-\frac_}}},j}(t)}_}}}}}+_(t)_(t).$$

(18)

Inhibitory inputs contributed to the plasticity model through a variable I(t). For the inhibitory trace, we used GABAA currents, which reflect influx of chloride, as the trigger of the process described by I(t). The inhibitory trace evolved as:

$$_}}}}\frac}}I(t)}}}t}=-I(t)+\mathop\limits__}}}}_}}}},k}(t)\left[u(t)-_}}}}_}}}}}\right],$$

(19)

where τI is the characteristic time of the inhibitory trace, and \(_}}}}_}}}},k}(t)\) is the conductance of the kth inhibitory synapse connected onto the postsynaptic neuron (or dendritic compartment) described as:

$$\frac}}_}}}}_}}}},k}(t)}}}t}=-\frac_}}}}_}}}},k}(t)}_}}}}_}}}}}}+_(t)_(t).$$

(20)

Notice that both Ej(t) and I(t) are in units of voltage because the conductance is unit free in our neuron model implementation (equation (3)).

Influence of distance between synapses

To incorporate distance-dependent influence of the activation of a synapse’s neighbors onto excitatory plasticity, we implemented the function \(_^}}}}(i,j)\) in equation (16). For simplicity, we considered that the amplitude of the distance-dependent influence decays with Gaussian-like shaped function of the synapses’ distance:

$$_^}}}}(i,j)=\exp \left[-\frac\right)}^\right]_}}}}}\mathop\limits_\exp \left[-\frac\right)}^\right]\right\}}^,$$

(21)

where NE is the number of excitatory synapses; i is the index of synapse undergoing plasticity; and j is the index of the its neighboring synapse, including j = i so that the strongest effect is the influx of the excitatory current by the synapse undergoing plasticity. In equation (21), the term Δx(i, j) is the electrotonic distance between synapses j and i, and the parameter σ is the characteristic distance (that is, standard deviation) of the contribution of excitatory synapses for the variable Ej(t). The term inside curly brackets on the right-hand side of equation (21) is a normalizing constant.

The sum of the co-dependent variables Ej(t) for a postsynaptic neuron based on the synapse-specific filtered NMDA currents, \(}_(t)\), can be written as:

$$\begin\mathop\limits__(t)&=&\mathop\limits_\mathop\limits_}_(t)\exp \left[-\frac\right)}^\right]_}}}}}\mathop\limits_\exp \left[-\frac\right)}^\right]\right\}}^\\ &=&\mathop\limits_}_(t)\mathop\limits_\exp \left[-\frac\right)}^\right]_}}}}}\mathop\limits_\exp \left[-\frac\right)}^\right]\right\}}^\\ &\approx &_}}}}\mathop\limits_}_(t),\,}}\,\,_}}}}\gg 1.\end$$

(22)

With the normalization used in equation (21), the average of the variable Ej(t) is approximately equal to the total synapse-specific filtered NMDA currents, \(}_(t)\) (equation (16)), which is independent of σ for a large number of synapses (NE ≫ 1). Notably, for very large σ values (σ ≫ NE), all synapses influence each other’s plasticity equally, so that its implementation can be simplified as:

$$_(t)=\mathop\limits_}_(t),\forall j.$$

(23)

Co-dependent excitatory synaptic plasticity

The co-dependent excitatory synaptic plasticity model is an STDP model regulated by excitatory and inhibitory inputs through Ej(t) and I(t). The weight of the jth synapse onto the the postsynaptic neuron (or dendritic compartment), wj(t), changed according to:

$$\begin\frac}}_(t)}}}t}&=&_}}}}(_(t),I(t);_(t),_}}}}(t))\\ &=&\left\_}}}}_^(t)_(t)-_}}}}_}}}}^(t)_(t)\right)}^\right]_}}}}(t)\right.\\ &&\left.-_}}}}_}}}}^(t)_(t)_(t)\right\}\exp \left[-^}\right)}^\right],\end$$

(24)

where ALTP, Ahet and ALTD are the learning rates of long-term potentiation, heterosynaptic plasticity and long-term depression, respectively. The additional parameter I* defines the level of control that inhibitory activity imposes onto excitatory synapses, with parameter γ defining the shape of the control. Variables Spost(t) and Sj(t) represent the postsynaptic and presynaptic spike trains, respectively, as described above for the neuron model (equations (6) and (8)). The trace of the presynaptic spike train is represented by \(_^(t)\), and the traces of the postsynaptic spike train (with different timescales) are represented by \(_}}}}^(t)\) and \(_}}}}^(t)\). They evolve in time according to:

$$\frac}}_^(t)}}}t}=-\frac_^(t)}_}+_(t),$$

(25)

$$\frac}}_}}}}^(t)}}}t}=-\frac_}}}}^(t)}_}}}}}+_}}}}(t),$$

(26)

and

$$\frac}}_}}}}^(t)}}}t}=-\frac_}}}}^(t)}_}+_}}}}(t).$$

(27)

For values of inhibitory trace larger than a threshold, I(t) > Ith, we effectively blocked excitatory plasticity to mimic complete shunting of backpropagating action potentials58 or additional blocking mechanisms that depend on inhibition23. We implemented maximum and minimum allowed values for excitatory weights, \(_^}}}}=10\) nS and \(_^}}}}=1^\) nS, respectively.

Co-dependent inhibitory synaptic plasticity

Similar to the excitatory learning rule, the co-dependent inhibitory synaptic plasticity is a function of spike times and synaptic currents. The weight of the jth inhibitory synapse onto the postsynaptic neuron (or dendritic compartment), wj(t), changed over time according to a differential equation given by:

$$\begin\frac}}_(t)}}}t}&=&_}}}}(_(t),I(t);_(t),_}}}}(t))\\ &=&_}}}}_(t)\left[_(t)-\alpha I(t)\right]\left[_}}}}(t)_(t)+_(t)_}}}}(t)\right].\end$$

(28)

Parameters AISP and α control the learning rate and the balance of excitatory and inhibitory currents, respectively. Variables xj(t) and ypost(t) are traces of presynaptic and postsynaptic spike trains, respectively, that create a symmetric STDP-like curve, with dynamics given by:

$$\frac}}_}}}}(t)}}}t}=-\frac_}}}}(t)}_}}}}}+_}}}}(t)$$

(29)

and

$$\frac}}_(t)}}}t}=-\frac_(t)}_}}}}}+_(t).$$

(30)

The STDP window is characterized by the time constant τiSTDP. The variable Ej(t) is given by equation (23). We implemented maximum and minimum allowed values for inhibitory weights, \(_^}}}}=70\) nS and \(_^}}}}=1^\) nS, respectively.

Experimental protocols: Fig. 2b,c,e,f,h–j and Extended Data Fig. 2d–k

We fitted three datasets with the co-dependent excitatory synaptic plasticity model to asses its dependency on voltage—that is, membrane potential (Fig. 2b)—on the frequency of presynaptic and postsynaptic spikes (Fig. 2c) and on the effect of co-induction of LTP at neighboring synapses (Fig. 2h–i).

Voltage-dependent STDP protocol. Following the original experiments16, we simulated five presynaptic and five postsynaptic spikes at 50 Hz, with 10 ms between presynaptic and postsynaptic spike times (pre-before-post; Δt = +10 ms), repeated 15 times with an interval of 10 s in between each pairing (Fig. 2b). The more depolarized the membrane potential, the bigger the effect of the NMDA currents, and, therefore, more LTP was induced. We combined three different ways to depolarize the postsynaptic neuron’s membrane potential: strength of synapse, current clamp and backpropagating action potential (see the Supplementary Modeling Note for details). Postsynaptic spike times were directly implemented in the co-dependent plasticity rule—that is, manually setting the spike times in equation (6), spike times that were also used to generate backpropagating action potentials (Supplementary Fig. 1; see the Supplementary Modeling Note for details). We implemented a parameter sweep on these three quantities (see the Supplementary Modeling Note for details), measuring the average depolarization during the pre-before-post interval of the simulation (200-ms interval starting at the first presynaptic spike in each burst). Due to the multiple ways to depolarize the postsynaptic membrane potential, we plotted a region (instead of a single line) in Fig. 2b indicating the possible weight changes for the same depolarization with the different depolarization methods.

Frequency-dependent STDP protocol. Following the protocol from the original experiments15, we simulated 60 presynaptic and postsynaptic spikes with either Δt = +10 ms (pre-before-post) or Δt = −10 ms interval (post-before-pre) with firing rates between 0.1 Hz and 50 Hz. In the simulations of the frequency-dependent protocol (Fig. 2c), postsynaptic spikes were induced by the injection of a current pulse, Iext(t) = 3 nA, for the duration of 2 ms. For a smooth curve, we incremented presynaptic and postsynaptic firing rates in steps of 0.1 Hz (500 simulations per pairing in total). The increase in presynaptic firing rate caused a bigger accumulation in NMDA currents, which increased LTP (Extended Data Fig. 2a). In the simulations with extra presynaptic partners (Fig. 2e,f and Extended Data Fig. 2d–k), we calculated the average synaptic change over 10 trials to account for the trial-to-trial variability due to the added external Poisson spike trains.

Distance-dependent STDP protocol. In the simulations of the distance-dependent protocol (Fig. 2h–i), postsynaptic spikes were induced by the injection of a current pulse, Iext(t) = 3 nA, for the duration of 2 ms. We simulated 60 presynaptic spikes with inter-spike interval of 500 ms, each followed by three postsynaptic spikes with inter-spike interval of 20 ms. For Fig. 2h, we varied the interval between the presynaptic spike and the first postsynaptic spike in a three-spike burst, defined as Δt. For Fig. 2i, we simulated the above protocol (pre-before-burst) with an interval Δt = 5 ms (‘strong LTP’) in a given synapse, followed by the same protocol with Δt = 35 ms (‘weak LTP’) in a neighboring synapse (Δx = 3 μm and σ = 3.16 μm in equation (21)), varying the interval between the strong and weak LTP inductions. For Fig. 2j, we simulated a similar protocol as the one in Fig. 2i, but we fixed the interval between the strong and weak LTP inductions (90 s) and varied the distance between the synapses.

Fitting. Fitting was done with brute force parameter sweep on four parameters for Fig. 2b,c (each fit with different values): ALTP, Ahet, ALTD and τE. For Fig. 2h–j, a similar brute force parameter sweep on five parameters was performed: ALTP, Ahet, ALTD, τE and σ, with the three plots having the same set of parameters.

Stability

The co-dependent plasticity model has a rich dynamics that involves changes in synaptic weights due to presynaptic and postsynaptic spike times as well as synaptic weight and input currents. In this section, we briefly analyze the fixed points for input currents and synaptic weights for general conditions of inputs and outputs.

Considering each synapse individually, we can write the average change in weights (from equation (24), ignoring inhibitory inputs) as:

$$\begin}}_(t)}}}t}\right\rangle}_&=&\left\langle _}}}}(t)\left[_}}}}_^(t)_(t)-_}}}}_}}}}^(t)_(t))}^\right]\right.\\ &&_}}}}_}}}}^(t)_(t)_(t)\right\rangle }_\\ \end$$

(31)

$$\begin}}_(t)}}}t}\right\rangle}_&=&_}}}}_^(t)_(t)_}}}}(t)\right\rangle }_\\ &&-_}}}}_}}}}(t)_}}}}^(t)_(t))}^\right\rangle }_\\ &&-_}}}}_(t)_}}}}^(t)_(t)\right\rangle }_\\ \end$$

(32)

$$\begin}}_(t)}}}t}\right\rangle}_&=&_}}}}_^(t)_(t)_}}}}(t)\right\rangle }_\\ &&-_}}}}_}}}}(t)\right\rangle }__}}}}^(t)\right\rangle }__(t))}^\right\rangle }_\\ &&-_}}}}_(t)\right\rangle }__}}}}^(t)\right\rangle }__(t)\right\rangle }_,\end$$

(33)

where 〈⋅〉t is the average over a time window bigger than the timescale of the quantities involved. In equation (33), we took into consideration that presynaptic spike times are not influenced by postsynaptic activity, and, thus, the average of the products in the last term on the right-hand side of equation (32) is the equal to the product of the averages. Additionally, we assumed no strong correlations between Ej(t) and Spost(t) due to the small fluctuations of the variable Ej(t). Correlations between presynaptic and postsynaptic spikes govern the LTP term and, thus, cannot be ignored. They also depend on the neuron model and amount of inhibition a neuron (or compartment) receives. We can conclude from equation (33) that the weights from silent presynaptic neurons will vanish due to the heterosynaptic term. In our model, these weights can vanish only in moments of disinhibition, when the inhibitory control over excitatory plasticity is minimum.

For our analysis, we consider that all neurons of the network have nearly stationary firing rates without strong fluctuations. Therefore, the spike trains can be rewritten as average firing rates:

$$_(t)\rangle }_=_,$$

(34)

and the traces from the spike trains become:

$$_^(t)\rangle }_=__,$$

(35)

where νj is the average firing rate of neuron j. The same is valid for the postsynaptic neuron’s firing rate as well as all other traces.

We consider the outcome of the excitatory plasticity rule when LTD is not present, ALTD = 0, which informs us on steady state for excitatory currents as a competition between LTP and heterosynaptic plasticity only. Moreover, we assume that the postsynaptic firing rate, νpost, is proportional to the total NMDA current:

$$_}}}}=^\left(\frac^}\mathop\limits_}_+1\right)-\frac_}}}}\rangle \langle _}}}}\rangle }_}}}}^},$$

(36)

where 〈νI〉 and 〈wI〉 are the population average firing rate and weight of inhibitory afferents, respectively, and ν *, E * and \(_}}}}^\) are parameters that depend on the neuron model (see the Supplementary Modeling Note for details). In this case, the steady state of the system is given by:

$$\begin\limits_}_\right\vert }__}}}} = 0}&=&\frac^}\left(1-\frac_}}}}\rangle \langle _}}}}\rangle }^_}}}}^}\right)\\ &&+\sqrt^}\left(1-\frac_}}}}\rangle \langle _}}}}\rangle }^_}}}}^}\right)\right]}^+\frac_}}}}\langle _\rangle _^}_}}}}_}}}}^}},\end$$

(37)

This is also the maximum value for excitatory currents for when LTD is present, as LTD can only decrease synaptic weights. To arrive in equation (37), we set equation (33) to zero and summed over j assuming weak correlations between presynaptic and postsynaptic spikes so that \(_^(t)_(t)_}}}}(t)\rangle }_=_^(t)\rangle }__(t)\rangle }__}}}}(t)\rangle }_\) (see the Supplementary Modeling Note for details). Notice that this fixed point depends on the presynaptic firing rates and the model parameters. For very low postsynaptic firing rates and weak excitatory weights, assuming two consecutive postsynaptic spikes and, thus, setting \(_}}}}^}}}}=1\) (rather than an average \(\langle _}}}}^}}}}\rangle =_}}}}_}}}}\ll 1\)), we find a threshold for which the learning rate of heterosynaptic plasticity induces vanishing of synapses:

$$_}}}}=\frac_}}}}^\langle _\rangle __}}}}}^\left[1+_}}}}\left(^-\frac_}}}}\rangle \langle _}}}}\rangle }_}}}}^}\right)\right]}.$$

(38)

For a recurrent network, we can assume that νj = νpost and thus:

$$}_^}}}}=\frac_}}}}_}_}}}}_}}}}},\forall j.$$

(39)

Notice that the maximum excitatory current onto a neuron embedded in a recurrent network is independent on firing rate of presynaptic and postsynaptic neurons.

In Fig. 3e–g, we simulated the co-dependent excitatory plasticity model with non-zero ALTP, Ahet and ALTD but without inhibitory control. Each excitatory input was simulated with a constant presynaptic firing rate, 0 < νj < 18 Hz, uniformly distributed, while the firing rate of all presynaptic inhibitory neurons was set to 18 Hz (details below). For each corresponding value in the x axis of Fig. 3e–g, we simulated 40 trials (one point per trial is plotted). We separated these 40 trials into four combinations of the parameters σ and τE (10 trials per parameter set) to confirm the independence of the steady state on these parameters: σ = 10 and τE = 1,000 ms; σ = 1,000 and τE = 10 ms; and σ = 1,000 and τE = 1,000 ms. In Fig. 3e–g, we plotted the theory as equation (37). In Fig. 3e, we plotted the learning rate for which weights may vanish as a dashed vertical line (equation (38)). The parameters from equation (36) were fitted by varying excitatory and inhibitory weights without any plasticity (see the Supplementary Modeling Note for details). Extra postsynaptic spikes were manually added to the plasticity rule implementation (equation (6)) at 1 Hz (Poisson process) to enforce plasticity when excitatory inputs were too weak (compared to inhibitory inputs) to elicit postsynaptic response. To test the effect of input firing rate and LTD with weight dependency, we also simulated a similar protocol (as in Fig. 3e) with different levels of excitatory input (all presynaptic neurons with the same firing rate), LTD and inhibitory gating (Supplementary Fig. 4). These simulations show that the excitatory input levels had minimal effect on the fixed point of excitatory currents.

Applying the same idea to the co-dependent inhibitory synaptic plasticity model, we get the following average dynamics for the jth inhibitory weight:

$$}}_(t)}}}t}\right\rangle }_=_}}}}_(t)\left[_(t)-\alpha I(t)\right]\left[_}}}}(t)_(t)+_(t)_}}}}(t)\right]\right\rangle}_$$

(40)

$$}}_(t)}}}t}\right\rangle }_\approx _}}}}\overline\left[\overline-\alpha \overline\right]\left[2_}}}}__}}}}\right],$$

(41)

where \(\overline=_\), and Ej(t) is the same for every inhibitory synapse connected onto the postsynaptic neuron (equation (23)) so that \(\overline=_(t)\rangle }_,_(t)=_(t),\forall j,k\). From equation (41), we can calculate the steady state for the inhibitory learning rule, which results in the balance between excitation and inhibition given by α:

$$\frac}}=\alpha .$$

(42)

Synaptic changes for simple spike patterns and fixed excitatory and inhibitory input levels

From equation (24) and equation (28), we calculated changes in excitatory and inhibitory synapses for simple spike patterns (Extended Data Fig. 1). We considered fixed excitatory and inhibitory inputs and calculated changes in a given excitatory synapse as:

$$\begin\Delta _}}}}&=&\left[_}}}}\exp \left(-\frac_}}}}}_}\right)E-_}}}}\exp \left(-\frac_}}}}}_}}}}}\right)^\right.\\ &&\left.-_}}}}\exp \left(-\frac_}}}}}_}\right)_\right]\exp \left[-^}\right)}^\right],\end$$

(43)

where ΔtLTP is the interval between presynaptic and postsynaptic spikes (pre-before-post); Δthet is the interval between two consecutive postsynaptic spikes; and ΔtLTD is the interval between postsynaptic and presynaptic spikes (post-before-pre). In a similar fashion, we calculated changes at a given inhibitory synapse as:

$$\Delta _}}}}=_}}}}E\left(E-\alpha I\right)\exp \left(-\frac_}}}}}\right),$$

(44)

where Δt is the interval between presynaptic and postsynaptic spikes, being positive for pre-before-post and negative for post-before-pre spike patterns.

InputsSingle output neuron (feedforward network)

Presynaptic spike trains for single neurons were implemented as follows. A spike of a presynaptic neuron j occurred in a given timestep of duration Δt with probability pj(t) if there was no spike elicited during the refractory period beforehand; \(_^\) for excitatory and \(_^\) for inhibitory inputs, respectively; and zero otherwise. Different simulation paradigms were defined by the input statistics, which are described below.

Constant firing rate. In Figs. 2e,f, 3 and 4, Extended Data Figs. 2d–k and 3 and Supplementary Figs. 2 and 4, presynaptic neurons fired spikes with a constant probability outside the refractory period. For a constant probability pj(t) = pj, the mean firing rate, νj, was therefore:

$$_=\frac__\right)}^_^/\Delta t}.$$

(45)

In Figs. 2e,f and 3c,d, Extended Data Figs. 2d–k and 3c,d and Supplementary Figs. 2 and 4, the firing rate for external neurons is indicated in the captions and legends. In Fig. 3c,d (colored points) and Fig. 3e–g, as well as Extended Data Fig. 3, the probability of external excitatory spikes was synapse specific, uniformly distributed: 0 < pj ⩽ 0.002, whereas the probability of external inhibitory spikes was pj = 0.002, resulting in 0 < νj⪅18.1 Hz and νj ≈ 18.1 Hz, respectively, considering a timestep Δt = 0.1 ms and refractory periods \(_}}}}^}}}}=5\) ms and \(_}}}}^}}}}=2.5\) ms. In Fig. 3c,d (gray points), the probability of external excitatory spikes was pj = 0.001, whereas the probability of external inhibitory spikes was pj = 0.002, resulting in νj ≈ 9 Hz and νj ≈ 18.1 Hz, respectively. In Fig. 4 and Supplementary Fig. 4, the probability of external excitatory and inhibitory spikes was pj = 5 × 10−4 and pj = 10−3 for excitatory and inhibitory afferents, resulting in νj ≈ 4.87 Hz and νj ≈ 9.75 Hz, respectively.

Variable firing rate (pathways). In Figs. 5 and 6, Extended Data Figs. 4–7 and Supplementary Fig. 3, presynaptic neurons fired spikes according to an inhomogeneous Poisson process.

For the receptive field plasticity simulations (Fig. 5, Extended Data Figs. 4 and 5 and Supplementary Fig. 3), we simulated eight input pathways. We defined a pathway as a group of 100 excitatory and 25 inhibitory afferents (spike trains of presynaptic neurons) with two components: a constant background firing rate and a fluctuating firing rate taken from an Ornstein–Uhlenbeck (OU) process as described below. The background firing rate for all 800 excitatory and 200 inhibitory afferents was given by a probability of \(_^}}}}=2\times 1^\) for excitatory and \(_^}}}}=4\times 1^\) for inhibitory afferents, with respective background firing rates of \(_^}}}}\approx 1.98\) Hz and \(_^}}}}\approx 3.96\) Hz for excitatory and inhibitory presynaptic neurons, respectively, considering a timestep Δt = 0.1 ms and refractory periods of \(_^=5\) ms and \(_^=2.5\) ms. The fluctuating firing rate of the pathway μ was created from an OU process. We used an auxiliary variable, yμ(t), that followed stochastic dynamics given by:

$$\frac}}_(t)}}}t}=-\frac_(t)}_}}}}}+_(t),$$

(46)

where τOU is the time constant of the OU process, and ξμ(t) is a random variable drawn from a Gaussian distribution with zero mean and unitary standard deviation. The fluctuating probability was then defined as:

$$_^(t)=^_(t)\right]}_,$$

(47)

where p* = 0.025 is the amplitude of the fluctuations, and [⋅]+ is a rectifying function. The probability of a presynaptic afferent j belonging to pathway μ to spike due to both background and fluctuating firing rate was given by:

$$_(t)=_^(t)+_^}}}}.$$

(48)

In Fig. 5 and Extended Data Fig. 5, we implemented two learning windows: first to learn the initial receptive field profile (Fig. 5b and Extended Data Figs. 5a,d; see Extended Data Fig. 4a) and later to learn the new configuration of the receptive field profile (Fig. 5c and Extended Data Fig. 5b,e; see Extended Data Fig. 4b). During both learning periods, which lasted 700 ms, we set the firing rate of all inhibitory neurons to background firing rate (constant) and the excitatory pathways as follows. During the first 500 ms, we set th