A Theoretical Framework to Derive Simple, Firing-Rate-Dependent Mathematical Models of Synaptic Plasticity

1. 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 spike-timings, 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 spike-dependency 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, firing-rate-dependent theoretical model from a biologically detailed spike-based model. For the detailed model, we consider two versions of a well-established calcium-based spike-timing-dependent 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 well-described by a compact, rate-dependent 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 spike-dependent pre- and postsynaptic processes. This calcium-dependent 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 calcium-dependent 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 rate-dependent 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, firing-rate-dependent mathematical formulation of the detailed dynamics of spike-timing 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 well-established calcium-based 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.

2.1. Neuronal Setups

Throughout this study we consider three basic neuronal motifs or setups: a presynaptic population with independent postsynaptic neuron (P₁; Figure 1, top left), a presynaptic population with dynamic postsynaptic neuron (P₂; Figure 1, top middle), and two presynaptic populations with dynamic postsynaptic neuron (P₃; Figure 1, top right).

Neuronal Population I (P₁): 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₂): This setup is similar to the P₁-setup. However, here, the activity of the postsynaptic neuron is not determined by a Poisson process (see P₁), instead, we consider a biologically detailed neuron model. For generality, throughout this study we use two different neuron models; the so-called Multi-timescale Adaptive Threshold model (MAT; Kobayashi et al., 2009) and the Adaptive Exponential Integrate-and-Fire 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 ($\dot{v}$) depends on u and w, i.e., $\dot{v}(u,w)$, which in turn influences the dynamics of synaptic plasticity (ẇ).

Neuronal Population II (P₃): This motif is an extension of the P₂-setup. Here, two presynaptic populations 1 and 2 are connected to the same postsynaptic neuron with each population consisting of independently firing neurons with population-specific frequency u₁ and u₂, respectively. Thus, the firing of the postsynaptic neuron depends on the presynaptic firing rates and strength of the connecting synapses of both presynaptic populations ($\dot{v}(u_1,u_2,w_1,w_2)$). 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. Calcium-Based 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 long-term potentiation (LTP) at a synapse with γ_p being the rate of synaptic increase, θ_p being the potentiation threshold. c_i(t) is the activity-dependent postsynaptic calcium-concentration, 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 long-term depression (LTD) at a synapse in accordance to the parameters γ_d and θ_d.

The noise term accounts for coincidental activity-dependent dynamics and non-plasticity related weight dynamics and is given by

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 $\dot{c_i}$, is given by:

where τ_Ca is the time constant of the calcium decay. The calcium concentration increases with pre- and postsynaptic spike times $t_k^i$ 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

and postsynaptically evoked transients of calcium given by

where the quantities τ_Ca, $t_k^i$, t_l, C_pre and C_post are analog to the linear calcium model. The nonlinearity factor

couples the pre- and postsynaptic calcium concentrations (Graupner et al., 2016).

2.3. Neuron Model

2.3.1. Multi-Timescale Adaptive Threshold Neuron (MAT)

The multi-timescale adaptive threshold (MAT) model (Kobayashi et al., 2009; Yamauchi et al., 2011) describes the postsynaptic membrane potential as a leaky integrator:

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 calcium-based plasticity model by summing the incoming spikes dependent on the synaptic weights:

N denotes the size of the presynaptic population, ρ_i the synaptic strength connecting the presynaptic neuron i to the postsynaptic neuron, $t_k^i$ the presynaptic spiking times, and α(N) is a population-size-dependent scaling constant for the current.

After each postsynaptic spike time t_j, the spiking threshold S gets adapted by

with S₀ being the resting threshold. The update of the threshold and its decay is described by multiple timescales, τ₁ and τ₂, and amplitudes, α₁ and α₂, and is given by

2.3.2. Adaptive Exponential Integrate-and-Fire Neuron (AEIF)

The adaptive exponential integrate-and-fire (AEIF) model (Brette and Gerstner, 2005) describes the postsynaptic membrane potential as

where ${\tau }_{\mathrm{\text{m}}}=\frac{C}{g_L}=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

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 spike-triggered 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 activity-dependent, 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).

Step I: Simulation of Calcium-Based 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₁) between 0 and 100 Hz (with a step size of 1 Hz) and the initialization of the synaptic strength w(t₀) 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 LC-model (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.

As the activities drive the change of synaptic weights, given each pair of pre- and postsynaptic spike train, we calculate the average synaptic strength ($\bar{w}$) of all synapses (ρ_i) across a period of two seconds:

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₁ or calculated as spike count over 500 ms for P₂ and P₃, denoted as v(u, w) for P₂ and $v(\underset{\_}{u},\underset{\_}{w})$ for the P₃-setup (here $\underset{\_}{x}$ indicates a vector). However, the resulting $\bar{w}(u,v,w,t)$ is analyzed further in Step II described next.

Step II: Transformation to Rate-Based Data for Synaptic Plasticity.

The mean synaptic strength $\bar{w}(t)$ from the numerical solution of the detailed synaptic plasticity model is further approximated by cubic splines to $\tilde{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

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 ${\int }_0^T{(\bar{w}(t)-\tilde{w}(t))}^2\mathrm{\text{d}}t\le s$ is satisfied.

By considering the synaptic dynamics $\tilde{w}(u,v,w,t)$ only around the first time point t = 0, on the one hand, the direct time-dependency is eliminated and, on the other hand, the weight-dependency of synaptic plasticity becomes controllable by neglecting its long-term development. Instead, we scan the weight-dependency of the weight change by considering different initial values for w. Thus, the change in synaptic strength of the population can be expressed as:

Please note that the resulting $\dot{ŵ}(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₁, P₂ and Z = 38 for P₃) of the same parameter set. By this, we obtain the final overall change in synaptic strength (from numerics)

In case of P₃, in which we consider two independent presynaptic populations, thus, v(u₁, u₂, w₁, w₂), we add indices in order to differentiate between both populations such that

describe the change of the average synaptic strength of population 1 and 2, respectively. As mentioned before, in the P₃-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 spike-dependencies to firing rate dependencies (Step II), in the following, we fit our rate-based 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

with i being the postsynaptic and j the presynaptic neuron.

To obtain such a function, we fit the rate-dependent 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 rate-based 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²) was determined from 5-fold cross-validation. 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 LC-model). 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.

1. 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).

2. 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₁) > ẇ(w₂) for w₁ < w₂).

3. 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*).

4. Invariability: If the pre- and postsynaptic activities are below 10 Hz, synaptic changes become negligible (ẇ ≈ 0 for u, v < 10).

5. Curvilinearity: For increasing pre- and/or postsynaptic activities (above 10 Hz), ẇ(u, v) switches from a convex function to a concave function.

6. Saturation: For high levels of neuronal activity, synaptic changes (ẇ) become independent of the actual frequency of neuronal activities.

7. w-Behavior: For the P₁-setup, the change of synaptic weights depends linearly on the actual synaptic weight (ẇ(w) ≈ w). For the P₂- and P₃-setup, ẇ(w) has one maximum.

8. Competition: (Only for the P₃-setup) If one population is highly active, the connections from the second to the postsynaptic neuron experience depression via competition.

9. Steadiness: (Only for the P₃-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, calcium-based synaptic plasticity model (Graupner and Brunel, 2012; Graupner et al., 2016). For this, we consider three different neuronal setups of one (P₁, P₂) or two (P₃) 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 spiking-based to rate-based and, finally, fit the activity-dependencies 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 (LC-model; 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 simulation-based dependencies evaluated by calculating the coefficient of determination between fit and simulation data from 5-fold cross-validation in the four-dimensional 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² = 0.981 for the P₁-setup (Figures 3 A.1,A.2,C.1), R² = 0.986 for the P₂-setup (Figures 3 A.3,A.4,C.2), and R² = 0.963 for the P₃-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₃) and, thus, the reference estimators basically describe the dynamics of the detailed, calcium-based 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 three-feature-estimators 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²) and qualitative small (compare Figure 3B with A).

In more detail, starting with the P₁-setup consisting of one presynaptic population connected to a postsynaptic neuron of clamped activity, the most accurate three-feature-estimator (Figures 3B.1,B.2) is given by

with an accuracy of R² = 0.810 (with weighting c_αβγ of feature u^αv^βw^γ; $c_{010}=0.007832\pm 5·10^{-6}$, $c_{011}=-0.009186\pm 8·10^{-6}$, $c_{102}=-0.000989\pm 4·10^{-6}$). This model fulfills all linear and nonlinear characteristics except the nonlinear characteristics of Invariability, Saturation, and Curvilinearity (Table 3).

Next, the most accurate three-feature-estimator for the P₂-setup (Figures 3B.3,B.4), which is a presynaptic population connected to a postsynaptic neuron with dynamic postsynaptic activity, is

with an accuracy of R² = 0.822 (with $c_{200}=0.0001900\pm 10^{-7}$, $c_{110}=0.0000527\pm 3·10^{-7}$, $c_{202}=-0.0001197\pm 4·10^{-7}$) and fulfilling the majority of the characteristics correctly except Saturation and Curvilinearity (Table 3). Finally, the best three-feature-estimator for the P₃-setup (Figures 3B.5,B.6), which is similar to the P₂-setup except that it consists of two presynaptic populations, is given by

with an accuracy of R² = 0.652 (with $c_{200}=0.0000227\pm 2·10^{-9}$, $c_{210}=0.0000004\pm 1·10^{-10}$, $c_{202}=-0.0000789\pm 7·10^{-9}$) and depicts five of nine characteristics correctly (all except Saturation, Curvilinearity, Steadiness, and Competition; Table 3).

Although the three-feature-estimators 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₁-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₂ and P₃, 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₁-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₂- and P₃-setup. Interestingly, besides the Hebbian terms, for all setups and models, the u²- and u²w-feature 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²-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²-terms (blue). Note that, if for instance all terms being linear in u are excluded (green), the R²-values of the resulting estimators are similar to the complete estimators emphasizing the importance of the u²-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.

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 LC-model with the MAT neuron model. We consider a unified accuracy measure $R_U^2$ 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

with $c_{200}=0.0000228\pm 4·10^{-9}$, $c_{210}=0.0000004\pm 2·10^{-10}$, $c_{201}=-0.0000789\pm 2·10^{-8}$. Please note that the c_αβγ-values for all unified estimators are derived from a 5-fold cross-validation across all setups and models.

The second estimator is

with $c_{200}=0.0000360\pm 4·10^{-9}$, $c_{210}=0.0000004\pm 2·10^{-10}$, $c_{202}=-0.0000786\pm 2·10^{-8}$. Both estimators or models contain the features u², u²v and only differ slightly in their stabilizing features u²w and u²w² (c₂₀₁, c₂₀₂ < 0). Quantitatively, they best match the neuronal setup P₂. 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² with c₁₀₂ < 0:

with $c_{110}=0.0001272\pm 3·10^{-8}$, $c_{111}=-0.0001444\pm 4·10^{-8}$, $c_{202}=-0.0011699\pm 2·10^{-7}$. Interestingly, this estimator matches best the P₁-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 activity-dependent synaptic plasticity dynamics of this estimator (compact model) and compared it with the detailed, calcium-dependent model for the different setups using the LC-model with MAT neuron model (Figure 4). Overall, the compact model matches the dynamics of the complex model especially in the long-term dynamics for LTP as well as LTD. Thus, the long-term dynamics of the rate-based model converge to the synaptic value of spike-timing-dependent plasticity.

However, the effect of competition between two neuronal populations is not well matched by the simple rate-based 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).

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, calcium-based model. In contrast to the detailed model, the resulting simple model can be used for further analytical and/or numeric analysis of, for instance, large-scale systems.

An alternative method is, for instance, to derive the cross-correlation function of neuronal firing (Brunel and Hakim, 1999; Ostojic et al., 2009) and, then, to incorporate this function into a mean-field 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 self-organization of neural circuits (Ocker et al., 2015; Tannenbaum and Burak, 2015). However, deriving the cross-correlation function of neuronal firing or the mean-field 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 calcium-based 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 here-developed 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₃ with 4.5 Mio. observations. Furthermore, we considered the estimators or models with the highest R² 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 rate-based 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²v (opposing to v²w from BCM theory) and u². Features observed to be correlated to LTD are vw, u²w (opposing to v²w from Oja's rule), and uvw. The u²-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²-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 ẇ-u-relation could indicate a u²-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.

Funding

The research was funded by the H2020-FETPROACT project Plan4Act (#732266) (JH and CT), by the Collaborative Research Center SFB-1286 (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).

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.

Acknowledgments

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#supplementary-material

References