Spike-timing prediction in cortical neurons with active dendrites

A complete single-neuron model must correctly reproduce the firing of spikes and bursts. We present a study of a simplified model of deep pyramidal cells of the cortex with active dendrites. We hypothesized that we can model the soma and its apical dendrite with only two compartments, without significant loss in the accuracy of spike-timing predictions. The model is based on experimentally measurable impulse-response functions, which transfer the effect of current injected in one compartment to current reaching the other. Each compartment was modeled with a pair of non-linear differential equations and a small number of parameters that approximate the Hodgkin-and-Huxley equations. The predictive power of this model was tested on electrophysiological experiments where noisy current was injected in both the soma and the apical dendrite simultaneously. We conclude that a simple two-compartment model can predict spike times of pyramidal cells stimulated in the soma and dendrites simultaneously. Our results support that regenerating activity in the apical dendritic is required to properly account for the dynamics of layer 5 pyramidal cells under in-vivo-like conditions.


INTRODUCTION
Partially neglected for a long time, dendrites have been recently shown to treat synaptic input in a surprising variety of modes (Stuart et al., 2007). Indeed, experiments have revealed that dendrites are excitable and that they can generate either sodium (Golding and Spruston, 1998), NMDA (Schiller et al., 2000) or calcium (Llinas and Sugimori, 1980) spikes. One particularly striking example is found in pyramidal cells of deep cortical layers. In these cells, a coincidence between a back-propagating action potential and dendritic input can trigger voltage-sensitive ion channels situated on the apical dendrite more than 300 μm from the soma (Larkum et al., 1999(Larkum et al., , 2001. The somatic membrane potential increases only after the activation of dendritic ion channels. This often resulting in a burst of action potentials. Bursts in these cells can therefore signal a coincidence of input from the soma (down) with inputs in the apical dendrites (top). Such top-down coincidence detection is one computation that is attributed to dendritic processes. Other allegedly dendritic computations include subtraction (Gabbiani et al., 2002), direction selectivity (Taylor et al., 2000), temporal sequence discrimination (Branco et al., 2010), binocular disparity (Archie and Mel, 2000), gain modulation (Larkum et al., 2004) and self-organization of neuron networks (Legenstein and Maass, 2011). These computations rely on the dendrite acting as an excitable subunit (Polsky et al., 2004;Stuart et al., 2007).
Models of large pyramidal neurons with active apical dendrites were first described by Traub et al. (1991) for the hippocampus. This model of the large CA3 pyramidal neurons included voltage-dependent conductances on the dendrites. It is a model based on the Hodgkin-Huxley description of ion channels. Cable properties of dendrites are taken into account by segmenting the dendrite into smaller compartments. The resulting set of equations is solved numerically. A simplified version of this model was advanced by Pinsky and Rinzel (1994). They have reduced the model to a dendritic compartment and a somatic compartment connected by an effective conductance. The model has a restricted set of five ion channels and accounts for bursting of CA3 pyramidal cells.
Models specific to deep cortical cells have been described by extending the approach of Traub et al. (1991);Schaefer et al. (2003) used morphological reconstruction to define compartments. This model could reproduce the top-down coincidence detection.
Using a simplified approach similar to Pinsky and Rinzel (1994), Larkum et al. (2004) have modeled dendrite-based gain modulation. The parameters in the model could be tuned to quantitatively reproduce the firing rate response of layer 5 pyramidal cells stimulated at the soma and the dendrites simultaneously. Larkum et al. (2004) concluded that a two-compartment model was sufficient to explain the time-averaged firing rate.
A more stringent requirement for neuron model validation, however, is to predict spike times (Keat et al., 2001;Pillow et al., 2005;Jolivet et al., 2006Jolivet et al., , 2008aGerstner and Naud, 2009). Given the low spike-time reliability of pyramidal neurons, spike time prediction is compared to the intrinsic reliability (Jolivet et al., 2006). This approach can be seen as predicting the instantaneous firing rate (Naud et al., 2011). Generalized integrate-and-fire models can predict instantaneous firing rate of layer 5 pyramidal neurons with substantial precision (Jolivet et al., 2008a;Gerstner and Naud, 2009;Naud et al., 2009) in the absence of dendritic stimulation. The question remains whether a neuron model can predict the spike times of layer 5 pyramidal neurons when both the dendrites and the soma are stimulated simultaneously.
We present a study of a simplified model of layer 5 pyramidal cells of the cortex with dendrites excitable with calcium spikes (Larkum et al., 2004(Larkum et al., , 2009). Following Larkum et al. (2004), we hypothesized that we can model the soma and its apical dendrite with two compartments, without significant loss in the accuracy of spike-timing predictions. We introduce experimentally measurable impulse-response functions (Segev et al., 1995), which transfer the effect of current injected in one compartment to current reaching the other. The impulse-response functions replace the instantaneous connection used in previous two-compartment models (Pinsky and Rinzel, 1994;Larkum et al., 2004) and acts as a third, passive, compartment. Each compartment was modeled with a pair of non-linear differential equations with a small number of parameters that approximate the Hodgkin-and-Huxley equations. The predictive power of this model was tested on electrophysiological experiments where noisy current was injected in both the soma and the apical dendrite simultaneously (Larkum et al., 2004).

METHODS
Methods are separated in four parts. First we present the model, second the experimental protocol, then fitting methods and finally the analysis methods. Figure 1 shows a schematic representation of the twocompartment model. In details, the model follows the system of differential equations: where I s is the current injected in the soma, I d the current injected in the dendrites, V s is the somatic voltage, V d is the dendritic voltage, m is the level of activation of a putative calcium current (I Ca = g 1m ), x is the level of activation of a putative calcium-activated potassium current (I K(Ca) = g 2x ), V T is the dynamic threshold for firing somatic spikes, I A is a spiketriggered current mediating adaptation, I BAP is the the current associated with the back-propagating action potential, sd is the filter relating the current injected in the soma to the current arriving in the dendrite and ds is the filter relating the current injected in the dendrite to the current arriving in the soma.

DESCRIPTION OF THE MODEL
The parameters are listed in Table 1.
As a control, we also consider an entirely passive model of dendritic integration. In this model, the current injected in the dendrite is filtered passively to reach the soma. The generalized passive model has an instantaneous firing rate: where λ 0 is a constant related to the reversal potential, κ s somatic membrane filter, κ ds is the filter relating the current injected in the dendrite to the voltage change in the soma, and η A is the effective spike-triggered adaptation.

EXPERIMENTAL PROTOCOL
Animal handling was in strict accordance with the guidelines given by the veterinary office of the canton Bern-Switzerland. Parasagittal brain slices of the somato-sensory cortex (300-350 m thick ) were prepared from 28-35 day-old Wistar rats. Slices were cut in ice-cold extracellular solution (ACSF), incubated at 34 • C for 20 min and stored at room temperature. During experiments, slices were superfused with in ACSF at 34 • C. The ACSF contained (in mM) 125 NaCl, 25 NaHCO3, 25 Glucose, 3 KCl, 1.25 NaH2PO4, 2 CaCl2 , 1 MgCl2 , pH 7.4, and was continuously bubbled with 5% CO2/95% O2. The intracellular solution contained (in mM) 115 K+-gluconate, 20 KCl, 2 Mg-ATP, 2 Na2-ATP, 10 Na2-phosphocreatine, 0.3 GTP, 10 HEPES, 0.1, 0.01 Alexa 594 and biocytin (0.2%), pH 7.2. Recording electrodes were pulled from thick-walled (0.25 mm) borosilicate gla-ss capillaries and used without further modification (pipette tip resistance 5-10 M for soma and 20-30 M for dendrites). Whole-cell voltage recordings were performed at the soma of a layer V pyramidal cell. After opening of the cellular membrane a fluorescent dye, Alexa 594 could diffuse in the entire neuron allowing to perform patch clamp recordings on the apical dendrite 600-700 μm from the soma. Both recordings were obtained using Axoclamp Dagan BVC-700A amplifiers (Dagan Corporation). Data was acquired with an ITC-16 board (Instrutech) at 10 kHz driven by routines written in the Igor software (Wavemetrics).
The injection waveform consisted of 6 blocks of 12 s. Each block is made of three parts: (1) one second of low-variance colored noise injected only in the soma, (2) one second of lowvariance colored noise injected only in the dendritic injection site, (3) ten seconds of high-variance colored noise whose injection site depends on the block: In the first block, the 10-s stimulus is injected only in the dendritic site, the second block delivers the 10-s stimulus in the soma only, and the four remaining blocks deliver simultaneous injections in the soma and the dendrites.
The colored noise was simulated with MATLAB as an Ornstein-Uhlenbeck process with a correlation time of 3 ms. The six blocks make a 72 s stimulus that was injected repeatedly without redrawing the colored noise (frozen-noise). Noise is frozen across repetitions to estimate intrinsic reliability, but not across blocks to ensure independent test set and training set. Twenty repetitions of the 72-s stimulus were carried out, separated by periods of 2-120 s. Out of the twenty repetitions, a set of seven successive repetitions were selected on the basis of high intrinsic reliability.

FITTING METHODS
Each kernel (κ s , κ ds , η A , ds , sd , I A , I BAP ) is expressed as a linear combination of non-linear basis (i.e., κ s (t) = i a i f i (t)). The rectangular function was chosen as the non-linear basis. The parameters weighting the contributions of the different rectangular functions are then linear in the derivative of the membrane potential for the two-compartment model and generalized linear for the passive model.
For the two-compartment model, we use a combination of regression methods and exhaustive search to maximize the mean square-error of the voltage derivative. The regression methods are similar to those previously used for estimating parameters with intracellular recordings. These methods are described in more details in Jolivet et al. (2006); Paninski et al. (2005); Mensi et al. (2012); Pozzorini et al. (2013). First, we distinguish two types of parameters, the parameters that can be expressed as a linear function of the observables and the parameters that cannot. For instance, the parameter g s is linear in the observable dV s /dt (Equation 1). Similarly, the amplitudes a i defining the filters are also linear parameters. There is a total of four nonlinear parameters in the two-compartment model, namely τ m , The fit of the somatic compartment essentially follows (Jolivet et al., 2006) but using multi-linear regression to fit the linear parameters. The fit of the dendritic compartment needs to iterate through the restricted set of non-linear parameters. All fits are performed only on the part of the data restricted for training the model. Each step in the fitting procedure uses the entire training set. 2a: Compute the first-order estimate of dV s /dt. 2b: Find the best estimates of the somatic parameters linear in dV s /dt given a set of non-linear parameters (D T , τ T , E T ). The best estimates are chosen through linear regression to minimize the mean square error of dV s /dt. 2c: Compute iteratively step 2b on a grid of the non-linear parameters and simulate the model with each set of nonlinear parameters in order to compute the coincidence rate (see Section 2.4). 2d: Take the parameters that yield the maximum coincidence factor.
For the generalized linear model, we use maximum likelihood methods (Paninski, 2004;Pillow et al., 2005). Expressing the kernels as a linear combination of rectangular bases we recover the generalized linear model. Here the link-function is exponential so that the likelihood is convex. We therefore performed a gradient ascent of the likelihood to arrive at the optimal parameters.

ANALYSIS METHODS
When one focuses on spike timing, one may want to apply methods that compare spike trains in terms of a spike-train metric (Victor and Purpura, 1996) or the coincidence rate (Kistler et al., 1997). Both measures can be used to compare a recorded spike train with a model spike train. A model which achieve an optimal match in terms of spike-train metrics will automatically account for global features of the spike train such as the interspike interval distribution.
Here we used the averaged coincidence rate (Kistler et al., 1997). The coincidence rate, like most other spike time metrics, can be related to the coefficient of correlation between the instantaneous firing rate of the model and the neuron (Naud et al., 2011). It can be seen as a similarity measure between pairs of spike trains, averaged on all possible pairs. To compute the pairwise coincidence rate, one first finds the number of spikes from the model that fall within an interval of = 4 ms after or before a spike from the real neuron. This is called the number of coincident events N nm between neuron repetition n and model repetition m. The coincidence rate is the ratio of the number of coincident events over the averaged number of events 0.5(N n + N m ), where N n is the number of spikes in the neuron spike train and N m is the number of spikes in the model spike train. This ratio is then scaled by the number of chance coincidences N Poisson = 2 N m N n /T. This formula comes from The pairwise coincidence rate nm is then averaged across all possible pairings of spike trains (trials) generated from the model with those from the neuron and gives the averaged coincidence rate . Averaging across all possible pairings of spike trains from the neuron with a distinct repetition of the same stimulus given to the same neuron gives the intrinsic reliability R.

RESULTS
Dual patch-clamp recordings were performed in L5 Pyramidal cells of Wistar rats (see Experimental Methods). A simplified two-compartment model (see Model Description) was fitted on the first 36 s of stimulation for all repetitions. The rest of the data (36 s) was reserved to evaluate the model's predictive power. The predictive power of the two-compartment model with active dendrites was then compared to a model without activity in the dendrites (see Section 2.1), the generalized linear passive model. Figure 2 summarizes the predictive power of the twocompartment model. The somatic and dendritic voltage traces are well captured (Figures 2A-D). The main cause for erroneous prediction of the somatic voltage trace is extra or missed spikes (Figures 2A,B lower panels). The dendritic voltage trace of the model follows the recorded trace both in a low dendritic-input regime ( Figure 2C) and in a high dendritic-input regime with dendritic "spikes" (Figure 2D). The greater spread of voltageprediction-error (Figure 2) is mainly explained by the larger range of voltages in the dendrites (somatic voltage prediction is strictly subthreshold whereas dendritic voltage prediction ranges from −70 to + 40 mV). The interspike interval distribution is well predicted by the model (Figure 2G).
The generalized passive model does not predict as many spike times ( Figure 2H). The intrinsic variability in the test set was 68% and the two-compartment model predicted 50%. The prediction falls to 36% in the absence of a dendritic non-linearity ( Figure 2H).
The fitted kernels show that spike triggered adaptation is a monotonically decaying current that starts very strongly and decays slowly for at least 500 ms ( Figure 3A). The backpropagating action potential is mediated by a strong pulse of current lasting 2-3 ms ( Figure 3B). The coupling ds from dendrite to soma has a maximal response after 2-3 ms and then decays so as to be slightly negative after 35 ms ( Figure 3C). The coupling sd from soma to dendrite follows qualitatively ds with smaller amplitudes and slightly larger delays for the maximum and minimum peaks (Figure 3D), consistent with the larger membrane time-constant in the soma than in the dendrites.
The two-compartment model can reproduce qualitative features associated with the dendritic non-linearity in the apical dendrite of L5 pyramidal neurons. We study two of these features: the critical frequency (Larkum et al., 1999) and the gain modulation (Larkum et al., 2004). The first relates to the critical somatic firing frequency above which a non-linear response is seen in the soma, reflecting calcium channel activation in the dendrites. To simulate the original experiment, we force 5 spikes in the soma at different frequencies and plot the integral of the dendritic voltage. The critical frequency for initiating a non-linear increase in summed dendritic voltage is 138 Hz (Figure 4A). Pérez-Garci et al. (2006) reported a critical frequency of 105 Hz while (Larkum et al., 1999) reported 85 Hz. This appears to vary across different cells and pharmacological conditions. The model also appears to perform gain modulation as in Larkum et al. (2004) (Figure 4B). The relation between somatic firing rate and mean somatic current depends on the dendritic excitability. The firing threshold but also the gain (or slope) of the somatic frequency vs. somatic current curve depend on the mean dendritic current. The gain modulation is attributed to a greater presence of bursts ( Figure 4B) caused by dendritic calcium-current activation at higher dendritic input. The link between burst and dendritic activity is reflected in the burst-and spike-triggered average injected current (Figures 4C,D) similar to Larkum et al. (2004). The burst-triggered current is greater for the dendritic injection, whereas the spike-triggered current is larger for somatic injection. Therefore bursts signal a higher dendritic current that was concomittant with an increased somatic current. This can be interpreted as a top-down coincidence detection.

DISCUSSION
A dendrite is said active when it sustains either sodium, calcium or NMDA spikes. Our model reflects calcium spikes in the dendritic compartment, but not dendritic sodium spikes or NMDA spikes. The parameters fitted (Table 1) are in agreement with voltage-activated calcium channels. An activation sensitivity D m of 5 mV is typical of many ion channels, and the time constant τ x of about 50 ms is slower than the high-voltage activated calcium channel which has a time constant of about 10 ms (Gerstner et al., 2014). The current injection in the apical dendrite presumably does not solicit NMDA spikes known to occur in the apical tuft (Larkum et al., 2009).
Even if the spike-time prediction is high, the fitted parameters may differ from the real biophyical parameters for various reasons. First, the fitting method we used avoids local minima combining convex fitting procedures with exhaustive search. Even if the steps in the procedure are convex, the sequence of such steps may not be convex. Therefore the fitted parameters may reflect a local minimum. Also, the drop in coincidence rate

(A)
Dendritic non-linearity is triggered by somatic spiking above a critical frequency. Somatic spike-trains of 5 spikes are forced in the soma of the mathematical model at different firing frequencies. The normalized integral of the dendritic voltage is shown as a function of the somatic spiking frequency. (B) Dendritic injection modulates the slope of the somatic spiking-frequency vs. current curve. The slope of the frequency vs mean somatic current as measured between 5 and 50 Hz is plotted as a function of the mean dendritic current. Both somatic and dendritic currents injected are Ornstein-Uhlenbeck processes with a correlation time of 3 ms and a standard deviation of 300 pA. (C) Spike-triggered average of the current injected in the soma (black) and in the dendrites (blue). (D) Burst-triggered average of the current injected in the soma (black) and in the dendrites (blue). The fact that the blue curve is higher than the black curve, and that this relation is inverted in (C), may be interpreted as a top-down coincidence detection by bursts. between the training set and the test set indicate that overfitting is present. This could be avoided by using a smaller number of non-linear bases for the filters, or the use of raised cosine functions instead of the rectangular ones. Lastly, the filters and reduced model parameters may lump together different biophysical processes. The particular shape of the filter will depend on the average membrane potential and the average firing rate. This is one reason why we did not estimate the filters empirically with a separate set of experiments, but instead we fitted the model on current injection designed to imitate the natural condition.

CONCLUSION
Using a two-compartment model interconnected with temporal filters, we were able to predict a substantial fraction of spike times. The predicted spike trains achieved an averaged coincidence rate of 50%. The scaled coincidence rate obtained by dividing by the intrinsic reliability (Jolivet et al., 2008a;Naud and Gerstner, 2012) was 72%, which is comparable to the state-of-the performance for purely somatic current injection which reaches up to 76% ). Comparing with a passive model for dendritic current integration, we found that the predictive power decreased to a scaled coincidence rate of 53%. Therefore we conclude that regenerating activity in the apical dendrite is required to properly account for the dynamics of layer 5 pyramidal cells under in-vivo-like conditions.