To introduce our theory, we consider a postsynaptic neuron computing some function, $f$, of the activity of its presynaptic partners, $\mathbf{u}$ (Figure 1A, top):

where $\dot{v}$ is the resultant temporal change of the activity of the postsynaptic neuron. We chose $\mathbf{u}$ and $v$ to be analog variables, rather than for example digital spike trains, in line with the vast bulk of theories of network computations (Hopfield, 1984; Dayan and Abbott, 2001; Pouget et al., 2003) and experimental results suggesting analog coding in the cortex (London et al., 2010; Shadlen and Newsome, 1998). In particular, we considered these variables to correspond to the coarse-grained (low-pass filtered) somatic membrane potentials of neurons (in particular, excluding the action potentials themselves, as often reported in experimental data; Carandini and Ferster, 2000), although the theory can equally be formalized in terms of instantaneous firing rates (Materials and methods, Figure 1-figure supplement 1).

Figure 1 with 3 supplements see all

Download asset Open asset

The standard description of neural circuit dynamics in Equation 1 hides an important informational bottleneck intrinsic to the operation of cortical circuits. While according to Equation 1, the postsynaptic neuron’s analog activity, $v$, is required to depend directly on the analog activities of its presynaptic partners, $\mathbf{u}$, in reality it only accesses these presynaptic activities through the spikes the presynaptic population transmits, $\mathbf{s}$, incurring a substantial loss of information (Alenda et al., 2010; Sengupta et al., 2014). Therefore, the function $g$ a neuron actually implements on its inputs can only depend directly on the presynaptic spikes, not the underlying activities (Figure 1A, bottom):

Importantly, while $f\left(\mathbf{u}\right)$ is dictated by the computational function of the circuit, the actual transformation of the synaptic input to the postsynaptic response, expressed by $g\left(\mathbf{s}\right)$, is determined by the morphological and biophysical properties of the cell. (For these purposes, we regard the presynaptic side of synapses, transforming presynaptic spike trains to synaptic transmission events, as conceptually being part of the postsynaptic cell’s $g\left(\mathbf{s}\right)$ function.) How can then the neuron integrate the incoming presynaptic spikes, as formalized by $g\left(\mathbf{s}\right)$, such that the resulting postsynaptic response best matches the required computational function, $f\left(\mathbf{u}\right)$, thereby alleviating the fundamental informational bottleneck of spiking-based communication?

Determining the best $g\left(\mathbf{s}\right)$ is nontrivial because the same presynaptic spike train may be the result of many different underlying presynaptic activities (Paninski, 2006), each potentially implying a different output of the computational function. This ambiguity is formalized mathematically as a posterior probability distribution, $\text{P}\left(\mathbf{u}\right|\mathbf{s})$, expressing the probability that the analog activities of the presynaptic cells might currently be $\mathbf{u}$ given their spike trains, $\mathbf{s}$ (Pfister et al., 2010; Ujfalussy et al., 2011). The optimal response, i.e. the $g\left(\mathbf{s}\right)$ that minimizes the average squared error relative to $f\left(\mathbf{u}\right)$, is the expectation of $f\left(\mathbf{u}\right)$ under the posterior:

Crucially, the expression for the posterior, given by Bayes’ rule, is:

Note that while Equations 3–4 do not reveal directly the specific biophysical properties the postsynaptic cell should have, they tell us phenomenologically what signal integration properties should result from its biophysical properties. In particular, they make it explicit that the optimal $g\left(\mathbf{s}\right)$ depends fundamentally on two factors (Figure 1B, top):

1. the function that needs to be computed, $f\left(\mathbf{u}\right)$, and

2. the statistics of presynaptic activities: $\text{P}\left(\mathbf{u}\right)$, the prior probability distribution characterizing the long-run statistics of multi-neural activity patterns in the presynaptic ensemble, and the likelihood $\text{P}\left(\mathbf{s}\right|\mathbf{u})$, expressing the potentially probabilistic relationship between analog activities (e.g. somatic membrane potential trajectories) and emitted spike trains.

In the following, we show that the outcome of the integration of presynaptic spike trains in cortical neurons approximates very closely the optimal response, and that dendritic nonlinearities are crucial for achieving this near-optimality. For this, 1) we make an assumption about the computational function of the postsynaptic cell, $f\left(\mathbf{u}\right)$ (Figure 1B, top right); 2) we constrain presynaptic statistics, $\text{P}\left(\mathbf{u}\right)$ and $\text{P}\left(\mathbf{s}\right|\mathbf{u})$, by in vivo data about cortical population activity patterns (Figure 1B, top left); and with these 3) we compute the optimal response they jointly determine for various stimulation patterns (Figure 1B, bottom left and right).

To specify our model, we considered the case when $f\left(\mathbf{u}\right)$ itself is linear. Although networks with purely linear dynamics can perform non-trivial computations already (Dayan and Abbott, 2001; Hennequin et al., 2014), in the general case, we do expect $f\left(\mathbf{u}\right)$ to be nonlinear, e.g. sigmoidal (Hopfield, 1984). Nevertheless, in typical electrophysiological experiments only a small fraction of the full dynamic range of a neuron’s total input is stimulated (Figure 1—figure supplement 2), and so we approximate the computational function, $f\left(\mathbf{u}\right)$, as being linear on this limited input range without loss of generality. (See Figure 1—figure supplement 3 for the application of the theory to the case of nonlinear $f$.) Yet, as we show below, for physiologically realistic statistics of presynaptic activity patterns, the optimal response combines input spike trains in highly nonlinear ways even in the case of linear computation, predicting experimentally characterized nonlinearities in dendritic input integration. In particular, second- and higher-order prior presynaptic correlations, represented by $\text{P}\left(\mathbf{u}\right)$, will have a major role in determining the form of the corresponding optimal response. The likelihood, $\text{P}\left(\mathbf{s}\right|\mathbf{u})$, also influences the optimal response, but only in its quantitative details, as it does not involve correlations across neurons: each neuron’s firing is independent from the others’, given its own somatic membrane potential (Materials and methods).

Previous suggestions for how postsynaptic neurons achieve reliable computation despite the substantial ambiguity about the individual presynaptic activities relied on the linear averaging of inputs arriving from a sufficiently large pool of presynaptic neurons (Dayan and Abbott, 2001; Pfister et al., 2010). However, linear averaging is only guaranteed to produce the correct output, as dictated by Equations 3-4, if the activities of presynaptic neurons are statistically independent under the prior distribution, i.e. $\text{P}\left(\mathbf{u}\right)={\prod }_i\text{P}\left(u_i\right)$ (Materials and methods). In contrast, the membrane potential (Crochet et al., 2011) and spiking (Cohen and Kohn, 2011) of cortical neural populations often show complex patterns of correlations, which include both ‘spatial’ (cross-correlations between different neurons) and temporal components (auto-correlations, i.e. the correlation of the activity of the same cell with itself at different moments in time). Thus, in this more general case, we expect the optimal response to involve a nonlinear integration of spike trains. While temporal correlations alone do not require nonlinear dendritic integration across synapses, only local nonlinearities within each synapse, as brought about e.g. by short term synaptic plasticity (Pfister et al., 2010), spatial correlations require the non-linear integration of spikes emitted by different presynaptic neurons.

To illustrate that presynaptic spatial correlations require nonlinear integration across synapses, we compared the best linear response to a given presynaptic spike pattern with the optimal response (Equation 3, as approximated by Equation 23) for two different input statistics that differed only in the correlations between the presynaptic cells but not in the activity dynamics or spiking of individual neurons (temporal correlations). To compute the postsynaptic response, we assumed that dendritic integration in the postsynaptic cell was linear but, in order to dissect the role of dendritic integration across synapses from the effects of nonlinearities in individual synapses, we allowed spikes from the same presynaptic neuron still to be integrated nonlinearly (Pfister et al., 2010). In the first case (Figure 1C), when the presynaptic neurons were independent, the best linear response was identical to the optimal response. However, if presynaptic neurons became correlated, the optimal response became nonlinear and the best linear response was unable to accurately follow the fluctuations in the input (Figure 1D).

Thus, inputs from presynaptic neurons whose activity tends to be correlated need to be nonlinearly integrated, while inputs from uncorrelated sources need to be integrated linearly. This could be naturally achieved in the same dendritic tree by clustering synapses of correlated inputs to efficiently engage dendritic nonlinearities, while distributing the synapses of uncorrelated inputs on different dendritic branches (Larkum and Nevian, 2008). Crucially, for correlated inputs it is also necessary that the dendritic nonlinearities have just the appropriate characteristics for the particular pattern of correlations in presynaptic activities.

In order to systematically study the nonlinearities in the optimal response in the face of naturalistic input patterns, we derived and analyzed its behavior for a flexible class of richly structured, correlated inputs. Our statistical model for presynaptic activities, specifying the parametric forms of $\text{P}\left(\mathbf{u}\right)$ and $\text{P}\left(\mathbf{s}\right|\mathbf{u})$ (Materials and methods and Figure 2—figure supplement 1), was able to generate a variety of multi-neural activity patterns resembling the statistical properties described in in vitro and in vivo multielectrode recordings of neuronal population activities (Figure 2A and D show two representative examples). Once we have specified the statistical model of presynaptic activities, it uniquely determined the optimal response to any given input pattern (Equations 3–4). Thus, we used the same statistical model in two fundamentally different ways: first, to generate “naturalistic” in vivo-like patterns of presynaptic membrane potential traces and spike trains; and second, to compute the optimal response pattern to any stimulation pattern, be it “naturalistic” or parametrically varying “artificial” as used in typical in vitro experiments.

Figure 2 with 1 supplement see all

Download asset Open asset

The optimal response determined by this statistical model, for essentially any setting of parameters, was inherently nonlinear because the additional effect of a presynaptic spike depended on the pattern of spikes that had been previously received from the presynaptic population. Temporal correlations in the presynaptic population caused the optimal response to depend on the spiking history of the same cell (Pfister et al., 2010), while crucially, the additional presence of spatial correlations introduced a dependency on the past spikes of other cells. Thus, the integrated effect of multiple spikes could not be computed as a simple linear sum of their individual effects in isolation. Specifically, a spike that was consistent with the information already gained from recent presynaptic spikes had only a small effect on the response (Figure 2B). Conversely, a spike that was unexpected based on the recent spiking history caused a larger change (Figure 2E).

As could be anticipated based on Equations 3-4, whether a spike counted as expected or unexpected relative to recently received spikes, and hence whether it had a small or large postsynaptic effect, depended on the long-run prior distribution of presynaptic activities, $\text{P}\left(\mathbf{u}\right)$. As a result, the same pattern of presynaptic spikes led to qualitatively different responses under different prior distributions. In particular, sublinear integration was optimal when presynaptic activities exhibited Gaussian random walks and thus they did not contain statistical dependencies beyond second order correlations (Figure 2A-C), as seen in the retina and cortical cultures (Schneidman et al., 2006). This was because with the activities of presynaptic neurons being positively correlated, successive spikes conveyed progressively less information about the presynaptic signal resulting in sublinear integration (Figure 2C) and the strength of the sublinearity depended on the magnitude of correlations (Ujfalussy and Lengyel, 2011). In contrast, supralinear integration was optimal when the presynaptic population exhibited coordinated switches between distinct states associated with large differences in the activity levels compared to activity-fluctuations within each state (Figure 2D–F). These switches led to higher order statistical dependencies as seen in the cortex in vivo, either due to population-wide modulation by cortical state (Gentet et al., 2010; Crochet et al., 2011), or due to stimulus-driven activation of particular cell ensembles (Harris et al., 2003; Ohiorhenuan et al., 2010; Miller et al., 2014). In this case, while observing a few spikes was consistent with random membrane potential fluctuations within the quiescent state, thus only warranting a small response, further spikes suggested that the presynaptic population was in the active state now and thus the response should be larger, leading to supralinear integration (Figure 2F).

Note, that nonlinearities in the optimal postsynaptic response needed not simply compensate for the nonlinearities in the presynaptic spike generation process, as captured by $\text{P}\left(\mathbf{s}\right|\mathbf{u})$, but they critically depended on the presynaptic correlations, as captured by $\text{P}\left(\mathbf{u}\right)$. Indeed, in Figures 1C, D and 2A–F, the same spiking nonlinearity was used and yet very different input integration was required depending on the form of the presynaptic statistics: linear integration for uncorrelated inputs (Figure 1B) and nonlinear integration for correlated inputs (Figure 1C), with sub- or supralinear integration being optimal depending on whether only second order (Figure 2A–C) or also higher order correlations were present in the presynaptic population (Figure 2D–F). Moreover, optimal input integration remained nonlinear even if the postsynaptic neuron computed a function of the presynaptic firing rates (rather than membrane potentials) which were linearly related to spikes (Figure 1—figure supplement 1).

The nonlinear input integration seen in the optimal response strongly resembled dendritic nonlinearities. Indeed, the basic biophysical mechanisms present in dendrites naturally yield nonlinearities that are qualitatively similar to those of the optimal response: purely passive properties lead to sublinear integration (Koch, 1999), whereas locally generated dendritic spikes endow dendrites with strong supralinearities (Nevian et al., 2007; Branco and Häusser, 2011). However, the full mathematical implementation of the optimal response is excessively complex (Materials and methods) and thus, there is unlikely to be a one-to-one mapping between the variables necessary for implementing it and the biophysical quantities available in dendrites. Therefore, we sought to establish a formal proof that dendritic-like dynamics can implement, even if approximately, the optimal response. For this, we considered two limiting cases of our statistical model of presynaptic activities, $\text{P}\left(\mathbf{u}\right)$ and $\text{P}\left(\mathbf{s}\text{|}\mathbf{u}\right)$, and compared the properties of the corresponding optimal response to a well-established simplified model of nonlinear dendritic integration, using a combination of analytical and numerical techniques.

First, we considered a limiting case in which the statistics of a large presynaptic population were strongly dominated by the simultaneous switching of presynaptic neurons between a quiescent and an active state (as shown in Figure 2D). In this limiting case we were able to show mathematically (see Materials and methods) that a simple, biophysically-motivated, canonical model of nonlinear dendritic integration (Poirazi and Mel, 2001) is able to produce responses that are near-identical to the optimal response for any sequence of presynaptic spikes (Figure 3A, see also Figures 4C). In this simple dendritic model (Figure 3A, inset; Equations 24–25), inputs within a branch are integrated linearly and the local dendritic response is then obtained by transforming this linear combination through a sigmoidal nonlinearity, which is a hallmark of supralinear behavior in dendrites (Poirazi et al., 2003b).

Figure 3 with 2 supplements see all

Download asset Open asset

Figure 4 with 2 supplements see all

Download asset Open asset

Second, we considered another limiting case in which the statistics of the presynaptic population were fully characterized by second-order correlations (as shown in Figure 2A). In this case, the same type of dendritic model, but with a sublinear input-output mapping, was able to approximate the optimal response very closely. Although a closed-form solution for the optimal nonlinear mapping could not be obtained in this case, it could be shown to be sublinear (Appendix), and was well approximated by a sigmoidal nonlinearity parametrized to be dominantly saturating (Figures 3B and Figure 3—figure supplement 1).

We also noted that it was the same type of sigmoidal nonlinearity which could implement supralinear and sublinear integration depending on the input regime (low background, synchronous spikes: supralinear; high background, asynchronous spikes: sublinear integration, compare Figure 3A and B, inset). This suggests that dendritic integration may adapt to systematic changes in presynaptic statistics, such as those brought about by transitioning between the desynchronized and synchronized states of the neocortex, or sharp waves and theta activity in the hippocampus, without having to change the parameters of its nonlinearity (Borst et al., 2005) (Figure 3—figure supplement 2). Indeed, Gasparini and Magee (2006) demonstrated that dendritic integration in hippocampal pyramidal cells was supralinear when inputs were highly synchronized (as they are during sharp waves), while integration was linear if the input was asynchronous (such as during theta activity).

While the foregoing analyzes proved that dendritic-like nonlinearities can closely approximate the optimal response in certain limiting cases, they do not address directly whether having such nonlinearities in input integration is crucial for attaining near-optimal computational performance for more realistic input statistics, or simpler forms of input integration could achieve similar computational power. To study this, we considered a scenario in which the presynaptic population consisted of four ensembles, such that neurons belonging to each ensemble underwent synchronized switches in their activity levels which were independent across the four ensembles, while there were also independent fluctuations in the activity of individual presynaptic neurons which were comparable in magnitude to those caused by these synchronized activity switches (Figure 4A). We then assessed the performance of four different variants of a simple dendritic model relative to that of the optimal response (Figure 4B): a model with linear dendrites and soma; a model in which only the soma was nonlinear, and two models in which nonlinearities resided in the dendrites with either random or clustered connectivity between the presynaptic assemblies and the dendritic branches.

We quantified the performance of each of the models based on how closely their output approximated the linear average of the analog presynaptic activities giving rise to the spike trains they were integrating (Figure 4—figure supplement 1, Materials and methods). For a fair comparison, we tuned the parameters of each variant of the dendritic model to obtain the best possible performance with it (Figure 4C). The model with nonlinear dendrites and clustered connectivity had near-optimal cross-validated performance (Figure 4D) while all other models performed significantly worse (n = 20 runs, t = 51, t = 35, t = 20, and P$<$10${}^{-15}$, P$<$10${}^{-15}$, P$<$10${}^{-13}$; respectively from left to right as shown in Figure 4D). This remained true when we varied the number and firing rate of presynaptic neurons over a wide range, and under a diverse set of qualitatively different population-level statistics, determining the dynamics of assembly switchings and within-assembly membrane potential correlations (Figure 4—figure supplement 2).

Taken together, these results demonstrate that the clustering of correlated inputs together with nonlinearities akin to those found in dendrites is necessary to achieve optimal estimation performance in the face of presynaptic correlations. However, in order to be tractable, our dendritic model was mathematically simplified and, as a result, only qualitatively reproduced the nonlinearities of real dendrites. Thus, we directly compared experimentally recorded responses in dendrites to the optimal response.

A crucial prediction of our theory is that dendritic nonlinearities act to achieve near-optimal responses in a way that the form of the nonlinearity is specifically matched to the long-run statistics of the presynaptic population. We tested this prediction in experiments in which two different types of cortical pyramidal neurons, from layer 2/3 of the neocortex (Figure 5A–E) and from area CA3 of the hippocampus (Figure 5F–J), received patterned dendritic stimulation using two-photon glutamate uncaging, and compared their subthreshold somatic responses with the optimal responses predicted by the theory.

Figure 5 with 2 supplements see all

Download asset Open asset

For generating our predictions of the optimal response in these two cell types, we fitted the parameters describing presynaptic statistics in our model, $\text{P}\left(\mathbf{u}\right)$ and $\text{P}\left(\mathbf{s}\right|\mathbf{u})$, to the statistical patterns in the activity of their respective presynaptic populations. For neocortical pyramidal cells, we fitted in vivo data available on the membrane potential fluctuations of layer 2/3 pyramidal cell-pairs in the barrel cortex during quiet wakefulness (Gentet et al., 2010; Crochet et al., 2011) (NC, Figure 5A, Table 2). For hippocampal pyramidal cells, we fitted presynaptic statistics to membrane potential fluctuations (Ylinen et al., 1995; English et al., 2014) and to multineuron spiking patterns of hippocampal pyramidal cells during sharp wave activity (Csicsvari et al., 1999; 2000) (HP, Figure 5F, Table 3). Due to the limitations of available hippocampal data sets, extracellular rather than intracellular data was used for fitting correlations. The motivation for our choice of the particular neocortical and hippocampal states used for fitting presynaptic statistics was two-fold. First, the general network state of the slice preparations in which we tested dendritic integration was likely most analogous to these states (A Gulyás, personal communication; see also Karlocai et al., 2014; Schlingloff et al., 2014), characterized by relatively suppressed neural excitability due to low levels of cholinergic modulation (Harris and Thiele, 2011; Eggermann et al., 2014). Second, the stimulation protocol used in our study (short bursts of synaptic stimuli following longer silent periods) was also most consistent with population activity during hippocampal sharp waves and quiet wakefulness in the cortex. In order to capture variability across the cells we recorded from, the parameters related to postsynaptic dendritic filtering (amplitude and decay of the response to a single stimulation, and the size of the dendritic subunit, Figure 5—figure supplement 1B-C) were tuned for the individual dendrites. Importantly, the parameters describing presynaptic statistics were fitted without regard to our dendritic experimental data, thus allowing a strong test of our predictions about dendritic integration (see Materials and methods).

We found that the non-linear integration of individual spike patterns in cortical neurons was remarkably well fit by the optimal response when it was tuned to the correct presynaptic statistics (Figure 5C,H). The systematic dependence of response amplitudes on the inter-stimulus interval (ISI) in individual cells (Figure 5D,I) was also well predicted by the optimal response. We quantified the quality of match between the predicted and experimentally recorded time course of responses across a population of n = 6 (neocortex) and n = 6 (hippocampus) dendrites under a range of conditions varying ISI or the number of stimuli, and found that the precision of our predictions was not statistically different from that expected from the inherent variability of responses in individual dendrites (Figure 5E,J; neocortex: t = 0.2, P = 0.85; hippocampus: t = 1.85, P = 0.12). In contrast, when the optimal response was tuned to unrealistic presynaptic statistics characterized purely by second-order correlations (cor2), or by a lack of any correlations implying statistically independent presynaptic firing (ind), the quality of fits became significantly worse (Figure 5E,J; neocortex: t = −4.6, P = 0.006 for cor2, and t = −4.9, P = 0.004 for ind; hippocampus: t = −4, P = 0.01 for cor2, and t = −4.9, P = 0.004 for ind).

Moreover, using realistic presynaptic statistics, but matching hippocampal rather than neocortical activities, also resulted in significantly worse fits for neocortical responses (Figures 5E; t= −3.6, P = 0.02). The converse was not observed in the case of hippocampal neurons (Figures 5J; t= 0.43, P = 0.68). This might be because hippocampal neurons also receive neocortical inputs (albeit on their apical not basal dendrites) that show similar population activity patterns to the ones we matched here for the neocortical cells (Isomura et al., 2006), while the primary sensory cortical pyramidal cells we recorded from do not receive direct input from the hippocampus. Nevertheless, when we analyzed the quality of fit between our predictions and recorded responses in hippocampal and neocortical data together, we found a small, but significant interaction between the source of the input statistics (neocortex or hippocampus) and the location of the postsynaptic neurons (ANOVA F = 5, P $<$ 0.05). This suggests that dendritic nonlinearities in cortical pyramidal neurons are specifically tuned to the dynamics of their presynaptic cortical ensembles. Furthermore, the blockade of NMDA receptor activation by AP5 resulted in dendritic responses that afforded substantially poorer fits by the model, even after refitting the postsynaptic parameters (Figure 5E,J, AP5). This indicated that the fine tuning of dendritic nonlinearities to input statistics relied on the action of NMDA receptors.

As dendrites in both of our cortical cell types integrated inputs supralinearly, as a further control, we analyzed similar data available from cerebellar stellate cell dendrites, which are known to integrate their inputs sublinearly (Abrahamsson et al., 2012) (Figure 5—figure supplement 2). In this case, we fitted the statistics of individual presynaptic cells to those of cerebellar granule cells. The correlations between these cells are less known, but we found that assuming simple second-order correlations made the optimal response a close match to dendritic responses. In contrast, the hippocampal- or neocortical-like statistics that were crucial for matching responses in cortical dendrites (Figure 5D,H) resulted in a substantially poorer fit in this cerebellar cell type. This demonstrates a double dissociation in the matching of cortical and subcortical neuron types to cortical and non-cortical input statistics.

The central insight of our theory is the relationship between presynaptic statistics and postsynaptic response, formalized as the optimal response. The optimal response can be expressed as a set of nonlinear differential equations that requires storing and continuously updating $~N^2$ variables within the dendritic tree, where $N$ is the number of synapses (Materials and methods). Thus, it is unlikely to be implemented by the postsynaptic neuron as such. Consequently, to demonstrate the biophysical feasibility of our theory, we derived a simple approximation to the optimal response that performs about equally well using just a few postsynaptic variables and that corresponds to a canonical descriptive model of dendritic integration (Poirazi et al., 2003; Poirazi et al., 2003b).

We found that simple second order correlations between presynaptic neurons imply sublinear integration which can be implemented by the saturating nonlinearity characteristic of passive dendrites. Conversely, the biophysical substrate for the type of supralinear integration that was optimal for state-switching dynamics likely involves NMDA receptors because the particular dendritic nonlinearites observed in the cortical cells in which we tested our theory are known to be mediated primarily through NMDA receptor activation (Branco et al., 2010; Makara and Magee, 2013; Major et al., 2013). Indeed, we found that pharmacological inactivation of NMDA receptors abolished the precise match between dendritic integration and presynaptic statistics in these neurons (Figure 5). Moreover, the local plateau potentials generated by NMDA currents have been shown to have graded response durations (Major et al., 2008), and the resulting nonlinearities could be continuously tuned between weaker and stronger forms (boosting and bistability, Major et al., 2013). These properties make NMDA receptor mediated dendritic nonlinearities ideally suited for being matched to presynaptic statistics, as the optimal response involves sustained dendritic depolarisations of varying duration (Figure 4) that depend parametrically on those statistics.

A central prediction of our theory that awaits confirmation is the existence of a tight relationship between the structure of correlations in the activity of presynaptic cells and the morphological clustering of their synapses on the postsynaptic dendrite. This is because our theory requires nonlinear integration of inputs from neurons with statistically dependent activity, while spikes from independent neurons need to be integrated linearly. Biophysical considerations suggest (Koch, 1999) and experimental data supports (Polsky et al., 2004; Losonczy and Magee, 2006) that, when synchronous, nearby inputs on the same dendritic branch are summed nonlinearly, whereas widely separated inputs are combined linearly. Consequently, our theory predicts that the correlation structure of the inputs will be mapped to the dendritic tree in a way that presynaptic neurons with strongly correlated activities target nearby locations while independent neurons innervate distinct dendritic subunits.

According to our theory, the kind of correlation relevant for determining synaptic clustering is the ‘marginal’ correlations between the membrane potentials of presynaptic neurons. Marginal correlations include both signal and noise correlations (Averbeck et al., 2006) and thus can reach substantial magnitudes even when noise correlations alone are small, as e.g. during desynchronized cortical states (Renart et al., 2010), especially for neurons with overlapping receptive fields (Froudarakis et al., 2014), and when measured between the membrane potentials of neurons rather than their spike counts (Dorn and Ringach, 2003; de la Rocha et al., 2007; Poulet and Petersen, 2008).

At the level of different dendritic regions, the segregation of different input pathways along the dendritic tree of hippocampal neurons supports this prediction (Witter et al., 1989; Druckmann et al., 2014). At the level of individual synapses, the degree and the existence of clustering among inputs showing correlated activity is currently debated. High resolution imaging revealed subcellular topography of sensory inputs in the tadpole visual system (Bollmann and Engert, 2009), clustered patterns of axonal activity in the parallel fibres that provide input to cerebellar Purkinje cells (Wilms and Häusser, 2015), and experience-driven synaptic clustering in the barn owl auditory localization pathway (McBride et al., 2008). Furthermore, it has been demonstrated that neighboring synapses are more likely to be coactive than synapses that are further away from each other in developing hippocampal pyramidal cells (Kleindienst et al., 2011) as well as in hippocampal cultures and in vivo in the barrel cortex during spontaneous activity (Takahashi et al., 2012). These results thus suggest clustering of correlated inputs.

In contrast, an interspersion of differently tuned orientation-, frequency- or whisker-specific synaptic inputs on the same dendritic segments was found in the mouse visual, auditory or somatosensory cortex, respectively, thus challenging the prevalence of synaptic clustering (Jia et al., 2010; Chen et al., 2011; Varga et al., 2011). However, in all these studies the stimuli used were non-naturalistic and varied along a single stimulus dimension only (direction of drifting gratings, pitch of pure tones, or the identity of the single whisker being stimulated), which may account for the apparent lack of clustering. In particular, our theory predicts clustering based on the long-term statistical dependencies between the responses of the presynaptic neurons for naturalistic inputs, which can be quite poorly predicted from their tuning properties for single stimulus dimensions (Harris et al., 2003; Fiser et al., 2004). In contrast, the statistical dependencies relevant for our theory are well represented by those found during spontaneous activity (Berkes et al., 2011). Indeed, studies finding evidence in favor of synaptic clustering analyzed the structure of synaptic input to dendritic branches during spontaneous network activity (McBride et al., 2008; Kleindienst et al., 2011; Makino and Malinow, 2011; Takahashi et al., 2012). Thus, presynaptic correlations for naturalistic stimulus sets may be predictive of synaptic clustering and providing more direct evidence for or against such clustering will offer a crucial test of our theory.

Although, in general, we expect single-neuron computations to be nonlinear (Zador, 2000), and our theory indeed applies to nonlinear computations (Figure 1—figure supplement 3), we assumed the postsynaptic computation to be linear for matching experimental data. This choice was justified by two reasons. First, it is difficult to determine, without making strong prior assumptions, what kind of nonlinear function the neuron actually computes; and so the choice of any particular such function would have been arbitrary. Note that even in relatively well-characterized cortical areas (such as the visual cortex) it is unknown how much of the computationally relevant output of individual neurons (such as orientation or direction selectivity) is brought about by specific nonlinearities in the input-output transformations of these neurons, or by multiple steps of feed-forward and recurrent processing carried out at various stages of the visual pathway between the retina and those neurons. Moreover, in some cases, even networks with linear dynamics can provide a remarkably close fit to experimentally observed cortical population dynamics (Hennequin et al., 2014). This issue may be best addressed in systems that are more specialized than the cortex so that there are well-supported hypotheses about the particular nonlinear computations individual neurons need to perform, such as the fly visual system (Single and Borst, 1998) or the mammalian and avian auditory brain stem (Agmon-Snir et al., 1998). In order to test our theory in these systems, in vivo multineural data will need to be collected from the afferent brain areas, preferably in the unanesthetized animal, for characterising the relevant statistical properties of the presynaptic population to which dendritic nonlinearities are adapted according to our prediction.

Second, any nonlinear function can be approximated to high precision by a linear function over a sufficiently limited input range. Currently available experimental techniques for systematically probing dendritic nonlinearities, including those used in our study, only provide data over such a very limited range ($\sim$0.1% of the number of excitatory inputs impinging a neuron, Megías et al., 2001). Inputs in this small range do not sufficiently engage global nonlinearities brought about by active somatic conductances or global events such as Ca${}^{\text{2+}}$ spikes. Thus, we could assume linear computation over this range without loss of generality (Figure 1—figure supplement 2). In fact, from this perspective, it is a non-trivial phenomenon to account for on its own right that stimulating such a small fraction of inputs already leads to observable nonlinearities in the postsynaptic dendrite. By defining the computation to be linear, we could demonstrate that such strong dendritic nonlinearities arise naturally in our theory, entirely due to the correlations in the prior input statistics, thus providing a functional account for this remarkable phenomenon.

Once patterned dendritic stimulation over a broader and more realistic range of inputs becomes feasible, our theory will provide a principled method for dissecting the roles of presynaptic correlations vs. genuine nonlinear computations in shaping dendritic nonlinearities. A sufficiently rich set of such data will allow the fitting of presynaptic parameters, as we did here, followed by fitting postsynaptic transfer functions to dendritic responses without having to make strong prior assumptions about their (linear) nature.

Our formalism was based on the assumption that cortical neurons only influence each other’s membrane potentials via the action potentials they emit. While there exist other, more analog forms of communication, such as the modulation of the effects of action potentials by subthreshold potentials (Clark and Häusser, 2006), the propagation of voltage signals through gap junctions (Vervaeke et al., 2012), and ephaptic interactions between nearby cells (Anastassiou et al., 2011), these either require slow membrane potential dynamics, small distances between interacting cells, or large degrees of population synchrony, and are thus generally believed to have a supplementary role beside spike-based communication (Sengupta et al., 2014). Note that our theory is self-consistent even though it considers spiking only in the presynaptic population and not in the postsynaptic neuron. This is because we assumed that the computationally relevant mapping is that between the membrane potentials of the presynaptic neurons and the postsynaptic cell (Figure 1A, Equation 1), and so, by induction, the spikes of the postsynaptic neuron will effect the mapping from its membrane potential to those of its postsynaptic partners.

We also assumed that presynaptic spikes deterministically and uniformly impact the postsynaptic response, and thus apparently neglected the stochasticity in synaptic transmission, and in particular systematic variations in synaptic efficacy brought about by short-term synaptic plasticity. Nevertheless, these presynaptic features are compatible with our theory. The stochasticity of synaptic transmission, due to a baseline level of synaptic failures, is straight-forward to incorporate in the model by reducing the effective presynaptic firing rate, which can thus be interpreted as a ‘transmission rate’ instead. In fact, we have already done this while matching hippocampal and neocortical presynaptic statistics (Table 1).

Short-term synaptic plasticity, resulting in dynamical changes in synaptic efficacy as a function of the recent spiking history of the presynaptic neuron, is not only a constraint in our framework, but as we have shown in related work, it can act itself as an optimal estimator of the membrane potentials of individual presynaptic neurons (Pfister et al., 2010). Thus, the effects of short-term plasticity can be regarded as a special case of what can be expected from our optimal response: when presynaptic neurons are statistically independent, spikes arriving at different synapses are integrated linearly, and local nonlinearities acting at individual synapses suffice (Figure 1C, see also Materials and methods). However, the importance of nonlinear interactions between inputs from different presynaptic neurons, brought about by dendritic nonlinearities, rapidly increases with the magnitude of presynaptic correlations, especially in large populations (Figure 1D, see also Ujfalussy et al., 2011).

These considerations suggest that short-term synaptic plasticity and dendritic nonlinearities have complementary roles in tuning the postsynaptic response to the statistics of the presynaptic population along the orthogonal dimensions of time and space. The former is useful in the face of temporal correlations private to individual presynaptic neurons (auto-correlations, e.g., brought about by spike frequency adaptation, Pfister and Surace, 2014), while the latter is matched to spatio-temporal correlation patterns present across the presynaptic population.

We focused on the nonlinear integration of excitatory inputs in the dendritic tree of cortical neurons that have been extensively studied and described over the past decades, giving rise to a strong body of converging evidence as to their characteristics and mechanisms (Spruston, 2008). Recent work studying the nonlinear interaction between inhibitory and excitatory inputs in active dendrites (Gidon and Segev, 2012; Jadi et al., 2012; Müller et al., 2012; Wilson et al., 2012; Lovett-Barron et al., 2012) demonstrated that local inhibition has a powerful control over the excitability of the dendritic tree.

However, it is not yet clear whether these inhibitory inputs are directly involved in the computation performed by the circuit, just as excitatory neurons but with negative signs (Koch et al., 1982), or, alternatively, they may have a more ancillary role in supporting computations carried out primarily by excitatory neurons (Vogels et al., 2011).

Our theory can be extended to include both possibilities, by allowing inhibitory inputs to contribute to the computational function, $f\left(\mathbf{u}\right)$, with negative weights, or by considering them as providing auxiliary information about the common state of the excitatory presynaptic ensemble, especially when this state is in the more suppressed regime. Indeed, our preliminary results suggest that such an extension of our theory (Ujfalussy and Lengyel, 2013) successfully accounts for the interaction of (excitatory) Schaffer collateral inputs with the feedforward inhibitory effects of the temporo-ammonic pathway (Remondes and Schuman, 2002), likely mediated by interneurons delivering dendritic inhibition (Dvorak-Carbone and Schuman, 1999).

In the present paper we focused on dendritic integration in pyramidal neurons because dendritic nonlinearities have traditionally been more extensively characterized in this cell type, but our theory equally applies to synaptic integration in other types of neurons, including inhibitory interneurons. Therefore, our theory predicts a qualitative similarity of dendritic integration in different neuron types (i.e. interneurons versus principal cells) when they receive inputs from overlapping presynaptic populations. Indeed, it has been recently found that inhibitory interneurons can exhibit dendritic NMDA spikes under certain experimental circumstances (Katona et al., 2011; Chiovini et al., 2014) in addition to standard sublinear integration. The differences between dendritic integration in excitatory and inhibitory neurons could be attributed to their different computational function, $f\left(\mathbf{u}\right)$, or differences in the specific presynaptic populations innervating them.

According to our theory, the optimal response depends on prior information about the input statistics. Consequently, for dendritic processing to approximate the optimal response, this prior information needs to be implicitly encoded in the form of the particular nonlinearity a dendrite expresses. Therefore, our theory predicts an ongoing adaptation of dendritic nonlinearities to presynaptic firing statistics across several time-scales.

First, there is a simple yet potent mechanism implicit in our theory that can ensure that a match of dendritic integration to presynaptic statistics is maintained as those statistics are changing over time. This is based on the observation that essentially instantaneous, albeit probably incomplete, adaptation can occur even without specific changes in the integrative properties of dendrites per se, simply due to the fact that a critical level of input synchrony is required to elicit dendritic spikes, and so the same sigmoid-looking dendritic transfer function can be used as superlinear, linear, or sublinear, depending on which part of its input range is being used (Figure 3—figure supplement 2).

Second, to match the more specific modulation of the statistics of presynaptic activities by global cortical states (Crochet et al., 2011; Mizuseki and Buzsaki, 2014), dendritic integration may also be modulated by these states. As different cortical states are typically accompanied by changes in the neuromodulatory milieu (Hasselmo, 1995; Harris and Thiele, 2011), neuromodulators may be the ideal substrates to ensure that dendritic integration also changes according to the current cortical activity mode. This may provide a functional account of changes in the excitability of the dendritic tree as dynamically regulated by acetylcholine and monoamines (Sjöström et al., 2008).

Third, experience-dependent synaptic plasticity can gradually change the statistics of the presynaptic population activity implying that the optimal form of input integration should also change as a function of experience. We propose that branch-specific forms of plasticity of dendritic excitability (Losonczy et al., 2008; Makara et al., 2009; Müller et al., 2012) may have a functional role in enabling dendrites to adjust the form of input integration to such slowly developing and long-lasting changes in the statistics of their inputs.

Finally, whether inputs from two presynaptic cells are integrated linearly or nonlinearly in a dendrite depends critically on the distance between their synapses within the dendritic tree (Polsky et al., 2004; Losonczy and Magee, 2006). Our theory requires nonlinear integration of inputs from neurons with statistically dependent activity, predicting a mapping of presynaptic correlations on the postsynaptic dendritic tree. Local electrical and biochemical signals can drive synaptic plasticity (Larkum and Nevian, 2008; Govindarajan et al., 2011; Winnubst et al., 2015) and rewiring (DeBello, 2008) leading to synaptic clustering of correlated inputs along the dendritic tree (Kleindienst et al., 2011; Takahashi et al., 2012).

A combination of all these mechanisms may be crucial for achieving and dynamically maintaining, at the level of individual neurons, a detailed matching of dendritic nonlinearities to presynaptic statistics. Thus, our theory provides a novel framework for studying a range of phenomena regarding the dynamical regulation of dendritic nonlinearities from the perspective of circuit-level computations.

In order to study the optimal form of input integration with realistic input statistics, we need to make two important assumptions. First, we need to assume a particular algebraic form for the computation that a neuron performs. Second, we need to define what the relevant presynaptic statistics are, that is, the membrane potential and spiking dynamics of the presynaptic population under naturalistic in vivo circumstances. Given these two assumptions, the theory uniquely defines the optimal response of a neuron to any input pattern. The optimal response has qualitatively similar properties whether computations are defined as mappings between pre- or postsynaptic membrane potentials or firing rates (Figure 1—figure supplement 1).

Throughout the paper the term input refers to the spatio-temporal spiking pattern impinging the neuron while the response of a neuron refers to its (subthreshold) somatic membrane potential (or firing rate, see below). All parameters used in the paper are given in Table 1 or in the caption of the corresponding Figure.

We assumed the postsynaptic computation to be linear, i.e. the dynamics of the postsynaptic membrane potential $v\left(t\right)$ evolves according to a weighted sum of the presynaptic membrane potential values, $\mathbf{u}\left(t\right)$ (cf. Equation 1):

where ${\tau }_{\text{post}}$ is the time constant of the postsynaptic neuron, $N$ is the number of presynaptic neurons, and $w_i$ is the computational weight assigned to presynaptic neuron $i$. As the postsynaptic neuron cannot access presynaptic membrane potentials, $\mathbf{u}$, directly only the spikes the presynaptic cells emit, $\mathbf{s}$, (Figure 1, Equation 2), the optimal response (that minimizes mean squared error) is the expectation of Equation 5 under the posterior distribution of the presynaptic membrane potential at time $t$, $\mathbf{u}\left(t\right)$ given the history of presynaptic spiking up to that time, $\mathbf{s}(0:t)$ (cf. Equation 3):

Throughout the paper we call the output of Equation 6 the optimal response and compare its behavior to input integration in the dendrites of cortical pyramidal cells.

Table 1 shows the values of the parameters in Equation 5 ($N$, ${\tau }_{\text{post}}$, and $w_i$) used in the simulations. In short, to illustrate the contributions of inference to Equation 6 (the term including the integral), we used ${\tau }_{\text{post}}=0$ in Figures 1, 3 and 4 as well as in all Supplemental Figures, unless otherwise stated. We used ${\tau }_{\text{post}}=10$ ms in Figure 2 to aid comparison with experimental data and fitted ${\tau }_{\text{post}}$ to data for Figure 5. Throughout the paper we used $w_i=1/N$, except in Figure 5 where we fit $N$ and $w_i=w$ jointly to the data.

Computing the posterior, $\text{P}\left(\mathbf{u}\left(t\right)\right|\mathbf{s}(0:t))$ in Equation 6 requires a model for the joint membrane potential and spiking statistics of the presynaptic population, $\text{P}(\mathbf{u},\mathbf{s})$ (see also Equation 4). For mathematical convenience, we present some of our results below in discrete time with time step size $\delta t$, which we will eventually take to zero to derive time-continuous equations. We distinguish discrete and continuous time results by using time as an index versus as an argument of the corresponding time-dependent quantities, e.g. ${\mathbf{u}}_t$ vs. $\mathbf{u}\left(t\right)$.

We describe the joint statistics of presynaptic membrane potentials and spikes by a hierarchical generative model that has three layers of variables (Figure 2—figure supplement 1A,B). The global state of the system is described by a single binary variable, $z$ that switches between a quiescent ($-$) and an active ($+$) state following first-order Markovian dynamics (see Appendix for the extension to an arbitrary number of states). The transition rates to the active and quiescent states are given by ${\mathrm{\Omega }}_{+}$ and ${\mathrm{\Omega }}_{-}$, respectively.

The dynamics of (subthreshold) membrane potentials $\mathbf{u}$ are modeled as a multivariate Ornstein-Uhlenbeck (mOU) process, which yields random walk-like behavior that (unlike simple Brownian motion) decays exponentially towards a baseline defined by the resting potential $\overline{\mathbf{u}}$, which in turn depends on the momentary global state of the system, $z_t$:

where $\tau$ is the presynaptic time constant of the exponential decay, and $\mathbf{Q}$ is the ‘process noise’ covariance matrix determining the variance of individual membrane potentials (together with $\tau$) and, importantly, also the correlations between presynaptic neurons. It is straightforward to extend the model by also making these parameters state- (or in fact, neuron-) dependent.

Note that both the state switching and mOU components of this model introduce both spatial and temporal statistical dependencies in the membrane potentials and spike trains of presynaptic cells. In the rest of this paper, we informally refer to any statistical dependency (second or higher order) as ‘correlation’, and we write ‘auto-correlation’ when we refer to the correlations between the membrane potential (firing rate) values of the same neuron at different times, and ‘cross-correlation’ when referring to the correlation between the activities of two different cells (at the same time, or at different times). Also note that temporal and spatial correlations can not be studied in complete isolation in the case of smoothly varying signals, such as membrane potentials, as the cross-correlation between the activity of two presynaptic neurons always has a characteristic temporal profile. While it is possible to consider a presynaptic neuronal population completely lacking spatial correlations (i.e. independent presynaptic neurons, as in Figure 1C), having a population with only spatial but not temporal correlations would require the membrane potentials of the individual neurons to be temporally white noise – which is so far removed from reality that we did not consider this case worth pursuing.

More specifically, the timescale of temporal correlations (auto-correlations) in the model depend on the transition rates of the switching component, ${\mathrm{\Omega }}_{+}$ and ${\mathrm{\Omega }}_{-}$, and the presynaptic time constant of the mOU component, $\tau$, such that cells are auto-correlated as long as $\tau$, ${\mathrm{\Omega }}_{+}^{-1}$, and ${\mathrm{\Omega }}_{-}^{-1}>0$. Spatial correlation (cross-correlations between different presynaptic neurons) also emerge from both components. First, the pairs of presynaptic neurons corresponding to the non-zero off-diagonal elements of $\mathbf{Q}$ matrix of the mOU component become correlated. Second, synchronous state transitions during state switching in the presynaptic ensemble introduce positive correlations. Importantly, in both the temporal and spatial domains, while the mOU process can only introduce second-order correlations (i.e. it makes membrane potentials be distributed according to a multivariate normal), the switching process introduces higher order correlations (such that membrane potentials are not normally distributed any more). These higher order correlations are stronger when the effect of state-switching is large relative to the membrane potential fluctuations within a single state.

Finally, instead of modeling the detailed dynamics of action potential generation, we model spiking phenomenologically by introducing a single discrete variable, $s_{i,t}$, for each presynaptic neuron that represents the number of spikes neuron $i$ fires in time step $t$. (Note that in the limit $\delta t\to 0$ this variable becomes binary, i.e. there can never be more than one spike fired in a $\delta t$ time window.) Spiking in each cell only depends on the membrane potential of that cell, and follows an inhomogeneous Poisson process with the firing rate, $r$, being an exponential function of the membrane potential (Gerstner and Kistler, 2002):

where $\beta$ describes how deterministically the cell switches to firing at threshold ($u=0$) and $g$ is the firing rate at that threshold. We modeled the absolute refractory period by not allowing the generation of spikes (i.e. setting $r_{i,t}=0$) within a time window of length ${\tau }_{\text{refr}}$ following each spike in a cell, regardless of its membrane potential.

The parameters of the presynaptic statistics used in the paper are given in Table 1. Examples of neural dynamics generated by the model are shown in Figures 1, 2, 4, and 5.

Our goal was to infer the posterior distribution of the membrane potential based on the spiking pattern observed up to time $t$, $\text{P}\left({\mathbf{u}}_t\right|{\mathbf{s}}_{0:t})$.

We first show that linear dendritic integration is sufficient when presynaptic neurons are statistically independent. We start by noting that by marginalising out the past membrane potential history of the presynaptic cells and using Bayes’ rule, the posterior can always be written as

and as the spikes of each neuron are independent from all other neurons conditioned on its own membrane potential history (Equation 8), this can be rewritten as

In the special case when we assume that presynaptic neurons are statistically independent, i.e. their prior factorizes $\text{P}\left({\mathbf{u}}_{0:t}\right)={\prod }_i\text{P}\left(u_{i,0:t}\right)$, the posterior also becomes factorised

which in continuous time reads simply as

Thus, taking our usual assumption that the postsynaptic computation is linear (Equation 5), the optimal response in Equation 6 can be written as

indicating that integration of inputs from different neurons is linear in this case (it is a weighted sum of terms each depending on just a single presynaptic neuron). However, even in this case, note that integration of input spikes from the same presynaptic neuron, i.e. the result of the integral over each $u_i\left(t\right)$ as a function of $s_i(0:t)$, is still nonlinear in general (Pfister et al., 2010). Indeed, Equation 14 including these local nonlinearties was used to compute the linear response in Figure 1C–D.

In the general case inference can be performed using filtering such that in each step we update the inferred state of the hidden variables, $z_t$ and ${\mathbf{u}}_t$, using information from two different sources: the likelihood of emitting a particular spiking pattern (observation) and the dynamics of the hidden variables combined with the previous estimate (innovation):

where the likelihood $\text{P}\left({\mathbf{s}}_t\right|{\mathbf{u}}_t)$ is defined by Equation 8, the dynamics of the global state variable $\text{P}\left(z_t\right|z_{t-\delta t})$ is first order, Markovian (see above) and the state-dependent membrane potential dynamics $\text{P}\left({\mathbf{u}}_t\right|{\mathbf{u}}_{t-\delta t},z_t)$ is given by Equation 7. Equation 15 thus defines a mapping between the posterior distribution of the hidden variables in the previous time step (last term on RHS) and their current distribution (LHS). The posterior over membrane potentials can then be obtained by simply marginalising out the state variable:

For the following, it is useful to represent the posterior as a product of two terms:

As the state variable is binary, its posterior is a Bernoulli distribution which we parametrize by $\zeta$, without loss of generality:

However, in general, the posterior of the membrane potentials conditioned on the current state, $z_t$, can be arbitrarily complex. To allow an analytical reduction of the inference process, we adopted an assumed density filtering approach in which this distribution is moment-matched in each time step by a multivariate normal distribution which is thus described by two (sets of) parameters, its mean, ${\mathbf{\mu }}_t^{\left(z\right)}$, and covariance, ${\mathbf{\text{Σ}}}_{\mathbf{t}}^{\mathbf{\left(}\mathbf{z}\mathbf{\right)}}$:

with the analogous equation for the posterior of ${\mathbf{u}}_t$ conditioned on $z_t$ being in the $-$ state.

One advantage of this parametric approach is that inference (filtering) can be implemented by updating only the parameters describing the (approximate) posterior distribution (Equations 18–19): ${\zeta }_t$, ${\mathbf{\mu }}_t^{\left(z\right)}$, and ${\mathbf{\text{Σ}}}_t^{\left(z\right)}$. In the Appendix we derive an analytical form for these parameter updates resulting in the following set of differential equations:

where ${\mathbf{\gamma }}^{\left(z\right)}$ (${\mathbf{\Gamma }}^{\left(z\right)}$) is a state-dependent vector (diagonal matrix) of which the elements ${\gamma }_i^{\left(z\right)}={\mathrm{\Gamma }}_{ii}^{\left(z\right)}=g{e}^{\beta {\mu }_i^{\left(z\right)}+\frac{1}{2}{\beta }^2{\mathrm{\Sigma }}_{ii}^{\left(z\right)}}$ are the expected firing rates of the neurons in a given state, ${\overline{\gamma }}^{\left(z\right)}={\sum }_i{\gamma }_i^{\left(z\right)}$ is the expected total population firing rate in state $z$, and $⟨\mathbf{\Gamma }⟩=\zeta {\mathbf{\Gamma }}^{+}+(1-\zeta ){\mathbf{\Gamma }}^{-}$ is the expected firing rate of the cells averaged across states. In these equations, the spike trains of the presynaptic neurons are represented by the sum of Dirac-delta functions in continuous time and are denoted by $\mathbf{s}\left(t\right)$, to be distinguished from its discrete time analog, ${\mathbf{s}}_t$, such that $\mathbf{s}\left(t\right)=\underset{\delta t\to 0}{\lim }{\mathbf{s}}_t/\delta t$ (see Equation A62 in the Appendix). The differential equations for the conditional mean and variance in the $-$ state, ${\dot{\mathbf{\mu }}}^{-}$ and ${\dot{\mathbf{\Sigma }}}^{-}$, are analogous to Equations 21–22. The absolute refractory period is taken into account by setting ${\gamma }_i={\mathrm{\Gamma }}_{ii}=0$ after each observed spike for the duration of the refractory period, ${\tau }_{\text{refr}}$, thus omitting the effect of the likelihood (terms containing ${\mathbf{\gamma }}^{\left(z\right)}$ or ${\mathbf{\Gamma }}^{\left(z\right)}$) from Equations 20–22.

The first term in Equation 20 captures the decay in ${\zeta }_t$ that is proportional to the difference in the state conditional firing rates in the absence of presynaptic spikes; the second term expresses the instantaneous change in ${\zeta }_t$ after observing a spike, proportional to both the state estimation uncertainty, $\zeta \left(1-\zeta \right)$, and the differences in the conditional firing rates (${\mathbf{\gamma }}^{+}-{\mathbf{\gamma }}^{-}$); and the last term captures the decay of ${\zeta }_t$ to its steady state in the absence of observations. The filtering equations for the conditional mean and covariance (Equations 21–22) are each composed of three terms: the first term expresses the decay of the variable towards its baseline in the absence of observations; the second term captures the effect of the current observation (i.e. the presence or absence of a spike) on the variable; and the third term describes the changes in the variable caused by potential state transitions. This can be viewed an extension and generalization of earlier work deriving the equivalents of Equations 21–22 for the special case of a single neuron without state-switching dynamics (Pfister et al., 2010).

Another advantage of the parametrization of the posterior we chose is that computing the optimal response, i.e. the posterior expectation of the simple linear functions that we consider in this paper, becomes straightforward (cf. Equation 6):

In order to verify the assumed density filtering approximations used above we numerically integrated the system of differential equations (Equations 20–22) using the software package R (R Core Team, 2012; Soetaert et al., 2010) and compared the results with those obtained using standard particle filters (Doucet et al., 2001). In these simulations we used 500,000 particles to evaluate Equation 15 with $1$ neuron and 2 states. Figure 2—figure supplement 1C shows that the results of assumed density filtering are essentially identical to those of particle filtering, confirming that the approximations we used were valid.

Here we first describe a simple canonical model of dendritic integration following Poirazi and Mel (Poirazi and Mel, 2001), and then show that it provides an approximation to the optimal response (Equations 20–23) in the limiting case in which presynaptic dynamics are dominated by simultaneous switching between a quiescent and an active state. In this simple dendritic model, inputs within a branch are integrated linearly:

where $v_{\mathrm{lin}}$ is the variable linearly integrating inputs with weight $\mathcal{ℬ}$, dendritic time constant $1/\mathcal{𝒜}$ and steady state value $-\mathcal{𝒞}/\mathcal{𝒜}$ (in the absence of spikes), and $s\left(t\right)={\sum }_i{s}_i\left(t\right)$ denotes the spike train of presynaptic neurons, collecting all spikes from the presynaptic population. The actual dendritic response, $v_{\mathrm{den}}$, is then given by mapping this linear response through a sigmoidal nonlinearity, scaled to be between $v_{\min }$ and $v_{\max }$ (Figure 3A, inset):

To demonstrate that this reduced model of dendritic integration closely approximates the optimal response, we first note that under appropriate conditions (${\tau }_{\text{post}}$ is small, $N$ is large, $\mathbf{Q}$ is diagonal, and $\beta$ is small relative to the diagonal elements of $\mathbf{Q}$) the dynamics of the optimal response are dominated by the state switching process (Equation 20; see Appendix). Thus, the optimal response essentially follows the inference about the global state variable, $\zeta$, up to linear rescaling and filtering:

with ${\overline{u}}^{+}$ and ${\overline{u}}^{-}$ respectively denoting the resting membrane potential in the active and quiescent states. As Equation 26 is linear, all nonlinear interactions, corresponding to dendritic nonlinearities, must be contained in the temporal dynamics of the posterior probability of this global state variable being in the active state, $\zeta \left(t\right)$, which can be expressed as

where constants $B$ and $C$ depend on the parameters of the presynaptic statistics (see Appendix). Note that the fact that $\zeta \left(1-\zeta \right)$ multiplies Equation 27 expresses the simple intuition that the size of the update to $\zeta$ (the posterior probability of $z=+$) in response to incoming information (presence, $B$ term, or absence of a spike, $C$ term) should be proportional to our current (posterior) uncertainty about $z$; and since the posterior is a Bernoulli distribution, the uncertainty associated with it is simply $\zeta \left(1-\zeta \right)$.

The solution of Equation 27 can be expressed in a form that is similar, albeit not identical (see below), to the canonical model for dendritic integration (Equations 24–25). This form requires the linear integration of incoming spikes

and the temporal evolution of $\zeta$ is expressed as a sigmoidal function of the linearly integrated inputs $\nu$:

Thus, dendrites with sigmoidal nonlinearity are near-optimal when their synaptic inputs switch between a quiescent and an active state.

The main difference between the dendritic integrator and the optimal response is that the dynamics of spike integration imply exponential decay towards a finite baseline in the former (there is a negative term in Equation 24 which is scaled by $v_{\mathrm{lin}}$) and steady decrease towards negative infinity in the latter (the only negative term in Equation 28 is a constant, independent of $\nu$). This is because the approximations we used for deriving Equation 27 were accurate only in the quasi-static case, when the state switching dynamics are infinitely slow. In this case, remote and more recent observations should have identical effects on the current value of $\zeta$ as they all correspond to the same underlying state. In the more general case, when state switching occurs with non-zero probability, more remote observations likely correspond to a state which has changed in the meantime, and should thus count for less, such that their effect on the current value of $\zeta$ should decay with time – leading to leaky integration of incoming spikes, similar to that in Equation 24.

To compare quantitatively the response of the linear and nonlinear dendrites to the optimal response in a computational task using realistic input statistics, we divided the presynaptic population into four groups (cell assemblies), where neurons within each group were statistically dependent (either through simple second-order correlations, Figure 4—figure supplement 2, or through sharing a common state variable, Figure 4 and Figure 4—figure supplement 2) while neurons from different groups were independent (Figure 4A). In this task we used 4 different versions of the simplified neuron model (Figure 4B):

• The linear model responded to all incoming spikes with an identical postsynaptic potential (PSPs) characterized by its amplitude ($w_{\mathrm{\ell }}$), decay time constant (${\tau }_{\mathrm{\ell }}$) and the resting potential (${\overline{v}}_{\mathrm{\ell }}$): 

where $s\left(t\right)$ is the total incoming spike train (as before). This model had three parameters.

• The model with a somatic nonlinearity had linear dendrites but a nonlinear soma. Motivated by our analytical calculations we used a sigmoidal nonlinearity: 

where we computed ${\stackrel{\mathit{~}}{\mathit{v}}}_{\mathrm{\ell }}\left(t\right)$ as defined in Equation 30 above with $w_{\mathrm{\ell }}=1$ and ${\overline{v}}_{\mathrm{\ell }}=0$ (as these parameters were interchangeable with ${\beta }_{\mathrm{g}}$ and ${\theta }_{\mathrm{g}}$). This model had five free parameters.

• The random dendrites model had four nonlinear dendritic subunits, each receiving inputs from a unique set of 10 neurons randomly selected from the four presynaptic assemblies, and integrating these inputs using a sigmoidal nonlinearity (Equations 30-31). Each subunit had its own set of five parameters, and the outputs of the subunits were simply averaged in the soma (without loss of generality), resulting in 20 parameters in total.

• The clustered dendrites model was similar to the random dendrites model, with the important difference that neurons in each cell assembly selectively targeted a single nonlinear dendritic branch. Since the presynaptic statistics was the same for all four branches, we constrained the parameters of the four dendritic subunits to be identical and this model had only five parameters (the same performance was achieved when we relaxed this constraint and optimised all 20 parameters).

To fit the models we generated 240 s-long samples of presynaptic activity and optimized the parameters of the models to minimize the squared error between the signal ($v$, the true average of the stimulated presynaptic potentials, cf. Equation 5 with $w_i=\frac{1}{N}$) and their estimates ($\widehat{v}$, the outputs of the models), averaged over the duration of the sample:

After training, we tested the different models in cross-validation, on a novel 120 s-long input sequence, and quantified their performance by the fraction of variance unexplained, i.e. the temporally averaged squared error, ${ϵ}_{\mathrm{estimation}}$, normalized by the variance of the signal (as a sensible upper limit on the error – achievable by an estimator that predicts the prior mean, ignoring incoming spikes altogether):

where $\mathrm{Var}\left[v\right]=\frac{1}{T-1}{\sum }_t{(v_t-\text{E}\left[v\right])}^2$ and $\text{E}\left[v\right]=\frac{1}{T}{\sum }_t{v}_t$. (Normalization was unnecessary during training because the parameters of the models that we were optimizing obviously did not influence the variance of the signal.) Figure 4 shows $1-{\overline{ϵ}}_{\mathrm{estimation}}$, i.e. the fraction of variance explained as ‘performance’, and Figure 1—figure supplement 3, Figure 3—figure supplement 2, and Figure 4—figure supplement 2 show ${\overline{ϵ}}_{\mathrm{estimation}}$ as ‘estimation error’.

To predict dendritic integration in hippocampal and neocortical neurons we fitted the parameters describing presynaptic statistics in our model, $\text{P}\left(\mathbf{u}\right)$ and $\text{P}\left(\mathbf{s}\right|\mathbf{u})$, to the statistical patterns in the activity of their respective presynaptic populations.

The basal dendrites of neocortical layer 2/3 pyramidal cells are targeted by neighbouring pyramidal neurons as well as by neurons from layer 4 (Douglas and Martin, 2004). We used in vivo intracellular paired recordings from layer 2/3 pyramidal neurons in the barrel cortex (Poulet and Petersen, 2008; Gentet et al., 2010; Crochet et al., 2011) to set the parameters of our model to reproduce the salient features of the presynaptic population dynamics during quiet wakefulness (Tables 1–2).

In the hippocampal experiments we stimulated synapses on the proximal dendrites of CA3 neurons targeted by recurrent collaterals of neighbouring pyramidal cells (Andersen et al., 2007). We fitted the presynaptic statistics to in vivo population activity patterns recorded from the hippocampus during quiet wakefulness, characterised by sharp wave (SPW) activity (Csicsvari et al., 2000). As intracellular recordings from CA3 pyramidal neurons during SPW activity in the awake animal are not available, we fitted the presynaptic statistics to awake extracellular data from CA3 (Csicsvari et al., 2000; Grosmark et al., 2012) and intracellular (Ylinen et al., 1995; English et al., 2014) data from CA1 pyramidal neurons (Tables 1 and 3).

We used four different parameter sets (models) to describe the activity of the presynaptic population (Figure 1). The parameters of the HP and NC models were fitted to in vivo recordings from the corresponding presynaptic populations as described above. As a control, we used two simpler models with no state switching dynamics. The cor2 model had correlated membrane potential fluctuations with all cross-correlations between presynaptic neurons being the same, $-\frac{1}{N-1}\le \rho \le 1$. Neurons in the last model, ind, had independent membrane potential fluctuations. The HP, NC, and ind models had no free parameters, while parameter $\rho$ of the cor2 model was left free and later tuned to fit dendritic data. Note that to fit supralinear cortical responses, $\rho$ had to be tuned to unnaturally large negative values in this model (Figure 5—figure supplement 1A) – and it still produced significantly poorer fits than the HP and NC models (Figure 5).

After setting the parameters of the presynaptic population, we computed the optimal response by numerically integrating Equations 20–22 in the software package R (R Core Team, 2012; Soetaert et al., 2010). When comparing the optimal response to experimental data, we assumed that each uncaging event corresponded to a single spike at the presynaptic axon terminal, and spines not showing measurable gluEPSP were considered to be non-stimulated. For all four presynaptic parameter sets, we varied two postsynaptic parameters to fit the responses of the optimal estimator to dendritic integration data, i.e. the somatic membrane potential traces recorded in our in vitro experiments. These two parameters were the weight $w$ of the presynaptic neurons and the time constant of the postsynaptic filtering, ${\tau }_{\text{post}}$ (Equation 5). To avoid overfitting, we assumed that all presynaptic neurons had equal weight, i.e. $\text{∀}i w_i=w$. A final free parameter that we had to consider was the number of synapses that were in the same functional cluster as the synapses we stimulated in our experiments – where the term ‘functional cluster’ refers to a set of synapses for which the presynaptic cells are correlated. This parameter was irrelevant for the ind model (by definition), it was fixed at 20 for the cor2 model (because its effects on the optimal response were largely indistinguishable from that of varying $\rho$, see above), and it was tuned to fit dendritic integration for the state-switching models (NC and HP). In sum, the number of free parameters used to fit dendritic integration data was 2 for the ind model and 3 for the cor2, NC, and HP models. We confirmed that the higher number of free parameters in the latter models did not result in an unfair advantage in fitting performance by using Bayesian Information Criterion (BIC), rather than squared error (see below), as our measure of performance. BIC includes an explicit term penalizing the number of parameters, and our results were not qualitatively affected by it: fitting the model using the relevant in vivo statistics resulted in 3500$\pm$940 (NC, mean$\pm$s.d.) and 2000$\pm$1100 (HP) higher BIC scores than when using independent statistics (where each unit of BIC difference corresponds to a likelihood that is higher by a factor of $e\simeq 2.71$).

We fitted each recorded neuron independently using these parameters by minimizing the mean squared error between the predicted, $\tilde{v}$, and the recorded postsynaptic membrane potential (averaged across repetitions of the same stimulation in the same cell), ${\overline{v}}^{*}$:

To be able to compare results across different neurons and different stimulation protocols, we normalized the error by the total variance of the data:

where $\mathrm{Var}\left[{\overline{v}}^{*}\right]=\frac{1}{T-1}{\sum }_t{({\overline{v}}_t^{*}-\text{E}\left[{\overline{v}}^{*}\right])}^2$ and $\text{E}\left[{\overline{v}}^{*}\right]=\frac{1}{T}{\sum }_t{\overline{v}}_t^{*}$. A natural lower bound of our fitting error was the intrinsic variability of the data, so we computed the mean of the variance of the experimental data across repetitions, normalized by the total variance of the data:

where $v_{t,l}^{*}$ is the raw data before averaging across repetitions, $L$ is the number of repetitions using the same stimulation protocol in the same cell and ${\overline{v}}_t^{*}=\frac{1}{L}{\sum }_L{v}_{t,l}^{*}$. Figures 5D,H and Figure 5—figure supplement 1G show ${\overline{ϵ}}_{\mathrm{fitting}}$ as the ‘fitting error’ and ${\overline{ϵ}}_{\min }$ as ‘var’. The best fitting parameter values for the postsynaptic time constant, ${\tau }_{\text{post}}$, and total number of neurons in a functional cluster, $N$, are shown in Figure 5—figure supplement 1B,C.

Here we derive the filtering equations for the switching multivariate OU (mOU) process introduced in the main text, such that we generalise the two-state model presented there to an arbitrary number of hidden states.