Abstract
The balance between excitatory and inhibitory inputs is a key feature of cortical dynamics. Such a balance is arguably preserved in dendritic branches, yet its underlying mechanism and functional roles remain unknown. In this study, we developed computational models of heterosynaptic spike-timing-dependent plasticity (STDP) to show that the excitatory/inhibitory balance in dendritic branches is robustly achieved through heterosynaptic interactions between excitatory and inhibitory synapses. The model reproduces key features of experimental heterosynaptic STDP well, and provides analytical insights. Furthermore, heterosynaptic STDP explains how the maturation of inhibitory neurons modulates the selectivity of excitatory neurons for binocular matching in the critical period plasticity. The model also provides an alternative explanation for the potential mechanism underlying the somatic detailed balance that is commonly associated with inhibitory STDP. Our results propose heterosynaptic STDP as a critical factor in synaptic organization and the resultant dendritic computation.
SIGNIFICANCE STATEMENT Recent experimental studies reveal that relative differences in spike timings experienced among neighboring glutamatergic and GABAergic synapses on a dendritic branch significantly influences changes in the efficiency of these synapses. This heterosynaptic form of spike-timing-dependent plasticity (STDP) is potentially important for shaping the synaptic organization and computation of neurons, but its functional role remains elusive. Through computational modeling at the parameter regime where previous experimental results are well reproduced, we show that heterosynaptic plasticity serves to finely balance excitatory and inhibitory inputs on the dendrite. Our results suggest a principle of GABA-driven neural circuit formation.
Introduction
Activity-dependent synaptic plasticity is essential for learning. Especially, the spike-timing difference between presynaptic and postsynaptic neurons is a crucial factor for synaptic learning (Bi and Poo, 1998; Caporale and Dan, 2008). Recent experimental results further reveal that the relative differences in spike timings at neighboring synapses on a dendritic branch have a significant influence on changes in synaptic efficiency at these synapses (Tsukada et al., 2005; Hayama et al., 2013; Paille et al., 2013; Bazelot et al., 2015; Oh et al., 2015). In particular, the timing of GABAergic input exerts a great impact on synaptic plasticity at nearby glutamatergic synapses. Similar phenomena have also been observed in biophysical simulations (Cutsuridis, 2011; Bar-Ilan et al., 2013). This heterosynaptic form of spike-timing-dependent plasticity (h-STDP) is potentially important for synaptic organization on the dendritic tree and the resultant dendritic computation (Mel and Schiller, 2004; Branco et al., 2010). However, the functional role of h-STDP remains elusive, partly due to the lack of a simple analytical model.
In the understanding of homosynaptic STDP, simple mathematical formulations of plasticity have played important roles (Gerstner et al., 1996; Song et al., 2000; Vogels et al., 2011). Motivated by these studies, we constructed a mathematical model of h-STDP based on calcium-based synaptic plasticity models (Shouval et al., 2002; Graupner and Brunel, 2012), and then considered the potential functional merits of the heterosynaptic plasticity. The model reproduces several effects of h-STDP that are observed in the hippocampal area CA1 and the striatum of rodents (Hayama et al., 2013; Paille et al., 2013), and provides analytical insights into underlying mechanisms. The model reveals that h-STDP causes a temporally precise balance (i.e., the detailed balance) between the timing of excitatory and inhibitory inputs on a dendritic branch, because of the inhibitory inputs that shunt long-term depression (LTD) at neighboring correlated excitatory synapses. This result suggests that, not only are the number and total current of excitatory/inhibitory synapses balanced at a branch (Liu, 2004; Wilson et al., 2007), but that the temporal input structure is also balanced, as observed in the soma (Dorrn et al., 2010; Froemke, 2015). Moreover, by considering dendritic computation, we demonstrate that the detailed balance is beneficial for detecting changes in input activity. The model also reconciles with the critical period plasticity for binocular matching observed in the V1 of mice (B. S. Wang et al., 2010, 2013), and provides an explanation for how GABA maturation modulates the selectivity of excitatory neurons during development.
Materials and Methods
In this study, we first constructed a model of a dendritic spine, and then based on that model, built models of a dendritic branch and a dendritic tree of a neuron. We also created an analytically tractable model of a spine by reducing the original spine model.
Spine model: dynamics.
Let us first consider the membrane dynamics of a dendritic spine. The membrane potential of a spine is mainly driven by activation of AMPA/NMDA receptors by presynaptic inputs, backpropagation of postsynaptic spikes, leaky currents, and current influx/outflux caused by excitatory/inhibitory synaptic inputs at nearby synapses. Hence, we modeled the membrane dynamics of a spine i with the following differential equation: where gN(ui) = αNui + βN, with αN and βN being constant coefficients. In the equation, ui is the membrane potential of the spine, and τm is the membrane time constant (for definitions of variables, see Table 1). Here, changes in conductance were approximated by current changes. The resting potential was normalized to zero for simplicity. The terms, xiA and xiN represent the glutamate concentration at AMPA and NMDA receptors, respectively. The function gN(ui) represents the voltage dependence of current influx through the NMDA receptors. This positive feedback is enhanced when additional current is provided through backpropagation. As a result, the model reproduces a large depolarization caused by coincident spikes between presynaptic and postsynaptic neurons. In addition, although the AMPA receptor also shows voltage dependence, here we neglected this dependence, as the relative change around the resting potential is small (Lüscher and Malenka, 2012). xiBP is the effect of backpropagation from the soma, and the last two terms of the equation represents heterosynaptic current, which is given as the sum of the inhibitory (excitatory) currents xjI (xjE) at nearby synapses. We defined the sets of nearby inhibitory and excitatory synapses as ΩiI and ΩiE, respectively, and their delays were denoted as dI and dE. The parameter for inhibitory heterosynaptic effect γI is not weight-dependent, because the inhibitory weight was kept constant throughout the paper. In addition, the parameter for excitatory heterosynaptic effect γE was approximated as a constant in Figure 2, E and F, as the synaptic weights changed only slowly, and was set as zero in the rest of simulations.
Each input xiQ (Q = A, N, BP, I, E) is given as the convoluted spikes: where sk represents the spike timing of the kth spike. Although convolution is calculated at the heterosynaptic synapse in the simulation, it does not influence the results, as the exponential decay is linear.
Spine model: plasticity.
We next consider the calcium influx to a spine through NMDA receptors and the voltage-dependent calcium channel (VDCC). For a given membrane potential ui, the calcium concentration at spine i can be written as follows: where gV(ui) = αVui represents calcium influx through VDCC, and gN(ui)xiN(t) is the influx from NMDA. Importantly, in this configuration, the hyperpolarization of the membrane potential through heterosynaptic inhibitory inputs suppresses the Ca2+ influx to the spine, because both gV(ui) and gN(ui) are modeled as monotonically increasing functions of the membrane potential. This is consistent with recent findings indicating that an inhibitory input significantly suppresses Ca2+ transients in both dendrites and spines (Chiu et al., 2013; Marlin and Carter, 2014; Müllner et al., 2015). Note that spine-projecting inhibitory synapses might mediate Ca2+ suppression in long-necked spines (Chiu et al., 2013).
The calcium concentration at the spine is a major indicator of synaptic plasticity, and many studies indicate that a high Ca2+ concentration on a spine typically induces LTP, whereas a low concentration often causes LTD, though the speed of Ca2+ rise is also known to affect the sign of plasticity (Lüscher and Malenka, 2012). Previous modeling studies showed that calcium-based synaptic plasticity models constructed on this principle well replicate various homosynaptic STDP time windows observed in in vitro experiments (Shouval et al., 2002; Graupner and Brunel, 2012). We therefore used this framework for modeling plasticity. In contrast to the previous calcium-based model where a synaptic weight was assumed to be binary (Graupner and Brunel, 2012), here we assumed that a synaptic weight is a continuous variable, and we additionally introduced an interim weight variable to ensure that the learning dynamics of synaptic weights is robust. As shown in Equation 5, the interim weight variable imposes an additional threshold mechanism to prevent minor synaptic modulation from affecting the synaptic weight (Petersen et al., 1998). This interim weight variable represents the approximate concentration of plasticity-related enzymes such as CaMKII or PP1 (Graupner and Brunel, 2007). In the proposed model, the interim weight yi and synaptic weight wi follow: [X]+ is a sign function that returns 1 if X ≥ 0, but returns 0 otherwise. In the model, the neural dynamics is defined in such a way that the somatic potential caused by a presynaptic spike linearly depends on its synaptic weight wi (Table 1). Thus, wi reflects the amplitude of EPSP. Note that, in this model setting, as observed in recent experiments (Gambino et al., 2014), backpropagation is not necessary for LTP if presynaptic inputs arrive when the membrane potential at the spine is well depolarized. The model reproduces various properties of homosynaptic STDP replicated by Graupner and Brunel (2012), because it is an extension of their model. For instance, it is known that the STDP time window depends on the frequency of the pre-post stimulation (Sjöström et al., 2001). We confirmed that this could indeed be observed in our model, by changing the interval of stimulation in the simulation of STDP (Fig. 1C). In addition, our model replicates the dendritic position dependence of STDP (Sjöström and Häusser, 2006) by mimicking dendritic attenuation with a reduced backpropagation (Fig. 1D, black line). Moreover, by increasing the amplitude of presynaptic stimulation, LTP is rescued (Fig. 1D, gray line) as observed in the experiment (Letzkus et al., 2006). Note that Figure 1, C and D, corresponds to Graupner and Brunel (2012), their Figures 4B and 5A, respectively.
As the heterosynaptic interaction in our model is essentially mediated by the voltage change, our model is also related to the voltage-dependent STDP model (Clopath et al., 2010). The key differences are that our model uses the local membrane potential instead of the somatic potential, and better approximates a previously proposed biophysical model of synaptic plasticity (Graupner and Brunel, 2007) than the phenomenological description used in the voltage-dependent STDP model.
Spine model: details.
In the simulation, we set the common parameters as τC = 18.0 ms, τM = 3.0 ms, τN = 15.0 ms, τA = 3.0 ms, τBP = 3.0 ms, τI = 3.0 ms, τE = 6.0 ms, τY = 50 s, dI = 0.0 ms, αN = 1.0, βN = 0.0, αV = 2.0, γA = 1.0, θp = 70, θd = 35, Cd = 1.0, Bp = 0.001, and Bd = 0.0005 (Table 2 shows the definitions and values of the parameters). Note that, due to positive feedback between Equations 1 and 3, the effective timescales of the calcium dynamics and NMDA channels become longer than the given values. In the model of STDP at the striatum, we additionally used γN = 0.05, γBP = 8.0, γI = 5.0, Cp = 2.3, and yth = 250, whereas for the model of Schaffer collateral synapses, we used γN = 0.2, γBP = 8.5, γI = 3.0, Cp = 2.2, yth = 750, dE = 1.0, and γE = 1.0. In the parameter search, the decay time constants were chosen within biologically reasonable ranges (Koch, 1998); αN, γA, Cd, and Bd were fixed at unitary values (i.e., at 1, except Bd, which was scaled to 0.0005), whereas the other parameters were manually tuned. The robustness of the parameter choices was subsequently confirmed numerically (see Fig. 3). Synaptic weight variables {w} were bounded to 0 < w < 500, and initialized at w = woE, which was defined as woE = 100. All other variables were initialized at zero in the simulation. Paired stimulation was applied every second for 100 s, and the synaptic weight changes were calculated from the values 400 s after the end of the stimulation. In the corticostriatal synapse model, the inhibitory spike was presented with the same timing as the presynaptic spike, whereas for the Schaffer collateral synapses, the inhibitory spikes were given 10 ms before the pre (post) spikes in the pre-post (post-pre) stimulation protocols. In the calculation of the interim weight variable y(t) in Figure 2, B, D, and F, we ignored the effect of exponential decay because of the difference in the timescale (τy ≫ 1 s). In the calculation of spike-timing difference, we subtracted 7.5 ms of axonal delay from the timing of presynaptic stimulation. Simulations of the differential equations were implemented using a Runge–Kutta method with a time step of 0.1 ms.
Dendritic hotspot model.
A dendritic hotspot model was constructed based on the Schaffer collateral synapse model described above. For simplicity, we hypothesized that the effect of dendritic geometry is negligible within a dendritic hotspot; hence the heterosynaptic current due to the inhibitory spike arrives at all the nearby excitatory spines at the same time. In addition, we also disregarded the excitatory-to-excitatory (E-to-E) interaction by setting γE = 0.0. Correlated spikes were generated using hidden variables as in previous studies (Vogels et al., 2011; Hiratani and Fukai, 2015). We generated five dynamic hidden variables and updated them at each time step by sμ(t + Δt) = (ζ − )(1 − αs) + sμ(t)αs, where αs = exp[− Δt/τS], τS = 10 ms, μ = 0, 1, …, 4, and ζ is a random variable uniformly chosen from [0, 1). In the simulation, the time step was set to Δt = 0.1 ms. The activities of the presynaptic neurons were generated by a rate-modulated Poisson process with riE(t) = rXE + rSEsμ(t) for excitatory neuron i modulated by the hidden variable μ (due to a non-negative constraint on riE(t), we set riE(t) = 0 when rXE + rSEsμ(t) < 0. Similarly, the presynaptic inhibitory neuron was described by a Poisson model with rI(t) = rXI + rSIs0(t). The activity of the postsynaptic neuron was given as a Poisson model with a fixed rate rpost. We set the parameters {rxE, rsE, rpost} such that all presynaptic and postsynaptic excitatory neurons show the same average firing rate at 5 Hz, to avoid the effect of firing-rate differences on synaptic plasticity.
In Figure 7, to explore the role of dendritic spiking in synaptic plasticity, we introduced the effect of dendritic spikes by changing Equation 1 as follows: where xds obeys , with sik being the kth spike of the ith presynaptic neuron. The sign function [X]+ returns 1 if X ≥ 0, but returns 0 otherwise. With these modifications, the membrane potential of each spine receives additional excitatory current when xds exceeds the threshold θds due to coincident inputs onto nearby spines.
We used γI = 1.2, βN = 1.0, γBP = 8.0, Cp = 2.11, and yth = 250 with other parameters being kept at the same values used in the original Schaffer collateral model (Table 2). The number of excitatory inputs to the branch NbE was set to NbE = 10. Except for Figure 5F, the mean delay of the inhibitory spikes was set to zero. Presynaptic activities were given by rXE = 1.0 Hz, rSE = 500.0, rXI = 2.0 Hz, and rSI = 1000.0 so that the average firing rate of presynaptic neurons became ∼5 Hz, whereas the postsynaptic firing rate was set to rpost = 5.0 Hz. The large values of rSE and rSI were chosen because the correlation factor sμ(t) was typically very small ( ≈ 0.02). In Figure 5C, the correlations were calculated between the dendritic membrane potential gb(ub) and hidden variables {sμ(t)}, where ub(t) ≡ ∑i=1NbEwiui(t)/(woENbE), and gb(u) were defined as gb(u) = u if u > ubo, otherwise gb(u) = ubo with ubo = −5.0. For the dendritic spike model shown in Figure 7, we used γsd = 1.0, τsd = 10.0 ms, and θds = 4woE. The other parameters were kept at the same values used above.
Two-layered neuron model.
Previous studies suggest that the complicated dendritic computation can be approximated by a two-layered single-cell model (Poirazi et al., 2003; London and Häusser, 2005). We therefore constructed a single cell model by assuming that each hotspot works as a subunit of a two-layered model. In this model, a dendritic subunit (i.e., a unit in the first layer of the two-layered model) corresponds to an electrically compartmentalized subregion of the dendritic tree, such as a thin terminal dendrite or a combination of oblique and terminal branches diverted from the main dendritic shaft. By contrast, in the dendritic hotspot model, a hotspot represents a group of excitatory synapses modulated by a common inhibitory input on a subregion of a dendritic branch. Thus, given that an inhibitory input can modulate excitatory synapses up to 10–15 μm away from the input site (Hayama et al., 2013), a dendritic subunit may contain multiple hotspots. However, for simplicity we supposed that a dendritic subunit corresponds to a hotspot in the two-layered model. We defined the mean potential of a dendritic subunit k as ubk(t) ≡ ∑i=1NbEwikuik(t)/(woENbE), and calculated the somatic membrane potential as usoma(t) ≡ ∑k=1Kgb(ubk(t)). Postsynaptic spikes were given as a rate-modulated Poisson model with the rate usoma(t)/Idv(t). The term Idv(t) is a divisive inhibition term introduced to keep the output firing rate at rpost. By using the mean somatic potential , Idv(t) was calculated as Idv(t) ≡ ūsoma(t)/rpost.
In the simulations described in Figure 6, we used Cp = 2.01, τv = 1 s, and K = 100, with the other parameters being kept at the same values used in the dendritic hotspot model. During the learning depicted in Figure 6B–E, we used the same input configuration as in the dendritic hotspot model. In Figure 6F–H, the activity level of the hidden variables {sμ(t)} was kept at a constant value sμ(t)= 0.25 during the 500 ms stimulation, whereas it was otherwise kept at zero. Additionally, the inhibitory presynaptic activities were set to rXI = 10 Hz and rSI = 2000. For the data presented in Figure 6, G and H, we modulated the firing rates of both the excitatory and inhibitory presynaptic neurons by changing the activity levels of the hidden variables {sμ(t)} from 0.1 to 0.5. The ratio of the change detecting spikes was defined as the ratio of the spikes occurring within 50 ms of the change to the total spike count.
In Figure 6, C and E, the standard STDP model (Song et al., 2000; Hiratani and Fukai, 2015) was implemented as follows: where Δt is the spike-timing difference between the post and presynaptic spikes (i.e., pre-post is LTP). Here, we introduced a 3 ms dendritic delay into the calculation of Δt. In addition, to induce branch-specific competition, we performed normalization wi = w̃i/∑i′∈branch(w̃i′/wo) at each dendritic branch at every time step. Neural dynamics was kept at the same dynamics used in the two-layered model described above. In the simulations presented, we used τp = 17 ms, τd = 34 ms, Ap = 1.0, Ad = 0.5, and ηstdp = 3.0.
The model of binocular matching.
For the model of the critical period plasticity of binocular matching depicted in Figure 8, we also used a two-layered single cell model. The neuron model has K = 100 dendritic branches, each receives NbE = 20 excitatory inputs and 1 inhibitory input. At each branch, half of the excitatory inputs are from the contralateral eye, and the other half are from the ipsilateral eye. Each excitatory input neuron has direction selectivity characterized by θk,iE, and shows rate-modulated Poisson firing with: where θ(t) is the direction of the visual stimulus at time t, Q is either contralateral or ipsilateral, and I0(βE) is the modified Bessel function of order 0. Similarly, the firing rate of an inhibitory neuron is given as rkI(t) = rxI exp [βI cos (θ(t) − θkI)]/I0(βI). For each excitatory input neuron, the mean direction selectivity {θk,iQ} was randomly chosen from a von Mises distribution exp [βS cos (θk,iQ − θQ)]/2πI0(βS), where Q = {contra, ipsi}. In the simulation, we used θcontra = −π/4, and θipsi = π/4. Correspondingly, the mean direction selectivity of an inhibitory neuron {θkI} was defined as the mean of its selectivity for ipsilateral and contralateral inputs (i.e., θkI = (θkI,ipsi + θkI,contra), where θkI,ipsi and θkI,contra were also randomly depicted from exp[βS cos(θkQ − θQ)]/2πI0(βS). The direction of the visual stimulus θ(t) changes randomly with θ(t + Δt) = θ(t) + σsrζG, where ζG is a Gaussian random variable and Δt is the time step of the simulation. To mimic monocular deprivation, in the shadowed area of Figure 8E, we replaced the contralateral-driven input neuron activity with a Poisson spiking having a constant firing rate of rmdE. In addition, to simulate the lack of contralateral-driven inputs to inhibitory neurons, we replaced the inhibitory activity with rkI(t) = rmdI + (rxI/2)exp[βI cos (θ(t) − θkI,ipsi)]/I0(βI). Similarly, in the firing response shown in Figure 8C, we measured direction selectivity by providing monocular inputs, while replacing the inputs from the other eye with homogeneous Poisson spikes with a firing rate of rmdE.
To evaluate the development of binocular matching we introduced three order parameters. First, the difference between the mean excitatory direction selectivity and the inhibitory selectivity at a branch k was evaluated by θb,kd = |arg(∑i wk,iEei(θk,iE−θkI))|. Similarly, the global direction selectivity difference between the inputs from the ipsilateral and contralateral eyes was defined by the following: where the function d̂[θ1, θ2] calculates the phase difference between the two angles. Finally, the direction selectivity index (DSI) for binocular input was calculated by: For the calculation of the monocular DSI, at each branch k, we took the sum over NbE/2 excitatory inputs corresponding to each eye, instead of all the NbE inputs.
In the simulation, we set γI = 2.5, Cp = 1.85, yth = 750.0, and ubo = 0.0, with the rest of parameters being kept at the values used in the dendritic hotspot model. The inputs parameters were set to βE = 4.0, βI = 2.0, βS = 1.0, θcontra = −π/4, θipsi = π/4, rXE = 5.0, rXI = 10.0, rmdE = 1.0, rmdI = 1.0, and σsr = 0.1 .
Reduced model.
If we shrink equations for membrane potential (Eq. 1) and calcium concentration (Eq. 3) into one, the reduced equation would be written as follows: where gc(X) = [X]+ηX captures the nonlinear effect caused by the pre-post coincidence [i.e., gc(X) returns ηX if X > 0, otherwise returns 0]. All inputs Xi, Xpost, XjI, and XjE were given as point processes, and dI and dE are heterosynaptic delays. The variable gc was calculated from the value of Ci at t = t − Δt to avoid pathological divergence due to the point processes. In the simulation, we simply used the value of Ci from the previous time step. For the interim weight y, we used the same equation as before. Note that Equation 11 is basically the same as the one by Graupner and Brunel (2012), except for the nonlinear term gc(C) and the heterosynaptic terms.
Let us consider the weight dynamics of an excitatory synapse that has only one inhibitory synapse in its neighbor. For analytical tractability, we consider the case when presynaptic, postsynaptic, and inhibitory neurons fire only one spikes at t = tpre, tpost, and tI, respectively. In the case of the CA1 experiment, because the GABA uncaging was always performed before the presynaptic and postsynaptic spike, the timing of the inhibitory spike is given as tI = min(tpre, tpost) − δI for δI > 0. In this setting, the change in the interim weight variable of the excitatory synapse is given as follows: where Similarly, in the case of the striatum experiment, by setting η = 0, the change in the interim weight variable is given as follows: where In the simulation, the parameters were set to τc = 30 ms, Cpost = 2.0, θp = 1.6, θd = 1.0, Bp = 2.25, and Bd = 1.0. Additionally, in the model of a Schaffer collateral synapse, we used δI = 1.0, Cpre = 1.0, CE = 0.30, and η = 2.0, whereas for the model of a corticostriatal synapse, we used δI = 5.0, Cpre = 0.75, CE = 0.0, and η = 0.0. In Figures 4, C and D, we used the parameter set for the model of a Schaffer collateral synapse.
As depicted in Figure 4D, the model also provides an analytical insight into the E-to-E interaction, in addition to the inhibitory-to-excitatory (I-to-E) interaction analyzed in the main result. In the E-to-E interaction, neighboring synapses receive small heterosynaptic calcium transient CE, instead of presynaptic input Cpre. We can therefore characterize the shapes of the STDP time windows by the heterosynaptic excitatory effect parameter CE, and the postsynaptic effect parameters Cpost (Fig. 4D). When the postsynaptic effect parameter Cpost satisfies θp < Cpost < θp + CIe−δI/τC, and the heterosynaptic effect parameter CE fulfills CIe−δI/τC < CE < θp, the STDP time window shows Hebbian-type timing dependency (Fig. 4D, top-middle orange region). However, if CE is smaller than CIe−δI/τC while satisfying θp + CIe−δI/τC − Cpost < CE, then the STDP curve becomes LTD dominant (Fig. 4D, top-left green region), as observed in previous experiments (Hayama et al., 2013; Oh et al., 2015). The excitatory heterosynaptic effect CE is expectedly smaller than the inhibitory effect CI, because the inhibitory potential is typically more localized (Gidon and Segev, 2012). Thus, CE < CIe−δI/τC is also expected to hold for small δI, suggesting robust heterosynaptic LTD at neighboring synapses.
Experimental design and statistical analysis.
Parameters used in the simulations are summarized in Table 2. The main simulation codes for the models are available at https://github.com/nhiratani/hstdp.
Results
Calcium-based synaptic plasticity model with current-based heterosynaptic interaction explains h-STDP
We constructed a model of a dendritic spine, as shown in Figure 1A (see Materials and Methods, Spine model). In the model, the membrane potential of the spine u(t) is modulated by the current influx/outflux via the AMPA and NMDA receptors (Fig. 1A, xA and gN(u)xN), spike backpropagation (xBP), and heterosynaptic currents from nearby excitatory and inhibitory synapses (xE and xI, respectively; see Table 1 for the definitions of variables). The calcium concentration in the spine c(t) is controlled through the NMDA receptors and the VDCCs gv(u) (Higley and Sabatini, 2012). Because both NMDA and VDCC are voltage-dependent (Lüscher and Malenka, 2012), the calcium level in the spine is indirectly controlled by presynaptic, postsynaptic, and heterosynaptic activities (Fig. 1B, top and middle). The voltage dependence of NMDA and VDCC [gN(u) and gv(u)] were assumed to be linear for simplicity. This linear assumption may overestimate the effect of a backpropagating action potential on the dendritic NMDA receptors in the overshooting phase, but the effect is expectedly insignificant because in the absence of Ca2+ spike, a backpropagated action potential tends to overshoot less at the dendrite compared with the soma due to dendritic attenuation (London and Häusser, 2005). In addition, although the synaptic input at an inhibitory synapse causes a positive Ca2+ influx by itself (Koch, 1998), activation of nearby inhibitory input hyperpolarizes the membrane potential at the local dendritic site, and may suppress the Ca2+ influx through VDCCs. Indeed, recent experimental results reveal that Ca2+ influx driven by backpropagating spikes or excitatory synaptic inputs is strongly reduced by an inhibitory input in a temporally and spatially precise manner (Hayama et al., 2013; Marlin and Carter, 2014; Müllner et al., 2015). On the basis of these observations, the Ca2+ level in our model is negatively regulated by the heterosynaptic inhibitory inputs through the hyperpolarization of the membrane potential. To model synaptic plasticity, we used a calcium-based model in which LTP/LTD are initiated if the Ca2+ level reaches above the LTP/LTD thresholds (Fig. 1B, middle, orange and cyan lines); this plasticity model reproduces the features of homosynaptic STDP very well (Shouval et al., 2002; Graupner and Brunel, 2012). We introduced the interim weight variable y(t) to capture the non-graded nature of synaptic weight change (Petersen et al., 1998). Thus, changes in Ca2+ level are first embodied in the interim weight y(t) (Fig. 1B, bottom), and are then reflected in the synaptic weight w(t) upon accumulation. The interim weight y(t) is expected to correspond with the concentration of active plasticity-related enzymes such as CaMKII or PP1 (Graupner and Brunel, 2007), and the synaptic weight w(t) reflects the somatic EPSP amplitude. Despite the modification, our model replicates the properties of homosynaptic STDP reproduced by Graupner and Brunel (2012) well (Fig. 1C,D).
We first consider the effect of inhibitory input on synaptic plasticity at nearby excitatory spines. A recent experimental study in a medium spiny neuron (Paille et al., 2013) revealed that synaptic connections from cortical excitatory neurons typically show anti-Hebbian type STDP under a pairwise stimulation protocol, but if the GABA-A receptor is blocked, the STDP time window flips to a Hebbian type STDP (Fig. 2A, circles). Our model can explain this phenomenon in the following way. Let us first consider the case when the presynaptic excitatory input arrives before the postsynaptic spike (the “pre-post” regime). If the GABAergic input is blocked, presynaptic and postsynaptic spikes jointly cause a large membrane depolarization at the excitatory spine. After repetitive stimulation, the calcium concentration rises above the LTP threshold (Fig. 2B, top-right, red line), hence inducing LTP (Fig. 2B, bottom-right, red line). By contrast, if the GABAergic input arrives coincidentally with the presynaptic input, depolarization at the excitatory spine is attenuated by a negative current influx though the inhibitory synapse. As a result, the calcium concentration cannot reach the LTP threshold although it is still high enough to eventually cause LTD (Fig. 2B, right, black lines). Similarly, when the postsynaptic spike arrives at the spine before the presynaptic spike (the post-pre regime) in the absence of GABAergic input, the delayed presynaptic spike causes a slowing of the decay in the level of calcium concentration that may induce LTD (Fig. 2B, left, red lines). To the contrary, if the GABAergic input is provided simultaneously with the presynaptic input, the decay in the calcium concentration is sped up because of the hyperpolarization of the membrane potential at the excitatory spine by the inhibitory input. As a result, LTP is more likely to occur (Fig. 2B, left, black lines). Therefore, when a GABAergic input arrives in coincidence with a presynaptic excitatory input, the STDP time window changes its sign in both the pre-post and the post-pre stimulation regimes (Fig. 2A, lines).
A GABAergic effect on excitatory synaptic plasticity is also observed in CA1 pyramidal neurons (Hayama et al., 2013). In this case, the post-pre stimulation does not induce LTD unless GABA uncaging is conducted near the excitatory spine immediately before the postsynaptic spike arrives at the spine, whereas LTP is induced by the pre-post stimulation regardless of GABA uncaging (Fig. 2C, squares). Our model can also replicate these results. In the pre-post stimulation, the membrane potential of the spine shows strong depolarization due to positive feedback through the NMDA receptor, even if inhibitory current is delivered through GABA (Fig. 2D, top-right, red and black lines). Thus, LTP occurs after repetitive stimulation (Fig. 2D, bottom-right, red and black lines). By contrast, in the post-pre protocol, the effects of triggering LTP and LTD tend to cancel each other in the absence of GABAergic input, whereas LTD becomes dominant under the influence of GABAergic input (Fig. 2D, left, red and black lines, respectively).
In addition to the I-to-E effect, the E-to-E effect is also observed in the case of CA1 pyramidal neurons (Hayama et al., 2013). If GABA uncaging is performed immediately before postsynaptic firing, LTD is also observed in neighboring excitatory spines (Fig. 2E, right, point). This E-to-E heterosynaptic effect is not observed in the absence of GABAergic input (Fig. 2E, left, points). Correspondingly, in the model, excitatory current influx from a nearby synapse causes mild potentiation of calcium concentration in cooperation with inhibitory current influx, eventually inducing LTD (Fig. 2F, left, gray lines). Note that, in this E-to-E effect, interactions of signaling molecules or competition for resources at a later stage of synaptic plasticity may also play a dominant role (Hayama et al., 2013).
To check the parametric robustness of the model, we uniformly sampled the values of all the main parameters from fixed ranges, and studied the sensitivity of the model performance to each parameter, by calculating the performance distribution over the distributions of the other parameters (Fig. 3A,B). Even if the values of the parameters were perturbed by 20–100% from the original values (i.e., the values used in Fig. 2A), we obtained similar sizes of fitting errors to Figure 2A in ∼2.5% of the simulations (Fig. 3A). Moreover, we found that anti-correlation between the VDCC coefficient and the backpropagation coefficient was crucial for replicating the experimental result (Fig. 3B, αV vs γBP and γBP vsαV), as their product determined the effective amplitude of the calcium transient caused by a postsynaptic spike (Eq. 1 + Eq. 3). A linear relationship between the LTP constant and LTD constant is also important for model fitting (Fig. 3B, Cp vs Cd and Cd vs Cp). Notably, fitting the experimental data from the striatum requires a larger coefficient of the heterosynaptic inhibitory effect than the value required for fitting the data from CA1, and thus the striatum model depends on stronger inhibition than the CA1 model (Fig. 3C, top). This is consistent with the observation of strong inhibition in the striatum (Mallet et al., 2005). We also found that faithful reproducing of the CA1 experimental data crucially depended on a high NMDA/AMPA ratio, whereas the striatum model was rather robust against this ratio (Fig. 3C, bottom).
Phase transitions underlying h-STDP
In the previous section, we introduced a biophysical model to establish its relevance to the corresponding biological processes and obtain insight into the underlying mechanism. However, not all of the components of the model are necessary to reproduce the observed properties of h-STDP. We next provide a simple analytically tractable model to investigate the generality of the proposed mechanism.
To this end, we simplify the model to one in which the calcium level at a spine is directly modulated by the presynaptic, postsynaptic, and heterosynaptic activities, as given below: Here, Ci(t) represents the Ca2+ concentration at spine i, Xi, and Xpost represent presynaptic and postsynaptic spikes respectively, dI and dE are heterosynaptic delays, and ΩiI and ΩiE are the sets of neighboring inhibitory and excitatory synapses (for details of the model, see Materials and Methods, Reduced model). Despite its simplicity, the model can qualitatively reproduce the heterosynaptic effects observed in striatal and CA1 neurons, although the quantitative accuracy is degraded (Fig. 4A and B, respectively). Importantly, the reduced model provides further analytical insights into the phenomena.
Let us consider how the inhibitory effect parameter CI controls the I-to-E heterosynaptic effect observed in the CA1 experiment. If we characterize the shape of the STDP time windows by the total number of their local minima/maxima, the parameter space can be divided into several different phases (Fig. 4C). If the LTP threshold θp satisfies Cpre < θp < Cpost, a Hebbian type STDP time window appears when the strength of heterosynaptic inhibitory effect CI satisfies (Cpost − θp)eδI/τC < CI < CpreeδI/τC (Fig. 4C, top, orange region; see Materials and Methods for details of the analysis). Here, we defined δI as the spike-timing difference between the inhibitory spike and the presynaptic or postsynaptic spikes in the pre-post or the post-pre stimulation protocols, respectively. If CI is larger than Cpreexp(δI/τC), a strong inhibitory effect causes LTD, even in the pre-post regime (Fig. 4C, green region), whereas LTD in the post-pre regime is suppressed when CI is smaller than (Cpre − θp) exp (δI/τC) (Fig. 4C, gray region). Thus, the heterosynaptic LTD observed in Figure 2C can be represented as the phase shift from the gray-colored region to the orange-colored region in Figure 4C, which is due to the change in the inhibitory effect CI. This analysis further confirms the condition for inducing heterosynaptic LTD where the heterosynaptic spike-timing difference δI should be smaller than the timescale of Ca2+ dynamics τC (Hayama et al., 2013). This is because δl < τC is necessary for a significant heterosynaptic LTD, and CI is typically smaller than Cpost and θp. In addition, heterosynaptic suppression of the pre-post LTP (green region) is very unlikely to occur, as it is necessary for CI to be larger than Cpreexp(δI/τC). This condition is difficult to satisfy even if δI = 0, because the heterosynaptic effect on Ca2+ dynamics in the spine is expected to be smaller than the homosynaptic effect (i.e., CI < Cpre). A similar analysis is possible for E-to-E interactions, although the phase diagram becomes complicated in this case (Fig. 4D; see Materials and Methods).
These analyses reveal that the heterosynaptic effects are always observable when the parameters of the calcium dynamics fall within a certain region of the parameter space, and underscore the robustness of h-STDP in our framework.
h-STDP induces the detailed dendritic E/I balance at dendritic hotspots
The results described so far suggest that the proposed model gives a good approximation of h-STDP. To investigate the possible functions of h-STDP, we next examined how this h-STDP rule shapes the synaptic organization on the dendrite of a simulated neuron. We first considered a model of a dendritic hotspot (Jia et al., 2010) that receives 10 excitatory inputs and one inhibitory input (Fig. 5A), because the heterosynaptic effect is typically confined within 10 μm of the synapse (Hayama et al., 2013). Excitatory inputs are organized into five pairs, with each pair of excitatory synapses receiving correlated inputs (Fig. 5B; see Materials and Methods, Dendritic hotspot model). In addition, the inhibitory input is correlated with one excitatory pair (Fig. 5A, A1 and A2). Here, we assumed that postsynaptic activity follows a Poisson process with a fixed rate, because the influence of a single hotspot on the soma is usually negligible. In addition, we neglected the effect of morphology and supposed that the heterosynaptic interaction occurs instantaneously within the hotspot. In this configuration, surprisingly, excitatory synapses correlated with the inhibitory input are potentiated, while other excitatory synapses experience minor depression (Fig. 5C, top). As a result, the dendritic membrane potential of the branch becomes less correlated with all the hidden signals, because the strong negative correlation with the blue signal is cancelled by the potentiated excitatory inputs, whereas weak positive correlations with other signals are diminished due to LTD at the corresponding excitatory synapses (Fig. 5C, bottom). Although input spikes to the hotspot is sparse and stochastic, the traces of excitatory and inhibitory currents shift toward the detailed balance after learning (Fig. 5D). A spike correlation between excitatory and inhibitory inputs is crucial for this potentiation of excitatory synapses, but in the model a small correlation is sufficient to produce a significant change in the synaptic weight (Fig. 5E). Moreover, this GABA-driven potentiation is only observable when inhibitory activity is precisely correlated with excitatory activities, and becomes larger when the inhibitory spike precedes excitatory spikes, rather than when it follow them (Fig. 5F). We also found that, when the heterosynaptic inhibitory effect γI is large enough to cause a strong hyperpolarization at nearby synapses, depression is observed at correlated excitatory synapses (Fig. 5G, blue area) instead of potentiation (Fig. 5G, red area). However, as can be seen in Figures 3 and 4, such a large inhibitory effect does not reproduce the STDP experiments, especially the data from the CA1 pyramidal neurons, and is thus unlikely to be observed in the actual brain. These results indicate that h-STDP induces a dendrite-specific temporally precise E/I balance by potentiating excitatory synapses that are correlated with inhibitory synapses.
To reveal the underlying mechanism of this E/I balance generation, we used the simulation data to calculate the probability of the calcium level reaching above the LTD/LTP thresholds after a presynaptic spike. The probabilities of LTP occurrences show similar trajectories after a presynaptic spike, regardless of whether the presynaptic activity is correlated with inhibitory input or not (Fig. 5H, blue and gray dotted lines, respectively). However, the peak probability of LTD occurrence is significantly lower for spines that are correlated with inhibitory inputs (Fig. 5H, blue vs gray solid lines), although in both cases the probability goes up after the presynaptic spike. This asymmetry between LTP and LTD is consistent with the following interpretation: LTD is mainly caused when the presynaptic neuron fires at a low firing rate and the postsynaptic neuron remains silent, both in the experiments (Malenka and Bear, 2004) and in our model (Fig. 5I, gray line). However, if an inhibitory input arrives at a nearby dendrite in coincidence with excitatory activity, the calcium boost caused by the excitatory presynaptic input is attenuated by the heterosynaptic inhibitory effect (Fig. 5I, black line). As a result, LTD is shunted by correlated inhibitory inputs. On the other hand, LTP is mainly caused by coincident presynaptic and postsynaptic spikes, which induce a large increase in calcium that overwhelms the heterosynaptic inhibitory effect. Thus, LTP at correlated excitatory synapses is not compromised by inhibitory activity at a nearby site (Fig. 5J). Therefore, correlated spines tend to be potentiated overall.
To evaluate the generality of the observed dendritic E/I balance, we extended the model to a two-layered single cell (Poirazi et al., 2003) by modeling each branch with one dendritic hotspot (Fig. 6A; see Materials and Methods, Two-layered neuron model), and investigated the dendritic organization of synaptic weight changes by h-STDP. In the simulation, we introduced a 10 ms delay between the excitatory and inhibitory stimulation (Froemke, 2015). Even in this case, when the dendritic branches of a postsynaptic neuron receive inputs from various neurons with different selectivity, each dendritic hotspot shapes its excitatory synaptic organization according to the selectivity of its inhibitory input (Fig. 6B,D; the frame colors in B represent the inhibitory selectivities). As a result, the excitatory synapses on the dendritic tree become clustered, as observed in previous experiments (Kleindienst et al., 2011; Takahashi et al., 2012). Note that, in our model, this clustering of excitatory synapses is caused by common inhibitory inputs, instead of direct interactions between excitatory spines.
To further clarify the importance of heterosynaptic interaction, we next compared the results of h-STDP with that of learning under the standard STDP (Song et al., 2000). In the model, we introduced a branch-specific homeostatic plasticity to induce competition among nearby synapses (see Materials and Methods for the details of the standard STDP model). In the standard STDP model, all dendritic branches developed similar synaptic distributions regardless of the differences in local inhibitory selectivity (Fig. 6C; the beige signal is learned in this example), because one of the hidden signals is captured by chance through self-organization (Song et al., 2000). As a result, under the standard STDP rule the synaptic organizations of branches remain akin to each other (Fig. 6E, gray line). By contrast, under the h-STDP rule, each dendritic branch acquires a synaptic structure according to its local inhibitory input, and individual dendritic branches become dissimilar from one another (Fig. 6E, black line).
We further investigated the possible function of this synaptic organization in information processing. To this end, we consecutively presented the five stimuli to the two-layered neuron model (Fig. 6F). Before learning, the neuron showed an almost constant response to the stimulation, with a small dip at the change points between each of the five stimulations (Fig. 6F, top). By contrast, after learning, the neuron showed transient bursting activity immediately after the onset of the new 500 ms stimulus window, and then rapidly returned to an almost silent state (Fig. 6F, middle). Hence, by h-STDP, a neuron can acquire sensitivity toward abrupt changes in stimuli (Fig. 6F, bottom, G). This sensitivity vanished if the selectivities of inhibitory synapses were randomly shuffled, suggesting the importance of the balance at each dendrite (Fig. 6H). This result indicates that although the detailed E/I balance at each dendrite has a small overall effect to the somatic membrane dynamics; the collective effect from all dendritic branches has a significant impact on the postsynaptic activity.
The effect of dendritic spikes
Previous studies on dendritic computation reveal the potential importance of dendritic spikes in synaptic plasticity (Smith et al., 2013; Kastellakis et al., 2016). Although we mainly consider a regime within a low firing rate (∼5 Hz), dendritic spikes may occur due to strong input spike correlation. We therefore extended the model discussed in the previous section by including dendritic spikes.
Here, we focus on Na+ spikes that are typically localized within a dendritic branch, not global Ca2+ spikes (London and Häusser, 2005), as we are interested in branch-specific synaptic organization. We modeled dendritic spikes as an excitatory heterosynaptic interaction with thresholding, and then added the interaction to the dendritic hotspot model (see Materials and Methods, Dendritic hotspot model). When the amplitude of a dendritic spike is comparable to the amplitude of backpropagating spike, the synapses set their weights differently, depending on the threshold of a dendritic spike. When the spike threshold is very low, the relative weight difference between the correlated and uncorrelated excitatory synapses converges to zero, because all excitatory synapses are potentiated in this regime. In contrast, under conditions of a very high threshold, dendritic spikes rarely occur, so the correlated synapses are only moderately potentiated, as in the control (Fig. 7A; the dotted line represents the control). Notably, when four to five spikes from synapses with the unit weight (w = wo) are sufficient to generate dendritic spikes, the relative weight differences becomes larger than the control, as the dendritic spikes selectively strengthen such synapses that are moderately potentiated by h-STDP (Fig. 7A). This effect is especially significant when inhibitory spikes are delayed by 5–10 ms on average (Fig. 7B). These results suggest that local dendritic spikes stabilize the synaptic weight structure generated through heterosynaptic STDP.
h-STDP explains the critical period plasticity of binocular matching
The results so far indicate that h-STDP induces GABA-driven reorganization of synaptic weights, which in turn enriches dendritic computation such as in the enhanced sensitivity to detect changes in input activity. To investigate its relationship with developmental plasticity, we next consider a model of critical period plasticity in binocular matching (B. S. Wang et al., 2010, 2013). In mice, 1 week after the eye opening, binocular neurons in V1 typically exhibit different orientation selectivity for inputs from the two eyes. Nevertheless, after another 2 weeks, the selective orientations for each eye become closer, and eventually they almost coincide with each other (B. S. Wang et al., 2010). Moreover, this binocular matching is disrupted by accelerating inhibitory maturation (B. S. Wang et al., 2013). Thus, the activity of inhibitory neurons plays a decisive role in shaping binocular matching, in addition to Hebbian plasticity at excitatory synapses.
We modeled this process with the two-layered single cell model introduced in Figure 6 (Fig. 8A, right; see Materials and Methods, The model of binocular matching). The input spike trains were modeled as rate-modulated Poisson processes driven by a circular variable θ, which corresponds to the direction of moving visual stimuli. We assumed the following: (1) inputs from ipsilateral and contralateral eyes already have some weak orientation selectivity at the eye-opening (B. S. Wang et al., 2010; Espinosa and Stryker, 2012), (2) inhibitory cells are driven by both ipsilateral and contralateral eyes (Yazaki-Sugiyama et al., 2009; Kuhlman et al., 2011), and (3) the average orientation selectivity of inhibitory inputs fall between the orientation selectivity for ipsilateral and contralateral excitatory inputs (Fig. 8A, left). This last assumption has not yet been supported by experimental evidence, but if inhibition is provided by neighboring interneurons, these inhibitory neurons are likely to be driven by similar sets of feedforward excitatory inputs to those driving the output neuron. For mathematical convenience, we consider direction selectivity instead of orientation selectivity, but the same argument holds for the latter.
In the simulation, we first ran the process without inhibition, and then introduced GABAergic inputs after a while (Fig. 8B,E, red vertical lines represent the starting points of inhibitory inputs), because maturation of the inhibitory neurons typically occurs in a later stage of the development (Hensch, 2005). Upon the introduction of inhibition, the mean preferred direction selectivity of excitatory synapses in each branch converges to that of the local inhibition, because of heterosynaptic plasticity (Fig. 8B, top; see Materials and Methods for details of evaluation methods), although the synaptic weight development was biased toward the overall direction selectivity of the postsynaptic neuron (Fig. 8D; the bias is toward the zero-degree direction). This dendritic E/I balancing reduces the difference between the direction selectivity of ipsilateral and contralateral inputs on average, because both become closer to the selectivity of the inhibitory input (Fig. 8B, middle). As a result, the binocular direction selectivity is strengthened (Fig. 8B, bottom), and the responses for monocular inputs approximately coincide with each other (Fig. 8C, right). Deprivation of the contralateral inputs immediately after the introduction of inhibition blocks binocular matching (Fig. 8E), in accordance with the experimental data (B. S. Wang et al., 2010).
Precocious GABA maturation has been reported to disrupt binocular matching (B. S. Wang et al., 2013). Our model suggests that the disruption is possibly related to the violation of the third assumption in the model. When the mean inhibitory direction selectivity is substantially different from the ipsilateral and the contralateral direction selectivity (Fig. 8F, at the parameter regions outside of the area surrounded by purple and green lines), h-STDP does not work effectively (Fig. 8F, top), and the difference between ipsilateral and contralateral inputs is not reduced (Fig. 8F, middle). As a result, the binocular direction selectivity is not improved by learning (Fig. 8F, bottom). These results indicate that the rate of maturation of GABA inputs and their effect on h-STDP are an important part of the underlying mechanisms of binocular matching in critical period plasticity.
Discussion
In this study, we first showed that a calcium-based plasticity model robustly captures several characteristics of the plasticity-related interactions between neighboring synapses that occur on a millisecond timescale; this was accomplished by the introduction of heterosynaptic interaction terms (Figs. 2–4). On the basis of this proposed model, we next investigated the possible functions of h-STDP. This study revealed that h-STDP causes the detailed dendritic E/I balance on dendritic hotspots (Figs. 5–7), which is beneficial for detecting changes in input activity (Fig. 6). Furthermore, we found that h-STDP can induce binocular matching upon GABA maturation, and can support an accurate input estimation (Fig. 8).
Experimental predictions
This study provides three experimentally testable predictions. First, our results provide a hypothesis for synaptic organization on the dendritic tree. Excitatory synaptic inputs to a dendritic hotspot often show correlated activities (Kleindienst et al., 2011; Takahashi et al., 2012). Our results indicate that an inhibitory input may also be correlated with excitatory inputs projecting to the nearby dendritic hotspot (Figs. 5, 6), especially on the dendritic tree of an excitatory neuron that is sensitive to changes in the external environment (Figs. 6, 8). Moreover, the model explains why the feature selectivity of these spines shows only a weak similarity, despite their correlations (Jia et al., 2010; Chen et al., 2011). When a synaptic cluster is carved by the heterosynaptic effect of common inhibitory inputs, and not by E-to-E interactions, the variability of feature selectivity within the cluster tends to be large, because inhibitory neurons typically have wider feature selectivity than excitatory neurons (Ma et al., 2010; Moore and Wehr, 2013). In addition, it should also be noted that, E-to-E heterosynaptic LTP is typically induced as a meta-plasticity over a timescale of minutes (Harvey and Svoboda, 2007), which by itself is insufficient to create a correlation-based synaptic cluster.
Second, the results in Figure 5 indicate that LTD at an excitatory synapse is offset by coincident inhibitory inputs to the nearby dendrite. Thus, LTD from low-frequency stimuli (Malenka and Bear, 2004) can be attenuated by coincident GABA uncaging near the stimulated spine. Note that this result would not contradict the previously reported GABA-driven heterosynaptic LTD by paired stimulation, because in that experiment, the excitatory spine was presumably too active to induce LTD in the absence of GABA (Hayama et al., 2013). Indeed, coincident GABAergic inputs may induce heterosynaptic LTD when combined with a moderately high-frequency presynaptic stimulation that in itself does not cause LTD (Blaise and Bronzino, 2003). The model also indicates that correlated inhibitory inputs are likely to suppress LTP at excitatory synapses if the heterosynaptic effect is sufficiently strong (Fig. 5G, blue area). This may be the case for spine-projecting inhibitory synapses.
The third implication of the model concerns the mechanism of binocular matching. Our model indicates that maturation of GABAergic inputs plays a critical role in binocular matching, and proposes a candidate mechanism for disruption of binocular matching by precocious GABA circuit maturation (B. S. Wang et al., 2013; Fig. 8). However, the phenomenon can also be explained by Hebbian plasticity plus some kind of meta-plasticity. If binocular matching is induced purely by Hebbian plasticity and not through a heterosynaptic mechanism, orientation selectivity after the matching should depend solely on the initial orientation selectivity of monocular inputs, assuming that the selectivity of presynaptic neurons remains the same. Especially when the contralateral input activity is larger than the ipsilateral input activity, the resultant orientation selectivity should approximately coincide with the original selectivity of the contralateral input. Alternatively, if the proposed mechanism is engaged during development, the refinement of orientation selectivity should also be influenced by the mean selectivity of the inhibitory input neurons. Thus, long-term imaging of monocular orientation selectivity of binocular neurons in V1 would reveal whether a covariance-based rule is sufficient to explain the phenomena, or whether some other mechanisms, including the one proposed here, also play a major role in the shift of orientation of selectivity.
Carrier of heterosynaptic interaction
Heterosynaptic plasticity has been observed over various spatial and temporal scales with different underlying molecular mechanisms (Nishiyama and Yasuda, 2015). In the case of heterosynaptic interactions involving milliseconds, single-atomic ions are strong candidates, as small molecules such as IP3 are too big to move rapidly from spine to spine (Santamaria et al., 2006). Under the assumption that changes in Ca2+ concentration at an unstimulated spine are crucial for heterosynaptic plasticity, Ca2+ influx/efflux from either intracellular or extracellular sources is necessary for induction of heterosynaptic plasticity. As inhibitory synaptic inputs often change the local Ca2+ concentration in the dendritic branch (Müllner et al., 2015), intracellular spreading of Ca2+ may be a major source of Ca2+ changes in nearby unstimulated spines. At the same time, because inhibitory inputs significantly modulate the membrane voltage of local dendrites (Gidon and Segev, 2012), a synaptic input should strongly drive Ca2+ influx/efflux through NMDA and VDCC from extracellular sources, even at nearby unstimulated spines. Additionally, most of the intracellular calcium ions are bound by calcium buffers (Higley and Sabatini, 2012), and the buffer concentration is also presumably important for induction of synaptic plasticity. In our model, both current-based interactions (Spine model) and calcium-based interactions (Reduced model) replicate the experimental results (Figs. 2 and 4, respectively). Nevertheless, our analytical result suggests that the heterosynaptic Ca2+ change typically needs to be comparable to the homosynaptic change to cause significant heterosynaptic plasticity through calcium-based interaction (Fig. 4C,D). Thus, this study highlights the possible importance of current-based interactions and a spine-specific influx/efflux of extracellular Ca2+ for inducing heterosynaptic plasticity.
Note that heterosynaptic interaction does not need to work on the order of milliseconds to influence the STDP time window. For example, E-to-E heterosynaptic LTD can be initiated by spreading of LTD-related molecules, not by messengers of neural activity (Hayama et al., 2013). Additionally, for a shift in the STDP time window, changes in the ratio of the calcium influx through NMDA and VDCC may play a crucial role (Paille et al., 2013).
Inhibitory cell types
Somatostatin-positive (SOM+) inhibitory neurons typically project to the apical dendrite, have a shorter membrane time constant than the typical timescale of calcium dynamics (Markram et al., 2004; Xu et al., 2013), and often show strong feature selectivity in comparison with other inhibitory neuron types (Ma et al., 2010). Thus, SOM+ is the likely candidate for heterosynaptic STDP. However, our results do not exclude parvalbumin-positive (PV+) inhibitory neurons, which usually have projections to proximal dendrites, and are typically fast-spiking (Markram et al., 2004). In particular, h-STDP through PV+ cells may play an important role in critical period plasticity (Takesian and Hensch, 2013).
Related theoretical studies
Previous theoretical studies show that excitatory heterosynaptic mechanisms such as dendritic spiking generate a functional synaptic clustering on the dendrite (Iannella and Tanaka, 2006; Legenstein and Maass, 2011; Kastellakis et al., 2016), and enriches computational capacity of the neuron (Poirazi and Mel, 2001; Legenstein and Maass, 2011). By contrast, we demonstrated in this study that an inhibitory synapse can induce clustering of nearby excitatory synapses (Figs. 5, 6). A characteristic of this inhibition-based clustering is the involvement of temporally precise activity to produce the detailed balance between local excitatory and inhibitory inputs. For instance, inhibition-based clustering is beneficial in the striatum, where temporally precise activity is crucial for motor coordination, and also in the hippocampus, where place cells exhibit temporally coordinated activity during spatial navigation. However, for modeling of contextual fear conditioning, excitation-based clustering is sufficient (Kastellakis et al., 2016), or potentially desirable, as temporally precise activity is not required for such a task.
For implementing the E/I balance at the soma, inhibitory STDP is a candidate underlying mechanism (Vogels et al., 2011; Kleberg et al., 2014). The proposed h-STDP model can be considered as an alternative explanation for the somatic detailed balance, because when all dendritic branches are balanced, the somatic membrane potential becomes naturally balanced. However, the model has further implications. First, unlike inhibitory STDP, in which inhibitory synapses passively counterbalance a pre-existing excitatory synaptic structure, in h-STDP the inhibitory synapses actively drive plasticity at nearby excitatory synapses (Fig. 5). Moreover, h-STDP enables a nonredundant synaptic weight organization on the dendrites by inducing the E/I balance locally at each dendritic branch (Fig. 6B). Previous synaptic plasticity models are unable to generate such dendritic synaptic weight distributions on their own, as the inhibitory STDP does not guarantee the balance at the dendrites, and the standard excitatory STDP does not support dendritic diversity (Fig. 6C). It is noted that, although h-STDP drives the local synaptic weights on a dendritic branch toward the detailed balance, the convergence to the exact balance is not guaranteed.
Recently, Yang and colleagues showed that an anti E/I balance on the dendritic tree can be beneficial for the gating of synaptic inputs, and suggested that heterosynaptic plasticity could be the underlying mechanism (Yang et al., 2016). In contrast, our model suggests that the anti-E/I balanced state is possible only under conditions of strong heterosynaptic inhibition, which is nonphysiological (Fig. 5G, blue area). This discrepancy may arise from the different definitions of presynaptic selectivity. In their work, the selectivity was defined based on the firing rate, whereas we used the spike correlation for defining the selectivity, as spike correlation presumably drives weight changes in STDP (Song et al., 2000).
Previous biophysical simulation studies reveal that synaptic plasticity at excitatory synapses critically depends on the inhibitory inputs at nearby dendrites (Cutsuridis, 2011; Bar-Ilan et al., 2013; Jedlicka et al., 2015; Wilmes et al., 2016), but these studies do not reveal much on the functional roles of the heterosynaptic plasticity. In particular, although the cancellation of plasticity by shunting inhibition was mentioned by Wilmes et al., (2016), their simple pairwise STDP model does not capture the differential effects of shunting inhibition on LTD and LTP depicted in our model (Fig. 5H–J). On the other hand, network modeling studies have found that heterosynaptic plasticity provides a homeostatic mechanism (Chen et al., 2013; Zenke et al., 2015), but in these models, heterosynaptic plasticity was modeled as a global homeostatic plasticity without any branch specificity, and its advantage over other homeostatic mechanisms was unclear. In this study, by considering intermediate abstraction with analytical but biologically plausible models, we proposed candidate mechanisms for experimental results that have not been modeled before, and revealed potential functions of h-STDP in neural circuit formation.
Future work
Although we fixed the weight of inhibitory synapses in our model to focus on the functions of h-STDP, inhibitory projections on excitatory neurons are known to show plasticity (Hennequin et al., 2017). In particular, a recent experimental study found that inhibitory and excitatory projections on the same postsynaptic cell show correlated weight changes, suggesting a heterosynaptic effect of inhibitory plasticity on excitatory plasticity (L. Wang and Maffei, 2014). This interaction among excitatory, inhibitory, and heterosynaptic plasticity should be studied in detail in the future.
In addition, an increasing number of recent studies indicates the importance of presynaptic changes for synaptic plasticity (Costa et al., 2015, 2017), suggesting the presence of active heterosynaptic plasticity at presynaptic axons. Correspondingly, previous experimental studies found vesicle superpools on axons that potentially regulate vesicle densities at neighboring boutons (Staras et al., 2010). On the other hand, computational studies on heterosynaptic plasticity including the present one are practically limited to changes in the postsynaptic dendrites. Hence, a theory on presynaptic heterosynaptic plasticity is awaited.
Footnotes
This work was partly supported by JSPS Doctorial Fellowship DC2 (N.H.), CREST JST (JPMJCR13W1 to T.F.), and KAKENHI No15H04265 and No16H01289 (T.F.). We thank Dr. Laurent Venance for kindly providing the experimental data and Dr. Yukiko Goda for comments on the paper.
The authors declare no competing financial interests.
- Correspondence should be addressed to Naoki Hiratani, Laboratory for Neural Circuit Theory, RIKEN Brain Science Institute, 2-1 Hirosawa, Wako, Saitama, Japan 351-0198. N.Hiratani{at}gmail.com