Naturalgradient learning for spiking neurons
Abstract
In many normative theories of synaptic plasticity, weight updates implicitly depend on the chosen parametrization of the weights. This problem relates, for example, to neuronal morphology: synapses which are functionally equivalent in terms of their impact on somatic firing can differ substantially in spine size due to their different positions along the dendritic tree. Classical theories based on Euclideangradient descent can easily lead to inconsistencies due to such parametrization dependence. The issues are solved in the framework of Riemannian geometry, in which we propose that plasticity instead follows naturalgradient descent. Under this hypothesis, we derive a synaptic learning rule for spiking neurons that couples functional efficiency with the explanation of several welldocumented biological phenomena such as dendritic democracy, multiplicative scaling, and heterosynaptic plasticity. We therefore suggest that in its search for functional synaptic plasticity, evolution might have come up with its own version of naturalgradient descent.
Editor's evaluation
The natural gradient has a long and rich history in machine learning. Here, the authors derive a biologically plausible implementation of naturalgradientbased plasticity for spiking neurons, which renders learning invariant under dendritic transformations. This new synaptic learning rule makes several very interesting experimental predictions with respect to the interplay of homo and heterosynaptic plasticity, and with regard to the scaling of plasticity by the presynaptic variance.
https://doi.org/10.7554/eLife.66526.sa0Introduction
Understanding the fundamental computational principles underlying synaptic plasticity represents a longstanding goal in neuroscience. To this end, a multitude of topdown computational paradigms have been developed, which derive plasticity rules as gradient descent on a particular objective function of the studied neural network (Rosenblatt, 1958; Rumelhart et al., 1986; Pfister et al., 2006; D’Souza et al., 2010; Friedrich et al., 2011).
However, the exact physical quantity to which these synaptic weights correspond often remains unspecified. What is frequently simply referred to as ${w}_{ij}$ (the synaptic weight from neuron $j$ to neuron i) might relate to different components of synaptic interaction, such as calcium concentration in the presynaptic axon terminal, neurotransmitter concentration in the synaptic cleft, receptor activation in the postsynaptic dendrite or the postsynaptic potential (PSP) amplitude in the spine, the dendritic shaft or at the soma of the postsynaptic cell. All these biological processes can be linked by transformation rules, but depending on which of them represents the variable with respect to which performance is optimized, the network behavior during training can be markedly different.
As an example we consider the parametrization of the synaptic strength either as PSP amplitude in the soma, ${w}^{\mathrm{s}}$, or as PSP amplitude in the dendrite, ${w}^{\mathrm{d}}$ (see also Figure 1 and Sec. ‘The naive Euclidean gradient is not parametrizationinvariant’). Reparametrizing the synaptic strength in this way implies an attenuation factor for each single synapse, but different factors are assigned across the positions on the dendritic tree. As a consequence, the weight vector will follow a different trajectory during learning depending on whether the somatic or dendritic parametrization of the PSP amplitude was chosen.
It certainly could be the case that evolution has favored a parametrizationdependent learning rule, along with one particular parametrization over all others, but this would necessarily imply suboptimal convergence for all but a narrow set of neuron morphologies and connectome configurations. An invariant learning rule on the other hand would not only be mathematically unambiguous and therefore more elegant, but could also improve learning, thus increasing fitness.
In some aspects, the question of invariant behavior is related to the principle of relativity in physics, which requires the laws of physics – in our case: the improvement of performance during learning – to be the same in all frames of reference. What if neurons would seek to conserve the way they adapt their behavior regardless of, for example, the specific positioning of synapses along their dendritic tree? Which equations of motion – in our case: synaptic learning rules – are able to fulfill this requirement?
The solution lies in following the path of steepest descent not in relation to a small change in the synaptic weights (Euclideangradient descent), but rather with respect to a small change in the inputoutput distribution (naturalgradient descent). This requires taking the gradient of the error function with respect to a metric defined directly on the space of possible inputoutput distributions, with coordinates defined by the synaptic weights. First proposed in Amari, 1998, but with earlier roots in information geometry (Amari, 1987; Amari and Nagaoka, 2000), naturalgradient methods (Yang and Amari, 1998; Rattray and Saad, 1999; Park et al., 2000; Kakade, 2001) have recently been rediscovered in the context of deep learning (Pascanu and Bengio, 2013; Martens, 2014; Ollivier, 2015; Amari et al., 2019; Bernacchia et al., 2018). Moreover, Pascanu and Bengio, 2013 showed that the naturalgradient learning rule is closely related to other machine learning algorithms. However, most of the applications focus on ratebased networks which are not inherently linked to a statistical manifold and have to be equipped with Gaussian noise or a probabilistic output layer interpretation in order to allow an application of the natural gradient. Furthermore, a biologically plausible synaptic plasticity rule needs to make all the required information accessible at the synapse itself, which is usually unnecessary and therefore largely ignored in machine learning.
The stochastic nature of neuronal outputs invivo (see, e.g. Softky and Koch, 1993) provides a natural setting for plasticity rules based on information geometry. As a model for biological synapses, natural gradient combines the elegance of invariance with the success of gradientdescentbased learning rules. In this manuscript, we derive a closedform synaptic learning rule based on naturalgradient descent for spiking neurons and explore its implications. Our learning rule equips the synapses with more functionality compared to classical error learning by enabling them to adjust their learning rate to their respective impact on the neuron’s output. It naturally takes into account relevant variables such as the statistics of the afferent input or their respective positions on the dendritic tree. This allows a set of predictions which are corroborated by both experimentally observed phenomena such as dendritic democracy and multiplicative weight dynamics and theoretically desirable properties such as Bayesian reasoning (MarceauCaron and Ollivier, 2007). Furthermore, and unlike classical errorlearning rules, plasticity based on the natural gradient is able to incorporate both homo and heterosynaptic phenomena into a unified framework. While theoretically derived heterosynaptic components of learning rules are notoriously difficult for synapses to implement due to their nonlocality, we show that in our learning rule they can be approximated by quantities accessible at the locus of plasticity. In line with results from machine learning, the combination of these features also enables faster convergence during supervised learning.
Results
The naive Euclidean gradient is not parametrizationinvariant
We consider a cost function $C$ on the neuronal level that, in the sense of cortical credit assignment (see e.g. Sacramento et al., 2017), can relate to some behavioral cost of the agent that it serves. The output of the neuron depends on the amplitudes of the somatic PSPs elicited by the presynaptic spikes. We denote these ‘somatic weights’ by $\mathit{w}}^{\mathrm{s}$, and may parametrize the neuronal cost as $C={C}^{\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}}[{\mathit{w}}^{\mathrm{s}}]$.
However, dendritic PSP amplitudes ${\mathit{w}}^{\mathrm{d}}$ can be argued to offer a more unmitigated representation of synaptic weights, so we might rather wish to express the cost as $C={C}^{\text{}\mathrm{d}}[{\mathit{w}}^{\mathrm{d}}]$. These two parametrizations are related by an attenuation factor $\mathit{\alpha}$ (between 0 and 1):
In general, this attenuation factor depends on the synaptic position and is therefore described by a vector that is multiplied componentwise with the weights. For clarity, we assume a simplified picture in which we neglect, for example, the spread of PSPs as they travel along dendritic cables. However, our observations regarding the effects of reparametrization hold in general.
It may now seem straightforward to switch between the somatic and dendritic representation of the cost by simply substituting variables, for example
To derive a plasticity rule for the somatic and dendritic weights, we might consider gradient descent on the cost:
At first glance, this relation seems reasonable: dendritic weight changes affect the cost more weakly then somatic weight changes, so their respective gradient is more shallow by the factor $\mathit{\alpha}$. However, from a functional perspective, the opposite should be true: dendritic weights should experience a larger change than somatic weights in order to elicit the same effect on the cost. This conflict can be made explicit by considering that somatic weight changes are, themselves, attenuated dendritic weight changes: $\mathrm{\Delta}{\mathit{w}}^{\mathrm{s}}=\mathit{\alpha}\mathrm{\Delta}{\mathit{w}}^{\mathrm{d}}$. Substituting this into Equation 3 leads to an inconsistency: $\mathrm{\Delta}{\mathit{w}}^{\mathrm{d}}={\mathit{\alpha}}^{2}\mathrm{\Delta}{\mathit{w}}^{\mathrm{d}}$. To solve this conundrum, we need to shift the focus from changing the synaptic input to changing the neuronal output, while at the same time considering a more rigorous treatment of gradient descent (see also Surace et al., 2020).
The naturalgradient plasticity rule
We consider a neuron with somatic potential $V$ (above a baseline potential ${V}_{\mathrm{rest}}$) evoked by the spikes $x}_{i$ of $n$ presynaptic afferents firing at rates $r}_{i$. The presynaptic spikes of afferent i cause a train of weighted dendritic potentials ${w}_{i}^{\mathrm{d}}{x}_{i}^{\u03f5}$ locally at the synaptic site. The $x}_{i$ denotes the unweighted synaptic potential (USP) train elicited by the lowpassfiltered spike train $x}_{i$. At the soma, each dendritic potential is attenuated by a potentially nonlinear function that depends on the synaptic location:
The somatic voltage above baseline thus reads as
We further assume that the neuron’s firing follows an inhomogeneous Poisson process whose rate
depends on the current membrane potential through a nonlinear transfer function $\varphi $. In this case, spiking in a sufficiently short interval $[t,t+\mathrm{d}t]$ is Bernoullidistributed. The probability of a spike occurring in this interval (denoted as ${y}_{t}=1$) is then given by
which defines our generalized linear neuron model (Gerstner and Kistler, 2002). Here, we used $\mathit{x}}^{\u03f5$ as a shorthand notation for the USP vector $({x}_{i}^{\u03f5})$. The full inputoutput distribution then reads
where the probability density of the input $\mathit{x}}_{t}^{\u03f5$ is independent of synaptic weights (see also Sec. ‘Detailed derivation of the naturalgradient learning rule’). In the following, we drop the time indices for better readability.
In view of this stochastic nature of neuronal responses, we consider a neuron that strives to reproduce a target firing distribution ${p}^{\ast}\phantom{\rule{negativethinmathspace}{0ex}}\left(y{\mathit{x}}^{\u03f5}\right)$. The required information is received as teacher spike train $Y}^{\ast$ that is sampled from ${p}^{*}$ (see Sec. ‘Sketch for the derivation of the somatic naturalgradient learning rule’ for details). In this context, plasticity may follow a supervisedlearning paradigm based on gradient descent. The KullbackLeibler divergence between the neuron’s current and its target firing distribution
represents a natural cost function which measures the error between the current and the desired output distribution in an informationtheoretic sense. Minimizing this cost function is equivalent to maximizing the loglikelihood of teacher spikes, since $p}^{\ast$ does not depend on $\mathit{w}$. Using naive Euclideangradient descent with respect to the synaptic weights (denoted by $\mathrm{\nabla}}_{\mathit{w}}^{\mathrm{E}$) results in the wellknown errorcorrecting rule (Pfister et al., 2006),
which is a spikebased version of the classical perceptron learning rule (Rosenblatt, 1958), whose multilayer version forms the basis of the errorbackpropagation algorithm (Rumelhart et al., 1986). On the singleneuron level, a possible biological implementation has been suggested by Urbanczik and Senn, 2014, who demonstrated how a neuron may exploit its morphology to store errors, an idea that was recently extended to multilayer networks (Sacramento et al., 2017; Haider et al., 2021).
However, as we argued above, learning based on Euclideangradient descent is not unproblematic. It cannot account for synaptic weight (re)parametrization, as caused, for example, by the diversity of synaptic loci on the dendritic tree. Equation 10 represents the learning rule for a somatic parametrization of synaptic weights. For a local computation of Euclidean gradients using a dendritic parametrization, weight updates decrease with increasing distance towards the soma (Equation 3), which is likely to harm the convergence speed toward an optimal weight configuration. With the multiplicative USP term $\mathit{x}}^{\u03f5$ in Equation 10 being the only manifestation of presynaptic activity, there is no mechanism by which to take into account input variability. This can, in turn, also impede learning, by implicitly assigning equal importance to reliable and unreliable inputs. Furthermore, when compared to experimental evidence, this learning rule cannot explain heterosynaptic plasticity, as it is purely presynaptically gated.
In general, Euclideangradient descent is wellknown to exhibit slow convergence in nonisotropic regions of the cost function (Ruder, 2016), with such nonisotropy frequently arising or being aggravated by an inadequate choice of parametrization (see Ollivier, 2015 and Figure 2). In contrast, naturalgradient descent is, by construction, immune to these problems. The key idea of natural gradient as outlined by Amari is to follow the (locally) shortest path in terms of the neuron’s firing distribution. Argued from a normative point of view, this is the only ‘correct’ path to consider, since plasticity aims to adapt a neuron’s behavior, that is, its inputoutput relationship, rather than some internal parameter (Figure 2). In the following, we therefore drop the index from the synaptic weights $\mathit{w}$ to emphasize the parametrizationinvariant nature of the natural gradient.
For the concept of a locally shortest path to make sense in terms of distributions, we require the choice of a distance measure for probability distributions. Since a parametric statistical model, such as the set of our neuron’s realizable output distributions, forms a Riemannian manifold, a local distance measure can be obtained in form of a Riemannian metric. The Fisher metric (Rao, 1945), an infinitesmial version of the ${D}_{\mathrm{KL}}$, represents a canonical choice on manifolds of probability distributions (Amari and Nagaoka, 2000). On a given parameter space, the Fisher metric may be expressed in terms of a bilinear product with the Fisher information matrix
The Fisher metric locally measures distances in the $p$manifold as a function of the chosen parametrization. We can then obtain the natural gradient (which intuitively may be thought of as ‘$\mathrm{\partial}C/\mathrm{\partial}p$’) by correcting the Euclidean gradient $\mathrm{\nabla}}_{\mathit{w}}^{\mathrm{E}}C:=\mathrm{\partial}C/\mathrm{\partial}\mathit{w$ with the distance measure above:
This correction guarantees invariance of the gradient under reparametrization (see also Sec. ‘Reparametrization and the general naturalgradient rule’). The naturalgradient learning rule is then given as $\dot{\mathit{w}}=\eta {\mathrm{\nabla}}_{\mathit{w}}^{\mathrm{N}}C$. Calculating the righthand expression for the case of Poissonspiking neurons (for details, see Detailed derivation of the naturalgradient learning rule and Inverse of the Fisher Information Matrix), this takes the form
where $\mathit{w}$ is an arbitrary weight parametrization that relates to the somatic amplitudes via a componentwise rescaling
Note that the attenuation function in Equation 4 represents a special case of Equation 14 for dendritic amplitudes.
We used the shorthand ${\gamma}_{\mathrm{s}}$, ${\gamma}_{\mathrm{u}}$, and ${\gamma}_{\mathrm{w}}$ for three scaling factors introduced by the natural gradient that we address in detail below. For easier reading, we use a shorthand notation in which multiplications, divisions and scalar functions of vectors apply componentwise.
Equation 13 represents the complete expression of our naturalgradient rule, which we discuss throughout the remainder of the manuscript. Note that, while having used a standard sigmoidal transfer function throughout the paper, Equation 13 holds for every sufficiently smooth $\varphi $.
Naturalgradient learning conserves both the error term $\left[{Y}^{\phantom{\rule{thinmathspace}{0ex}}\ast}\varphi (V)\right]$ and the USP contribution $\mathit{x}}^{\u03f5$ from classical gradientdescent plasticity. However, by including the relationship between the parametrization of interest $\mathit{w}$ and the somatic PSP amplitudes $\mathit{f}\phantom{\rule{thinmathspace}{0ex}}(\mathit{w})$, naturalgradientbased plasticity explicitly accounts for reparametrization distortions, such as those arising from PSP attenuation during propagation along the dendritic tree. Furthermore, naturalgradient learning introduces multiple scaling factors and new plasticity components, whose characteristics will be further explored in dedicated sections below (see also Sec. ‘Global scaling factor’ and ‘Empirical Analysis of $\gamma}_{\mathrm{u}$ and $\gamma}_{\mathrm{w}$’ for more details).
First of all, we note the appearance of two scaling factors (more details in ‘Input and outputspecific scaling’). On one hand, the size of the synaptic adjustment is modulated by a global scaling factor ${\gamma}_{\mathrm{s}}$, which adjusts synaptic weight updates to the characteristics of the output nonlinearity, similarly to the synapsespecific scaling by the inverse of $\mathit{f}}^{\phantom{\rule{thinmathspace}{0ex}}\mathrm{\prime}$. Furthermore ${\gamma}_{\mathrm{s}}$ also depends on the output statistics of the neuron, harmonizing plasticity across different states in the output distribution (see Sec. ‘Global scaling factor’). On the other hand, a second, synapsespecific learning rate scaling accounts for the statistics of the input at the respective synapse, in the form of a normalization by the afferent input rate $c}_{\u03f5}/\mathit{r$, where ${c}_{\u03f5}$ is a constant that depends on the PSP kernel (see Sec. ‘Neuron model’). Unlike the global modulation introduced by ${\gamma}_{\mathrm{s}}$, this scaling only affects the USPdependent plasticity component. Just as for Euclideangradientbased learning, the latter is directly evoked by the spike trains arriving at the synapse. Therefore, the resulting plasticity is homosynaptic, affecting only synapses which receive afferent input.
However, in the case of naturalgradient learning, this inputspecific adaptation is complemented by two additional forms of heterosynaptic plasticity (Sec. ‘Interplay of homosynaptic and heterosynaptic plasticity’). First, the learning rule has a bias term ${\gamma}_{\mathrm{u}}$ which uniformly adjusts all synapses and may be considered homeostatic, as it usually opposes the USPdependent plasticity contribution. The amplitude of this bias does not exclusively depend on the afferent input at the respective synapse, but is rather determined by the overall input to the neuron. Thus, unlike the USPdependent component, this heterosynaptic plasticity component equally affects both active and inactive synaptic connections. Furthermore, naturalgradient descent implies the presence of another plasticity component ${\gamma}_{\mathrm{w}}\phantom{\rule{thinmathspace}{0ex}}\mathit{f}\left(\mathit{w}\right)$ which adapts the synapses depending on their current weight. More specifically, connections that are already strong are subject to larger changes compared to weaker ones. Since the proportionality factor ${\gamma}_{\mathrm{w}}$ only depends on global variables such as the membrane potential, this component also affects both active and inactive synapses.
The full expressions for ${\gamma}_{\mathrm{s}}$, ${\gamma}_{\mathrm{u}}$, and ${\gamma}_{\mathrm{w}}$ are functions of the membrane potential, its mean and its variance, which represent synapselocal quantities. In addition, ${\gamma}_{\mathrm{u}}$ and ${\gamma}_{\mathrm{w}}$ also depend on the total input $\sum _{i=1}^{n}{x}_{i}^{\u03f5}$ and the total instantaneous presynaptic rate $\sum _{i=1}^{n}{r}_{i}$. However, under reasonable assumptions such as a high number of presynaptic partners and for a large, diverse set of empirically tested scenarios, we have shown that these factors can be reduced to simple functions of variables that are fully accessible at the locus of individual synapses:
where ${c}_{\u03f5}$, ${c}_{\mathrm{w}}$ are constants (Sec. ‘Global scaling factor’ and ‘Empirical Analysis of $\gamma}_{\mathrm{u}$ and $\gamma}_{\mathrm{w}$’). The above learning rule along with closedform expressions for these factors represent the main analytical findings of this paper.
In the following, we demonstrate that the additional terms introduced in naturalgradientbased plasticity confer important advantages compared to Euclideangradient descent, both in terms of of convergence as well as with respect to biological plausibility. More precisely, we show that our plasticity rule improves convergence in a supervised learning task involving an anisotropic cost function, a situation which is notoriously hard to deal with for Euclideangradientbased learning rules (Ruder, 2016). We then proceed to investigate naturalgradient learning from a biological point of view, deriving a number of predictions that can be experimentally tested, with some of them related to in vivo observations that are otherwise difficult to explain with classical gradientbased learning rules.
Naturalgradient plasticity speeds up learning
Nonisotropic cost landscapes can easily be provoked by nonhomogeneous input conditions. In nature, these can arise under a wide range of circumstances, for elementary reasons that boil down to morphology (the position of a synapse along a dendrite can affect the attenuation of its afferent input) and function (different afferents perform different computations and thus behave differently). To evaluate the convergence behavior of our learning rule and compare it to Euclideangradient descent, we considered a very generic situation in which a neuron is required to map a diverse set of inputs onto a target output.
In order to induce a simple and intuitive anisotropy of the error landscape, we divided the afferent population into two equally sized groups of neurons with different firing rates (Figure 3A, firing rates were chosen as 10 Hz for group one and 50 Hz for group 2). The input spikes were lowpassfiltered with a difference of exponentials (see Sec. ‘Neuron model for details’). This resulted in an asymmetric cost function (KLdivergence between the student and the target firing distribution, Equation 9), as visible from the elongated contour lines (Figure 3D and E). We further chose a realizable teacher by simulating a different neuron with the same input populations connected via a predefined set of target weights $\mathit{w}}^{\ast$. For the weight path plots in Figure 3D and E fixed target weight ${\mathit{w}}^{*}={(\frac{0.3}{n},\frac{0.5}{n})}^{T}$ was chosen, whereas for the learning curve in Figure 3F, the target weight components were randomly sampled from a uniform distribution on $[1/n,1/n]$ (see Sec. ‘Supervised Learning Task’ for further details). Figure 3B and C shows that our naturalgradient rule enables the student neuron to adapt its weights to reproduce the teacher voltage $V}^{\phantom{\rule{thinmathspace}{0ex}}\ast$ and thereby its output distribution.
In the following, we compare learning in two student neurons, one endowed with Euclideangradient plasticity (Equation 10, Figure 3D) and one with our naturalgradient rule (Equation 13, Figure 3E). To better visualize the difference between the two rules, we used a twodimensional input weight space, that is, one neuron per afferent population. While the negative Euclideangradient vectors stand, by definition, perpendicular to the contour lines of $C$, the negative naturalgradient vectors point directly towards the target weight configuration ${\mathit{w}}^{*}$. Due to the anisotropy of $C$ induced by the different input rates (see also Figure 2B), Euclideangradient learning starts out by mostly adapting the highrate afferent weight and only gradually begins learning the lowrate afferent. In contrast, natural gradient adapts both synaptic weights homogeneously. This is clearly reflected by paths traced by the synaptic weights during learning.
Overall, this lead to faster convergence of the naturalgradient plasticity rule compared to Euclideangradient descent. In order to enable a meaningful comparison, learning rates were tuned separately for each plasticity rule in order to optimize their respective convergence speed. The faster convergence of naturalgradient plasticity is a robust effect, as evidenced in Figure 3F by the average learning curves over 1,000 trials.
In addition to the functional advantages described above, naturalgradient learning also makes some interesting predictions about biology, which we address below.
Democratic plasticity
As discussed in the introduction, classical gradientbased learning rules do not usually account for neuron morphology. Since attenuation of PSPs is equivalent to weight reparametrization and our learning rule is, by construction, parametrizationinvariant, it naturally compensates for the distance between synapse and soma. In Equation 13, this is reflected by a componentwise rescaling of the synaptic changes with the inverse of the attenuation function $\mathit{f}}^{\phantom{\rule{thinmathspace}{0ex}}\mathrm{\prime}$, which is induced by the Fisher information metric (see also Figure 8 and the corresponding section in the Materials and methods). Under the assumption of passive attenuation along the dendritic tree, we have
where d_{i} denotes the distance of the ith synapse from the soma. More specifically, $\mathit{\alpha}(\mathit{d})={e}^{\mathit{d}/\lambda}$, where λ represents the electrotonic length scale. We can write the naturalgradient rule as
For functionally equivalent synapses (i.e. with identical input statistics), synaptic changes in distal dendrites are scaled up compared to proximal synapses. As a result, the effect of synaptic plasticity on the neuron’s output is independent of the synapse location, since dendritic attenuation is precisely counterbalanced by weight update amplification.
We illustrate this effect with simulations of synaptic weight updates at different locations along a dendritic tree in Figure 4. Such ‘democratic plasticity’, which enables distal synapses to contribute just as effectively to changes in the output as proximal synapses, is reminiscent of the concept of ‘dendritic democracy’ (Magee and Cook, 2000). These experiments show increased synaptic amplitudes in the distal dendritic tree of multiple cell types, such as rat hippocampal CA1 neurons; dendritic democracy has therefore been presumed to serve the purpose of giving distal inputs a ‘vote’ on the neuronal output. Still, experiments show highly diverse PSP amplitudes in neuronal somata (Williams and Stuart, 2002). Our plasticity rule refines the notion of democracy by asserting that learning itself rather than its end result is rescaled in accordance with the neuronal morphology.
Whether such democratic plasticity ultimately leads to distal and proximal synapses having the same effective vote at the soma depends on their respective importance towards reaching the target output. In particular, if synapses from multiple afferents that encode the same information are randomly distributed along the dendritic tree, then democratic plasticity also predicts dendritic democracy, as the scaling of weight changes implies a similar scaling of the final learned weights. However, the absence of dendritic democracy does not contradict the presence of democratic plasticity, as afferents from different cortical regions might target specific positions on the dendritic tree (see, e.g., Markram et al., 2004). Furthermore, also with democratic plasticity in place, the functional efficacy of synapses at different locations on the dendritic tree will change differently based on different initial conditions. The experimental findings by Froemke et al., 2005, who report an inverse relationship between the amount of plasticity at a synapse and its distance from soma, are therefore completely consistent with the predictions of our learning rule, since they compared synapses that were not initially functionally equivalent.
Note also that dendritic democracy could, in principle, be achieved without democratic plasticity. However, this would be much slower than with naturalgradient learning, especially for distal synapses, as discussed in Sec. ‘Naturalgradient plasticity speeds up learning’.
Input and outputspecific scaling
In addition to undoing distortions induced by, for example, attenuation, the naturalgradient rule predicts further modulations of the homosynaptic learning rate. The factor ${\gamma}_{\mathrm{s}}$ in Equation 13 represents an outputdependent global scaling factor (for both homo and heterosynaptic plasticity):
Together with the ${\varphi}^{\prime}/\varphi $ factor in Equation 13, the effective scaling of the learning rule is approximately $1/{\varphi}^{\prime}$. In other words, it increases the learning rate in regions where the sigmoidal transfer function is flat (see also Sec. ‘Global scaling factor’). This represents an unmediated reflection of the philosophy of naturalgradient descent, which finds the steepest path for a small change in output, rather than in the numeric value of some parameter. The desired change in the output requires scaling the corresponding input change by the inverse slope of the transfer function.
Furthermore, synaptic learning rates are inversely correlated to the USP variance ${\sigma}^{2}({\mathit{x}}^{\u03f5})$ (Figure 5). In particular, for the homosynaptic component, the scaling is exactly equal to
(see Equation 13 and Sec. ‘Neuron model’). In other words, naturalgradient learning explicitly scales synaptic updates with the (un)reliability – more specifically, with the inverse variance – of their input. To demonstrate this effect in isolation, we simulated the effects of changing the USP variance while conserving its mean. Moreover, to demonstrate its robustness, we independently varied two contributors to the input reliability, namely input rates (which enter ${\sigma}^{2}({\mathit{x}}^{\u03f5})$ directly) and synaptic time constants (which affect the PSPkerneldependent scaling constant ${c}_{\u03f5}$). Figure 5 shows how unreliable input leads to slower learning, with an inverse dependence of synaptic weight changes on the USP variance. We note that this observation also makes intuitive sense from a Bayesian point of view, under which any information needs to be weighted by the reliability of its source. Furthermore, the approximately inverse scaling with the presynaptic firing rate is qualitatively in line with observations from Aitchison and Latham, 2014 and Aitchison et al., 2021, although our interpretation and precise relationship is different.
Interplay of homosynaptic and heterosynaptic plasticity
One elementary property of update rules based on Euclideangradient descent is their presynaptic gating, that is, all weight updates are scaled with their respective synaptic input $\mathit{x}}^{\u03f5$. Therefore, they are necessarily restricted to homosynaptic plasticity, as studied in classical LTP and LTD experiments (Bliss and Lomo, 1973; Dudek and Bear, 1992). As discussed above, naturalgradient learning retains a rescaled version of this homosynaptic contribution, but at the same time predicts the presence of two additional plasticity components. Contrary to homosynaptic plasticity, these components also adapt synapses to currently nonactive afferents, given a sufficient level of global input. Due to their lack of input specificity, they give rise to heterosynaptic weight changes, a form of plasticity that has been observed in hippocampus (Chen et al., 2013; Lynch et al., 1977), cerebellum (Ito and Kano, 1982), and neocortex (Chistiakova and Volgushev, 2009), mostly in combination with homosynaptic plasticity. A functional interpretation of heterosynaptic plasticity, to which our learning rule also alludes, is as a prospective adaptation mechanism for temporarily inactive synapses such that, upon activation, they are already useful for the neuronal output.
Our naturalgradient learning rule Equation 13 can be more summarily rewritten as
where the three additive terms represent the variancenormalized homosynaptic plasticity, the uniform heterosynaptic plasticity and the weightdependent heterosynaptic plasticity:
with the common proportionality factor $\eta {\gamma}_{\mathrm{s}}\left[{Y}^{\phantom{\rule{thinmathspace}{0ex}}\ast}\varphi (V)\right]\frac{{\varphi}^{\mathrm{\prime}}(V)}{\varphi (V)}{\mathit{f}}^{\phantom{\rule{thinmathspace}{0ex}}\mathrm{\prime}}{\left(\mathit{w}\right)}^{1}$ composed of the learning rate, the outputdependent global scaling factor, the postsynaptic error, a sensitivity factor and the inverse attenuation function, in order of their appearance. The effect of these three components is visualized in Figure 6B. The homosynaptic term $\mathrm{\Delta}{\mathit{w}}^{\mathrm{h}\mathrm{o}\mathrm{m}}$ is experienced only by stimulated synapses, while the two heterosynaptic terms act on all synapses. The first heterosynaptic term $\mathrm{\Delta}{\mathit{w}}^{\mathrm{h}\mathrm{e}{\mathrm{t}}_{\mathrm{u}}}$ introduces a uniform adjustment to all components by the same amount, depending on the global activity level. For a large number of presynaptic inputs, it can be approximated by a constant (see Sec. ‘Empirical analysis of $\gamma}_{\mathrm{u}$ and $\gamma}_{\mathrm{w}$’). Furthermore, it usually opposes the homosynaptic change, which we address in more detail below.
In contrast, the contribution of the second heterosynaptic term $\mathrm{\Delta}{\mathit{w}}^{\mathrm{h}\mathrm{e}{\mathrm{t}}_{\mathrm{w}}}$ is weightdependent, adapting all synapses in proportion to their current strength. This corresponds to experimental reports such as Loewenstein et al., 2011, which found in vivo weight changes in the neocortex to be proportional to the spine size, which itself is correlated with synaptic strength (Asrican et al., 2007). Our simulations in Sec. ‘Empirical Analysis of $\gamma}_{\mathrm{u}$ and $\gamma}_{\mathrm{w}$’ show that $\mathrm{\Delta}{\mathit{w}}^{\mathrm{h}\mathrm{e}{\mathrm{t}}_{\mathrm{w}}}$ is roughly a linear function of the membrane potential (more specifically, its deviation with respect to its baseline), which is reflected by the last approximation in Equation 21. Since the latter can be interpreted as a scalar product between the afferent input vector and the synaptic weight vector, it implies that spikes transmitted by strong synapses have a larger impact on this heterosynaptic plasticity component compared to spikes arriving via weak synapses. Thus, weaker synapses require more persistent and strong stimulation to induce significant changes and ‘override the status quo’ of the neuron. Since, following a period of learning, afferents connected via weak synapses can be considered uninformative for the neuron’s target output, this mechanism ensures a form of heterosynaptic robustness towards noise.
The homo and heterosynaptic terms exhibit an interesting relationship. To illustrate the nature of their interplay, we simulated a simple experiment (Figure 7A) with varying initial synaptic weights for both active and inactive presynaptic afferents. Stimulated synapses (Figure 7B) are seen to undergo strong potentiation (LTP) for very small initial weights; the magnitude of weight changes decreases for larger initial amplitudes until the neuron’s output matches its teacher, at which point the sign of the postsynaptic error term flips. For even larger initial weights, potentiation at stimulated synapses therefore turns into depression (LTD), which becoms stronger for higher initial values of the stimulated synapses’ weights. This is in line with the error learning paradigm, in which changes in synaptic weights seek to reduce the difference between a neuron’s target and its output.
For unstimulated synapses (Figure 7C), we observe a reversed behavior. For small weights, the negative uniform term $\mathrm{\Delta}{\mathit{w}}^{\mathrm{h}\mathrm{e}{\mathrm{t}}_{\mathrm{u}}}$ dominates and plasticity is depressing. As for the homosynaptic case, the sign of plasticity switches when the weights become large enough for the error to switch sign. Therefore, in the regime where stimulated synapses experienced potentiation, unstimulated synapses are depressed and viceversa. This reproduces various experimental observations: on one hand, potentiation of stimulated synapses has often been found to be accompanied by depression of unstimulated synapses (Lynch et al., 1977), such as in the amygdala (Royer and Paré, 2003) or the visual cortex (Arami et al., 2013); on the other hand, when the postsynaptic error term switches sign, depression at unstimulated synapses transforms into potentiation (Wöhrl et al., 2007; Royer and Paré, 2003).
While plasticity at stimulated synapses is unaffected by the initial state of the unstimulated synapses, plasticity at unstimulated synaptic connections depends on both the stimulated and unstimulated weights. In particular, when either of these grow large enough, the proportional term $\mathrm{\Delta}{\mathit{w}}^{\mathrm{h}\mathrm{e}{\mathrm{t}}_{\mathrm{w}}}$ overtakes the uniform term $\mathrm{\Delta}{\mathit{w}}^{\mathrm{h}\mathrm{e}{\mathrm{t}}_{\mathrm{u}}}$ and heterosynaptic plasticity switches sign again. Thus, for very large weights (top right corner of Figure 7C), heterosynaptic potentiation transforms back into depression, in order to more quickly quench excessive output activity. This behavior is useful for both supervised and unsupervised learning scenarios (Zenke and Gerstner, 2017), where it was shown that pairing Hebbian terms with heterosynaptic and homeostatic plasticity is crucial for stability.
In summary, we can distinguish three plasticity regimes for naturalgradient learning (Figure 7D–G). In two of these regimes, heterosynaptic and homosynaptic plasticity are opposed (O1, O2), whereas in the third, they are aligned and lead to depression (S). The two opposing regimes are separated by the zeroerror equilibrium line, at which plasticity switches sign.
Discussion
As a consequence of the fundamentally stochastic nature of evolution, it is no surprise that biology withstands confinement to strict laws. Still, physicsinspired arguments from symmetry and invariance can help uncover abstract principles that evolution may have gradually discovered and implemented into our brains. Here, we have considered parametrization invariance in the context of learning, which, in biological terms, translates to the fundamental ability of neurons to deal with diversity in their morphology and inputoutput characteristics. This requirement ultimately leads to various forms of scaling and heterosynaptic plasticity that are experimentally welldocumented, but can not be accounted for by classical paradigms that regard plasticity as Euclideangradient descent. In turn, these biological phenomena can now be seen as a means to jointly improve and accelerate errorcorrecting learning.
Inspired by insights from information geometry, we applied the framework of natural gradient descent to biologically realistic neurons with extended morphology and spiking output. Compared to classical errorcorrecting learning rules, our plasticity paradigm requires the presence of several additional ingredients. First, a global factor adapts the learning rate to the particular shape of the voltagetospike transfer function and to the desired statistics of the output, thus addressing the diversity of neuronal response functions observed in vivo (Markram et al., 2004). Second, the homosynaptic component of plasticity is normalized by the variance of presynaptic inputs, which provides a direct link to Bayesian frameworks of neuronal computation (Aitchison and Latham, 2014; Jordan et al., 2020). Third, our rule contains a uniform heterosynaptic term that opposes homosynaptic changes, downregulating plasticity and thus acting as a homeostatic mechanism (Chen et al., 2013; Chistiakova et al., 2015). Fourth, we find a weightdependent heterosynaptic term that also accounts for the shape of the neuron’s activation function, while increasing its robustness towards noise. Finally, our naturalgradientbased plasticity correctly accounts for the somatodendritic reparametrization of synaptic strengths.
These features enable faster convergence on nonisotropic error landscapes, in line with results for multilayer perceptrons (Yang and Amari, 1998; Rattray and Saad, 1999) and ratebased deep neural networks (Pascanu and Bengio, 2013; Ollivier, 2015; Bernacchia et al., 2018). Importantly, our learning rule can be formulated as a simple, fully local expression, only requiring information that is available at the locus of plasticity.
We further note an interesting property of our learning rule, which it inherits directly from the Fisher information metric that underlies natural gradient descent, namely invariance under sufficient statistics (Cencov, 1972). This is especially relevant for biological neurons, whose stochastic firing effectively communicates information samples rather than explicit distributions. Thus, downstream computation is likely to require a reliable samplebased, that is, statistically sufficient, estimation of the afferent distribution’s parameters, such as the sample mean and variance. This singles out our naturalgradient approach from other secondorderlike methods as a particularly appealing framework for biological learning.
Many of the biological phenomena predicted by our invariant learning rule are reflected in existing experimental results. Our democratic plasticity can give rise to dendritic democracy, as observed by Magee and Cook, 2000. Moreover, it also relates to results by Letzkus et al., 2006 and Sjöström and Häusser, 2006 which describe a sign switch of synaptic plasticity between proximal and distal sites that is not easily reconcilable with naive gradient descent. As the authors themselves speculate, the most likely reason is the (partial) failure of action potential backpropagation through the dendritic tree. In some sense, this is a mirror problem to the issue we discuss here: while we consider the attenuation of input signals, these experimental findings could be explained by an attenuation of the output signal. A simple way of incorporating this aspect into our learning rule (but also into Euclidean rules) is to apply a (nonlinear) attenuation function to the teacher signal $Y}^{\phantom{\rule{thinmathspace}{0ex}}\ast$ (Equations 10 and 13). At a certain distance from the soma, this attenuation would become strong enough to switch the sign of the error term ${Y}^{\phantom{\rule{thinmathspace}{0ex}}\ast}\varphi (V)$ and thereby also, for example, LTP to LTD. However, addressing this at the normative level of natural gradient descent as opposed to a direct tweak of the learning rule is less straightforward.
Our rule also requires heterosynaptic plasticity, which has been observed in neocortex, as well as in deeper brain regions such as amygdala and hippocampus (Lynch et al., 1977; Engert and Bonhoeffer, 1997; White et al., 1990; Royer and Paré, 2003; Wöhrl et al., 2007; Chistiakova and Volgushev, 2009; Arami et al., 2013; Chen et al., 2013; Chistiakova et al., 2015), often in combination with homosynaptic weight changes. Moreover, we find that heterosynaptic plasticity generally opposes homosynaptic plasticity, which qualitatively matches many experimental findings (Lynch et al., 1977; White et al., 1990; Royer and Paré, 2003; Wöhrl et al., 2007) and can be functionally interpreted as an enhancement of competition. For very large weights, heterosynaptic plasticity aligns with homosynaptic changes, pushing the synaptic weights back to a sensible range (Chistiakova et al., 2015), as shown to be necessary for unsupervised learning (Zenke and Gerstner, 2017). In supervised learning it helps speed up convergence by keeping the weights in the operating range.
These qualitative matches of experimental data to the predictions of our plasticity rule provide a welldefined foundation for future, more targeted experiments that will allow a quantitative exploration of the relationship between heterosynaptic plasticity and natural gradient learning. In contrast to the mostly deterministic protocols used in the referenced literature, such experiments would require a Poissonlike, stochastic stimulation of the student neuron. To further facilitate a quantitative analysis, the experimental setup should, in particular, allow a clear separation between student and teacher, for example via a pairing protocol (see e.g. Royer and Paré, 2003, but note that their protocol was deterministic).
A further prediction that follows from our plasticity rule is the normalization of weight changes by the presynaptic variance. We would thus anticipate that increasing the variance in presynaptic spike trains (through, e.g. temporal correlations) should reduce LTP in standard plasticity induction protocols. Also, we expect to observe a significant dependence of synaptic plasticity on neuronal response functions and output statistics. For example, flatter response functions should correlate with faster learning, in contrast to the inverse correlation predicted by classical learning rules derived from Euclideangradient descent. These propositions remain to be tested experimentally.
By following gradients with respect to the neuronal output rather than the synaptic weights themselves, we were able to derive a parametrizationinvariant errorcorrecting plasticity rule on the singleneuron level. Errorcorrecting learning rules are an important ingredient in understanding biological forms of error backpropagation (Sacramento et al., 2017). In principle, our learning rule can be directly incorporated as a building block into spikebased frameworks of error backpropagation such as (Sporea and Grüning, 2013; Schiess et al., 2016). Based on these models, topdown feedback can provide a target for the somatic spiking of individual neurons, toward which our learning rule could be used to speed up convergence.
Explicitly and exactly applying natural gradient at the network level does not appear biologically feasible due to the existence of crossunit terms in the Fisher information matrix $G$. These terms introduce nonlocalities in the naturalgradient learning rule, in the sense that synaptic changes at one neuron depend on the state of other neurons in the network. However, methods such as the unitwise naturalgradient approach (Ollivier, 2015) could be employed to approximate the natural gradient using a blockdiagonal form of $G$. The resulting learning rule would not exhibit biologically implausible codependencies between different neurons in the network, while still retaining many of the advantages of full naturalgradient learning. For spiking networks, this would reduce global naturalgradient descent to our local rule for single neurons.
Materials and methods
Neuron model
Request a detailed protocolWe chose a Poisson neuron model whose firing rate depends on the somatic membrane potential above the resting potential $V}_{\text{rest}}=70\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{V$. They relate via a sigmoidal activation function
with $\beta =0.3,\phantom{\rule{thinmathspace}{0ex}}\theta =10\mathrm{m}\mathrm{V}$, and a maximal firing rate $\varphi}_{\mathrm{m}\mathrm{a}\mathrm{x}}=100\phantom{\rule{thinmathspace}{0ex}}\mathrm{H}\mathrm{z$. This means the activation function is centered at −60 mV, saturating around −50 mV. Note that the derivation of our learning rule does not depend on the explicit choice of activation function, but holds for arbitrary monotonically increasing, positive functions that are sufficiently smooth. For the sake of simplicity, refractoriness was neglected. For the same reason, we assumed synaptic input to be currentbased, such that incoming spikes elicit a somatic membrane potential above baseline given by
where
denotes the unweighted synaptic potential (USP) evoked by a spike train
of afferent i, and w_{i} is the corresponding synaptic weight. Here, ${t}_{i}^{\mathrm{f}}$ denote the firing times of afferent i, and the synaptic response kernel $\u03f5$ is modeled as
where ${\u03f5}_{0}$ is a scaling factor with units $\text{}\mathrm{mV}\text{}\mathrm{ms}$, and unless specified otherwise, we chose ${\u03f5}_{0}=1\text{}\mathrm{mV}\text{}\mathrm{ms}$. For slowly changing input rates, mean and variance of the stationary unweighted synaptic potential are then given as (Petrovici, 2016)
and
with
Equations 29 and 30 are exact for constant input rates. Note that it is Equation 30 that induces the inverse dependence of the homosynaptic term on the input rate discussed around Equation 19.
Unless indicated otherwise, simulations were performed with a membrane time constant ${\tau}_{\mathrm{m}}=10\text{}\mathrm{ms}$ and a synaptic time constant ${\tau}_{\mathrm{s}}=3\text{}\mathrm{ms}$. Hence, USPs had an amplitude of 60 mV and were normalized with respect to area under the curve and multiplied by the synaptic weights. Initial and target weights were chosen such that the resulting average membrane potential was within operating range of the activation function. As an example, in Figure 3F, the average initial excitatory weight was 0.005, corresponding to an EPSP amplitude of 300V.
In Figure 5, the scaling factor ${\u03f5}_{0}$ was additionally normalized proportionally to the input rate r_{i} at the synapse in order to keep the mean USP constant and allow a comparison based solely on the variance.
Sketch for the derivation of the somatic naturalgradient learning rule
Request a detailed protocolThe choice of a Poisson neuron model implies that spiking in a small interval $[t,t+\mathrm{d}t]$ is governed by a Poisson distribution. For sufficiently small interval lengths $\mathrm{d}t$, the probability of having a single spike in $[t,t+\mathrm{d}t]$ becomes Bernoulli with parameter $\varphi \left({V}_{t}\right)\mathrm{d}t$. The aim of supervised learning is to bring this distribution closer to a given target distribution with density ${p}^{*}$. The latter is delivered in form of a teacher spike train
where the probability of having a teacher spike (denoted as ${Y}^{*}=1$) in $[t,t+\mathrm{d}t]$ is Bernoulli with parameter
Note that to facilitate the convergence analysis, in Figure 3 we chose a realizable teacher ${\varphi}_{t}^{\ast}=\varphi \left(\sum _{i=1}^{n}{w}_{i}^{\ast}{x}_{i}^{\u03f5}\right)$ generated by a given set of teacher weights ${\mathit{w}}^{*}$.
We measure the error between the desired and the current inputoutput spike distribution in terms of the KullbackLeibler divergence given in Equation 9, which attains its single minimum when the two distributions are equal. Note that while the ${D}_{\mathrm{KL}}$ is a standard measure to characterize ‘how far’ two distributions are apart, its behavior can sometimes be slightly unintuitive since it is not a metric. In particular, it is not symmetric and does not satisfy the triangle inequality.
Classical error learning follows the Euclidean gradient of this cost function, given as the vector of partial derivatives with respect to the synaptic weights. A short calculation (Sec. ‘Detailed derivation of the naturalgradient learning rule’) shows that the resulting Euclideangradientdescent learning rule is given by Equation 10. By correcting the vector of partial derivatives for the distance distortion between the manifold of inputoutput distributions and the synaptic weight space, given in terms of the Fisher information matrix $G\left(\mathit{w}\right)$, we obtain the natural gradient (Equation 12). We then followed an approach by Amari, 1998 to derive an explicit formula for the product on the right hand side of Equation 12.
In Sec. ‘Detailed derivation of the naturalgradient learning rule’, we show that given independent input spike trains, the Fisher information matrix defined in Equation 11 can be decomposed with respect to the vector of input rates $\mathit{r}$ and the current weight vector $\mathit{w}$ as
Here,
is the covariance matrix of the unweighted synaptic potentials, and ${c}_{1},{c}_{2},{c}_{3}$ are coefficients (see Equation 82 for their definition) depending on the mean ${\mu}_{V}$ and variance ${\sigma}_{V}^{2}$ of the membrane potential, and on the total rate
Through repeated application of the ShermanMorrisonFormula, the inverse of $G\left(\mathit{w}\right)$ can be obtained as
Here the coefficients ${\gamma}_{\mathrm{s}},{g}_{1},{g}_{2},{g}_{3},{g}_{4}$, which are defined in Equation 102, are again functions of mean and variance of the membrane potential, and of the total rate. Consequently, the naturalgradient rule in terms of somatic amplitudes is given by
Note that the formulas for,
arise from the product of the inverse Fisher information matrix and the Euclidean gradient, using $\mathbf{1}}^{T}{\mathit{x}}^{\u03f5}=\sum _{i=1}^{n}{x}_{i}^{\u03f5$ and ${{\mathit{w}}^{\mathrm{s}}}^{T}{\mathit{x}}^{\u03f5}=V$. Due to the complicated expressions for ${g}_{1},\mathrm{\dots},{g}_{4}$ (Equation 82), Equation 39 only provides limited information about the behavior of ${\gamma}_{\mathrm{u}}$ and ${\gamma}_{\mathrm{w}}$. Therefore, we performed an empirical analysis based on simulation data (Sec. ‘Empirical analysis of $\gamma}_{\mathrm{u}$ and $\gamma}_{\mathrm{w}$’; Figure 11). In a stepwise manner we first evaluated ${g}_{1},\mathrm{\dots},{g}_{4}$ under various conditions, which revealed that the products with g_{2} and g_{3} in Equation 39 are neglible in most cases compared to the other terms, hence
Furthermore, for a sufficient number of input afferents, we can approximate $g}_{1}\approx {q}^{1$. Since $q={c}_{\u03f5}{\u03f5}_{0}\mathbb{E}\left[\sum _{i=1}^{n}{x}_{i}^{\u03f5}\right]$, by the central limit theorem we have
for large $n$ and ${c}_{\mathrm{u}}={\u03f5}_{0}=1$. Moreover, while the variance of g_{4} across weight samples increases with the number and firing rate of input afferents, its mean stays approximately constant across conditions. This lead to the approximation
where ${c}_{\mathrm{w}}$ is constants across weights, input rates and the number of input afferents.
To evaluate the quality of our approximations, we tested the performance of learning in the setting of Figure 3 when ${\gamma}_{\mathrm{u}}$ and ${\gamma}_{\mathrm{w}}$ were replaced by their approximations (Equations 112; 113). The test was performed for several input patterns (Sec. ‘Evaluation of the approximated naturalgradient rule’, Sec. ‘Performance of the approximated learning rule’). It turned out that a convergence behavior very similar to naturalgradient descent could be achieved with ${c}_{\mathrm{u}}=0.95$, which worked much better in practice than ${c}_{\mathrm{u}}=1$. For ${c}_{\mathrm{w}}$ a choice of 0.05 which was close to the mean of ${c}_{\mathrm{w}}$ worked well.
For these input rate configurations and choices of constants, we additionally sampled the negative gradient vectors for random initial weights and USPs (Sec. ‘Evaluation of the approximated naturalgradient rule’, Sec. ‘Performance of the approximated learning rule’) and compared the angles and length difference between naturalgradient vectors and the approximation to the ones between natural and Euclidean gradient.
Reparametrization and the general naturalgradient rule
Request a detailed protocolTo arrive at a more general form of the naturalgradient learning rule, we consider a parametrization $\mathit{w}$ of the synaptic weights which is connected to the somatic amplitudes via a smooth componentwise coordinate change $\mathit{f}=\left(\phantom{\rule{thinmathspace}{0ex}}{f}_{1},\dots ,{f}_{n}\right)$, such that ${w}_{i}^{\mathrm{s}}={f}_{i}\left({w}_{i}\right)$. A Taylor expansion shows that small weight changes then relate via the derivative of $\mathit{f}$
On the other hand, we can also express the cost function in terms of $\mathit{w}$, with $C[\mathit{w}]={C}^{\phantom{\rule{thinmathspace}{0ex}}\mathrm{s}}[\phantom{\rule{thinmathspace}{0ex}}\mathit{f}({\mathit{w}}^{\mathrm{s}})]$, and directly calculate the Euclidean gradient of $C$ in terms of $\mathit{w}$. By the chain rule, we then have
Plugging this into Equation 43, we obtain an inconsistency: $\mathrm{\Delta}{w}_{i}^{\mathrm{s}}={f}_{i}^{\phantom{\rule{thinmathspace}{0ex}}\mathrm{\prime}}{\left({w}_{i}\right)}^{2}\mathrm{\Delta}{w}_{i}^{\mathrm{s}}$. Hence the predictions of Euclideangradient learning depend on our choice of synaptic weight parametrization (Figure 1).
In order to obtain the naturalgradient learning rule in terms of $\mathit{w}$, we first express the Fisher Information Matrix in the new parametrization, starting with Equation 60
Inserting both Equation 44 and Equation 47 into Equation 12, we obtain the naturalgradient rule in terms of $\mathit{w}$ (Equation 13). As illustrated in Figure 8, unlike for Euclideangradient descent, the result is consistent with Equation 43.
Simulation details
All simulations were performed in python and used the numpy and scipy packages. Differential equations were integrated using a forward Euler method with a time step of 0.5 ms.
Supervised learning task
Request a detailed protocolA single output neuron was trained to spike according to a given target distribution in response to incoming spike trains from n independently firing afferents. To create an asymmetry in the input, we chose one half of the afferents’ firing rates as 10 Hz, while the remaining afferents fired at 50 Hz. The supervision signal consisted of spike trains from a teacher that received the same input spikes. To allow an easy interpretation of the results, we chose a realizable teacher, firing with rate $\varphi \left({V}^{\phantom{\rule{thinmathspace}{0ex}}\ast}\right)$, where $V}^{\phantom{\rule{thinmathspace}{0ex}}\ast}={\mathit{w}}^{\ast T}{\mathit{x}}^{\u03f5$ for some optimal set of weights $\mathit{w}}^{\ast$. However, our theory itself does not include assumptions about the origin and exact form of the teacher spike train.
For the learning curves in Figure 3F, initial and target weight components were chosen randomly (but identically for the two curves) from a uniform distribution on $\mathcal{U}\left(1/n,1/n\right)$, corresponding to maximal PSP amplitudes between −600 μV and 600 μV for the simulation with $n=100$ input neurons. Learning curves were averaged over 1000 initial and target weight configurations. In addition, the minimum and maximum values are shown in Figure 9, as well as the mean Euclidean distance between student and teacher weights and between student and teacher firing rates. We did not enforce Dales’ law, thus about half of the synaptic input was inhibitory at the beginning but sign changes were permitted. This means that the mean membrane potential above rest covered a maximal range [−30 mV, 30 mV]. Learning rates were optimized as $\eta}_{\mathrm{n}}=6\ast {10}^{4$ for the naturalgradientdescent algorithm and $\eta}_{\mathrm{e}}=4.5\ast {10}^{7$ for Euclideangradient descent, providing the fastest possible convergence to a residual root mean squared error in output rates of 0.8 Hz. To confirm that the convergence behavior of both naturalgradient and Euclideangradient learning is robust, we varied the learning rate between $0.5{\eta}_{0}$ and $2{\eta}_{0}$ and measured the time until the $D}_{\mathrm{K}\mathrm{L}$ first reached a value of $5*{10}^{5}$ (Figure 9). Here $\eta}_{0}={\eta}_{\mathrm{e}$ for Euclideangradient descent and $\eta}_{0}={\eta}_{\mathrm{n}$ for naturalgradient descent.
Per trial, the expectation over USPs in the cost function was evaluated on a randomly sampled test set of 50 USPs that resulted from input spike trains of 250 ms. The expectation over output spikes was calculated analytically.
For the weight path simulation with two neurons (Figure 3D–E), we chose a fixed initial weight $\mathit{w}}_{0}={\left(\frac{0.3}{n},\frac{0.5}{n}\right)}^{T$, a fixed target weight $\mathit{w}}^{\ast}={\left(\frac{0.15}{n},\frac{0.15}{n}\right)}^{T$, and learning rates $\eta}_{\mathrm{n}}=2.5\ast {10}^{4$ and $\eta}_{\mathrm{e}}=7.\ast {10}^{7$. Weight paths were averaged over 500 trials of 6000s duration each.
The vector plots in Figure 3D–E display the average negative normalized natural and Euclideangradient vectors across 2000 USP samples per synapse ($n=2$) on a grid of weight positions on $\left[\frac{0.4}{n},\frac{0.6}{n}\right]}^{2$, with the first coordinate of the gridpoints in $\left\{\frac{0.28}{n},\frac{0.1}{n},\frac{0.08}{n},\frac{0.26}{n},\frac{0.44}{n}\right\}$ and the second in $\left\{\frac{0.22}{n},\frac{0.02}{n},\frac{0.26}{n},\frac{0,5}{n}\right\}$. Each USP sample was the result of a $1\text{}\mathrm{s}$ spike train at rate ${r}_{1}=10\text{}\mathrm{Hz}$ and ${r}_{2}=50\text{}\mathrm{Hz}$ respectively. The contour lines were obtained from 2000 samples of the ${D}_{\mathrm{KL}}$ along a grid on $\left[\frac{0.4}{n},\frac{0.6}{n}\right]}^{2$ (distance between two grid points in one dimension: $\frac{0.006}{n}$) and displayed at the levels 0.001, 0.003, 0.005, 0.009, 0.015, 0.02, 0.03, 0.04.
For the plots in Figure 3B–C, we used initial, final, and target weights from a sample of the learning curve simulation. We then randomly sampled input spike trains of 250 ms length and calculated the resulting USPs and voltages according to Equation 26 and Equation 25. The output spikes shown in the raster plot were then sampled from a discretized Poisson process with $\mathrm{d}t=5.*{10}^{4}$. We then calculated the PSTH with a bin size of 12.5 ms.
Distance dependence of amplitude changes
Request a detailed protocolA single excitatory synapse received Poisson spikes at 5 Hz, paired with Poisson teacher spikes at 20 Hz. The distance from the soma was varied between 0 μm and 460 μm. Learning was switched on for 5 s with an initial weight corresponding to 0.05 at the soma, corresponding to a PSP amplitude of 3 mV. Initial dendritic weights were scaled up with the proportionality factor $\alpha {(d)}^{1}$ depending on the distance from the soma, in order for input spikes to result in the same somatic amplitude independent of the synaptic position. Example traces are shown for $\alpha {(d)}^{1}=3$ and $\alpha {(d)}^{1}=7$.
Variance dependence of amplitude changes
Request a detailed protocolWe stimulated a single excitatory synapse with Poisson spikes, while at the same time providing Poisson teacher spike trains at 80 Hz. To change USP variance independently from mean, unlike in the other exercises, the input kernel in Equation 28 was additionally normalized by the input rate. USP variance was varied by either keeping the input rate at 10 Hz while varying the synaptic time constant ${\tau}_{s}$ between 1 ms and 20 ms, or fixing ${\tau}_{s}$ at 20 ms and varying the input rate between 10 Hz and 50 Hz.
Comparison of homo and heterosynaptic plasticity
Out of $n=10$ excitatory synapses of a neuron, we stimulated 5 by Poisson spike trains at 5 Hz, together with teacher spikes at 20 Hz, and measured weight changes after 60 s of learning. To avoid singularities in the homosynaptic term of learning rule, we assumed in the learning rule that unstimulated synapses received input at some infinitesimal rate, thus effectively setting $\mathrm{\Delta}{\mathit{w}}^{\mathrm{hom}}=0$. Initial weights for both unstimulated and stimulated synapses were varied between $\frac{1}{n}$ and $\frac{5}{n}$. For reasons of simplicity, all stimulated weights were assumed to be equal, and tonic inhibition was assumed by a constant shift in baseline membrane potential of 5 mV. Example weight traces are shown for initial weights of $\frac{1.6}{n},\frac{2.5}{n},$ and $\frac{4.5}{n}$ for both stimulated and unstimulated weights. The learning rate was chosen as $\eta =0.01$.
For the plots in Figure 7B C, we chose the diverging colormap matplotlib.cm.seismic in matplotlib. The colormap was inverted such that red indicates negative values representing synaptic depression and blue indicates positive values representing synaptic potentiation. The colormap was linearly discretized into 500 steps. To avoid too bright regions where colors cannot be clearly distinguished, the indices 243–257 were excluded. The colorbar range was manually set to a symmetric range that includes the max/min values of the data. In Figure 7D, we only distinguish between regions where plasticity at stimulated and at unstimulated synapses have opposite signs (light green), and regions where they have the same sign (dark green).
Approximation of learningrule coefficents
Request a detailed protocolWe sampled the values for ${g}_{1},\mathrm{\dots},{g}_{4}$ from Equation 82 for different afferent input rates. The input rate $r$ was varied between 5 Hz and 55 Hz for $n=100$ neurons. The coefficients were evaluated for randomly sampled input weights (20 weight samples of dimension $n$, each component sampled from a uniform distribution $\mathcal{U}\left(5/n,5/n\right)$).
In a second simulation, we varied the number $n$ of afferents between 10 and 200 for a fixed input rate of 20 Hz, again for randomly sampled input weights (20 weight samples of dimension $n$, each component sampled from a uniform distribution $\mathcal{U}\left(5/n,5/n\right)$).
In a next step, we compared the sampled values of g_{1} as a function of the total input rate $n*r$ to the values of the approximation given by $g}_{1}\approx {q}^{1$ ($r$ between 5 Hz and 55 Hz, $n$ between 10 and 200 neurons, 20 weight samples of dimension $n$, each component sampled from a uniform distribution $\mathcal{U}\left(5/n,5/n\right)$).
Afterwards, we plotted the sampled values of ${\gamma}_{\mathrm{u}}$ as a function of the approximation $s$ (Equation 111, $r$ between 5 Hz and 55 Hz, $n=100$, 20 weight samples of dimension $n$, each component sampled from a uniform distribution $\mathcal{U}\left(5/n,5/n\right)$, 20 USPsamples of dimension $n$ for each rate/weightcombination).
Next, we investigated the behavior of $\gamma}_{\mathrm{w}$ as a function of ${g}_{4}V$ ($r$ between 5 Hz and 55 Hz, $n=100$, 20 weight samples of dimension $n$, each component sampled from a uniform distribution $\mathcal{U}\left(5/n,5/n\right)$, 20 USPsamples of dimension $n$ for each rate/weightcombination), and in last step, as a function of ${c}_{w}V$ with a constant ${c}_{w}=0.05$.
Evaluation of the approximated naturalgradient rule
Request a detailed protocolWe evaluated the performance of the approximated naturalgradient rule in Equation 114 (with ${c}_{\mathrm{u}}=0.95$ and ${c}_{\mathrm{w}}=0.05$) compared to Euclideangradient descent and the full rule in Equation 13 in the learning task of Figure 3 under different input conditions (n=100, Group 1: $10\text{}\mathrm{Hz}$/ Group 2: 30 Hz, Group 1: 10 Hz/ Group 2: 50 Hz, Group 1: 20 Hz/ Group 2: 20 Hz, Group 1: 20 Hz/ Group 2: 40 Hz). The learning curves were averaged over 1000 trials with input and target weight components randomly chosen from a uniform distribution on $\mathcal{U}\left(1/n,1/n\right)$. Learning rate parameters were tuned individually for each learning rule and scenario according to Table 1. All other parameters were the same as for Figure 3F.
For the angle histograms in Figure 12AB, we simulated the natural, Euclidean and approximated natural weight updates for several input and initial weight conditions. Similar to the setup in Figure 3 we separated the $n=100$ input afferents in two groups firing at different rates (Group1/Group2: 10 Hz/10 Hz, 10 Hz/30 Hz, 10 Hz/50 Hz, 20 Hz/20 Hz, 20 Hz/40 Hz). For each input pattern, 100 Initial weight components were sampled randomly from a uniform distribution $\mathcal{U}\left(5/n,5/n\right)$, while the target weight was fixed at ${\mathit{w}}^{*}={(\frac{0.15}{n},\frac{0.15}{n})}^{T}$. For each initial weight, 100long input spike trains were sampled and the average angle between the naturalgradient weight update and the approximated naturalgradient weight update at $t=1\text{}\mathrm{s}$ was calculated. The same was done for the average angle between the natural and the Euclidean weight update.
Detailed derivation of the naturalgradient learning rule
Request a detailed protocolHere, we summarize the mathematical derivations underlying our naturalgradient learning rule (Equation 13). While all derivations in Sec. ‘Detailed derivation of the naturalgradient learning rule’ and Sec. ‘Inverse of the Fisher Information Matrix’ are made for the somatic parametrization and can then be extended to other weight coordinates as described in Sec. ‘Reparametrization and the general naturalgradient rule’, we drop the index $\mathrm{s}$ in ${\mathit{w}}^{\mathrm{s}}$ for the sake of readability.
Supervised learning requires the neuron to adapt its synapses in such a way that its inputoutput distribution approaches a given target distribution with density ${p}^{*}$. For a given input spike pattern $\mathit{x}$, at each point in time, the probability for a Poisson neuron to fire a spike during the interval $[t,t+\mathrm{d}t]$ (denoted as ${y}_{t}=1$) follows a Bernoulli distribution with a parameter ${\varphi}_{t}\mathrm{d}t=\varphi \left({V}_{t}\right)\mathrm{d}t$, depending on the current membrane potential. The probability density of the binary variable y_{t} on {0, 1}, describing whether or not a spike occurred in the interval $[t,t+\mathrm{d}t]$, is therefore given by
and we have
where $p}_{\mathrm{u}\mathrm{s}\mathrm{p}$ denotes the probability density of the unweighted synaptic potentials $\mathit{x}}_{t}^{\u03f5$. Measuring the distance to the target distribution in terms of the KullbackLeibler divergence, we arrive at
Since the target distribution does not depend on the synaptic weights, the negative Euclidean gradient of the $D}_{\mathrm{K}\mathrm{L}$ equals
We may then calculate
where Equation 55 follows from the fact that $y}_{t}\in \{0,1\$ and for Equation 56 we neglected the term of order $\mathrm{d}{t}^{2}$ which is small compared to the remainder. Plugging Equation 56 into Equation 52 leads to the Euclideangradient descent online learning rule, given by
Here,
is the teacher spike train. We obtain the negative natural gradient by multiplying Equation 57 with the inverse Fisher information matrix, since
with the Fisher information matrix $G\left(\mathit{w}\right)$ at $\mathit{w}$ being defined as
Exploiting that, just like for the target density, the density of the USPs does not depend on $\mathit{w}$, so
we can insert the previously derived formula Equation 56 for the partial derivative of the loglikelihood. Hence, using the tower property for expectation values and the definition of $p}_{\mathit{w}$ (Equation 49, Equation 50), Equation 60 transforms to
In order to arrive at an explicit expression for the naturalgradient learning rule, we further decompose the Fisher information matrix, which will then enable us to find a closed expression for its inverse.
Inspired by the approach in Amari, 1998, we exploit the fact that a positive semidefinite matrix is uniquely defined by its values as a bivariate form on any basis of ${\mathbb{R}}^{n}$. Choosing a basis for which the bilinear products with $G\left(\mathit{w}\right)$ are of a particularly simple form, we are able to decompose the Fisher Information Matrix by constructing a sum of matrices whose values as a bivariate form on the basis equal are equal to those of $G\left(\mathit{w}\right)$. Due to the structure of this particular decomposition, we may then apply wellknown formulas for matrix inversion to obtain $G{\left(\mathit{w}\right)}^{1}$.
Consider the basis $\mathcal{B}=\{\mathit{w},{\mathit{b}}_{1},\dots ,{\mathit{b}}_{n1}\}$ such that the vectors $\sqrt{{\mathrm{\Sigma}}_{\mathrm{u}\mathrm{s}\mathrm{p}}}\phantom{\rule{thinmathspace}{0ex}}\mathit{w},\sqrt{{\mathrm{\Sigma}}_{\mathrm{u}\mathrm{s}\mathrm{p}}}\phantom{\rule{thinmathspace}{0ex}}{\mathit{b}}_{1},\dots ,\sqrt{{\mathrm{\Sigma}}_{\mathrm{u}\mathrm{s}\mathrm{p}}}\phantom{\rule{thinmathspace}{0ex}}{\mathit{b}}_{n1}$ are orthogonal to each other. Here, $\mathrm{\Sigma}}_{\mathrm{u}\mathrm{s}\mathrm{p}$ denotes the covariance matrix of the USPs which in the case of independent Poisson input spike trains is given as
In this case, the matrix square root reduces to the componentwise square root. Note that for any $\mathit{b},{\mathit{b}}^{\mathbf{\prime}}\in \mathcal{B}$ with $\mathit{b}\ne {\mathit{b}}^{\mathbf{\prime}}$, the random variables $\mathit{b}}^{T}{\mathit{x}}_{t}^{\u03f5$ and $\mathit{b}}^{\mathrm{\prime}T}{\mathit{x}}_{t}^{\u03f5$ are uncorrelated, since
We make the mild assumptions of having small afferent populations firing at the same input rate, and that the basis $\mathcal{B}$ is constructed in such way that the basis vectors are not too close to the coordinate axes, such that the products $\mathit{b}}^{T}{\mathit{x}}_{t}^{\u03f5$ are not dominated by a single component. Then, for sufficiently large $n$, every linear combination of the random variables $\mathit{b}}^{T}{\mathit{x}}_{t}^{\u03f5$ and $\mathit{b}}^{\mathrm{\prime}T}{\mathit{x}}_{t}^{\u03f5$ is approximately normally distributed, thus, the two random variables follow a joint bivariate normal distribution. Furthermore, uncorrelated random variables that are jointly normally distributed are independent. Since functions of independent random variables are also independent, this allows us to calculate all products of the form $\mathit{b}}^{T}G\left(\mathit{w}\right){\mathit{b}}^{\mathbf{\prime}$ for $\mathit{b},{\mathit{b}}^{\mathbf{\prime}}\in \mathcal{B}\setminus \{\mathit{w}\}$, as we can transform the expectation of products
into products of expectations. Taking into account that $V(t)={\mathit{w}}^{T}{\mathit{x}}_{t}^{\u03f5}$ and $\mathbb{E}\left[{\mathit{x}}^{\u03f5}\right]=\mathit{r}$, we arrive at
Here, $I}_{1},{I}_{2},{I}_{3$ denote the generalized voltage moments given as
The integral formulas follow from the fact that for a large number of input afferents, and under the mild assumption of all synaptic weights roughly being of the same order of magnitude, the membrane potential approximately follows a normal distribution with mean $\mu}_{\mathrm{v}}={\u03f5}_{0}\sum _{i=1}^{n}{w}_{i}{r}_{i$ and variance $\sigma}_{\mathrm{v}}^{2}=\frac{1}{{c}_{\u03f5}}\sum _{i=1}^{n}{w}_{i}^{2}{r}_{i$. Here, $\u03f5}_{0}=1\text{mV ms$ and $c}_{\u03f5$ is given by Equation 31, see Sec. ‘Neuron model’.
Based on the above calculations, we construct a candidate $\stackrel{~}{G}(\mathit{w})$ for a decomposition of $G\left(\mathit{w}\right)$. We start with the matrix ${c}_{1}\left({\u03f5}_{0}^{2}\mathit{r}{\mathit{r}}^{T}+{\mathrm{\Sigma}}_{\mathrm{u}\mathrm{s}\mathrm{p}}\right)$, since its easy to see that
Exploiting the orthogonal properties according to which we constructed $\mathcal{B}$ to carefully add more terms such that also the other identities in Equation 73 hold, we arrive at
Here,
To check that indeed $\stackrel{~}{G}(\mathit{w})=G\left(\mathit{w}\right)$, it suffices to check the values of $\stackrel{~}{G}$ as a bilinear form on the basis $\mathcal{B}$.
Inverse of the Fisher information matrix
Request a detailed protocolAs the expectation of the outer product of a vector with itself, $G\left(\mathit{w}\right)$ is per construction symmetric and positive semidefinite. From the previous calculations, it follows that for elements $b$ of a basis $\mathcal{B}$ of ${\mathbb{R}}^{n}$ the products $\mathit{b}}^{T}G\left(\mathit{w}\right)\mathit{b$ are strictly positive. Hence $G\left(\mathit{w}\right)$ is positive definite and thus invertible. We showed that,
We introduce the notation.
Then,
For the following calculations, we will repeatedly use the identities
By the ShermanMorrisonWoodbury formula, the inverse of an invertible rank one correction $(A+\mathit{u}{\mathit{v}}^{T})$ of an invertible matrix A is given by
Applying this to invert $G\left(\mathit{w}\right)$, as a first step, we consider the term $M}_{1}+{M}_{2$ and identify ${M}_{1}=A$ and $M}_{2}=\mathit{u}{\mathit{v}}^{T$. Its inverse is given by
with
Applying the ShermanMorrisonWoodbury formula a second time, this time with ${M}_{1}+{M}_{2}=A$ and $M}_{3}=\mathit{u}{\mathit{v}}^{T$, we obtain
Here,
Plugging in the definitions of M_{3} and M_{1}, and using that
we arrive at
After resorting the terms and grouping, we obtain the inverse of the Fisher Information Matrix as
with
Note that instead of applying the ShermanMorrisonWoodbury formula twice one could also directly use the version for rank2 corrections, known as the Woodbury identity (Max, 1950). This can be seen by rewriting Equation 84 as
with
and
By the Woodbury identity for rank2 corrections, we then have
with $\stackrel{~}{U}={c}_{1}^{1}{\mathrm{\Sigma}}_{\mathrm{u}\mathrm{s}\mathrm{p}}^{1}U$.
Analysis of the learning rule coefficients
In order to gain an intuitive understanding of Equation 13 and to judge its suitability as an invivo plasticity rule, we require insights into the behavior of the coefficients ${\gamma}_{\mathrm{s}},{\gamma}_{\mathrm{u}},$ and $\gamma}_{\mathrm{w}$ under various circumstances.
Global scaling factor
Request a detailed protocolThe global scaling factor $\gamma}_{\mathrm{s}$ is given as
The above formula reveals that $\gamma}_{\mathrm{s}$ is closely tied to the firing nonlinearity of the neuron, as well as the statistics of output. The scaling by the inverse slope of the output nonlinearity amplifies the synaptic update in regions where a small change in weight and thus in the membrane potential would not lead to a noticeable change in output distribution. A further scaling with $\varphi$ additionally amplifies the synaptic change in highoutput regimes (Figure 10). This is in line with the spirit of natural gradient that ensures that the size of the synaptic weight update is homogeneous in terms of $D}_{\mathrm{K}\mathrm{L}$ change rather than in absolute weight terms. Furthermore, the rescaling is based on the average output statistics, which the synapse might have access to via backpropagating action potentials (Stuart et al., 1997), rather than an instantaneous value.
Empirical analysis of $\gamma}_{\mathrm{u}$ and $\gamma}_{\mathrm{w}$
Request a detailed protocolWhile the formula for $\gamma}_{\mathrm{s}$ provides a rather good intuition of this coefficients’ behavior, from the derivations in the previous sections, it becomes clear that such a straightforward interpretation is not readily available from the formulas for $\gamma}_{\mathrm{u}$ and $\gamma}_{\mathrm{w}$. Defined as
and
where $g}_{1},\dots {g}_{4$ are given in Equation 102, these coefficients depend, apart from the membrane potential and its first and second moments, on the total input and its mean the total rate. This raises the question whether the synapse, which in general only has limited access to global quantities, can implement $\gamma}_{\mathrm{u}$ and $\gamma}_{\mathrm{w}$. We therefore used an empirical analysis through simulations to obtain a more detailed insight.
As a starting point, we sampled $g}_{1},\dots ,{g}_{4$ for various input conditions (Figure 11A–D, refer to Sec. ‘Approximation of learningrule coefficents’ for simulation details). Here, we varied the afferent input rate $r$ between $5\text{}\mathrm{Hz}$ and $55\text{}\mathrm{Hz}$ for $n=100$ neurons and evaluated the value of the respective coefficient for a randomly sampled input weight. In a second simulation (Figure 11E–H, Sec. ‘Approximation of learningrule coefficents’), we varied the number $n$ of afferent inputs between 10 and 200 neurons for fixed input rate $20\text{}\mathrm{Hz}$, again with randomly chosen input weights. This revealed an approximately inverse proportional relationship between the average input rate $r$ and $g}_{1},{g}_{2$ and g_{3} respectively. Furthermore, these coefficients seemed to be approximately inversely proportional to the number $n$ of afferent inputs. However, this was not true for g_{4}, whose mean seemed to stay approximately constants across input rates, although the scattering across the mean value (for different weight samples) increased.
While the value of the membrane potential stays bounded also for a large number of inputs (due to the normalization of the synaptic weight range with $\frac{1}{n}$), the total sum of USPs increases with an increasing number of inputs. Therefore, for large enough $n$, the term ${g}_{3}V$ can be neglected, so
A closer look at the behavior of g_{1} shows that, for a sufficiently large number of neurons or high input rates, $g}_{1}\approx {q}^{1$, which we verified by simulations (Figure 11I, Sec. ‘Approximation of learningrule coefficents’). In consequence,
To verify this approximation, we sampled the values of $\gamma}_{\mathrm{u}$ for various conditions (Figure 11J, Sec. ‘Approximation of learningrule coefficents’) and compared them against the approximation in Equation 111, which confirmed our approximation. Since the input rate is the mean value of the USP, assuming large enough populations with the same input rate and a sufficient number of input afferents, by the central limit theorem we have
with $c}_{\mathrm{u}}={\u03f5}_{0$. However, in practice, for ${\u03f5}_{0}=1\text{}\mathrm{mV}\text{}\mathrm{ms}$, a learning behavior much closer to natural gradient was obtained when $c}_{\mathrm{u}$ was slightly smaller than 1 such as ${c}_{\mathrm{u}}=0.95$ (Sec. ‘Quadratic transfer function’).
As a starting point to approximate $\gamma}_{\mathrm{w}$, we noticed that the mean of g_{4} stayed approximately constant when varying the input rate or the number of input afferents. On the other hand, g_{2} rapidly tends to zero in those cases, so we assumed that $g}_{2}\sum _{i=1}^{n}{x}_{i}^{\u03f5$ stays either constant or goes to zero in the limit of large n. Since $c}_{\u03f5}{\u03f5}_{0}{g}_{2$ seemed to be rather small compared to g_{4}, we hypothesized ${\gamma}_{\mathrm{w}}\approx {g}_{4}V$, which was confirmed by simulations (Figure 11K, Sec. ‘Approximation of learningrule coefficents’). As a second step (Figure 11L, Sec. ‘Approximation of learningrule coefficents’), since g_{4} seemed to be constant in the mean, we approximated
where simulations with ${c}_{\mathrm{w}}=0.05$, close to the mean of g_{4} showed a learning behavior close to naturalgradient learning (Figure 12C–F). Replacing $\gamma}_{\mathrm{u}$ and $\gamma}_{\mathrm{w}$ in Equation 13 by the expressions in Equation 112 and Equation 113, we obtain the approximated naturalgradient rule
Performance of the approximated learning rule
Request a detailed protocolSimulations of natural, Euclidean and approximated naturalgradient weight updates for several input patterns and randomly sampled initial conditions (Sec. ‘Evaluation of the approximated naturalgradient rule’) showed that the average angles (both in the Euclidean metric and in the Fisher metric) between the true and approximated naturalgradient weight update were small compared to the average angle between Euclidean and naturalgradient weight update (Figure 12A–B). This was confirmed by the learning curves for several tested input conditions in the setting of Figure 3, since the performance of the approximation lay in between the natural and the Euclidean gradient’s performance (Figure 12C–F, simulation details: Sec. ‘Evaluation of the approximated naturalgradient rule’). It can hence be regarded as a tradeoff between optimal learning speed, parameter invariance and biological implementability. Note that plugging the above calculations into Equation 13 still keeps invariance under coordinatewise smooth parameter changes.
Quadratic transfer function
Request a detailed protocolWhile for all our simulations, we used the sigmoidal transfer function given in Equation 24, the derivations that lead to Equation 100 also hold for other choices of transfer function. A particularly simple and thereby instructive result can be obtained for a transfer function $\varphi $ that satisfies
This is the case for a rectified quadratic transfer function
where $\mathrm{\Theta}$ is the Heaviside step function. To prevent issues with division by zero and to ensure that Equation 115 holds for all relevant membrane voltages $V$, we assume $\theta \ll {\mu}_{\mathrm{v}}{\sigma}_{\mathrm{v}}$. Then, Equation 67 reduces to
where the right side no longer depends on $\mathit{w}$, thus yielding ${\gamma}_{\mathrm{w}}=0$. Furthermore, we have
and therefore ${c}_{2}={c}_{3}=0$ (see Equation 82). Plugging this into Equation 81, we obtain
with $\mathrm{\Sigma}}_{\mathrm{u}\mathrm{s}\mathrm{p}$ defined in Equation 35. Inverting this using the ShermanMorrisonWoodbury Formula, we arrive at
Inserting this version of the Fisher information matrix into Equation 12, we get ${\gamma}_{\mathrm{s}}=1$ and ${\gamma}_{\mathrm{w}}=0$ (see Equation 107 and 109), thus obtaining a simplified version of the naturalgradient learning rule as
with $q$ as in Equation 36. This is in line with our empirical approximation for $\gamma}_{\mathrm{u}$ for the case of the sigmoidal transfer function (Equation 111). It is interesting to note that, for a large number of inputs, we have $\gamma}_{\mathrm{u}}\to {c}_{\u03f5$ and the average of the parenthesis in Equation 123, taken in isolation, reduces to zero. However, this does not mean that, on average, there is no plasticity. Instead, we can say that, for a quadratic activation regime, the natural gradient is purely driven by the correlation between the postsynaptic error and the (normalized and centered) presynaptic input.
Continuoustime limit
Request a detailed protocolUnder the Poissonprocess assumption, firing a spike in one of the time bins is independent from the spikes in other bins, therefore the probability for firing in two disjunct time intervals $[{t}_{1}+\mathrm{d}t]$ and $[{t}_{2}+\mathrm{d}t]$ is given as
Since $\mathbb{E}{\left[\frac{\mathrm{\partial}\mathrm{log}{p}_{\mathit{w}}}{\mathrm{\partial}\mathit{w}}\right]}_{{p}_{\mathit{w}}}=0$, we have
so the Fisher information matrix is additive.
In the continuoustime limit, where the interval $[0,T]$ is decomposed as $[0,T]={\cup}_{i=0}^{k1}[{t}_{i},{t}_{i+1}]$, with ${t}_{0}=0,{t}_{k}=T$ and $k\u27f6\mathrm{\infty}$, we have $\mathrm{d}t=\frac{T}{k}\u27f60$. Therefore,
Then, under the assumption that $T$ is small and firing rates are approximately constant on $[0,T]$, for the Fisher information matrix, we have
Data availability
The data used to generate the figures and the simulation code are available on GitHub, (copy archived at swh:1:rev:4da5611348ea969aff461886f85842fe36a8571a).
References

Synaptic plasticity as Bayesian inferenceNature Neuroscience 24:565–571.https://doi.org/10.1038/s41593021008095

BookDifferential geometrical theory of StatisticsIn: Amari S, editors. In Differential Geometry in Statistical Inference. Hayward: Institute of Mathematical Statistics. pp. 19–94.

Natural Gradient Works Efficiently in LearningNeural Computation 10:251–276.https://doi.org/10.1162/089976698300017746

ConferenceFisher information and natural gradient learning in random deep networksIn The 22nd International Conference on Artificial Intelligence and Statistics. pp. 694–702.

BookMethods of Information Geometry. Translations of Mathematical MonographsAmerican Mathematical Society.

Synaptic strength of individual spines correlates with bound Ca2+calmodulindependent kinase IIThe Journal of Neuroscience 27:14007–14011.https://doi.org/10.1523/JNEUROSCI.358707.2007

ConferenceExact natural gradient in deep linear networks and its application to the nonlinear caseAdvances in Neural Information Processing Systems. pp. 5941–5950.

BookOptimal Decision Rules and Optimal InferenceAmerican Mathematical Society, Rhode Island. Translation from Russian.

Heterosynaptic plasticity prevents runaway synaptic dynamicsThe Journal of Neuroscience 33:15915–15929.https://doi.org/10.1523/JNEUROSCI.508812.2013

Heterosynaptic plasticity in the neocortexExperimental Brain Research 199:377–390.https://doi.org/10.1007/s0022100918595

Homeostatic role of heterosynaptic plasticity: models and experimentsFrontiers in Computational Neuroscience 9:89.https://doi.org/10.3389/fncom.2015.00089

Spatiotemporal credit assignment in neuronal population learningPLOS Computational Biology 7:e1002092.https://doi.org/10.1371/journal.pcbi.1002092

BookSpiking Neuron Models: Single Neurons, Populations, PlasticityCambridge University Press.https://doi.org/10.1017/CBO9780511815706

ConferenceLatent equilibrium: A unified learning theory for arbitrarily fast computation with arbitrarily slow neuronsAdvances in Neural Information Processing Systems.

ConferenceA natural policy gradientAdvances in Neural Information Processing Systems. pp. 1531–1538.

Learning rules for spike timingdependent plasticity depend on dendritic synapse locationThe Journal of Neuroscience 26:10420–10429.https://doi.org/10.1523/JNEUROSCI.265006.2006

ConferenceNatural Langevin dynamics for neural networksIn International Conference on Geometric Science of Information. pp. 451–459.https://doi.org/10.1007/9783319684451

Interneurons of the neocortical inhibitory systemNature Reviews. Neuroscience 5:793–807.https://doi.org/10.1038/nrn1519

Riemannian metrics for neural networks I: feedforward networksInformation and Inference 4:108–153.https://doi.org/10.1093/imaiai/iav006

Information and the accuracy attainable in the estimation of statistical parametersBulletin of Calcutta Mathematical Society 37:81–91.

Analysis of natural gradient descent for multilayer neural networksPhysical Review E 59:4523–4532.https://doi.org/10.1103/PhysRevE.59.4523

The perceptron: A probabilistic model for information storage and organization in the brainPsychological Review 65:386–408.https://doi.org/10.1037/h0042519

Somatodendritic Synaptic Plasticity and Errorbackpropagation in Active DendritesPLOS Computational Biology 12:e1004638.https://doi.org/10.1371/journal.pcbi.1004638

The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPsThe Journal of Neuroscience 13:334–350.

Supervised learning in multilayer spiking neural networksNeural Computation 25:473–509.https://doi.org/10.1162/NECO_a_00396

Action potential initiation and backpropagation in neurons of the mammalian CNSTrends in Neurosciences 20:125–131.https://doi.org/10.1016/s01662236(96)100758

On the choice of metric in gradientbased theories of brain functionPLOS Computational Biology 16:e1007640.https://doi.org/10.1371/journal.pcbi.1007640

Acute and longterm effects of MK801 on direct cortical input evoked homosynaptic and heterosynaptic plasticity in the CA1 region of the female ratThe European Journal of Neuroscience 26:2873–2883.https://doi.org/10.1111/j.14609568.2007.05899.x

Hebbian plasticity requires compensatory processes on multiple timescalesPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 372:20160259.https://doi.org/10.1098/rstb.2016.0259
Article and author information
Author details
Funding
European Union 7th Framework Programme (604102)
 Walter Senn
Horizon 2020 Framework Programme (720270)
 Walter Senn
 Mihai A Petrovici
Horizon 2020 Framework Programme (785907 and 945539)
 Walter Senn
 Mihai A Petrovici
Swiss National Fonds (personal grant 310030L_156863 (WS))
 Walter Senn
Swiss National Fonds (Sinergia grant CRSII5_180316)
 Walter Senn
Manfred Stärk Foundation (personal grant)
 Mihai A Petrovici
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work has received funding from the European Union 7th Framework Programme under grant agreement 604102 (HBP), the Horizon 2020 Framework Programme under grant agreements 720270, 785907 and 945539 (HBP) and the SNF under the personal grant 310030L_156863 (WS) and the Sinergia grant CRSII5_180316. Calculations were performed on UBELIX (http://www.id.unibe.ch/hpc), the HPC cluster at the University of Bern. We thank Oliver Breitwieser for the figure template, and Simone Surace and João Sacramento, as well as the other members of the Senn and Petrovici groups for many insightful discussions. We are also deeply grateful to the late Robert Urbanczik for many enlightening conversations on gradient descent and its mathematical foundations. Furthermore, we wish to thank the Manfred Stärk Foundation for their ongoing support.
Version history
 Received: January 13, 2021
 Accepted: February 21, 2022
 Version of Record published: April 25, 2022 (version 1)
 Version of Record updated: April 17, 2023 (version 2)
Copyright
© 2022, Kreutzer et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,411
 views

 213
 downloads

 6
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Medicine
Background:
Preterm birth is the leading cause of neonatal morbidity and mortality worldwide. Most cases of preterm birth occur spontaneously and result from preterm labor with intact (spontaneous preterm labor [sPTL]) or ruptured (preterm prelabor rupture of membranes [PPROM]) membranes. The prediction of spontaneous preterm birth (sPTB) remains underpowered due to its syndromic nature and the dearth of independent analyses of the vaginal host immune response. Thus, we conducted the largest longitudinal investigation targeting vaginal immune mediators, referred to herein as the immunoproteome, in a population at high risk for sPTB.
Methods:
Vaginal swabs were collected across gestation from pregnant women who ultimately underwent term birth, sPTL, or PPROM. Cytokines, chemokines, growth factors, and antimicrobial peptides in the samples were quantified via specific and sensitive immunoassays. Predictive models were constructed from immune mediator concentrations.
Results:
Throughout uncomplicated gestation, the vaginal immunoproteome harbors a cytokine network with a homeostatic profile. Yet, the vaginal immunoproteome is skewed toward a proinflammatory state in pregnant women who ultimately experience sPTL and PPROM. Such an inflammatory profile includes increased monocyte chemoattractants, cytokines indicative of macrophage and Tcell activation, and reduced antimicrobial proteins/peptides. The vaginal immunoproteome has improved predictive value over maternal characteristics alone for identifying women at risk for early (<34 weeks) sPTB.
Conclusions:
The vaginal immunoproteome undergoes homeostatic changes throughout gestation and deviations from this shift are associated with sPTB. Furthermore, the vaginal immunoproteome can be leveraged as a potential biomarker for early sPTB, a subset of sPTB associated with extremely adverse neonatal outcomes.
Funding:
This research was conducted by the Perinatology Research Branch, Division of Obstetrics and MaternalFetal Medicine, Division of Intramural Research, Eunice Kennedy Shriver National Institute of Child Health and Human Development, National Institutes of Health, U.S. Department of Health and Human Services (NICHD/NIH/DHHS) under contract HHSN275201300006C. ALT, KRT, and NGL were supported by the Wayne State University Perinatal Initiative in Maternal, Perinatal and Child Health.

 Computational and Systems Biology
 Genetics and Genomics
Runs of homozygosity (ROH) segments, contiguous homozygous regions in a genome were traditionally linked to families and inbred populations. However, a growing literature suggests that ROHs are ubiquitous in outbred populations. Still, most existing genetic studies of ROH in populations are limited to aggregated ROH content across the genome, which does not offer the resolution for mapping causal loci. This limitation is mainly due to a lack of methods for the efficient identification of shared ROH diplotypes. Here, we present a new method, ROHDICE, to find large ROH diplotype clusters, sufficiently long ROHs shared by a sufficient number of individuals, in large cohorts. ROHDICE identified over 1 million ROH diplotypes that span over 100 SNPs and are shared by more than 100 UK Biobank participants. Moreover, we found significant associations of clustered ROH diplotypes across the genome with various selfreported diseases, with the strongest associations found between the extended HLA region and autoimmune disorders. We found an association between a diplotype covering the HFE gene and hemochromatosis, even though the wellknown causal SNP was not directly genotyped or imputed. Using a genomewide scan, we identified a putative association between carriers of an ROH diplotype in chromosome 4 and an increase in mortality among COVID19 patients (Pvalue=1.82×10^{11}). In summary, our ROHDICE method, by calling out large ROH diplotypes in a large outbred population, enables further population genetics into the demographic history of large populations. More importantly, our method enables a new genomewide mapping approach for finding diseasecausing loci with multimarker recessive effects at a population scale.