Front. Comput. Neurosci.Frontiers in Computational NeuroscienceFront. Comput. Neurosci.16625188Frontiers Media S.A.10.3389/fncom.2019.00026NeuroscienceOriginal ResearchA Theoretical Framework to Derive Simple, FiringRateDependent Mathematical Models of Synaptic PlasticityLappalainenJanne^{1}HerpichJuliane^{1}^{2}TetzlaffChristian^{1}^{2}^{*}^{1}Department of Computational Neuroscience, Third Institute of PhysicsBiophysics, GeorgAugustUniversity, Göttingen, Germany^{2}Bernstein Center for Computational Neuroscience, GeorgAugustUniversity, Göttingen, Germany
Edited by: Nicolas Brunel, Duke University, United States
Reviewed by: Johnatan Aljadeff, Imperial College London, United Kingdom; Joaquín J. Torres, University of Granada, Spain
This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Synaptic plasticity serves as an essential mechanism underlying cognitive processes as learning and memory. For a better understanding detailed theoretical models combine experimental underpinnings of synaptic plasticity and match experimental results. However, these models are mathematically complex impeding the comprehensive investigation of their link to cognitive processes generally executed on the neuronal network level. Here, we derive a mathematical framework enabling the simplification of such detailed models of synaptic plasticity facilitating further mathematical analyses. By this framework we obtain a compact, firingratedependent mathematical formulation, which includes the essential dynamics of the detailed model and, thus, of experimentally verified properties of synaptic plasticity. Amongst others, by testing our framework by abstracting the dynamics of two wellestablished calciumdependent synaptic plasticity models, we derived that the synaptic changes depend on the square of the presynaptic firing rate, which is in contrast to previous assumptions. Thus, the herepresented framework enables the derivation of biologically plausible but simple mathematical models of synaptic plasticity allowing to analyze the underlying dependencies of synaptic dynamics from neuronal properties such as the firing rate and to investigate their implications in complex neuronal networks.
STDPsynaptic plasticitycalciumregression analysislearningactivitydependencyH2020 Future and Emerging Technologies10.13039/100010664Deutsche Forschungsgemeinschaft10.13039/5011000016591. Introduction
Synaptic plasticity serves as an essential mechanism of neuronal networks being linked to diverse functional properties such as the cognitive mechanisms of learning and memory (Hebb, 1949; Martin et al., 2000). Amongst others, to understand the details of this link between synaptic plasticity and functional properties of a neuronal system, theoretical or mathematical models of synaptic plasticity are formulated and investigated (Dayan and Abbott, 2001; Gerstner et al., 2012, 2014). On the one hand, detailed biological models of synaptic plasticity are formulated closely related to experiments, which provide molecular details or synaptic dynamics given diverse stimulation protocols (Earnshaw and Bressloff, 2006; Graupner and Brunel, 2010; Graupner and Brunel, 2012; Antunes et al., 2016; Gallimore et al., 2018). However, to capture the richness of the underlying molecular processes and to match a wide repertoire of experimental findings, the resulting theoretical models become mathematically rather complex and often depend on detailed spiketimings, which impedes further analysis especially on the neuronal network level. On the other hand, to investigate the link between synaptic plasticity and functional properties of neuronal systems, compact, simple theoretical models of synaptic plasticity are developed neglecting the underlying biological details (Tsodyks and Feigelman, 1988; Gerstner and Kistler, 2002; Tetzlaff et al., 2011, 2013; Knoblauch and Sommer, 2016). In contrast to the spikedependency of the detailed models, these compact models usually depend on the average firing rates of neurons (Gerstner and Kistler, 2002). However, the detailed formulation of these compact models are rather phenomenological and they are only loosely linked to the detailed biological models.
Here, we will provide a mathematical method enabling us to derive a compact, firingratedependent theoretical model from a biologically detailed spikebased model. For the detailed model, we consider two versions of a wellestablished calciumbased spiketimingdependent plasticity model (Graupner and Brunel, 2012; Graupner et al., 2016; Li et al., 2016). We will show that this model can be simplified such that the general dynamics induced by it can be welldescribed by a compact, ratedependent model.
Many experimental studies have been performed to investigate the multiplicity of molecular processes, which are involved in synaptic plasticity. Thereby, calcium currents are identified as key players of at least two molecular processes that mediate the change in synaptic strength by synaptic plasticity (Kandel et al., 2000; Lüscher and Malenka, 2012; Korte and Schmitz, 2016). First, inside the cell, calcium triggers different pathways resulting to changes in the number and phosphorylation of AMPA receptors yielding the potentiation (strengthening) and depression (weakening) of synapses (Choquet and Triller, 2013; Huganir and Nicoll, 2013). The second process is triggered by endocannabinoids and induces depression of synapses (Sjöström et al., 2001; Sjöström et al., 2003; Nevian and Sakmann, 2006). The endocannabinoids are transmitted retrogradely from the postsynaptic neuron to the presynaptic neuron changing the neurotransmitter release capability. For this, calcium plays a crucial role (Hashimotodani et al., 2005; Maejima, 2005). These molecular findings are integrated into a theoretical model, in which the postsynaptic calcium concentration within the dendritic spine directly accounts for synaptic plasticity (Graupner and Brunel, 2012). The calcium concentration, in turn, is modulated by spikedependent pre and postsynaptic processes. This calciumdependent model matches several experimental findings and represents a biologically detailed theoretical model of synaptic plasticity.
After introducing our mathematical method to derive compact models from detailed models (see section 2), we simplify the calciumdependent model of synaptic plasticity to derive a compact model with synaptic changes dependent on pre and postsynaptic firing rates. We repeat this procedure for different versions of the considered synaptic plasticity model as well as neuron model. For all cases, the resulting compact model reliably matches the dynamics of the detailed model indicating that, despite the mathematical simplification, the compact model captures the essentials of synaptic dynamics. Further analyses of the resulting compact models reveal, amongst others, that the dynamics of synaptic plasticity are dominated by the square of the presynaptic firing rate, which is in contrast to phenomenologically derived compact theoretical models. Thus, the here developed method enables the derivation of simple, compact, and ratedependent models from detailed ones allowing further investigations on the network level without loss of essential details of the synaptic dynamics.
2. Materials and Methods
In this study, we derive a simplified, compact, firingratedependent mathematical formulation of the detailed dynamics of spiketiming dependent plasticity. For this, we analyze the firing rate dependency of the detailed model and, based on this, estimate the functional relations between activity and synaptic weight changes by regression analysis. As detailed model, we consider a wellestablished calciumbased synaptic plasticity model, which describes a multitude of plasticity effects based on the underlying molecular dynamics of calcium currents in the postsynaptic dendritic spine (Graupner and Brunel, 2012; Graupner et al., 2016). To take the influence of network effects and postsynaptic activity dynamics into account, we focus in our analysis on three fundamental network motifs. These motifs consist of one or two different presynaptic neuronal populations connected to one postsynaptic neuron. In the following, first, we introduce the different considered motifs and the detailed models of synaptic plasticity as well as used neuron models (please see Table 1 for used parameter values). Then, we describe our generic approach to investigate the firing rate dependencies of synaptic plasticity to derive a simplified, compact mathematical description.
Used values for model parameters.
Synaptic plasticity
MATneuron
AEIFneuron
Parameter
Unit
Linear calcium
Nonlinear calcium
Parameter
Unit
Value
Parameter
Unit
Value
τ_{Ca}
ms
22.27212
18.93044
τ_{m}
ms
5
τ_{m}
ms
9.367
C_{pre}
0.84410
0.86467
R
MΩ
50
E_{L}
mV
70.6
C_{post}
1.62138
2.30815
S_{0}
mV
20
Δ_{T}
mV
2
θ_{d}
1
1
α_{1}
mV
30
V_{T}
mV
50.4
θ_{p}
2.009289
4.99780
α_{2}
mV
2
R
MΩ
33.33
γ_{d}
137.7586
111.82515
τ_{1}
ms
10
τ_{z}
ms
144
γ_{p}
597.08922
894.23695
τ_{2}
ms
200
a
nS
4
τ
s
520.76129
707.02258
refr. time
ms
2
b
nA
0.0805
α(N)
mV/s
400 (P2) / 200 (P3)
α(N)
mV/s
170 (P2) / 80.5 (P3)
The values for the different calciumbased synaptic plasticity models are taken from Graupner et al. (2016), for the MATmodel from Kobayashi et al. (2009), and for the AEIFmodel from Brette and Gerstner (2005). Please note that those values are derived from experimental data. The scaling constant α(N) was empirically determined by simulation, such that the maximal postsynaptic firing frequency stays within plausible ranges of activity (maxv(u, w) ≈ 130 Hz).
2.1. Neuronal Setups
Throughout this study we consider three basic neuronal motifs or setups: a presynaptic population with independent postsynaptic neuron (P_{1}; Figure 1, top left), a presynaptic population with dynamic postsynaptic neuron (P_{2}; Figure 1, top middle), and two presynaptic populations with dynamic postsynaptic neuron (P_{3}; Figure 1, top right).
Dynamics of the detailed model of calciumbased synaptic plasticity for three different neuronal setups P_{1}, P_{2}, and P_{3}. The different setups are shown schematically in the first row. The following rows show (second row; A) the average firing rates, (third row; B) the postsynaptic input current, (fourth row; C) the calcium concentration of an individual synapse using the LCmodel (red) and population average (blue) with standard deviation (blue shading), and (fifth row; D) the synaptic weight of an individual synapse (red) and population average (blue) with standard deviation (blue shading). Columns 1 and 2: Dynamics of the P_{1}setup for initial synaptic strengths of w(t_{0}) = 0.6 and pre and postsynaptic spikefiring of u = v = 40Hz (column 1; LTP) and u = 40Hz and v = 5Hz (column 2; LTD). Columns 3 and 4: Same as in columns 1 and 2 for the P_{2}setup with presynaptic spikefiring of u = 60Hz (column 3; LTP) and u = 30Hz (column 4; LTD). Columns 5 and 6: Same as in previous columns for P_{3}setup of two, competing presynaptic populations 1 (blue) and 2 (green) with initial synaptic weights of w_{1}(t_{0}) = 0.9, w_{2}(t_{0}) = 0.6 and presynaptic spikefiring of u_{1} = 5Hz, u_{2} = 60Hz (column 5; no competition) and u_{1} = 5Hz, u_{2} = 80Hz (column 6; strong competition). In all cases the populations consist of N = 1, 000 synapses. The continuous presynaptic activities were obtained by averaging all N spiketrains and temporally filtering with an alpha function window of 100ms width. The postsynaptic activities were obtained with an alpha function window of 500ms width and using the MATmodel. The calcium and current traces were filtered with a moving average of 25ms width.
Neuronal Population I (P_{1}): This motif consists of one population of N = 1, 000 presynaptic neurons (single lines abstracted to one arrow) connected to one postsynaptic neuron (circle). All presynaptic neurons fire independently from each other by a Poisson process with frequency u. The postsynaptic neuron is active by a Poisson process with frequency v. Note, in this setup the firing of the postsynaptic neuron is independent of the activity of the presynaptic population, since all spike trains are generated independently by probabilistic Poisson processes. These processes can yield, at a specific point in time, different synaptic weights for different presynaptic neurons. Thus, we consider the average synaptic strength of the presynaptic population denoted by w.
Neuronal Population I with neuron model (P_{2}): This setup is similar to the P_{1}setup. However, here, the activity of the postsynaptic neuron is not determined by a Poisson process (see P_{1}), instead, we consider a biologically detailed neuron model. For generality, throughout this study we use two different neuron models; the socalled Multitimescale Adaptive Threshold model (MAT; Kobayashi et al., 2009) and the Adaptive Exponential IntegrateandFire model (AEIF; Brette and Gerstner, 2005). Thus, the activity of the presynaptic population can influence the firing of the postsynaptic neuron. Consequently, the change in the postsynaptic firing rate (v˙) depends on u and w, i.e., v˙(u,w), which in turn influences the dynamics of synaptic plasticity (ẇ).
Neuronal Population II (P_{3}): This motif is an extension of the P_{2}setup. Here, two presynaptic populations 1 and 2 are connected to the same postsynaptic neuron with each population consisting of independently firing neurons with populationspecific frequency u_{1} and u_{2}, respectively. Thus, the firing of the postsynaptic neuron depends on the presynaptic firing rates and strength of the connecting synapses of both presynaptic populations (v˙(u1,u2,w1,w2)). Note, the only difference between both populations is the firing frequency, as all other parameters are set to be the same. Consequently, we only investigate the influence of the firing rate of the first population as results also hold for the second population.
2.2. CalciumBased Synaptic Plasticity Model
The strength of an individual synapse connecting a presynaptic neuron i with the postsynaptic one, denoted by ρ_{i}, underlies the following dynamic (Graupner and Brunel, 2012;Graupner et al., 2016):
The first term describes the processes of longterm potentiation (LTP) at a synapse with γ_{p} being the rate of synaptic increase, θ_{p} being the potentiation threshold. c_{i}(t) is the activitydependent postsynaptic calciumconcentration, which is described either by a linear calcium model (LC; see section 2.2.1) or a more complex nonlinear calcium model (NLC; see section 2.2.2). H(x) denotes the Heaviside function. Analogously, the second term describes the processes of longterm depression (LTD) at a synapse in accordance to the parameters γ_{d} and θ_{d}.
The noise term accounts for coincidental activitydependent dynamics and nonplasticity related weight dynamics and is given by
Noisei(t)=στH(ci(t)θp)+H(ci(t)θd)ηi(t),
where τ is the same temporal constant as in Equation 1, σ is a constant amplitude parameter, and η_{i}(t) represents Gaussian white noise with mean zero and unit variance drawn for each synapse individually.
2.2.1. Linear Calcium Dynamics (LC)
In the linear calcium (LC) model, the dynamics of the intracellular calcium concentration of the postsynaptic spine for each synapse, denoted as ci˙, is given by:
c˙i=−ciτCa+CpreΣkδ(t−tki)+CpostΣlδ(t−tl),
where τ_{Ca} is the time constant of the calcium decay. The calcium concentration increases with pre and postsynaptic spike times tki and t_{l} of amplitudes C_{pre} and C_{post}, respectively. Delays between spiking events and the flow of ions are neglected as such delays only shift the independent, probabilistic spike events in the time domain. Please note, instead of Equation 3, we also consider a more complex set of equations to model the calcium dynamics as described in the following.
2.2.2. Nonlinear Calcium Dynamics (NLC)
In the nonlinear calcium (NLC) model, the calcium concentration c_{i}(t) = c_{i, pre}(t) + c_{i, post}(t) is determined by presynaptically evoked transients of calcium, described by
c˙i,pre=−ci,preτCa+CpreΣkδ(t−tki),
and postsynaptically evoked transients of calcium given by
The multitimescale adaptive threshold (MAT) model (Kobayashi et al., 2009; Yamauchi et al., 2011) describes the postsynaptic membrane potential as a leaky integrator:
τmV˙=V(t)+RI(t),
where τ_{m} is the membrane time constant, R the membrane resistance, and I the excitatory postsynaptic current. The excitatory postsynaptic current provides the linkage to the calciumbased plasticity model by summing the incoming spikes dependent on the synaptic weights:
I(t)=α(N)τmR∑i=1Nρiδ(ttki).
N denotes the size of the presynaptic population, ρ_{i} the synaptic strength connecting the presynaptic neuron i to the postsynaptic neuron, tki the presynaptic spiking times, and α(N) is a populationsizedependent scaling constant for the current.
After each postsynaptic spike time t_{j}, the spiking threshold S gets adapted by
S(t)=∑jΔ(ttj)+S0,
with S_{0} being the resting threshold. The update of the threshold and its decay is described by multiple timescales, τ_{1} and τ_{2}, and amplitudes, α_{1} and α_{2}, and is given by
The adaptive exponential integrateandfire (AEIF) model (Brette and Gerstner, 2005) describes the postsynaptic membrane potential as
τmV˙=V+EL+ΔTexp(VVTΔT)Rz+RI,
where τm=CgL=CR is the membrane time constant composed of the membrane capacitance C and the leak conductance g_{L}, E_{L} is the leak reversal potential, V_{T} the spiking threshold, Δ_{T} a slope factor, I the excitatory postsynaptic current as given in Equation 8, and z is the adaptation current determined by
τzz˙=a(VEL)z.
Here, τ_{z} denotes the adaptation time constant and a the subthreshold adaptation. Once a spike is triggered, the variable z is increased by an amount b which denotes the spiketriggered adaptation and V is reset to the reset potential V_{r} = E_{L}.
2.4. Acquisition of Data, Transformation, and Analysis
The acquisition of the data resulting from the detailed model and its transformation to an activitydependent, compact model of synaptic plasticity is carried out in a similar manner for all three different neuronal setups, two plasticity, and two neuron models. In the following, the corresponding steps taken are explained (Figure 2).
Generic approach to derive the firing rate dependencies of calciumbased STDP. Step I involves the simulation of the detailed calciumbased STDP model. Step II transforms the data obtained in Step I to ratedependent data of synaptic plasticity. In Step III the ratedependent data is analyzed using a regression algorithm and the results are evaluated.
Step I: Simulation of CalciumBased STDP.
We solve the differential equations of synaptic (and neuronal) dynamics numerically by the Euler method with time step Δt = 0.5 ms for different presynaptic input frequencies. The parameter space of the system is given by the presynaptic firing frequency u (and postsynaptic activity v for P_{1}) between 0 and 100 Hz (with a step size of 1 Hz) and the initialization of the synaptic strength w(t_{0}) between 0 and 1 (in arbitrary units) with a step size of 0.05. Please see Table 2 for an overview of the analyzed parameter spaces. An exemplary time evolution of the involved quantities is shown in Figure 1 for the LCmodel (section 2.2.1) with the MAT neuron model (section 2.3.1). Due to the fast time constant for the calcium concentration (τ_{Ca}), the calcium concentration is initialized at zero.
Mathematical notations for the different neuronal setups.
Postsynaptic firing
Change in synaptic strength
Parameter space
P_{1}
v is preset
ẇ^{num}(u, v, w)
Ω_{1} = {(u, v, w):u, v ∈ [0, 100], w ∈ [0, 1]}
P_{2}
v = v(u, w)
ẇ^{num}(u, v(u, w), w)
Ω_{2} = {(u, v, w):u ∈ [0, 100], v = v(u, w), w ∈ [0, 1]}
P_{3}
v=v(u_,w_)
w˙1num(u1,v(u_,w_),w1)
Ω3={(u_,v,w_):u1,2∈[0,100],v=v(u_,w_),w1,2∈[0,1]}
As the activities drive the change of synaptic weights, given each pair of pre and postsynaptic spike train, we calculate the average synaptic strength (w̄) of all synapses (ρ_{i}) across a period of two seconds:
w̄(u,v,w,t)=1N∑i=1Nρi(t),
where N is the size of the presynaptic population and ρ_{i} is given by solving Equation 1 with corresponding calcium dynamics numerically. Alongside, the numerical simulation yields the activity of the postsynaptic neuron. The latter is preset for P_{1} or calculated as spike count over 500 ms for P_{2} and P_{3}, denoted as v(u, w) for P_{2} and v(u_,w_) for the P_{3}setup (here x_ indicates a vector). However, the resulting w̄(u,v,w,t) is analyzed further in Step II described next.
Step II: Transformation to RateBased Data for Synaptic Plasticity.
The mean synaptic strength w̄(t) from the numerical solution of the detailed synaptic plasticity model is further approximated by cubic splines to w~(t). By this, the result is regularized from random fluctuations and enables us to calculate a reliable numerical derivative at the initial point in time. This approximation is described by
w~(u,v,w,t)=Splines=0.1k=1(w̄(u,v,w,t)),
where Spline is the UnivariateSpline function in the python library scipy.interpolate with parameters k and s. Here, k is the degree of the smoothing spline and s specifies the number of knots. The routine increases the number of knots until the smoothing condition ∫0T(w̄(t)w~(t))2dt≤s is satisfied.
By considering the synaptic dynamics w~(u,v,w,t) only around the first time point t = 0, on the one hand, the direct timedependency is eliminated and, on the other hand, the weightdependency of synaptic plasticity becomes controllable by neglecting its longterm development. Instead, we scan the weightdependency of the weight change by considering different initial values for w. Thus, the change in synaptic strength of the population can be expressed as:
ŵ˙(u,v,w)=w~˙(u,v,w,t=0).
Please note that the resulting ŵ˙(u,v,w) does not significantly depend on the initial distribution of ρ_{i} (see Figure S1). We average the change in synaptic strength of the population over several statistically independent initializations z (total number of intializations Z = 100 for P_{1}, P_{2} and Z = 38 for P_{3}) of the same parameter set. By this, we obtain the final overall change in synaptic strength (from numerics)
w˙num(u,v,w)=1Z∑z=1Zŵ˙z(u,v,w).
In case of P_{3}, in which we consider two independent presynaptic populations, thus, v(u_{1}, u_{2}, w_{1}, w_{2}), we add indices in order to differentiate between both populations such that
w˙1num(u1,v,w1) and w˙2num(u2,v,w2)
describe the change of the average synaptic strength of population 1 and 2, respectively. As mentioned before, in the P_{3}setup both populations dynamics are similar. Therefore, we will only consider the weight changes and variables related to the first population, if not stated otherwise.
In summary, based on these definitions for all neuronal setups (Table 2), the (average) synaptic weight changes of the detailed model ẇ commonly depends on the average population firing rate u, the postsynaptic firing rate v, and the average synaptic strength w, which are all local variables of synaptic plasticity (Gerstner and Kistler, 2002).
Step III: Regression Analysis.
Given the transformation of the numerical data from spikedependencies to firing rate dependencies (Step II), in the following, we fit our ratebased data of synaptic changes ẇ^{num}(u, v, w) by regression analysis to obtain compact mathematical formulations of the underlying weight dynamics. Thus, we will obtain a function of the form of
w˙ij=f(uj,vi,wij),
with i being the postsynaptic and j the presynaptic neuron.
To obtain such a function, we fit the ratedependent data (Step II) by a Taylor series expansion containing combinations of the involved state variables u, v, w (for better readability we omit the neuronal indices) with coefficients c_{αβγ}. The indices α, β, and γ indicate the order of the state variables u, v, and w in the corresponding feature u^{α}v^{β}w^{γ}:
Here, we limit the amount of considered features (individual terms) to the 27 listed in Table 4 to reduce the number of possible combinations for the regression analysis. These are all features with α, β, γ ∈ {0, 1, 2}; thus, individual state variables of these features are maximally of the second order. The choice of the second order is related to general ratebased equations used in literature, such as the Hebb rule (Hebb, 1949), the Oja rule (Oja, 1982), or the BCM rule (Bienenstock et al., 1982). In all these equations, no quantity is higher than the second order (Tetzlaff et al., 2011). Of course, higher orders or rather more features can be considered if required.
The regression algorithm, we use, is a weighted linear least squares regression and thus minimizes the weighted sum of squared residuals. The weights are inferred from the variance across the results from the Z independent initializations. As a measure of accuracy, the weighted coefficient of determination (R^{2}) was determined from 5fold crossvalidation. Note that it was further unbiased with respect to the number of considered features (Theil, 1958; Willett and Singer, 1988).
Qualitative Analysis of Resulting Compact Models
From the procedure described above we obtain a compact model consisting of a simplified differential equation, which resembles the dynamics of the corresponding complex model. However, to provide a better understanding of the underlying properties of the resulting compact models, we will assess whether these fulfill several qualitative characteristics of synaptic plasticity identified from numerical results (here we focus on the LCmodel). For this, we define nine different qualitative characteristics, which are provided in the following. Please note that these characteristics are defined such that a compact model can either fulfill them or not.
LTP Area: If the average synaptic weight is below a certain value, essentially for all combinations of pre and postsynaptic frequencies, LTP dominates the synaptic changes. The LTP Area weight value is about 0.3 (ẇ(w) > 0 for w > 0.3).
LTD Increase: If the synaptic weight is between the LTP Area weight value and the equilibrium value of the synaptic dynamics (w^{*}), the influence of LTD on the synaptic dynamics increases with larger weight values which, in turn, decreases the overall influence of LTP with LTP still being the dominant process (ẇ > 0 with ẇ(w_{1}) > ẇ(w_{2}) for w_{1} < w_{2}).
LTD Area: If the synaptic weight is above the equilibrium value w^{*}, LTD dominates the synaptic weight dynamics for all combinations of pre and postsynaptic frequencies (ẇ(w) < 0 for w > w^{*}).
Invariability: If the pre and postsynaptic activities are below 10 Hz, synaptic changes become negligible (ẇ ≈ 0 for u, v < 10).
Curvilinearity: For increasing pre and/or postsynaptic activities (above 10 Hz), ẇ(u, v) switches from a convex function to a concave function.
Saturation: For high levels of neuronal activity, synaptic changes (ẇ) become independent of the actual frequency of neuronal activities.
wBehavior: For the P_{1}setup, the change of synaptic weights depends linearly on the actual synaptic weight (ẇ(w) ≈ w). For the P_{2} and P_{3}setup, ẇ(w) has one maximum.
Competition: (Only for the P_{3}setup) If one population is highly active, the connections from the second to the postsynaptic neuron experience depression via competition.
Steadiness: (Only for the P_{3}setup) If the activity of one population approaches a maximum level, the influence of the other population on the system dynamics is negligible.
3. Results
In this study, we develop a simple, compact mathematical model by inferring the neuronal firing rate dependency of synaptic plasticity from a detailed, calciumbased synaptic plasticity model (Graupner and Brunel, 2012; Graupner et al., 2016). For this, we consider three different neuronal setups of one (P_{1}, P_{2}) or two (P_{3}) presynaptic neuronal population(s) connected onto one postsynaptic neuron and analyze the dependencies of synaptic changes (ẇ) on the (average) pre and postsynaptic neuronal activities (u, v) and the corresponding synaptic weight (w). For generality, we consider two different models of synaptic plasticity and two different neuron models (see section 2). In other words, for a given neuron model, plasticity model, and neuronal setup, first, we simulate the dynamics of the detailed plasticity model given different levels of pre and postsynaptic activities. Then, we transfer the data from spikingbased to ratebased and, finally, fit the activitydependencies by regression analysis. This analysis provides the synaptic changes given different terms or features of pre, postsynaptic activity, and the synaptic weight and is repeated for other combinations of models and setups. For details please see section 2. Please note that, in the following, we will focus on synaptic plasticity with linear calcium dynamics (LCmodel; section 2.2.1) with the MAT neuron model (section 2.3.1) to explain the basics of our approach and results.
We restricted the regression analysis by considering maximal second order terms resulting in 27 different features of u, v, w each individually weighted (see Table 4 for all considered features). Thus, for all neuronal setups, we searched for the most relevant combinations of 27 features that fit best the simulationbased dependencies evaluated by calculating the coefficient of determination between fit and simulation data from 5fold crossvalidation in the fourdimensional space of ẇ, w, u, v (see Figure 3A for different intersections). The resulting best fits serve as reference estimator of the corresponding setup. The reference estimators reached accuracies of R^{2} = 0.981 for the P_{1}setup (Figures 3 A.1,A.2,C.1), R^{2} = 0.986 for the P_{2}setup (Figures 3 A.3,A.4,C.2), and R^{2} = 0.963 for the P_{3}setup (Figures 3A.5,A.6,C.3). To compare the resulting estimators more qualitatively with the simulation data, we defined nine characteristics of the simulated synaptic weight changes and analyzed whether the resulting estimators or compact models match these (see Material and Methods and Table 3). The reference estimators meet all characteristics (except Saturation and, additionally, Steadiness in P_{3}) and, thus, the reference estimators basically describe the dynamics of the detailed, calciumbased synaptic dynamics. However, common simple models of synaptic plasticity consist of a small number of individual terms (Gerstner and Kistler, 2002; Tetzlaff et al., 2011; Gerstner et al., 2014). Thus, we repeated the regression analysis searching for the combination of three features fitting best the data (Figure 3B) and compared the most accurate threefeatureestimators with the corresponding reference estimator. As expected, the lower number of considered features fits the data less accurate (see Figure 3C), however, the resulting differences remain quantitative (large values of R^{2}) and qualitative small (compare Figure 3B with A).
Firing ratedependent synaptic plasticity dynamics for the three different neuronal setups P_{1} (left) P_{2} (middle), and P_{3} (right). Data for the LCmodel with the MAT neuron model are shown. (A) Activitybased change in synaptic strength for the reference fit with 27 features (F = 27) for predefined initial synaptic weights (odd panels) and neuronal activities (even panels) for each neuronal setup, respectively. Deviations from simulated data are negligible. (B) Same as in (A) for the best model with three features (F = 3). (C) For a wide variety of number of features the corresponding models (indicated by the coefficient of determination R^{2}) match the simulated data describing the dependencies of synaptic weight changes on activities and current weight.
Assessment of best fit models with nine qualitative characteristics.
Setups
P_{1}
P_{2}
P_{3}
P_{1}
P_{2}
P_{3}
R^{2}
0.810
0.822
0.652
0.981
0.986
0.963
RU2



0.967
Estimator
(Equation 19)
(Equation 20)
(Equation 21)
Reference (F= 27)
Quality
1 LTP Area
✓
✓
✓
✓
✓
✓
2 LTD Increase
✓
✓
✓
✓
✓
✓
3 LTD Area
✓
✓
✓
✓
✓
✓
4 Invariability
✗
✓
✓
✓
✓
✓
5 Saturation
✗
✗
✗
✗
✗
✗
6 Curvilinearity
✗
✗
✗
✓
✓
✓
7 wBehavior
✓
✓
✓
✓
✓
✓
8 Competition


✗


✓
9 Steadiness


✗


✗
(Left) Assessment of the best individual model for each neuronal setup and (Right) the reference fit.
In more detail, starting with the P_{1}setup consisting of one presynaptic population connected to a postsynaptic neuron of clamped activity, the most accurate threefeatureestimator (Figures 3B.1,B.2) is given by
P1setup: w˙=c010v+c011vw+c102uw2
with an accuracy of R^{2} = 0.810 (with weighting c_{αβγ} of feature u^{α}v^{β}w^{γ}; c010=0.007832±5·106, c011=0.009186±8·106, c102=0.000989±4·106). This model fulfills all linear and nonlinear characteristics except the nonlinear characteristics of Invariability, Saturation, and Curvilinearity (Table 3).
Next, the most accurate threefeatureestimator for the P_{2}setup (Figures 3B.3,B.4), which is a presynaptic population connected to a postsynaptic neuron with dynamic postsynaptic activity, is
P2setup: w˙=c200u2+c110uv+c202u2w2
with an accuracy of R^{2} = 0.822 (with c200=0.0001900±107, c110=0.0000527±3·107, c202=0.0001197±4·107) and fulfilling the majority of the characteristics correctly except Saturation and Curvilinearity (Table 3). Finally, the best threefeatureestimator for the P_{3}setup (Figures 3B.5,B.6), which is similar to the P_{2}setup except that it consists of two presynaptic populations, is given by
P3setup: w˙=c200u2+c210u2v+c202u2w2
with an accuracy of R^{2} = 0.652 (with c200=0.0000227±2·109, c210=0.0000004±1·1010, c202=0.0000789±7·109) and depicts five of nine characteristics correctly (all except Saturation, Curvilinearity, Steadiness, and Competition; Table 3).
Although the threefeatureestimators do not match the nonlinear aspects (e.g., Curvilinearity) of the detailed synaptic plasticity dynamics, the resulting rules (Equations 19–21) reliably match the basics of synaptic plasticity, for instance the regimes of LTP and LTD. Thus, they are representing simple, compact mathematical models of synaptic plasticity specifically for each setup.
Given the different models with numbers of features ranging from three to ten (Figure 3C), we analyzed the consistency of each feature within all the eight best estimators per setup to derive which features are generally powerful in describing synaptic plasticity. For the P_{1}setup, the accuracies of these estimators range from 0.810 for three features to 0.975 for ten features (Figure 3C.1). For the setups P_{2} and P_{3}, the accuracies lie in [0.822, 0.979] and [0.652, 0.956], respectively (Figures 3C.2,C.3).
Next, we repeated these analyses for all models and derived the absolute frequency for each feature occurring in the eight best estimators (Table 4). For the P_{1}setup the most frequently occurring features are uv, uw and vw. The Hebbian correlation terms (uv and uvw) are also often present in the P_{2} and P_{3}setup. Interestingly, besides the Hebbian terms, for all setups and models, the u^{2} and u^{2}wfeature occur frequently. This suggests that the square of the presynaptic firing rate is an important component influencing synaptic weight dynamics. To verify this finding, we excluded all features having a u^{2}term and derived new estimators. These estimators match the data from the detailed model less accurately (red in Figure 3C) than the estimators with u^{2}terms (blue). Note that, if for instance all terms being linear in u are excluded (green), the R^{2}values of the resulting estimators are similar to the complete estimators emphasizing the importance of the u^{2}terms. Thus, we conclude a general trend toward a significant correlation between synaptic plasticity and the square of presynaptic activity, which is not often considered in literature.
Number of occurrences of each feature for the eight best estimators for each setup and different models.
P_{1}
P_{2}
P_{3}
LC
NLC
LC + MAT
NLC + MAT
LC + AEIF
LC + MAT
NLC + MAT
LC + AEIF
1
1








2
u
2

4

5
3

6
3
v
8

4
4
1



4
w








5
u^{2}
7

8

8
8

8
6
uv
2
8
4
8
7
7
3
7
7
uw
4
6
2

2



8
v^{2}

6

3
2
5


9
vw
8
7
4
1
5
6
6
3
10
w^{2}








11
u^{2}v

2
2
2

1
5
1
12
u^{2}w
5
6
7
6
6
7
2
7
13
uv^{2}
1
1

3
3
3
5
2
14
uvw
1
6
5
6
1


3
15
uw^{2}
2







16
v^{2}w


3
3
3



17
vw^{2}


4


1
2
4
18
u^{2}v^{2}
3
3
2


3
4
4
19
u^{2}vw


1
3


3

20
u^{2}w^{2}


1

2
1
6
2
21
uv^{2}w
1
2

3
3
2
5

22
uvw^{2}
3
1


2



23
v^{2}w^{2}
3

1
2
1

3

24
u^{2}v^{2}w

1

3

1
4
1
25
u^{2}vw^{2}
1
1




1
2
26
uv^{2}w^{2}

1

1
1
1


27
u^{2}v^{2}w^{2}
1
1

4

3
3
2
For each setupmodel combination, the three most occurring features are highlighted in gray. LC: linear calcium model (section 2.2.1); NLC: nonlinear calcium model (section 2.2.2); MAT: multitimescale adaptitive threshold neuron model (section 2.3.1); AEIF: adaptive exponential integrateandfire neuron model (section 2.3.2).
After deriving the best compact model for each setup (Equations 19–21), next, we will derive the model which most precisely describes synaptic plasticity regarding all three different neuronal setups. For that purpose we are focusing on the LCmodel with the MAT neuron model. We consider a unified accuracy measure RU2 being the mean over all individual accuracies minus their standard deviation. Subtracting the standard deviation is meant to avoid a bias with respect to a certain setup.
In the following, we will discuss the three most accurate unified estimators resulting from this analysis (Table 5): The first most accurate unified estimator is
w˙=c200u2+c210u2v+c201u2w
with c200=0.0000228±4·109, c210=0.0000004±2·1010, c201=0.0000789±2·108. Please note that the c_{αβγ}values for all unified estimators are derived from a 5fold crossvalidation across all setups and models.
Assessment of the best three unified models describing the collective dynamics of all three neuronal setups with nine qualitative characteristics.
Setups
P_{1}
P_{2}
P_{3}
P_{1}
P_{2}
P_{3}
P_{1}
P_{2}
P_{3}
R^{2}
0.629
0.677
0.626
0.604
0.739
0.652
0.798
0.601
0.630
RU2
0.620
0.610
0.599
Estimator
(Equation22)
(Equation23)
(Equation24)
Quality
1 LTP Area
✓
✓
✓
✓
✓
✓
✓
✓
✓
2 LTD Increase
✗
✓
✗
✗
✗
✗
✓
✓
✓
3 LTD Area
✓
✗
✗
✓
✓
✗
✓
✓
✓
4 Invariability
✓
✓
✓
✓
✓
✓
✓
✓
✓
5 Saturation
✗
✗
✗
✗
✗
✗
✗
✗
✗
6 Curvilinearity
✗
✗
✗
✗
✗
✗
✗
✗
✗
7 wBehavior
✓
✓
✗
✗
✓
✓
✓
✓
✓
8 Competition


✗


✗


✗
9 Steadiness


✗


✗


✗
Only data for the LCmodel with the MAT neuron model are shown.
The second estimator is
w˙=c200u2+c210u2v+c202u2w2
with c200=0.0000360±4·109, c210=0.0000004±2·1010, c202=0.0000786±2·108. Both estimators or models contain the features u^{2}, u^{2}v and only differ slightly in their stabilizing features u^{2}w and u^{2}w^{2} (c_{201}, c_{202} < 0). Quantitatively, they best match the neuronal setup P_{2}. Qualitatively, they fail in correctly describing several characteristics of synaptic plasticity (Table 5).
On the other hand, the third most accurate unified estimator contains Hebbian terms (uv and uvw) as well as a stabilizing one uw^{2} with c_{102} < 0:
w˙=c110uv+c111uvw+c102uw2
with c110=0.0001272±3·108, c111=0.0001444±4·108, c202=0.0011699±2·107. Interestingly, this estimator matches best the P_{1}setup and it qualitatively generalizes better than the two most accurate unified estimators. It also correctly depicts diverse characteristics of synaptic plasticity. In addition, this estimator best captures the dynamics considering the nonlinear calcium dynamics (Table 6). We simulated the activitydependent synaptic plasticity dynamics of this estimator (compact model) and compared it with the detailed, calciumdependent model for the different setups using the LCmodel with MAT neuron model (Figure 4). Overall, the compact model matches the dynamics of the complex model especially in the longterm dynamics for LTP as well as LTD. Thus, the longterm dynamics of the ratebased model converge to the synaptic value of spiketimingdependent plasticity.
R_{2}values for the unified estimators across all models.
c200u2+c210u2v+c201u2w (Equation 22)
c200u2+c210u2v+c202u2w2 (Equation 23)
c110uv+c111uvw+c102uw2 (Equation 24)
P1 LC
0.629
0.604
0.798
P2 LC + MAT
0.677
0.739
0.601
P3 LC + MAT
0.626
0.652
0.630
P1 NLC
0.485
0.517
0.800
P2 NLC + MAT
0.261
0.511
0.665
P3 NLC + MAT
0.317
0.427
0.644
P2 LC + AEIF
0.743
0.799
0.434
P3 LC + AEIF
0.700
0.715
0.456
Comparison of spiketiming dependent synaptic plasticity with ratebased synaptic plasticity according to the third unified estimator (Equation 24). Here, we show the results using the LCmodel with the MAT neuron model. (A,B): Dynamics of P_{1} for global initial synaptic strengths of w = 0.3 and pre and postsynaptic spikefiring of u = v = 35Hz (A; LTP) and w = 0.7 and u = v = 20Hz (B; LTD). The ratebased plasticity (RBP) is illustrated in red whereas the average spiketiming dependent plasticity (STDP) is illustrated in blue. (C,D) Same as in (A,B) for P_{2} with global initial synaptic strengths of w = 0.3 and presynaptic spikefiring of u = 65Hz (C; LTP) and w = 0.9 u = 36Hz (D; LTD). The ratebased plasticity (RBP) is illustrated in red whereas the average spiketiming dependent plasticity (STDP) is illustrated in blue. (E,F) Same as in C,D for competitive behavior in P_{3} of two presynaptic populations 1 (blue) and 2 (green) with initial synaptic weights of w_{1} = 0.9, w_{2} = 0.6 and presynaptic spikefiring of u_{1} = 5Hz, u_{2} = 75Hz (E; no competition) and u_{1} = 5Hz, u_{2} = 90Hz (F; strong competition). The ratebased plasticities for populations 1 and 2 are illustrated in yellow and red respectively. The corresponding average spiketiming dependent plasticity is illustrated in blue and green. In all cases the populations consist of N = 1, 000 synapses.
However, the effect of competition between two neuronal populations is not well matched by the simple ratebased model (Figures 4E,F). If we allow a fourth feature in the estimator, the resulting compact model describes competition, too (Figure 5). This is also an example indicating that, if a specific property of the detailed model is desired, one has to consider more features to obtain a compact model which also includes this property. Similarly, if a compact model is desired which matches more accurate the detailed model, more features have to be incorporated (see e.g., Figure 3C).
Comparison of spiketiming dependent synaptic plasticity with ratebased synaptic plasticity according to the third unified estimator (Equation 24) with the additional feature vw^{2}. As in Figure 4, we show the results using the LCmodel with the MAT neuron model. (A,B): Dynamics of P_{1} for global initial synaptic strengths of w = 0.3 and pre and postsynaptic spikefiring of u = v = 35Hz (A; LTP) and w = 0.7 and u = v = 20Hz (B; LTD). The ratebased plasticity (RBP) is illustrated in red whereas the average spiketiming dependent plasticity (STDP) is illustrated in blue. (C,D) Same as in (A,B) for P_{2} with global initial synaptic strengths of w = 0.3 and presynaptic spikefiring of u = 65Hz (C; LTP) and w = 0.9 u = 36Hz (D; LTD). The ratebased plasticity (RBP) is illustrated in red whereas the average spiketiming dependent plasticity (STDP) is illustrated in blue. (E,F) Same as in (C,D) for competitive behavior in P_{3} of two presynaptic populations 1 (blue) and 2 (green) with initial synaptic weights of w_{1} = 0.9, w_{2} = 0.6 and presynaptic spikefiring of u_{1} = 5Hz, u_{2} = 75Hz (E; no competition) and u_{1} = 5Hz, u_{2} = 90Hz (F; strong competition). The ratebased plasticities for populations 1 and 2 are illustrated in yellow and red respectively. The corresponding average spiketiming dependent plasticity is illustrated in blue and green. In all cases the populations consist of N = 1, 000 synapses.
In summary, our results lead to the conclusion that features of squared presynaptic activity are well describing and likely be associated to synaptic plasticity. Nevertheless, the third unified estimator (Equation 24) provides evidence that models of other features exist that better describe the qualitative linear characteristics of synaptic plasticity throughout different neuronal setups while they quantitatively perform worth.
4. Discussion
We developed a method to derive a simple, compact mathematical model of synaptic plasticity based on a biologically detailed, calciumbased model. In contrast to the detailed model, the resulting simple model can be used for further analytical and/or numeric analysis of, for instance, largescale systems.
An alternative method is, for instance, to derive the crosscorrelation function of neuronal firing (Brunel and Hakim, 1999; Ostojic et al., 2009) and, then, to incorporate this function into a meanfield model of synaptic plasticity (Kempter et al., 1999). This method is mathematically more accurate as our method and it can also be used to investigate the selforganization of neural circuits (Ocker et al., 2015; Tannenbaum and Burak, 2015). However, deriving the crosscorrelation function of neuronal firing or the meanfield model of synaptic plasticity for diverse detailed neuron and plasticity models is very complex and, up to now, not always possible. On the other hand, the complexity of our method depends less on the considered neuron and synaptic plasticity models.
We considered two different versions of the calciumbased STDP model with parameters which where determined by Graupner et al. (2016) on slices of the visual cortex from macaque monkeys. Of course, these parameters (and also from the neuron model) could vary for different brain regions. However, with the heredeveloped method one can easily link these parameters to simpler models of plasticity. Given that the function of a brain region is mainly determined on the network level, by using the resulting simple model one can investigate the influence of the parameters on the functional properties of the brain regions.
Overall, we cannot guarantee that we applied the optimal weighting of the observations and thus found the best linear unbiased estimator. Technically, it could be checked if the whole covariance matrix of the N observations in Z initializations would provide a better weighting. However, this is computationally not feasible, especially for the case of P_{3} with 4.5 Mio. observations. Furthermore, we considered the estimators or models with the highest R^{2} while taking the next best model into account could yield other reasonable learning rules matching, for instance, more qualitative characteristics (see for instance Equation 24). If a compact model is required, which matches the detailed model with a higher accuracy, instead of three features, one can easily allow more features being present in the final estimator (Figure 3C). Furthermore, by considering features of higher order than two in Step III, the resulting final estimators could describe additional aspects of the detailed model. However, given the increased number of features, the regression procedure would require more computational resources.
The results for individual learning rules of the considered neuronal populations did not significantly match with, for instance, the BCM (Bienenstock et al., 1982) or the Oja's rule (Oja, 1982). Although the results significantly match Hebbian correlation learning, the stabilization mechanisms are different to commonly used ones. Overall, regarding the treated ratebased models, our applied methods on three different neuronal setups suggest that features from the BCM rule and Oja's rule are not the best estimators for the underlying plasticity data. From our methods, features observed to be correlated to LTP are the Hebbain term uv and terms, which depend on the square of the presynaptic activity such as u^{2}v (opposing to v^{2}w from BCM theory) and u^{2}. Features observed to be correlated to LTD are vw, u^{2}w (opposing to v^{2}w from Oja's rule), and uvw. The u^{2}terms could imply that, on the network level, synaptic dynamics are more sensitive to variations in the inputs (u) than in the network activity (v) providing a stronger influence of the actual inputs on the overall network dynamics. However, such dynamics and the u^{2}dependency require further theoretical and experimental studies. For instance, by varying systematically the rate of stimulation triggering synaptic plasticity (similar to Sjöström et al., 2001), the resulting ẇurelation could indicate a u^{2}dependency. Overall, the plausibility of certain terms remain unclear and corresponding molecular pathways of synaptic plasticity still have to be investigated to find pieces of evidence for these terms. However, on the one hand, our method provides a way to link complex synaptic plasticity dynamics with dynamics on the network level and, on the other hand, the individual terms from our method allow to identify important ingredients of synaptic plasticity enabling an (functional) ordering of diverse molecular pathways.
Author Contributions
CT study design. JL data acquisition. JL, JH, and CT data analysis and manuscript written.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We acknowledge support by the German Research Foundation and the Open Access Publication Funds of the Göttingen University.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2019.00026/full#supplementarymaterial
ReferencesAntunesG.RoqueA. C.Simoesde SouzaF. M. (2016). Stochastic induction of longterm potentiation and longterm depression. BienenstockE. L.CooperL. N.MunroP. W. (1982). Theory for the development of neuron selectivity: orientation specificity and binocular interaction in visual cortex. BretteR.GerstnerW. (2005). Adaptive exponential integrateandfire model as an effective description of neuronal activity. BrunelN.HakimV. (1999). Fast global oscillations in networks of integrateandfire neurons with low firing rates. ChoquetD.TrillerA. (2013). The dynamic synapse. DayanP.AbbottL. F. (2001). EarnshawB. A.BressloffP. C. (2006). Biophysical model of AMPA receptor trafficking and its regulation during longterm potentiation/longterm depression. GallimoreA. R.KimT.TanakaYamamotoK.De SchutterE. (2018). Switching on depression and potentiation in the cerebellum. GerstnerW.KistlerW. M. (2002). Mathematical formulations of Hebbian learning. GerstnerW.KistlerW. M.NaudR.PaninskiL. (2014). GerstnerW.SprekelerH.DecoG. (2012). Theory and simulation in neuroscience. GraupnerM.BrunelN. (2010). Mechanisms of induction and maintenance of spiketiming dependent plasticity in biophysical synapse models. GraupnerM.BrunelN. (2012). Calciumbased plasticity model explains sensitivity of synaptic changes to spike pattern, rate, and dendritic location. GraupnerM.WallischP.OstojicS. (2016). Natural firing patterns imply low sensitivity of synaptic plasticity to spike timing compared with firing rate. HashimotodaniY.OhnoShosakuT.TsubokawaH.OgataH.EmotoK.MaejimaT.. (2005). Phospholipase Cβ serves as a coincidence detector through its Ca^{2+} dependency for triggering retrograde endocannabinoid signal. HebbD. O. (1949). HuganirR. L.NicollR. A. (2013). AMPARs and synaptic plasticity: the last 25 years. KandelE. R.SchwartzJ. H.JessellT. M.SiegelbaumS. A.HudspethA. J.. (2000). KempterR.GerstnerW.Van HemmenJ. L. (1999). Hebbian learning and spiking neurons. KnoblauchA.SommerF. T. (2016). Structural plasticity, effectual connectivity, and memory in cortex. KobayashiR.TsuboY.ShinomotoS. (2009). Madetoorder spiking neuron model equipped with a multitimescale adaptive threshold. KorteM.SchmitzD. (2016). Cellular and system biology of memory: timing, molecules, and beyond. LiY.KulviciusT.TetzlaffC. (2016). Induction and consolidation of calciumbased homo and heterosynaptic potentiation and depression. LüscherC.MalenkaR. C. (2012). NMDA receptordependent longterm potentiation and longterm depression (LTP/LTD). MaejimaT. (2005). Synaptically driven endocannabinoid release requires Ca^{2+}assisted metabotropic glutamate receptor subtype 1 to phospholipase Cβ4 signaling cascade in the cerebellum. MartinS. J.GrimwoodP. D.MorrisR. G. M. (2000). Synaptic plasticity and memory: an evaluation of the hypothesis. NevianT.SakmannB. (2006). Spine Ca^{2+} signaling in spiketimingdependent plasticity. OckerG. K.LitwinKumarA.DoironB. (2015). Selforganization of microcircuits in networks of spiking neurons with plastic synapses. OjaE. (1982). Simplified neuron model as a principal component analyzer. OstojicS.BrunelN.HakimV. (2009). How connectivity, background activity, and synaptic properties shape the crosscorrelation between spike trains. SjöströmP. J.TurrigianoG. G.NelsonS. B. (2001). Rate, timing, and cooperativity jointly determine cortical synaptic plasticity. SjöströmP. J.TurrigianoG. G.NelsonS. B. (2003). Neocortical LTD via coincident activation of presynaptic NMDA and cannabinoid receptors. TannenbaumN. R.BurakY. (2015). Shaping neural circuits by high order synaptic interactions. TetzlaffC.KolodziejskiC.TimmeM.TsodyksM.WörgötterF. (2013). Synaptic scaling enables dynamically distinct short and longterm memory formation. TetzlaffC.KolodziejskiC.TimmeM.WörgötterF. (2011). Synaptic scaling in combination with many generic plasticity mechanisms stabilizes circuit connectivity. TheilH. (1958). TsodyksM.FeigelmanM. (1988). Enhanced storage capacity in neural networks with low level of activity. WillettJ. B.SingerJ. D. (1988). Another cautionary note about R 2: its use in weighted leastsquares regression analysis. YamauchiS.KimH.ShinomotoS. (2011). Elemental spiking neuron model for reproducing diverse firing patterns and predicting precise firing times.
Funding. The research was funded by the H2020FETPROACT project Plan4Act (#732266) (JH and CT), by the Collaborative Research Center SFB1286 (Project C1) (CT), and by the International Max Planck Research School for Physics of Biological and Complex Systems by stipends of the country of Lower Saxony with funds from the initiative Niedersächsisches Vorab and of the University of Göttingen (JH).