A synaptic learning rule for exploiting nonlinear dendritic computation.

Information processing in the brain depends on the integration of synaptic input distributed throughout neuronal dendrites. Dendritic integration is a hierarchical process, proposed to be equivalent to integration by a multilayer network, potentially endowing single neurons with substantial computational power. However, whether neurons can learn to harness dendritic properties to realize this potential is unknown. Here, we develop a learning rule from dendritic cable theory and use it to investigate the processing capacity of a detailed pyramidal neuron model. We show that computations using spatial or temporal features of synaptic input patterns can be learned, and even synergistically combined, to solve a canonical nonlinear feature-binding problem. The voltage dependence of the learning rule drives coactive synapses to engage dendritic nonlinearities, whereas spike-timing dependence shapes the time course of subthreshold potentials. Dendritic input-output relationships can therefore be flexibly tuned through synaptic plasticity, allowing optimal implementation of nonlinear functions by single neurons.

Introduction

An essential role of each neuron in a circuit is to transform a barrage of synaptic input into a meaningful stream of action potential output. The majority of input to a neuron is received on the dendritic tree, potentially allowing transformations that extend far beyond the simple sum-and-threshold operation that is commonly assumed (London and Häusser, 2005; Silver, 2010; Poirazi and Papoutsi, 2020). Understanding the cellular basis of brain function therefore demands an understanding of how single-neuron input-output functions are constrained by physiology and how they can be tuned for processing information. Dendritic morphology controls the interactions between inputs distributed in space and time and determines the attenuation and shape of synaptic potentials propagated to the soma (Rall, 1967; Vetter et al., 2001; Williams and Stuart, 2002; Koch, 2002). Voltage-dependent conductances modulate these dynamics in a highly nonlinear manner (Stuart et al., 1997; Larkum et al., 1999; Magee, 1999; Schiller et al., 2000; Polsky et al., 2004; Nevian et al., 2007; Major et al., 2013), yielding modes of synaptic integration spanning sublinear to supralinear regimes (Cash and Yuste, 1999; Polsky et al., 2004; Tran-Van-Minh et al., 2015). In pyramidal neurons, it has been proposed that the combination of dendritic morphology and local NMDA receptor-dependent nonlinearities form a hierarchical processing structure with substantial computational power (Mel, 1992b; Archie and Mel, 2000; Poirazi and Mel, 2001). In the first stage of processing, synaptic input is integrated nonlinearly within individual dendrites, followed by a subsequent stage in which current flowing from dendrites is integrated at the soma, elegantly summarized as an equivalence of single neurons to multilayer neural networks (Poirazi et al., 2003a, 2003b; Jadi et al., 2014; Ujfalussy et al., 2018; Beniaguev et al., 2021; Jones and Kording, 2021a). However, this influential theory is incomplete, as it is yet to be comprehensively determined how or whether computations that capitalize on dendritic physiology can be learned.

Studies of dendritic processing have focused largely on how the responses of neurons differ when the rates, spatial distribution or timing of synaptic inputs are varied (London and Häusser, 2005; Silver, 2010). Computation in this context, such as selective responses to clustered input (Mel, 1992a; Cazé et al., 2013) or ordered sequences (Rall, 1964; Branco et al., 2010; Bhalla, 2017), is thus inherited from presynaptic firing patterns and connectivity. Although these computations may be accessible through precise programs of development or rearrangement of axonal connections (Mel, 1992a; Poirazi and Mel, 2001), a reliance on targeted wiring for exploiting dendritic mechanisms neglects synaptic plasticity as a dominant form of learning in the brain. Indeed, the concept of a neuron as a multilayer network instead evokes an image of a highly flexible device, whereby appropriate responses are learned through tuning of synaptic weights. Although many biological factors underlying plasticity are understood in molecular detail (Zucker and Regehr, 2002; Sjöström et al., 2008), at the computational level, the extent to which changes in synaptic strength can be harnessed to control the input-output function of a neuron is unknown.

Foundational theoretical studies have addressed this problem in simplified cell models, using supervised learning algorithms to investigate the limits of computation in single-compartment “point neurons” (Rosenblatt, 1958; Brunel et al., 2004; Legenstein et al., 2005; Gütig and Sompolinsky, 2006). More recent models that incorporate additional coarse-grained compartments to represent a dendritic tree have begun to reveal a variety of enhancements that may be accessible with a spatially extended morphology (Wu and Mel, 2009; Legenstein and Maass, 2011; O’Donnell and Sejnowski, 2014; Urbanczik and Senn, 2014; Schiess et al., 2016; Hawkins and Ahmad, 2016; Kastellakis et al., 2016; Guerguiev et al., 2017; Sacramento et al., 2018; Ujfalussy and Makara, 2020; Moldwin et al., 2021; Sezener et al., 2021; Jones and Kording, 2021a, 2021b). However, as model detail is often traded off against mathematical tractability in these approaches, it is possible that crucial features of neuronal function may have been lost, or artificially gained, by excluding important biological constraints. A realistic dendritic tree, for instance, presents not only opportunities for computation but also significant impediments that must be overcome (Häusser and Mel, 2003; Moldwin and Segev, 2020). Studies using detailed biophysical models have instead focused on exploring the consequences of experimentally motivated unsupervised rules (Bono and Clopath, 2017) or plasticity rules that implicitly assume only linear integration (Steuber et al., 2007; Moldwin and Segev, 2020), but whether these forms of learning realize the full potential of single neuron computation is unclear.

Here, we show that the computational power of dendritic processing can be flexibly exploited through synaptic plasticity. Although our approach is general, we focus on understanding how dendritic morphology and NMDA receptor-dependent excitability can be recruited without requiring structured connectivity. We analyze a detailed biophysical model of a layer 2/3 pyramidal neuron using a set of variational equations that reveal how the somatic membrane potential depends on the history of input to the dendritic tree. Using this to construct a local learning rule, we train the model to perform a nonlinear classification task. With rate-coded input, we find that synaptic weights evolve to selectively engage dendritic nonlinearities, consistent with the proposed equivalence to multilayer networks. We then extend this theory to encompass dendritic integration of temporal signals, revealing an effective spatiotemporal strategy for processing bursts of synaptic input. The trained models predict observable signatures of optimal processing by single pyramidal neurons, reflecting computations that are neither hand-tuned nor inherited but arise naturally when plasticity is governed by the constraints of physiology.

Results

To understand the dendritic computations that can be learned through synaptic plasticity, we developed a compartmental model of a pyramidal neuron and a learning rule to tune its synaptic weights. The main technical contribution of this paper is that we show how a synaptic credit assignment problem can be addressed by using cable theory to account for complex nonlinear interactions in the dendritic tree and that this can be approximated at each synapse using only local signals. By training the model to perform a task, we aim to uncover the computational strategies accessible to single neurons, without presupposing the ultimate implementation. Specifically, after introducing the model and learning rule, we consider a nonlinear feature-binding problem similar to that studied by Legenstein and Maass (2011) and Cazé et al. (2013). This requires the model neuron to classify conjunctions of synaptic input patterns by producing spikes in response to a prescribed set of inputs but remaining silent for others. As this biologically motivated task demands nonlinear processing of synaptic input, it is ideal for examining the ability of single neurons to learn network-level computations.

Synaptic integration in a layer 2/3 pyramidal neuron

We constructed the model using a detailed morphology of a layer 2/3 pyramidal neuron from mouse primary visual cortex (Allen Institute for Brain Science, 2015; ID 502359001) and fitted a somatic spiking mechanism with characteristic linear f-I curve and adaptive firing rate (Pospischil et al., 2008). Synapses were randomly distributed across basal and apical dendrites, comprising 800 excitatory synapses with AMPA and voltage-dependent NMDA conductances and 200 inhibitory synapses with GABA conductances. The voltage dependence of the NMDA conductance results in supralinear integration of excitatory synaptic input to single dendritic branches, once a threshold level of input is reached (Figure 1). This well-described nonlinearity (Schiller et al., 2000; Polsky et al., 2004; Nevian et al., 2007; Branco et al., 2010; Branco and Häusser, 2011) is the foundation of the multilayer-network analogy, with the responses of individual dendrites resembling the activation functions of artificial network units (Mel, 1992b; Poirazi et al., 2003b; London and Häusser, 2005). We refer to the model with excitable NMDA synapses as the “active model,” as the response to synaptic input is dominated by the voltage dependence of the NMDA conductance. In a subset of simulations, we extend this model to consider a more excitable regime in which voltage-dependent intrinsic conductances are also present throughout the dendritic tree (Figure S2). To isolate the contributions of NMDA-dependent excitability and morphology, we also consider two alternative models: a passive model in which we omit the voltage dependence of the NMDA receptors and a “point neuron” model that is biophysically identical to the active model but has all synapses located at the soma. In qualitative contrast to the active model, the passive model integrates sublinearly within dendritic branches because the dominant local effect of depolarization due to synaptic input is a reduction in synaptic driving force. The point neuron model, whose synapses are not subject to the high local voltages experienced in dendrites (Spruston, 2008), integrates approximately linearly at the soma (Figure 1).

Figure 1

Synaptic integration depends on dendritic morphology and local nonlinearity

Layer 2/3 pyramidal cell morphology with distinct basal and apical dendritic domains. Plots show the simulated peak somatic response to increasing numbers of excitatory synaptic inputs at the indicated locations, compared with the peak of the linear sum of the same number of unitary EPSPs. The basal and apical inputs are located at path distances of 95 and 275 μm from the soma, respectively. Voltage-dependent NMDA receptors yield supralinear integration within dendritic branches (active model, red lines), whereas integration in a purely passive model is sublinear (passive model, blue lines). Integration is approximately linear when the synapses of the active model are relocated to the soma (point neuron model, black line).

For a given pattern of input and synaptic weights (defined as peak synaptic conductance), the dynamics of the models are described by discretized dendritic cable equations, coupled to equations for the dynamics of active somatic conductances:In Equation 1, denotes the vector of voltages for each segment of the morphology. The nonlinear function describes the activation of synaptic conductances, generation of somatic spikes, and the flow of axial and transmembrane currents. The variables in Equation 2 represent Hodgkin-Huxley gating variables for the spiking mechanism (m, h: fast transient Na+; n: delayed rectifier K+; p: slow voltage-dependent ). The functions describe the gating kinetics, which depend on the somatic voltage . We provide the full details of the model dynamics in STAR Methods, (Equation 8), (Equation 9), (Equation 10), (Equation 11), (Equation 12).

For efficient learning via synaptic plasticity, a fundamental question must be addressed: during a barrage of synaptic input, what is the influence of each synapse on the spiking output of a neuron? When the solution to this credit assignment problem is known, powerful learning strategies can be used in which synaptic weights are selectively modified relative to their ability to control the output. A compelling demonstration of this principle, based on an integrate-and-fire point neuron, is provided by the tempotron learning rule (Gütig and Sompolinsky, 2006, 2009; Urbanczik and Senn, 2009). In the tempotron, synapses are modified in proportion to the gradient of the membrane potential with respect to the weights, allowing the voltage to be efficiently pushed above or below spiking threshold as demanded by a given computational task. Although proven to be highly effective for learning in simple point neurons, to our knowledge, fuller extensions that are able to exploit the complex dynamics of dendritic integration have not yet been developed. With this goal in mind, we augment Equations 1 and 2 with a system of variational equations for credit assignment that computes gradients of the dynamically evolving voltage with respect to the synaptic weights.

For the weight of synapse i, the variational equations are given byIntegrated numerically for each , these equations track the dependence over time of the voltage throughout the morphology on the synaptic weights and thereby the relative contributions of each synapse within the population. In particular, we are interested in the dependence that underlies the spiking output at the soma, . When the specific identity of a synapse is not important, we suppress the index and express this general quantity as . Because the derivatives are propagated through the model dynamics, we capture the influence of dendritic morphology, interactions between excitatory and inhibitory synapses, dendritic NMDA spikes, and active somatic integration. We provide the full details of the variational equations in STAR Methods, (Equation 13), (Equation 14), (Equation 15), (Equation 16), (Equation 17).

Using the numerical solution of Equations 3 and 4 to visualize at times preceding action potentials, we find the influence of a synapse varies considerably with location, time of activation, and local interactions within the dendritic tree (Figure 2A). For instance, if excitatory synapses are active in spatial clusters, the impacts of weight changes are amplified because of voltage-dependent coupling of NMDA currents within the local branch (Figure 2A, example i). Potent inhibitory control is seen for inhibitory synapses that are proximal to active excitatory synapses, relative to the soma (Figure 2A, example ii, morphology lower right), and also on distal dendritic tips (Figure 2A, example ii, morphology upper left), in agreement with previous work (Gidon and Segev, 2012). By extending the standard compartmental equations to track internal model dynamics, we thus access an explicit readout of the interactions between input statistics and dendritic biophysics, and the influence of every synapse on the soma.

Figure 2

Two local variables determine the impact of synaptic plasticity on somatic output

(A) Example simulation of the active model stimulated with Poisson input into excitatory (black) and inhibitory (magenta) synapses. Bottom left: somatic voltage trace. Top left: raster plots of synaptic input preceding two somatic spikes. Synapses located on the same dendrite are grouped together on the y axis. Markers are scaled by the magnitude of influence on the somatic voltage, , immediately prior to the spike, normalized by the maximum within excitatory and inhibitory groups. In this example, the variational equations were solved numerically for each individual synaptic activation by making dummy copies of synapses that were active more than once. Right: spatial distribution of activated synapses from example (ii).

(B) Polynomial fits of somatic spike-triggered average in the active model, to be used as plasticity kernels in the learning algorithm.

(C) The approximations in (B) accurately predict the voltage gradients computed from numerical integration of (Equation 13), (Equation 14), (Equation 15), (Equation 16), (Equation 17) (fitted on 75% of the simulated data and tested on the remaining 25%). For visibility, scatterplot shows randomly sampled points from bins of 0.1 mV nS−1 width along the x axis (up to 100 points per bin). R2 values are computed from the correlation between actual and approximated values over all held-out data.

(D) The voltage at a synapse at the time of somatic spikes depends on multiple factors, allowing their implicit representation in the learning rule. Shown is the semipartial correlation computed from a linear model fitted on 75% of the data and tested on the remaining 25%.

See also Figure S1.

Local variables can assign synaptic credit for somatic spikes

Equations 3 and 4 provide a general means to solve the problem of synaptic credit assignment in single neurons, without sacrificing the rigorous constraints and predictive power of biophysical simulations. However, there are two limitations to constructing a learning rule directly from this system. First, the inherent complexity means that numerical solution over potentially many training iterations is computationally expensive. Second, the equations are non-local, as the computed influence of a synapse formally depends on activity in even very distant regions of the dendritic tree, information unlikely to be generally available to synapses of a real neuron. To overcome these issues and gain further insight into the underlying principles, we sought an approximation using variables local to each synapse.

To construct the approximation, we performed an extensive set of simulations in which input rates, weights, and synaptic locations were varied. Pooling over synapses in either basal or apical dendrites, we computed a postsynaptic spike-triggered average of as a function of the time that a given synapse was activated before the postsynaptic spike, and the local dendritic voltage at the time of the postsynaptic spike. These functions could be well described by polynomial fits (Figure 2B, active model shown). Approximating values of at the time of postsynaptic spikes using the fitted spike-triggered averages accounted for ∼90% of the variance of held-out data (Figure 2C). The temporal components of the fits resemble time-reversed and cable-filtered synaptic conductances; the window of excitatory integration reflects a mix of AMPA and NMDA kinetics, and the peak influence of apical synapses is shifted to earlier activation times and reduced in magnitude relative to basal synapses (Figures S1A–S1C). In the active model, the peaks in the voltage components for both excitatory and inhibitory synapses coincide with the voltage-dependent activation of NMDA receptors between ∼−30 and −40 mV, revealing a small time and voltage window in which synaptic weight changes are maximally effective at modifying somatic output. The results for the passive model are similar, but with the peak in the voltage components for the excitatory synapses shifted to the somatic spiking threshold of ∼−50 mV (Figure S1D). We approximated for the point neuron using only synaptic activation time, as in this case the local voltage is the somatic membrane potential itself (Figure S1C).

The output of Equations 3 and 4 can be accurately captured in terms of just two local variables, synaptic activation time and dendritic voltage, which is surprising given the complex morphology and impact of features such as spatiotemporal clustering of inputs (Figure 2A). This important sensitivity is conferred through the local voltage dependence, as we find that such features are represented indirectly through the voltage signal within each branch. We used linear regression to quantify the unique contributions of key structural and dynamic variables to the local voltage by their semipartial correlation (i.e., ability to predict the residual voltage after removing components explained by other variables). At the time of somatic spikes, the voltage at a synapse depends strongly on input resistance (related to its finer-scale location), and the degree of synchronous activation of other synapses on the same branch (within a 100 ms window) (Figure 2D). The voltage depends comparatively weakly on the weight or kinetic state (value of the temporal component of conductance) of a synapse itself, once other variables are accounted for (Figure 2D). This implies that sets of both excitatory and inhibitory synapses can be strongly coupled by the voltage within branches, and in principle, voltage-dependent plasticity can thereby act with sensitivity to the spatiotemporal structure of the input. Learning controlled by local variables can therefore harness the morphology and dynamics of a pyramidal neuron without an explicit representation of either.

We use these results to construct a supervised synaptic learning rule, inspired by the tempotron (Gütig and Sompolinsky, 2006), that accounts for nonlinear processing in the dendritic tree. In the following, we consider an application where a neuron should produce somatic spikes on presentation of noisy input patterns from a preferred class (+) but remain silent during presentation of patterns from a nonpreferred class (−). The goal of learning is to minimize the expected number of incorrect classifications across the full set of input patterns.

We use a greedy algorithm in which input patterns are presented sequentially and synaptic plasticity acts to either suppress or reinforce somatic spiking activity. We assume that a supervisory system maintains a running average of the classification error for each input pattern and uses this to guide the learning process. On presentation of (+) patterns, the supervisor encourages spiking by depolarizing the soma. On both (+) and (−) patterns, the supervisor controls the sign and globally scales the magnitude of plasticity, such as via a neuromodulator (Seol et al., 2007). Synaptic weights are modified whenever a somatic action potential is fired, which we assume is signaled throughout the dendrites by action potential backpropagation (Stuart and Sakmann, 1994; Stuart et al., 1997). We use the local approximations of as plasticity kernels that assign the specific weight change for each synapse as a function of synaptic activation time and local dendritic voltage. Altogether, the weight update rule is given byThe first factor in Equation 5 is composed of a learning rate α and a running average error provided by the supervisory system ( for (+) patterns and for (−) patterns). This term adaptively modulates the magnitude of plasticity, such that weight updates are largest when the neuron has been performing poorly, and learning ceases once the noisy input patterns can be classified without error. In Figure S5 we show that this term can also be computed online from the spiking output of the neuron, in which case the supervisor need only provide a binary classification label. The second factor in Equation 5 acts locally to assign synaptic credit for the somatic spike; denotes the appropriate plasticity kernel for synapse i (basal, apical; excitatory, inhibitory), is the time of the kth synaptic input relative to the somatic spike, and is the local dendritic voltage at the time of the somatic spike. Thus, for example, an excitatory synapse in the basal dendrites that was activated 40 ms before a somatic spike, and subject to a local voltage of −40 mV at the time of the somatic spike, would be updated in proportion to the value read from the top left panel of Figure 2B: , after unit conversions. Because the plasticity kernels approximate the gradient of somatic voltage with respect to the weights, synapses are modified in proportion to their relative ability to reinforce or suppress somatic spikes.

Can a single neuron learn nonlinear functions?

We challenged the active, passive, and point neuron models with a nonlinear feature-binding task, framed as binary classification of combinations of synaptic input patterns (Legenstein and Maass, 2011; Cazé et al., 2013; Tran-Van-Minh et al., 2015). This task is an abstraction of the tendency of cortical neurons to respond to selective conjunctions of variables, such as spatial orientation and frequency, or multi-modal stimuli (Leventhal et al., 1995). We first focused on the rate-coded input regime on which previous theories of pyramidal cell computation are based (Poirazi et al., 2003b; Jadi et al., 2014) and chose classification labels such that the task could not be solved by linear integration (Legenstein and Maass, 2011). We consider streams of input from two presynaptic populations representing different classes of stimulus features, for concreteness, shape, and color. The neuron must learn to spike on specific associations of features from each class, such as “green triangle” and “orange square,” but no other combinations, such as “green square,” even though similar groups of synapses are active on both preferred and nonpreferred associations (Figure 3A). Note that this task is similar to learning an exclusive OR (XOR), as historically studied in a machine learning context (Minsky and Papert, 2017), but differs by requiring selective responses to combinations of different input patterns, rather than input magnitudes. Biologically, this corresponds to the presynaptic populations encoding stimulus identity through patterns of activity rather than the binary presence or absence of stimuli as in the XOR. Previously proposed solutions to this problem require strong assumptions about plasticity of the coupling of dendritic branches to the soma (Legenstein and Maass, 2011), motivated by observations in hippocampal neurons (Losonczy et al., 2008), or structured connectivity that organizes specific inputs into clusters (Cazé et al., 2013; Tran-Van-Minh et al., 2015). Here, we neither assume a biophysical implementation nor hard-wire a solution through connectivity but let synaptic plasticity act on randomly distributed inputs to learn the task.

Figure 3

A single neuron can learn nonlinear functions

(A) Nonlinear feature-binding problem. Synapses representing different stimulus features were randomly distributed throughout basal and apical dendrites. In this example, the neuron should only spike in response to the associations “green triangle” and “orange square” as indicated by the classification labels (bottom).

(B) Example simulations of a model before (gray) and after (black) training on the task defined in (A). Each combination of features is presented in turn via rate-coded Poisson input, interspersed with background noise. For clarity, only input to excitatory synapses is shown.

(C) Performance (fraction correct) of models trained on ten random instantiations of the task (left bars). In the somatic inhibition condition (middle bars), models were trained with all inhibitory synapses placed at the soma. Performance collapsed when dendritic voltage dependence was omitted from the learning rule (right bars).

(D) Classification of associations is made by differential supralinear or sublinear integration. Input was presented to the indicated domains of trained models with somatic spiking blocked. The peak somatic depolarization measured when features were presented together was compared with the sum of responses when presented independently (averaged over 20 presentations of each association pair, then over label types).

All bars denote means; p values are from two-tailed Wilcoxon signed-rank tests between groups for n = 10 independent replications. See also Figures S2 and S3.

To map the task to patterns of synaptic input, we defined each feature component by a sparse, randomly generated vector of input rates taking values of either 0 or 40 Hz (see STAR Methods, (Equation 18), (Equation 19), (Equation 20)). Input patterns were realized in simulations as Poisson spike trains comprising an initial 100 ms of background activity to randomize the initial conditions, followed by a 400 ms “stimulus presentation” of pairs of features. We used a sparseness of 1/16 (fraction of active synapses) to construct the rate vectors, giving an average of 1,000 presynaptic spikes per stimulus presentation (i.e., one spike per synapse, facilitating later comparison with candidate temporal coding schemes). We confirmed that our results do not depend qualitatively on these specific parameter choices, which are comparable with related work (Poirazi et al., 2003b; Legenstein and Maass, 2011). Pairs of features were presented in random order over 1,000 epochs of training and then tested 20 times without input from the supervisor. After training, the active model could robustly discriminate preferred associations from nonpreferred associations or background noise, demonstrating the ability of our learning rule to tune the synaptic weights of a detailed biophysical model for nonlinear computation (Figure 3B).

Which mechanisms are recruited to solve this task? Comparing classification performance across the three models confirms the advantage of dendritic morphology, with the point neuron performing at chance levels (Figure 3C; active, 81% ± 6%; passive, 78% ± 6%; point neuron, 53% ± 4%; mean ± SD). Dendritic excitability, however, is not essential, as the passive model performed almost as well as the active model. The performance of both models remained high when trained with all inhibitory synapses placed at the soma (Figure 3C; active, 81% ± 4%; passive, 79% ± 5%; mean ± SD), pointing to the engagement of predominantly excitatory nonlinearities. The voltage dependence of the learning rule was crucial for recruiting these nonlinearities, as performance collapsed when restricting it to depend only on synaptic activation time (Figure 3C; active, 56% ± 7%; passive, 55% ± 5%; mean ± SD).

To explore the possible roles of supralinear and sublinear modes of integration, we simulated the response of the trained models to different feature components, individually and in pairs, while blocking somatic spiking. Summation of pairs of features was supralinear in basal dendrites of the active model for all associations, but substantially more so for the preferred class, relative to nonpreferred (Figure 3D). By contrast, summation in passive basal dendrites was sublinear, with a small but consistent bias toward stronger sublinearity on nonpreferred associations, and similarly for passive apical dendrites. Unexpectedly, apical dendrites of the active model also exhibited a pronounced sublinear response, reflecting a level of depolarization such that the saturating phase of the sigmoid nonlinearity was engaged. We found similar results when the model was extended to include dendritic Na+, K+, and HCN channels to increase the intrinsic excitability of the dendritic tree (Figure S2). Therefore, in agreement with the central premise of Poirazi et al. (2003b), selective engagement of active mechanisms in basal dendrites can indeed be harnessed for nonlinear computation. However, voltage-dependent plasticity can also recruit sublinear forms of processing, and with this flexibility, the optimal apical strategy predicted for this morphology in fact resembles passive integration.

Examining the weights of the trained models, we find that the learning rule rediscovers the clustering principles proposed in terms of structured connectivity in the binary neuron model of Cazé et al. (2013). Cazé et al. (2013) proposed that with an expansive nonlinearity, excitatory inputs from two features forming a preferred association should be spatially clustered to enable supralinear integration. Inputs from two features forming a nonpreferred association should instead be dispersed, thereby ensuring only linear integration and a smaller somatic response (Figure 4A, left). Conversely, with a suppressive nonlinearity, inputs forming preferred associations should be dispersed, with linear integration yielding a larger somatic response than the sublinear integration of clustered inputs (Figure 4A, right). Translated to our setting, in which connectivity is random, we find these strategies can be realized functionally through learning by tuning the spatial distribution of input strength (Figure 4B). To quantify this, we constructed a spatial profile of excitation strength for each feature in an association pair by summing the weighted excitatory input rates on each dendritic branch and then computed the correlation between profiles. Comparing profiles in basal and apical domains separately, differences in functional clustering, as defined by the overlap in patterns of excitation strength, reflect the pattern of supra- or sublinear summation observed in the simulations (Figure 4C). Comparing profiles of excitation and inhibition in a similar manner, we find that they overlap selectively in basal dendrites of the active model on nonpreferred associations, serving to further suppress the response to those association pairs (Figure 4D).

Figure 4

Learning tunes the spatial distribution of input strengths

(A) Schematic illustrating how the feature-binding task can be solved via the structured connectivity mechanism proposed by Cazé et al. (2013). Arrows denote the targeting of excitatory input from a given feature to a compartment. With supralinear integration, responses to clustered input are enhanced relative to responses to dispersed input, and conversely with sublinear integration.

(B) When connectivity is random, the strategies of (A) can be realized functionally by tuning synaptic weights through learning. Left: spatial distribution of excitatory input strengths (weight × input rate) in a trained model. Inputs are color-coded by the features they represent and the classification labels defined by the matrix below. Right: profiles of excitatory input strength for the model depicted on the left. The height of each point is proportional to the sum of weighted input rates on a branch. Preferred associations (e.g., and ; blue and green) have strong inputs to common basal dendrites but separate apical dendrites. Conversely, strong inputs of nonpreferred associations (e.g., and ; blue and orange) are dispersed in basal dendrites and clustered in apical dendrites.

(C) Left: functionally clustered or dispersed configurations are reflected in the spatial correlation between weighted input to dendritic branches. and represent excitatory input from two features. Right: correlation between spatial profiles of excitation (weighted input rates, summed within branches) from association pairs after learning.

(D) As in (C) but comparing excitation and inhibition. In this case, the excitatory and inhibitory contributions from both features are summed before computing the correlation. In basal dendrites, spatially selective inhibition serves to suppress the response to (−) patterns.

All bars denote means; p values are from two-tailed Wilcoxon signed-rank tests between groups for n = 10 independent replications.

Although the performance of the active and passive models was surprisingly similar, there are several notable differences. Although the degree of functional clustering was comparable in both cases, this was translated into much larger differences in nonlinear summation in the active model (Figure 3D). This effect may be due to the presence of steeper dendritic nonlinearities in the active case (Figure 1), allowing a larger response difference between classes. Moreover, we find that the active model is more robust to input noise, tested by varying the rates of stimulus-dependent and background synaptic activity and the total number of synapses (Figures S3A–S3C). We also performed additional simulations in which the connectivity was highly structured, such that the solutions learned with random synapse placement in Figure 3 were either increasingly “hard-wired” or prohibited. Strikingly, in the latter case, whereas the passive model failed to learn the task, the active model learned to use an alternative strategy of sublinear processing in basal dendrites (Figures S3D and S3E). Active dendritic processing therefore enables more robust and flexible solution of the task. Generally, however, our results demonstrate that the power of multilayer integration is broadly accessible to single neurons through active and passive dendritic mechanisms, even with random connectivity. Spatial processing of rate-coded input can be learned in either case through selective enhancement of co-localized synaptic weights.

Computation with temporal signals

We have shown how a single pyramidal neuron can learn to compute nonlinear functions of input rates. A long-standing question is whether analogous computational principles govern temporal processing within a branched dendritic tree (Poirazi et al., 2003b; London and Häusser, 2005; Jadi et al., 2014). One challenge in addressing this question is choosing a form of temporal input, as there are many possibilities for encoding information in spike timing, and the organization of in vivo synaptic input is unknown. We therefore proceed in two stages. Using the biophysical model as a guiding constraint, we first identify a form of input that maximizes performance on the feature-binding task, while respecting the stochasticity of in vivo spike generation and transmission. We then explore this regime as a hypothesis of temporal coding that is optimally suited to single neuron computation.

Extending the approach described above, we represented input features by randomly generated patterns of time-dependent presynaptic rates. We parameterized a broad space of candidate input schemes by the time-averaged firing rate of active synapses and the number of precisely timed events per active synapse (implemented as Gaussian bumps of elevated activity), while holding the total population rate constant (Figure 5A). This parameterization defines a space of possible input schemes ranging from a sparse rate code (as implemented in Figure 3) to a dense temporal code in which every synapse receives one presynaptic spike on average within a precise window of time (similar to Gütig and Sompolinsky, 2006, but here with stochastic spike generation). Interpolating between these extremes is a mixed regime in which inputs carry both rate and temporal information, communicated through single or multiple bumps of activity of varying size.

Figure 5

Sparse, precisely timed bursts of input maximize classification performance

(A) Left: parameterization of candidate rate and temporal coding schemes by the time-averaged input rate to active synapses and the number of precisely timed elevations of the input rate (implemented as Gaussian bumps). The total presynaptic population rate is constrained to be the same for all parameters. To enforce this constraint, and with scaling to maintain physiological instantaneous input rates, patterns also differ in sparseness (fraction of active synapses) and temporal precision (width of rate elevations). Note that we use a decreasing order for the x axis in the temporal region; having multiple precisely timed events per synapse more closely resembles a rate code than a single event per synapse as the input spikes are more uniformly distributed in time. Right: example rate functions generated for an association pair in the feature-binding task by different parameter choices. (i) sparse rate code, (ii) dense temporal code, (iii) and (iv) mixed regimes comprising bursts of activity. For clarity, only 20 synapses are shown.

(B) Example 7 × 7 matrix of associations to be classified. Classification labels are randomly assigned for each replication.

(C) Performance (fraction correct) of trained active, passive, and point neuron models, averaged over 10 replications for each input condition. The optimal form of input (asterisk) was the same for all models. R and T denote the rate and temporal schemes used for comparison in (D) and (E).

(D) Example realizations of Poisson input to a synapse for the rate, optimal, and temporal conditions. Bursts in the optimal condition are temporally localized but do not suffer from the transmission failures of the temporal condition.

(E) Detailed comparison of performance across the three models from an independent set of simulations.

Bars denote means; p values are from two-tailed Wilcoxon signed-rank tests between groups for n = 10 independent replications.

Using the same learning rule and training procedure as before, we found the introduction of temporal signals resulted in near-perfect performance on the 2 × 2 association task of Figure 3 (not shown). We therefore increased the computational load by extending the feature-binding task to classification of 7 × 7 randomly labeled associations (Figure 5B). We trained each of the three models on the 7 × 7 task using input patterns generated under the various parameterized input regimes. We find that the form of input that maximizes classification performance is the same for all models (Figure 5C). In the optimal input regime, sparsely activated synapses communicate in single ∼50 ms bursts, each comprising an expected eight presynaptic spikes (Figure 5D). Notably, across all models, synaptic input of this form resulted in a pronounced performance enhancement over the more commonly assumed purely rate or temporal schemes, as assessed with an independent set of simulations (Figure 5E; active, 73% ± 5%/90% ± 3%/70% ± 3%; passive, 72% ± 5%/91% ± 2%/71% ± 5%; point neuron, 62% ± 5%/89% ± 4%/72% ± 2%; mean ± SD for rate/optimal/temporal input). This advantage was observed for computational loads tested up to 10 × 10 associations and was robust to background noise and trial-by-trial perturbations to synaptic weights and burst timing (Figure S4).

In Figure 6A, we show an example simulation of the active model in the optimal input condition solving a nonlinear subset of the task, analogous to Figure 3A. To understand the implementation, we examined the subthreshold potentials and synaptic weights of the trained models as before. Simulating the subthreshold response to individual feature components revealed a common temporal strategy. After training, the time courses of somatic responses are shaped such that only preferred pairs of features will sum constructively (Figure 6B). For all models, the somatic responses to features forming preferred associations are therefore temporally correlated, whereas responses to features forming nonpreferred associations are temporally anticorrelated (Figure 6C). Analogous to the rate-coded case (Figure 4), the implementation can be understood at the synaptic level in terms of the profiles of input that have been shaped through plasticity. For each feature component and dendritic domain we computed a temporal profile of excitation and inhibition as a time-dependent sum over synapses, scaled by synaptic weights. In the point neuron model, temporal profiles of excitation are positively correlated for preferred associations and negatively correlated for nonpreferred associations (Figure 6D). The temporal correlation between inhibition and excitation also differs with association type, but with a reversal of sign (Figure 6E). In the active and passive models, in which the influence of inhibition is spatially restricted within the tree, we find that this strategy is implemented within each dendritic domain (Figures 6D and 6E). Simulating the response to basal and apical input separately and comparing somatic membrane potentials reveals that these locally computed signals are also coordinated globally. Subthreshold potentials arising from input to each domain are preferentially aligned on preferred associations (active, 0.38 ± 0.1/0.04 ± 0.09; passive, 0.27 ± 0.12/0.21 ± 0.09; mean ± SD of temporal correlation between somatic membrane potentials from basal and apical input for preferred/nonpreferred associations), meaning that responses can be learned that are shaped locally by dendritic inhibition before appropriate summation at the soma.

Figure 6

Synaptic plasticity can shape subthreshold potentials to implement a temporal feature-binding strategy

(A) Example simulation of a model trained on the 7 × 7 association task in the optimal precisely timed burst input condition. The response on a nonlinear 2 × 2 subset of the task is shown, as defined by the classification labels (left). Each combination of features is presented in turn, interspersed with background noise. The markers in the raster are scaled in proportion to synaptic weight. For clarity, only input to excitatory synapses is shown.

(B) Example of average subthreshold somatic membrane potentials arising from presentation of input features in isolation, with somatic spiking blocked. Shaded area is SD from 20 presentations. The components and correspond to those simulated in (A). Note that and will sum constructively to produce a spike at ∼200 ms, but and will not.

(C) Across all models, after training the subthreshold potentials arising from pairs forming preferred associations are temporally correlated, whereas those arising from pairs forming nonpreferred associations are anticorrelated.

(D) Left: analogous to the spatial clustering strategy of Figure 4, over learning, synaptic weights evolve to temporally align patterns of excitation to bind preferred associations. and represent excitatory input from two features. Right: correlations between temporal profiles of excitation (weighted input rates, summed over synapses) from pairs forming preferred or nonpreferred associations.

(E) As in (D) but for the alignment of excitation and inhibition. and represent inhibitory input from two features. Excitation-inhibition correlations are calculated after summing the excitatory and inhibitory contributions of each feature component in a pair. Weighted excitatory and inhibitory input is aligned on nonpreferred associations, serving to suppress somatic output.

Bars denote means (averaged over 20 presentations of each association, then over label types); p values are from two-tailed Wilcoxon signed-rank tests between groups for n = 10 independent replications.

In summary, these results show that from a random basis of synaptic input bursts of suitable density and precision, reliable analog signals can be constructed and selectively combined for computation. In contrast to the rate-coded scheme, this optimal representation allows binding of input features through integration in either single or multiple compartments. However, as performance was consistently high across all models in our simulations, further contributions of dendritic processing were unresolved. We now show that as this input regime contains information in both rates and spike timing, spatial processing can also be synergistically recruited as the capacity of the temporal processing strategy is reached.

Synergistic recruitment of spatial and temporal processing

To understand what can be learned by a single neuron beyond computations involving purely spatial or temporal processing, we challenged the models with the 7 × 7 task, though under progressively shorter durations of stimulus presentation to gradually saturate the temporal capacity (Figure 7A). The performance of all models decreased with decreasing stimulus duration, reaching that of an equivalent rate code as the duration was reduced below the width of a single input burst (Figure 7B). The performance of the point neuron, however, which represents the limit of a purely temporal strategy, fell off the most rapidly. At intermediate durations, the active and passive models performed well above the lower bound represented by the equivalent rate code, even as the point neuron dipped below (active, 80% ± 3%; passive, 77% ± 4%; point neuron, 67% ± 4%; mean ± SD for 75 ms duration). This suggests a concurrence of spatial and temporal processing, enhancing performance beyond that of either strategy alone.

Figure 7

Synergistic recruitment of spatial and temporal processing

(A) Schematic of presynaptic rates underlying the precisely timed burst input scheme with compressed stimulus presentation time.

(B) Average model performance after training as a function of stimulus duration for ten independent replications per condition. Dashed lines are average performance under a rate code of the same sparseness and time-averaged maximum rate (20 Hz), presented for 400 ms. Note that the x axis is not a linear scale. Shaded areas are SEM.

(C) Analysis of the relative contribution of spatial and temporal processing as a function of stimulus duration. The signature of each strategy is imprinted on the synaptic weights through learning, allowing the classification label of a given association pair to be predicted from knowledge of the weights and input rates. Plots show the prediction accuracy of logistic regression models fitted to predict the classification labels on the basis of correlations between spatial (, squares), temporal (, circles), and spatiotemporal (, diamonds) input profiles after training for the active (red) and passive (blue) models. Spatiotemporal correlations in the active model are more predictive of class labels than spatial correlations alone, implying a local organization of temporal signals within individual branches.

(D) Schematic of spatiotemporal feature-binding strategies. Traces represent the excitation of dendritic branches with weighted input from two stimulus features. With supralinear integration, the response to preferred associations (+) can be synergistically enhanced by tuning weights such that excitation is both clustered and synchronous (denoted by red arrow). Input from nonpreferred associations (−) should instead be dispersed and asynchronous. Feature binding with sublinear integration demands dispersed, synchronous input from preferred associations and clustered, asynchronous input to suppress responses to nonpreferred associations. Although less compatible than supralinear processing, local sublinear integration could compensate on branches where temporal segregation is incomplete (blue arrow). Inhibitory input can also act in both cases to sharpen temporal responses and aid suppression of (−) pairs (not shown).

See also Figure S6.

The analysis of trained models in Figures 4C and 4D and Figures 6D and 6E shows that distinct signatures of the spatial and temporal feature-binding strategies are imprinted on the synaptic weights through learning. We repeated this analysis to determine the contributions of spatial and temporal processing across stimulus durations. As above, we computed weighted spatial and temporal input correlations for all association pairs, separately for each dendritic domain, and for excitation and inhibition. We use the shorthand (spatial correlations) and (temporal correlations) to refer to these groups of measurements. We find that differences in between pairs of preferred and nonpreferred associations decrease with stimulus duration for all models (Figures S6A and S6B). In the active and passive models, this apparent reduction in temporal processing power is accompanied by a graded increase in spatial organization, revealed by differences in , such that both signatures are present where the advantage over the point neuron is greatest (Figure S6C). We quantified this by fitting regression models to predict the classification label of each association pair from and (Figure 7C). This analysis reveals two contributions to the superior performance of the active and passive models. Sensitivity to both spatial and temporal structure gives the opportunity to use one strategy when the other would fail because of the statistics of a particular input pattern. Spatial and temporal processing can also be jointly recruited to act within a single pattern presentation.

Regression prediction accuracy is highest when regressing classification labels on for long durations and for short durations, as expected (Figure 7C). For intermediate durations, both are approximately equally effective. Pooling over the 75 ms and 100 ms conditions, we find a fraction of association labels can be correctly predicted from , but not from (active, 25%; passive, 22% of labels). On this set (as defined by the active model) the average task performance of the point neuron model was 77%. A similar fraction can be correctly predicted from , but not (active, 22%; passive, 23% of labels). On this set the task performance of the point neuron was reduced to 63%, close to the rate-coded lower bound (Figure 7B). These two groups reflect different subsets of associations where either the spatial or temporal strategy is implemented in isolation. However, the labels of a large fraction of associations can be correctly predicted from and also separately from (active, 44%; passive, 45% of labels), reflecting another subset in which both forms of processing are present. We explored this further by predicting labels from learned spatiotemporal input correlations, , computed by concatenating the temporal profiles of input to individual branches. For the active model, is uniformly more predictive of class label than alone (Figure 7C), driven mostly by weights in basal dendrites. This implies that learning has organized temporal signals even within single dendritic branches. Clustered and synchronous excitation can thereby strongly engage supralinear integration to bind preferred associations (Figure 7D). By contrast, prediction accuracy with in the passive model was degraded. Here, the spatial and temporal strategies are in opposition, as preferred associations demand dispersed yet synchronous input, whereas nonpreferred associations demand clustered yet asynchronous input (limiting engagement of dendritic nonlinearities). Although local sublinear processing could serve to compensate on nonpreferred associations where temporal segregation is incomplete (Figure 7D), the passive model lacks the inherent compatibility of the active model for spatiotemporal integration on this task, which may explain the consistent performance gap between the two (Figure 7B).

We conclude that with synaptic plasticity acting on spatiotemporal input to dendrites, complementary processing strategies can be flexibly selected, and even combined, to solve a challenging nonlinear computational task. When spatial and temporal features of synaptic input patterns both carry information about a stimulus, supralinear dendritic integration provides an ideal foundation for exploiting their interaction.

Discussion

Nonlinear dendritic integration enables single neurons to perform sophisticated computations on their inputs. For neural circuits to capitalize on this processing capacity beyond what can be hard-wired during development requires dendritic computations to be learned. We have shown how this can be efficiently accomplished through synaptic plasticity. Through an extension of dendritic cable theory, we quantified for the first time how the spiking output of a neuron is influenced by small changes in the strengths of its many interacting synaptic inputs. We used this analysis to develop a supervised spike-timing and voltage-dependent learning rule that optimally adjusts synaptic weights to control the input-output function of a detailed pyramidal neuron model. Training the model to perform a nonlinear classification task, we demonstrated that the biophysical properties of dendrites can be harnessed without structured connectivity, which greatly expands the repertoire of dendritic computations that can be exploited by the brain. Computations that use high-order spatial or temporal features of random synaptic input patterns can be learned through the same form of plasticity, yielding synaptic weight distributions that engage supralinear or sublinear dendritic integration or precisely control the time course of subthreshold potentials. When information is encoded in both the rates and spike timing of the input, spatial and temporal processing strategies can be combined to maximize the computational power of single neurons.

The foundation of our approach is a precise mathematical description of the influence of synaptic weight changes on somatic output. Despite the inherent complexity of dendritic integration, we found this information is available to individual synapses through a simple functional dependence on synaptic activation time and dendritic voltage. A neuron can therefore optimally update its synaptic weights using local signals, while still accounting for nonlinear interactions throughout the dendritic tree. Encapsulating this insight in a supervised learning rule allowed us to explore principles of single neuron computation under a range of input conditions and model assumptions. That very different computations were learned with rate or temporal input, and even in distinct dendritic domains, demonstrates the sensitivity of our learning rule and its general applicability. The algorithmic perspective brought by our approach complements previous theoretical studies that have invoked spike-timing and voltage-dependent plasticity on the basis of experimental observations (Legenstein and Maass, 2011; Clopath et al., 2010; Bono and Clopath, 2017; Ebner et al., 2019). Although we have not sought to reproduce the results of specific experiments, our rule is biologically grounded insofar as it operates within the biophysical constraints imposed by our detailed model and weight updates are computed from variables that are plausibly accessible to dendritic synapses. We acknowledge, however, that several biological details of the implementation have been left unspecified. In particular, we have modeled the supervisory signal at an abstract level; to make the underlying assumptions explicit, a more detailed treatment in future work could include a specific circuit for neuromodulatory control (Gütig, 2014). Similarly, the implementation of the local plasticity kernels could be developed by modeling signaling cascades at the synapse that process pre- and postsynaptic spike times and the local dendritic voltage (Ebner et al., 2019). To further elucidate connections between the algorithmic and biological aspects of synaptic plasticity, another important extension of our analysis will be to characterize how the influence of synaptic weight changes depends on the wider range of dendritic mechanisms found throughout the brain. We hypothesize that, among other factors, experimentally observed differences in plasticity induction within and across cell types (Sjöström et al., 2008) may reflect specializations to allow optimal assignment of synaptic credit from local dendritic signals.

Training the model on a rate-coded feature-binding task, we showed that a single neuron can learn to implement nonlinear functions of randomly distributed input by engaging dendritic nonlinearities. Therefore, not only do pyramidal neurons share structural similarities with multilayer artificial networks (Poirazi et al., 2003b; Jadi et al., 2014), but they can capitalize on this foundation without requiring special connectivity rules or structural plasticity. Modulo the remaining mechanistic questions detailed above, this flexibility would reduce the developmental burden of precise axonal targeting and allow the output of a neuron to easily adapt to changes in input tuning or network demands. The learned implementation differs from the branch-strength potentiation modeled by Legenstein and Maass (2011), as this mechanism, observed in hippocampal neurons, is not a feature of our model. Instead, our model predicts a form of in-branch functional clustering (Kastellakis et al., 2015), similar to that imposed in the abstract Boolean model of Cazé et al. (2013), in which it was shown that distinct wiring schemes could take advantage of either supralinear or sublinear dendritic integration. Our results demonstrate that both possibilities can in fact emerge naturally from identical input connectivity and the same learning rule. As a consequence, we find opposite implementations in apical and basal dendrites, whereas the rules of multistage integration in these domains have previously been assumed to be the same (Jadi et al., 2014). Why sublinear processing in apical dendrites is optimal for the task is unclear; it may be that saturating voltages are required to overcome attenuation from the tuft (Moldwin and Segev, 2020). This intriguing result reinforces the importance of analyzing models built with realistic morphological constraints. More broadly, these observations and the comparable performance of the purely passive model underscore the potential of diverse forms of dendritic nonlinearity for computation (Poirazi and Mel, 2001; Ujfalussy et al., 2015). The power of hierarchical dendritic processing may extend far beyond the excitable dendrites of pyramidal cells to include dendrites exhibiting largely passive integration (Abrahamsson et al., 2012; Vervaeke et al., 2012), or where a mixture of both supralinear and sublinear integration is found within inhibitory interneurons (Tzilivaki et al., 2019).

Solving the same task under a range of temporal input schemes revealed an optimal regime of input comprising sparsely distributed bursts. Although it was feasible in this study to survey only a fraction of possible temporal coding schemes, our results highlight several general principles. First, as argued previously (Lisman, 1997; Krahe and Gabbiani, 2004), burst firing of presynaptic cells provides reliable units of information to postsynaptic cells. Earlier models trained on temporal pattern discrimination showed impressive performance when presented with precisely timed and reliable single spikes (Gütig and Sompolinsky, 2006, 2009; Schiess et al., 2016). We found that using a similar dense temporal scheme, although accounting for the stochasticity of spike generation, failed to improve upon the performance of a much simpler rate code (Figure 5E). However, performance improved dramatically when synapses were activated in bursts, ensuring robust signaling while retaining spike-timing information. Short-term synaptic dynamics, not considered here, could further accentuate this advantage (Lisman, 1997) and even allow multiple streams of information to be multiplexed within single cells (Körding and König, 2001; Naud and Sprekeler, 2018; Payeur et al., 2021). Second, as with connectivity, plasticity acting on globally unstructured temporal signals provides a general and flexible basis for computation. Well-studied structured alternatives include highly synchronous or sequential activation of synapses, which, through intrinsic biophysical mechanisms, can elicit differential responses to specific sets of inputs (Rall, 1964; Ariav et al., 2003; Branco et al., 2010; Branco and Häusser, 2011; Bhalla, 2017). We found, however, that with patterns defined by randomly timed events, our model could learn arbitrary responses by shaping the time course of subthreshold potentials. In this case, dendritic processing tuned by plasticity can decide the stimuli that are most relevant to a cell. Finally, our results reveal the remarkable computational potential that exists beyond a strict rate or temporal coding dichotomy. Dendrites confer exquisite sensitivity to both representations, and through spatiotemporal integration, processing of rate and temporal information can be synergistically combined.

Key predictions of our model can be tested experimentally without the need to measure synaptic weights or input rates. We propose experiments that, although challenging, could be performed in layer 2/3 pyramidal cells with currently available techniques. Performing chronic in vivo Ca2+ imaging of pyramidal cell populations while training animals to discriminate stimulus conjunctions as in Figure 3, we predict the emergence of individual neurons with nonlinear selectivity. This would confirm that such computations can be learned, even in early sensory areas, rather than being hard-wired through connectivity. In neurons exhibiting this selectivity, the signature of the optimal spatial strategy could then be identified through in vivo Ca2+ imaging of dendritic branches (Kerlin et al., 2019; Palmer et al., 2014). As schematized in Figure 7D, in basal dendrites, preferred associations are predicted to elicit a strong local Ca2+ signal from activating NMDA receptors in a small number of branches, whereas nonpreferred associations are predicted to elicit comparatively weaker responses in a larger number of branches (conversely for apical dendrites). The predicted signature of the temporal strategy could be identified through two-photon targeted patch-clamp recordings (Margrie et al., 2003; Kitamura et al., 2008) from the imaged population. In conjunction-selective cells, preferred associations are predicted to evoke a consistent peak in somatic voltage underlying a spiking response, whereas multiple subthreshold peaks are predicted for nonpreferred associations, reflecting the temporal alignment or misalignment of excitation arising from individual stimulus features.

Beyond the present study, the strength of our general approach is that it can be applied to investigate a wide range of possible computations, as well as other biophysical mechanisms and cell types. The governing Equations 1, 2, 3, and 4 are indeed easily extended to any number of active conductances distributed throughout a given dendritic morphology. Thus, although we have focused on voltage-dependent NMDA receptors as the major driver of supralinear synaptic integration (Branco and Häusser, 2011; Major et al., 2013), the computational roles of additional dendritic conductances can also be quantitatively assessed. We have taken a first step toward this goal, demonstrating an increase in classification performance in the presence of currents and dendritic Na+ spikes (Figure S2). However, a remaining challenge in this more excitable regime is finding an accurate local approximation of , as used for learning in the strictly NMDA-dependent active model (Figure 2B). This will likely require a dependence on additional local variables, such as high-pass-filtered voltage signals that can detect fast dendritic Na+ spikes. Further investigation of the impact of Na+ spikes on stimulus tuning (Smith et al., 2013; Goetz et al., 2021) and temporal processing (Ariav et al., 2003) is a priority for future work. An intriguing nonmonotonic dendritic nonlinearity was also recently found in human layer 2/3 neurons (Gidon et al., 2020). In principle, this inverted U-shaped response to input magnitude could allow an XOR to be implemented within a single compartment (Zador et al., 1991; Gidon et al., 2020), whereas implementation with the more common sigmoid nonlinearities of our model would require two (in a similar fashion to feature binding but with stronger recruitment of inhibition). Whether human neuron electrophysiology might confer additional benefits beyond this specific computation remains to be elucidated and would be of great interest to explore through task optimization of a detailed model, as we have done here. Applying our methods to study dendritic integration in neurons of cortical layer 5, the hippocampus, and cerebellum could also provide fresh insights into other active phenomena, such as associative interactions that produce apical tuft Ca2+ plateaus (Larkum et al., 1999; Bittner et al., 2015) and temporal regulation by currents (Magee, 1999; Angelo et al., 2007). Although we found that all of our major results were reproduced in a second layer 2/3 pyramidal cell model (Figure S7), it will be very revealing to explore how local rules of plasticity and computation may differ across cell types with distinct dendritic architecture and physiology. We suspect, for instance, that the segregated apical tuft of larger layer 5 neurons would preclude a direct generalization of the learning rule presented here. In this case, a two-stage rule may be more applicable, in which credit is assigned locally to synapses for producing apical Ca2+ plateaus, and also globally in proportion to the ability of the plateau to drive somatic output. Although many important questions about the implications of dendritic physiology for learning and computation remain open, the tools we have developed will allow these to now be systematically addressed.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Michael Häusser (m.hausser@ucl.ac.uk).

Materials availability

This study did not generate any new reagents.

Method details

Biophysical model

Morphology and passive properties

We used two adult mouse V1 layer 2/3 pyramidal cell morphologies from the Allen Cell Types Database (Allen Institute for Brain Science, 2015). The morphology shown in Figure 1 (ID 502359001) was used for the majority of simulations. The second morphology (ID 521409057) was used to confirm the reproducibility of the results (Figure S7). Possible reconstruction errors in the form of pinched sections of dendrite were manually corrected by linear interpolation between adjacent segments, and the diameter profile of each branch was then smoothed with a moving average filter. Dendritic sections were subsampled such that each branch comprised at least two compartments, and the maximum compartment length was less than 10 μm (462 and 510 total dendritic compartments for the first and second morphology). For computational efficiency and to constrain the number of free parameters, we removed the detailed reconstructed axon and used a lumped axo-somatic compartment of radius 10 μm for action potential generation. The specific membrane capacitance , specific membrane resistance , axial resistivity , and leak conductance reversal potential were set as in the layer 2/3 pyramidal neuron model of Branco et al. (2010), giving a membrane time constant . We did not include dendritic spine compartments in the model due to the computational cost and uncertainty surrounding parameters such as spine-neck resistance. We confirmed that our main results did not differ qualitatively when modeling the influence of spines on membrane surface area using the common approach of scaling the dendritic capacitance and leak conductance by a factor of two (not shown) (Holmes and Rall, 1992; Hay et al., 2011).

Active conductances

For a biophysical spiking mechanism, we included fast transient channels (reversal potential , maximum conductance ) and fast delayed-rectifier channels (reversal potential , maximum conductance ) in the axo-somatic compartment. Slow persistent channels (reversal potential , maximum conductance , adaptation time-constant ) were included for spike-rate adaptation (Pospischil et al., 2008). Channel kinetics were implemented with standard models using the Hodgkin-Huxley formalism, as described by Pospischil et al. (2008) (see also (Equation 8), (Equation 9), (Equation 10), (Equation 11), (Equation 12)).

In the set of simulations exploring the role of intrinsic dendritic excitability (Figure S2), the and conductances were extended throughout the dendritic tree (dendritic maximum conductances: , , ). A uniform density of HCN channels was also included in both the dendritic and axo-somatic compartments (reversal potential , maximum conductance ) using the model of Kole et al. (2006). Dendritic conductance parameters were tuned such that the model produced fast dendritic spikes in response to synaptic input (Smith et al., 2013), and a voltage sag of in response to hyperpolarizing current injection, consistent with deeper mouse layer 2/3 pyramidal neurons (Kalmbach et al., 2018).

Synaptic conductances

Excitatory synapses with AMPA and NMDA conductances, and inhibitory synapses with GABA conductances were modeled with double-exponential activation kinetics. That is, for a presynaptic spike arriving at time , the time course of the activation for is given bywhere and are rise and decay time constants specific to each conductance type, and is a normalization factor that ensures peaks at a maximum of 1. The time-dependent activation function is scaled by a synaptic weight to determine the conductance. Weights were specific to each synapse and variable within the simulations. For excitatory synapses, the weight describes a combined NMDA and AMPA conductance, with a fixed NMDA/AMPA ratio . Rise and decay time constants for each of the conductance types were as follows. AMPA: , ; NMDA: , ; GABA: , , similar to previous work (Doron et al., 2017).

In the active and point neuron models, the NMDA voltage dependence was modeled by multiplying the time-dependent conductance by a local-voltage-dependent sigmoidas described by Jahr and Stevens (1990). We used parameters and to define the shape of the nonlinearity. While some recent studies have used steeper nonlinearities (Poleg-Polsky, 2015; Doron et al., 2017) or more complex multi-state receptor models (Branco et al., 2010), in our morphologies, we found that Equation 7 with the stated parameters was sufficient to produce the characteristic sigmoid response to synaptic input (Figure 1). The passive model was constructed by setting and leaving all other biophysical parameters unchanged. The reversal potentials for excitatory and inhibitory synapses were and , respectively.

Dynamics and variational equations

To simplify the notation, we provide the model equations for the general case in which and conductances are extended throughout the entire dendritic tree. The model simulated in the main text in which these conductances are restricted to the axo-somatic compartment (summarized by Equations 1, 2, 3, and 4), is recovered by setting the dendritic conductance parameters (, , for ) to zero. For a given pattern of presynaptic input, the dynamics of the voltage in each compartment evolve by the coupled system of equationsThe first two sums on the right hand side of Equation 8 describe the excitatory and inhibitory synaptic currents. The matrix with elements projects excitatory conductances with weights into the compartments in which the synapses are located, indexed by i, and normalizes the resulting currents by compartment surface area to give the appropriate change in voltage. Similarly, the inhibitory currents are described by weights and projection . The functions , , and denote the time-dependent activation of AMPA, NMDA and GABA conductances by input spike trains, constructed by linear summation of the responses to individual spikes described by Equation 6. The NMDA/AMPA ratio γ controls the relative contributions of excitatory conductance types. The subsequent three terms describe the active and conductances, expressed here in generality for both the somatic and dendritic compartments. The final two terms describe leak and axial conductances, with the matrix encoding the branching structure and geometry of the dendritic tree (see Dayan and Abbott (2001); Chapter 6). Finally, (Equation 9), (Equation 10), (Equation 11), (Equation 12) describe the kinetics of the Hodgkin-Huxley gating variables for the active conductances. Variables and control the transient conductance in compartment i, controls the fast conductance, and controls the slow conductance. The voltage-dependent rates for each case, and for , are as described in Pospischil et al. (2008), section 2.2.

Variational equations for the system are derived by taking partial derivatives of (Equation 8), (Equation 9), (Equation 10), (Equation 11), (Equation 12) with respect to any parameters of interest (similar to the point-neuron parameter-fitting technique of Doya et al. (1994)). When solved numerically in parallel with the model dynamics described above, these equations track how small perturbations to parameters are propagated through time and the morphology to influence the voltage throughout the tree. Here, we consider the excitatory and inhibitory weight parameters, with the aim of using the gradient of somatic voltage with respect to synaptic weights to guide learning.

We denote by the partial derivative of the voltage in compartment i by synaptic weight (similarly for the gating variables). The component with index corresponds to the axo-somatic compartment, which we often express using the notation , suppressing the index for the weight when the specific identity of a synapse is not important.

For the excitatory weights,where , and denote derivatives with respect to voltage. The expression for inhibitory weights is similar.

In a subset of the simulations presented in Figure S2 we also include voltage-dependent HCN channels in the axo-somatic and dendritic compartments. For this we used the Hodgkin-Huxley type channel model of Kole et al. (2006), which describes the channel kinetics in terms of a single gating variable that we denote by q. The HCN conductance can be incorporated into the model in the same manner as the other active conductances, by adding the additional term to Equation 8 and an additional first order equation for the gating variable analogous to (Equation 9), (Equation 10), (Equation 11), (Equation 12). The associated terms that appear in the variational equations are computed by taking partial derivatives. Note that for channel models constructed using the Hodgkin-Huxley formalism these operations are formulaic, allowing our general approach to be extended to additional conductance mechanisms.

We solve the system numerically using a custom implicit Euler solver written in Python 3 (Van Rossum and Drake, 2009), using Numpy (Harris et al., 2020) and Numba (Lam et al., 2015). Theoretical background for the implementation was drawn from Hines (1984); Dayan and Abbott (2001); Hines and Carnevale (2006), and the accuracy of solutions to (Equation 8), (Equation 9), (Equation 10), (Equation 11), (Equation 12) was validated during development using NEURON 7.5 (Hines and Carnevale, 2006). After using this custom simulator to derive the local plasticity kernels as described below, simulations of the model were performed using NEURON 7.5 for computational efficiency.

Synaptic input patterns

To explore the space of possible input regimes, we parameterized a model for generating random input patterns for the feature-binding task. A total number of synapses were randomly divided into two equally sized subsets, representing presynaptic populations associated with two feature classes (described as feature X and feature Y). Input patterns for specific features within a class were defined as vectors of time dependent rates for the synapses assigned to the class. Input patterns were realized in simulations as non-homogeneous Poisson spike trains, generated by rejection sampling. Unless stated otherwise, each input pattern consisted of 100 ms of background activity at constant rate for all synapses, followed by a stimulus presentation between times and . To constrain the parameter space, the time- and population-averaged firing rate of the stimulus component was fixed at per synapse to give an average of one presynaptic spike per synapse over the 400 ms stimulus presentation.

Two parameters define the characteristics of a pattern: the time-averaged firing rate of active synapses , and the number of precisely timed events K. To satisfy the imposed constraint that ensemble-averaged activity is maintained at a constant level, also determines the sparseness of activity – few synapses active with high firing rates, or many active with low firing rates. Synapses were thus assigned time-averaged rates with probability , or otherwise. Precisely timed events were implemented as Gaussian bumps of elevated firing rate centered at uniformly distributed times. To ensure physiological instantaneous firing rates, Gaussian widths were scaled with the number of presynaptic spikes expected to occur during an event (which depends on ) from an initial value of . Altogether, the stimulus-dependent component of the rate function for the ith synapse, , was initialized asThe multiplicative factor in the case in Equation 18 scales the time-dependent rates such that the integral over , ignoring small boundary effects, is equivalent to that when . When pairs of features are presented together, for instance and , the expected number of presynaptic spikes over all synapses is thus for all choices of parameters.

Plasticity kernels

We constructed plasticity kernels to be used in the learning algorithm (Figure 2B; Figures S1A–S1D) by approximating the output of (Equation 13), (Equation 14), (Equation 15), (Equation 16), (Equation 17) in terms of the local dendritic voltage and synaptic activation time. Specifically, we approximated the influence of synaptic weight changes on the somatic membrane potential at times immediately preceding somatic spikes, . In this way, when the kernels are used to update synaptic weights at the time of somatic spikes, synapses are selectively modified relative to their ability to enhance or suppress somatic output in response to a given input pattern.

To construct the approximation, for each of the active, passive and point neuron models we ran 5000 10 s simulations using randomly generated Poisson input. In each simulation, input rates were drawn independently for each synapse from a lognormal distribution , yielding a bombardment of synaptic input featuring both isolated single spikes and high frequency spike trains. Synaptic weights were set for excitatory and inhibitory synapses by randomly perturbing the initial values and , giving a range of weights . For the passive model, the excitatory weights were reduced by a factor of five to give a similar postsynaptic firing rate.

Guided by preliminary simulations, we focused our analysis on windows of 150 ms preceding somatic spikes. We excluded windows in which additional output spikes had occurred, to ensure we collected independent measurements (leaving somatic spikes for analysis in each condition). For each somatic spike we computed the influence of synaptic weight changes on the somatic voltage at the upswing of the action potential (defining as 2 ms before the voltage crossed 0 mV). We used this time point as a functional definition of spike threshold crossing; once the autonomous voltage-dependent spiking mechanism is engaged during the upswing, the somatic voltage becomes far less susceptible to synaptic control (i.e., ‘the horse has bolted’). Integrating (Equation 13), (Equation 14), (Equation 15), (Equation 16), (Equation 17) for a detailed model is computationally expensive, so was performed over the selected 150 ms windows, rather than the entire 10 s simulation run. This was accomplished by simulating the model again over the selected window using initial conditions and input set from the same time point in the original run. We also made two small modifications in these sub-simulations, which were found during development to enable more reliable predictions about the influence of weight changes on spiking output. First, we made dummy copies of any synapses that were activated more than once, such that every synapse received a single presynaptic spike, while the input to the neuron remained functionally identical. This allowed us to compute and then approximate for each individual input, rather than the value aggregated over a sequence of randomly timed inputs. Second, we blocked the action potential at the end of the sub-simulation window by setting the fast axo-somatic and conductances to zero. This avoids occasional spuriously high values arising from the steep nonlinearity of the action potential itself, but has negligible influence on the preceding sub-threshold integration that is the underlying cause. Along with the computed values of for each somatic spike, we recorded the associated presynaptic spike times, synaptic weights and dendritic voltages.

For the active and passive models, we pooled the simulated data across synapse locations in basal and apical dendrites. We computed averages of as a function of the difference in timing between pre- and post-synaptic spikes, and the local dendritic voltage at the time of the post-synaptic spike. In Figure 2D we show that the local voltage at the time of the post-synaptic spike encodes the level of activity of neighboring excitatory and inhibitory synapses over the past 100 ms. This is what allows our rule to selectively enhance or suppress coactive groups of synapses within a branch. Note that if plasticity were instead to depend on the voltage at the time of a presynaptic spike, much of this crucial dependence would be lost because the rule would be blind to any subsequent activity in the branch. Averages were computed after discretizing time and voltage in 1 ms and 1 mV bins. The plasticity kernels were then constructed by fitting the averaged two-dimensional functions with polynomials of degree 8 in both variables. Because the simulated data did not cover all points of the parameter space, the boundaries at , and were first padded with zeros to allow a smooth interpolation over the whole domain. We also restricted the temporal dimension to a limit of 100 ms time difference between pre- and post-synaptic spikes, beyond which the values of the spike-triggered averages are effectively zero. After fitting, the kernels were thresholded to remove negative values from the excitatory kernels and positive values from the inhibitory kernels, which was found to improve convergence in pilot simulations. For the point neuron model, and simulations using temporal kernels (Figure 3C), we fitted the temporal spike-triggered average (Figure S1A–S1C) with a degree 10 polynomial after padding the boundary at . The accuracy of the fitted kernel approximations was assessed by fitting on of the data and testing on the remaining (Figure 2C), and was used to inform the choices of polynomial degrees used in the fits.

Learning rule

Our aim was to develop a supervised spike-based plasticity rule that could take advantage of the complex morphology and biophysics of pyramidal cell dendritic trees. Importantly, we sought a formulation that could be applied under diverse forms of noisy synaptic input and that did not presuppose the implementation of a given computation (for instance, by including an explicit mechanism to encourage synaptic clustering). Inspired by the tempotron (Gütig and Sompolinsky, 2006), we used a greedy algorithm that updates synaptic weights in proportion to their influence on the somatic membrane potential, . Although the tempotron rule is based on a single-compartment leaky integrate-and-fire neuron model, we found that the underlying principles provide a powerful heuristic for implementing spike-based learning in a detailed biophysical model. We first outline the computational task and tempotron rule, and then detail the modifications that lead to the rule represented by Equation 5 of the main text.

We consider a general binary classification problem, in which a neuron must produce at least one postsynaptic spike in response to preferred patterns of input, and remain silent for nonpreferred patterns. For the specific feature-binding task that we study, input patterns are defined by different pairs of stimulus associations. The preferred and nonpreferred classes of input patterns are denoted by the symbols (+) and (−). Upon presentation of an input pattern p, realized as stochastic presynaptic spiking activity, a binary variable records the output of the neuron, taking the value if at least one spike is emitted, and if the neuron remains silent. The target output is given by a label , taking the value for (+) patterns and for (−) patterns. We define a signed classification error aswhich takes the value for correct trials, for (+) pattern errors, and for (−) pattern errors. The goal of learning is to minimize the expected absolute value of for all input patterns. As a proxy for the expected value, we use an average over a sample of trials, defining the total classification error for the whole set of patterns bywhere bars denote an average over pattern presentations. Unless stated otherwise, we set .

Tempotron learning

For binary classification with tempotron learning, input patterns are presented sequentially over many trials and synaptic weights are modified when the neuron makes an error. Plasticity during presentation of a pattern p is guided by gradient descent on a single-pattern loss function that measures how far the somatic voltage deviates from spiking threshold (assumed to be fixed). In our notation, the loss function is given byIn Equation 23, denotes the spiking threshold and denotes the voltage at the time of the first post-synaptic spike, or the maximum subthreshold voltage attained if no spikes were emitted. By definition of , and assuming a spike is fired when crosses , on correct trials, and for both (+) and (−) pattern error trials.

After presentation of an input pattern, synaptic weights are updated in proportion to the negative gradient of the loss function with respect to the weights,We write Equation 24 as an approximation because it ignores the possibly discontinuous dependence of the value of itself on the weights (see Urbanczik and Senn (2009) for a detailed discussion). The first factor on the right-hand side of Equation 24 acts as a corrective supervisory signal – plasticity is gated on only on error trials, with the sign determined by the class of pattern (positive for (+) pattern errors and negative for (−) pattern errors). The second factor implements synaptic credit assignment by scaling the magnitude of weight update for each synapse by its contribution to the voltage at time . Intuitively, iterating over repeated presentations of input patterns, locally optimal weight updates are applied to push peaks in the somatic voltage below threshold on (−) patterns, and to push peaks above threshold on (+) patterns.

Adaptive learning with dendritic synapses

The tempotron rule can be directly generalized to our biophysical model by using values of computed from Equation 13. However, with dendritic synapses, this rule becomes implausibly nonlocal. In (+) pattern error trials, where no spike was emitted, must somehow be computed and transmitted to a synapse from the subthreshold activity at the soma. Furthermore, with an extended morphology, the value of for a given synapse formally depends on activity in even distant dendrites. To address these issues and better satisfy the requirements of our noisy task, we make three modifications to the tempotron rule. First, we introduce an explicit mechanism that allows the timing of somatic events to be communicated to dendritic synapses. Second, we use an approximation that allows accurate synaptic credit assignment using local dendritic signals. Third, in pursuit of minimizing (Equation 22), we replace the gating factor with modulation by the running average error . The first and third modifications are achieved using an adaptive supervisory signal that guides the learning process. Subsequently, we show that these dynamics permit a straightforward extension to an online implementation, where learning is controlled by feedback from the spiking output of the neuron.

In our approach, we assume that an external supervisory system maintains a running average of past errors, . The supervisory system guides learning by encouraging spiking when it is desired, and modulating the sign and magnitude of plasticity. On presentation of (+) patterns during training, a depolarizing teaching current is activated at the somaThe teaching current could be provided through regulation of synaptic input or active conductances. In simulations, for simplicity, we use a direct somatic current injection. We set the parameter , yielding a maximum current that is too small to elicit spiking on its own, but sufficient in the presence of synaptic input. When the neuron has been performing poorly, the teaching current adds a constant positive offset to the somatic potential, , where denotes the somatic input resistance measured in the absence of synaptic activity. This has the effect of raising otherwise-subthreshold peaks in the somatic voltage above the spiking threshold. The time of a peak, now marked by a backpropagating action potential (Stuart and Sakmann, 1994; Stuart et al., 1997), can thus be broadcast to the dendrites, and does not require secondary computation or transmission after the pattern has been presented. As performance improves, the teaching current gradually decreases until positive classifications are made without assistance.

During presentation of either (+) of (−) patterns, synaptic weights are modified whenever a somatic action potential is fired. The specific contribution of each synapse to the somatic voltage is computed by approximating the output of Equation 13 with a local plasticity kernel (Figure 2B), summed over presynaptic spikes,In Equation 26, the index k runs over presynaptic spikes that arrived at synapse i in the past 100 ms, is the appropriate plasticity kernel for the synapse (basal, apical; excitatory, inhibitory), is the time of the kth synaptic input relative to the somatic spike, and is the local dendritic voltage at the time of the somatic spike. Intuitively, to implement Equation 26, each synapse requires a plasticity kernel that is common among all synapses in basal and apical domains, to know that a postsynaptic spike has occurred, to know the local voltage, and to know how far in the past presynaptic spikes arrived.

Analogous to Equation 24, weight updates are then determined by the product of the locally approximated gradient of the voltage and a global scalar error term, scaled by the learning rate α (Equation 5, in the main text),Using the average in Equation 27 introduces adaptive dynamics to the learning process. When past performance has been consistently poor, updates are large and favor exploration of the weight space. As performance improves, the magnitudes of updates decrease, which favors refinement of the solution. Learning ceases upon minimization of . For comparison, we also tested an alternative rule in which plasticity was instead gated only by the error on the current pattern (for this purpose, positive classifications were considered as errors if assisted by the teaching current). We found that while both rate and temporally coded tasks could still be learned, performance was inconsistent and highly sensitive to the learning rate (Figure S5).

Training and testing procedure

Unless otherwise stated, all models were trained using the same procedure and learning rule. input patterns, divided into (+) and (−) classes, were presented in random order in each of 1000 total epochs of training (500 epochs for the parameter sweep in Figure 5C). For computational efficiency, the simulations were interrupted at the time of the first somatic spike during a pattern presentation, weights were updated using Equation 27, and the next pattern was presented. Training ended when the maximum number of epochs was reached, or after ten successive epochs without errors . With synaptic weights defined in the simulations in units of , we used an initial learning rate , which decayed as a function of training epoch as . The learning rate parameters were selected based on pilot simulations using the association task. In case of runaway growth of the weights, we enforced a maximum weight by clipping at , though in practice this maximum was never reached. Model training required of CPU time on a single core for learning the task, and for the task.

After training, all patterns were presented 20 times with both plasticity and the depolarizing teaching current turned off. Testing performance was quantified from these simulations as , using .

Online implementation

The history dependence of our learning rule provides a basis for an online implementation, in which errors are integrated over time within a single pattern presentation. In this case we assume that the supervisor provides only a binary classification label when a pattern is presented. All other quantities are computed online from spiking output and fed back to the neuron to guide plasticity.

The input patterns to be classified were defined as above by time series of presynaptic rates. During training, to allow for temporal integration and feedback, patterns were presented for a duration , ranging from . For the precisely timed burst input condition, this was achieved by periodic extension of the input patterns over the training presentation time (Figure S5C). All patterns were superimposed on a constant background rate of to every synapse, which remained active throughout the whole presentation period (equivalent to a background noise level of 0.5, in Figure S3B).

We augmented the model with two first-order equations that compute the classification error for the current input pattern as a function of recent spiking activity and the label ,Equation 28 computes an exponentially weighted average of the spike count over an averaging time . The first term denotes an impulse applied whenever the somatic voltage cross 0 mV from below. In practice, when the system is discretized for simulations, the variable r is instantaneously incremented by in any time step in which a spike occurs. Equation 29 uses the sampled rate to compute an exponentially weighted average of the classification error over averaging time . In the first term of Equation 29, denotes the Heaviside step function ( for and for ), which thresholds the rate r, analogous to the variable in Equation 21. To enforce the spikes/silence binary classification we use two separate thresholds. For (+) patterns, when , the rate is compared to an upper threshold . For (−) patterns, when , the rate is compared to a lower threshold . The time-averaged error defined by Equation 29 converges to zero when the rate is stably maintained above on (+) patterns, and below on (−) patterns. In the simulations presented in Figure S5, we used parameters , , , and .

For each presentation of an input pattern, variables r and are initialized at 0, and numerically integrated in parallel with the model voltage dynamics. The variable – now recording a temporal rather than trial average – determines the magnitude of the teaching current activated during (+) patterns as in Equation 25. For the first ms, the model is simulated without plasticity to allow the system to come to a steady state. At time , plasticity is turned on and weight updates are applied using Equation 27 (substituting the time-averaged ) at the time of every somatic spike. The time at which plasticity is turned on is largely arbitrary; generally, we found that learning was more stable when adding a delay that avoids the transient dynamics that arise from initializing Equations 28 and 29 at 0, and the model at rest, for each new pattern presentation. In Figure S5D we show examples of the online learning dynamics during single pattern presentations.

Models were trained using this procedure for 500 epochs. After training, using identical conditions to the offline learning simulations, models were tested on pattern presentations of duration of , comprising 100 ms of background input, followed by 400 ms of background and stimulus-dependent input.

Spatial and temporal processing analysis

Subthreshold integration

After training, models were simulated with the fast axo-somatic Na+ and K+ conductance parameters set to zero to block somatic spiking ( and in Equation 8).

To quantify the supralinearity and sublinearity or responses in Figures 3D, S2C, and S7G, input patterns representing individual feature components were presented individually and then together in association pairs. This was performed separately for basal and apical dendrites by setting synaptic weights to zero in one or the other domain. Somatic voltage traces were recorded and averaged over 20 replications. The nonlinearity of summation for each association was defined as the peak of the average membrane potential when both features were presented together divided by the peak of the sum of the traces from features presented separately. Peaks were computed from the second half of the 400 ms stimulus presentation time to capture the steady-state response.

To quantify the temporal alignment of responses in Figures 6C and S7H, somatic voltage traces were recorded from presentation of individual feature components and averaged over 20 replications (including both basal and apical input together for the active and passive models). The Pearson correlation was computed between average traces from each feature in an association pair, using the entire stimulus presentation time. To compare the output of the basal and apical domains for each association, pairs of features were presented together with input restricted to one or the other domain, and the Pearson correlation was computed between the average basal-input and apical-input traces.

Weighted input profiles

To understand the computational strategies learned by the models, we examined how the synaptic weights determine the spatial and temporal distribution of input strength. We defined a spatial profile of input by a vector with each element computed as the sum of time-averaged excitatory or inhibitory input rates to a single dendritic branch, scaled by synaptic weights . Here, and denote the weights and rates associated with branch b. Similarly, replacing branches with 1 ms time steps, we defined a temporal profile of input by a time-dependent sum of rates over all synapses in a given domain (basal, apical or soma), . Spatiotemporal profiles were constructed by computing a temporal profile of input into each branch, , and then concatenating the branch time series. We use these definitions to quantify the manner in which learning has tuned the interactions between pairs of features forming associations. In Figure 4C, we quantify the spatial clustering of weighted excitatory input patterns by computing the Pearson correlation between their respective spatial input profiles. In Figure 4D, we quantify the spatial clustering of weighted excitatory and inhibitory input patterns. In this case, we first sum the excitatory profiles of each pair, separately sum the inhibitory profiles, and then compute the Pearson correlation between the total excitatory and inhibitory profiles. Similarly, in Figures 6D and 6E, we quantify the temporal alignment of weighted input patterns by computing the Pearson correlation between their respective temporal input profiles.

Regression of classification labels

In Figure 7 we simulate a regime of input in which both spatial and temporal processing strategies can be learned. To unpack the relative contributions, we use regression models to predict the classification label of a given association pair from the profiles of input that have been shaped by plasticity. Distinct signatures of spatial and temporal processing are shown in Figures 4 and 6. For the spatial strategy, preferred patterns are associated with positive correlations between excitatory spatial profiles (spatial clustering) and nonpreferred patterns with negative correlations (spatial dispersion). Analogously, for the temporal strategy, preferred patterns are associated with positive correlations between excitatory temporal profiles (temporal alignment) and nonpreferred patterns with negative correlations (temporal misalignment). The reverse contingencies tend to hold when comparing excitation and inhibition. In Figure 7, pooling across all simulated data for each stimulus duration, we fit separate logistic regression models to predict association label, (+) or (−), from correlations between the spatial, temporal and full spatiotemporal profiles of association pairs (using LogisticRegression from the Scikit-learn library in Python). Prediction accuracy was determined using leave-one-out cross validation. Regressors included correlations between profiles of excitation, and between profiles of excitation and inhibition, computed separately for basal and apical domains. In the temporal model, we also included correlations between global temporal profiles, computed as a sum over all synapses.

Noise robustness

Models were trained and tested under a variety of conditions to assess the robustness to multiple sources of noise: trial-by-trial spike count variability (Figure S3A), background synaptic activity (Figures S3B and S4B), mislabelling during training (Figure S3C), perturbations to synaptic weights (Figure S4C), and trial-by-trial burst timing variability (Figure S4D).

Robustness to spike count variability was tested by varying the rates of presynaptic input between via parameter in Equation 19. Although our model for generating input patterns keeps the total population rate fixed for different values (by adapting the sparseness of activation), smaller increases the standard deviation of the number of spikes arriving at each synapse, relative to the mean. To test the dependence on spike-count variability at the presynaptic population level, for each value of , we also varied the total number of synapses from . To maintain similar output responses to the main model at the onset of training, initial values of synaptic weights were scaled by a multiplicative factor of .

Robustness to background synaptic activity was tested by varying the background rate parameter . In the simulations presented in Figures 3, 4, 5, 6, and 7, was set to for all synapses for the first 100 ms of a simulation, and set to zero during the stimulus presentation. For the simulations in Figures S3B and S4B, was varied between , and remained at that level during the stimulus presentation. We quantify the background noise level by the ratio of the background and stimulus-dependent population rates, (see Equations 19 and 20). When the noise level is 1, during the stimulus presentation, as many background spikes arrive across the population of synapses as do ‘signal’ spikes that depend on the specific input pattern being presented. To maintain similar output responses at the onset of training, initial values of synaptic weights we scaled by a multiplicative factor of to account for the additional input spikes. For the association task (Figure S3B), these simulations were repeated with , and initial weight values scaled by an additional factor of .

Robustness to label noise was tested by varying the probability that models were given an incorrect classification label during training. For each presentation of an input pattern, the label was switched with probability ranging between . During training, the label was used as in the noise-free case to determine both the sign of plasticity and the presence of a depolarizing current at the soma. The error term , which scales the magnitude of plasticity, was computed with reference to the noisy labels.

Robustness to synaptic weight jitter was tested by perturbing the weights of trained models with random multiplicative noise. Models were tested on 20 presentations of each input pattern used in training. Each time a pattern was presented, a noise term was drawn independently for every synapse from a uniform distribution over the interval , with the scaling factor ranging from . Weights were multiplied by , and thresholded at zero to avoid negative weights. In Figure S4C, model performance is quantified as a function of the standard deviation of the scaling factor; corresponds to a weight jitter SD of .

Robustness to timing variability was tested for the precisely timed burst input condition by perturbing the times of input bursts. Models were tested on 20 presentations of each input pattern used in training. We considered two types of perturbation – shared and independent. For the shared case, each time a pattern was presented, a single noise term was drawn from a uniform distribution over the interval and all burst times were periodically shifted within the 400 ms stimulus presentation window by this value. In the independent case, noise terms were drawn for every synapse, and burst times were shifted independently. The parameter was varied between , and the noise level quantified by the standard deviation of the noise term. Note that for a timing jitter standard deviation of 100 ms in the independent noise case, corresponding to , the burst times are completely scrambled.

Structured connectivity

We assume throughout the majority of this study that synapses are distributed randomly, without reference to the particular classification problem that needs to be solved. For comparison, we also assessed the performance of models in which the solution was increasingly ‘hard-wired’ by the placement of synapses on specific sets of branches (Figures S3D and S3E). Guided by the results of Figure 4C, for the active model, when two features and formed a preferred association, their excitatory synapses were placed on a common set of basal dendrites, and separate sets of apical dendrites. Conversely, when two features formed a non-preferred association, their synapses were placed on separate basal dendrites and common apical dendrites. This arrangement should allow the model to take maximal advantage of supralinear integration in basal dendrites and sublinear integration in apical dendrites to solve the task. For the passive model, for both basal and apical dendrites, we separated synapses representing features forming preferred associations and clustered synapses forming nonpreferred associations.

To implement the clustering we randomly selected two equally sized and possibly overlapping sets of basal branches (from a possible 16), and two equally sized and possibly overlapping sets of apical branches (from a possible 8). Synapses to be clustered within a dendritic domain were randomly distributed across the same set of branches, whereas synapses to be dispersed were distributed across different sets. We parameterized the degree of structured connectivity in terms of the overlap between the two selected branch sets within each domain (1 minus the fraction of shared branches). For completely unstructured connectivity (0 in the parameterization), the two sets selected within each domain were identical, each consisting of all 16 basal or 8 apical branches. In this condition, as in the simulations in the main text, synapse placement does not depend on the classification labels. For maximally structured connectivity (1 in the parameterization), the two sets selected within each domain were disjoint, consisting of 8 basal or 4 apical branches each, such that convergence of input on common dendrites was completely prescribed by the classification labels. At intermediate levels the two sets selected within each domain partially overlapped. For instance, at 0.5 in the parameterization, the two sets of basal dendrites consisted of 12 branches each, half of which were common to both. For each condition we trained models using the classification labels that were favored by the imposed connectivity (‘favored labels’), and also with the contingencies switched (‘reversed labels’).

Active intrinsic conductances

To assess the possible role of additional active dendritic conductances beyond NMDA receptors, we extended the model to include a low density of dendritic voltage-dependent and channels, as well as axo-somatic and dendritic HCN channels (see Dynamics and variational equations). Dendritic conductances are responsible for fast dendritic spikes observed in pyramidal neurons, and serve to sharpen stimulus tuning and action potential timing (Ariav et al., 2003; Smith et al., 2013). HCN conductances are responsible for depolarizing currents, regulating the coupling between apical and somatic compartments and counterbalancing the effects of dendritic filtering (Magee, 1999; Harnett et al., 2015). We chose parameters such that dendritic spikes were reliably triggered by synaptic input and the HCN-dependent voltage sag response to hyperpolarizing input was consistent with experimental observations (Kalmbach et al., 2018). We did not seek to precisely reproduce the experimental data, but rather to capture these effects qualitatively in a regime of dendritic excitability beyond that of the main NMDA-dependent active model.

We initially proceeded as before by performing many simulations with random input to fit plasticity kernels, as in Figure 2B. However, we found that using the same spike-triggered average approach as for the main models, we could not account for a comparable amount of the variance in . Specifically, whereas the local approximation accounts for of the variance in the active-NMDA model (Figure 2C), in the presence of dendritic spikes this was reduced to in basal dendrites and in apical dendrites. As this discrepancy would make a comparison between models unbalanced, we instead performed simulations in which the plasticity rule of Equation 27 used values of computed directly for each synapse via numerical integration of (Equation 13), (Equation 14), (Equation 15), (Equation 16), (Equation 17) (i.e., we did not apply the approximation of Equation 26). We used the same sub-simulation procedure as described in Plasticity kernels, with the difference that the equations were integrated directly for synapses receiving multiple inputs, rather than copying synapses to solve separately for each individual synaptic activation. We applied this strategy in simulations with the original passive and active models, an active model with additional HCN conductances, an active model with additional dendritic and conductances, and an active model with additional HCN, and conductances. All models were trained and tested on the rate-coded association task of Figure 3. To test whether the models with additional dendritic mechanisms obeyed the same rules of integration, we simulated and analyzed the subthreshold response to input features as in Figure 3D, setting the fast axo-somatic and conductances to zero, but leaving all other conductances intact.

Quantification and statistical analysis

Two-tailed Wilcoxon signed-rank tests were used to compare the performance and properties of models that were trained and tested on the same input patterns. Sample sizes and p values are provided in the figures and figure legends.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Code for performing simulations	This paper	https://github.com/babicknell/Dendrites; https://doi.org/10.5281/zenodo.5524314
Reconstructed cell morphologies	Allen Institute for Brain Science, 2015	https://celltypes.brain-map.org; ID: 502359001, 521409057
Software and algorithms
Python 3	Van Rossum and Drake, 2009	https://www.python.org
NEURON	Hines and Carnevale, 2006	https://www.neuron.yale.edu/neuron
Numpy	Harris et al., 2020	https://numpy.org
Numba	Lam et al., 2015	http://numba.pydata.org