Data-driven reduction of dendritic morphologies with preserved dendro-somatic responses

  1. Willem AM Wybo  Is a corresponding author
  2. Jakob Jordan
  3. Benjamin Ellenberger
  4. Ulisses Marti Mengual
  5. Thomas Nevian
  6. Walter Senn
  1. Department of Physiology, University of Bern, Switzerland


Dendrites shape information flow in neurons. Yet, there is little consensus on the level of spatial complexity at which they operate. Through carefully chosen parameter fits, solvable in the least-squares sense, we obtain accurate reduced compartmental models at any level of complexity. We show that (back-propagating) action potentials, Ca2+ spikes, and N-methyl-D-aspartate spikes can all be reproduced with few compartments. We also investigate whether afferent spatial connectivity motifs admit simplification by ablating targeted branches and grouping affected synapses onto the next proximal dendrite. We find that voltage in the remaining branches is reproduced if temporal conductance fluctuations stay below a limit that depends on the average difference in input resistance between the ablated branches and the next proximal dendrite. Furthermore, our methodology fits reduced models directly from experimental data, without requiring morphological reconstructions. We provide software that automatizes the simplification, eliminating a common hurdle toward including dendritic computations in network models.


Morphological neuron models have been instrumental in neuroscience (Segev and London, 2000). Major experimental discoveries, for instance that N-methyl-D-aspartate (NMDA) channels (MacDonald and Wojtowicz, 1982) produce local dendritic all or none responses (Schiller et al., 2000; Major et al., 2008), or that dendritic Ca2+ spikes mediate coincidence detection between distal inputs and somatic action potentials (APs) (Larkum et al., 1999), have been combined in morphological models to arrive at a consistent picture of dendritic integration: the dendrite is an intricate system of semi-independent subunits (Mel, 1993; Poirazi et al., 2003b; Poirazi et al., 2003a), amenable to dynamic regulation (Poleg-Polsky et al., 2018; Wybo et al., 2019), and able to distinguish specific input patterns (Branco and Häusser, 2010; Laudanski et al., 2014).

Nevertheless, morphological models are not without shortcomings. They are highly complex and consist of thousands of coupled compartments, each receiving multiple non-linear currents. The parameters of the models, typically fitted with evolutionary algorithms to electro-physiological recordings (Hay et al., 2011; Almog and Korngreen, 2014; Van Geit et al., 2016), number in the tens of thousands. Since recordings can only be obtained at a few dendritic sites, these fits are under-constrained, and thus susceptible to over-fitting. Accordingly, the single-neuron fitting challenge, where model performance was measured on unseen spike trains, was not won by a biophysical model, but by an abstract spiking model (Gerstner and Naud, 2009; DiLorenzo and Victor, 2013). Finally, many network-level observations can be explained without morphological models (Gerstner et al., 2012).

These shortcomings underline the need to find the essential computational repertoire of a neuron: the set of computations needed to understand brain functions such as learning and memory. Model simplification is crucial in this endeavor, as it elucidates the lowest level of complexity at which computational features are preserved. Conceptually, the simplification thus extracts the essential elements required for the computation from the underlying biophysics. Simulation-wise, the reduced model requires few resources and can thus be integrated in large-scale networks. Experimentally, the reduced model likely leads to a well-constrained fit.

Past simplification efforts can be grouped in two categories: approaches that use traditional compartments, but with adapted parameters (Pinsky and Rinzel, 1994; Destexhe, 2001; Tobin et al., 2006; Bahl et al., 2012; Marasco et al., 2013; Amsalem et al., 2020), and approaches that rely on more advanced mathematical techniques (Kellems et al., 2009; Kellems et al., 2010; Wybo et al., 2013; Wybo et al., 2015). Many authors apply ad hoc morphological simplifications and then adjust model parameters to conserve geometrical quantities, such as surface area (Davison et al., 2000; Hendrickson et al., 2011; Marasco et al., 2013), or electrical quantities, such as attenuation (Destexhe, 2001) and transfer impedance (Amsalem et al., 2020). Other authors propose two-compartment models whose parameters are adjusted to reproduce certain response properties (Pinsky and Rinzel, 1994; Naud et al., 2014). Nevertheless, these approaches often lack the flexibility to be adapted to a wide range of dendritic computations. Furthermore, particular afferent spatial connectivity motifs might be lost, so that essential computations are not captured by such a reduction. Hence, in a network model, the computational relevance of these motifs might unknowingly be ignored. Advanced mathematical techniques on the other hand are very flexible, as they can incorporate the response properties of the morphology implicitly. (Wybo et al., 2013 and Wybo et al., 2015 propose a system of convolutions whose kernels capture the passive properties of the compartmental model, whereas Kellems et al., 2009 and Kellems et al., 2010 linearly transform the compartmental model into a low-dimensional basis.) However, these techniques are not supported by standard simulation software, such as NEURON (Carnevale and Hines, 2004), a considerable hurdle toward their integration in the canonical neuroscience toolset.

Here, we introduce a simplification method with the flexibility of advanced mathematical techniques, but using traditional compartments. The approach accommodates any morphological model in public repositories (McDougal et al., 2017) and commonly used ion channels (Podlaski et al., 2017). The reduced models represent the optimal approximation to the dendritic resistance matrix, evaluated at any spatial resolution of choice. We fit ion channels by approximating the quasi-active resistance matrices (Koch and Poggio, 1985). All our fits are uniquely solvable linear least-squares problems. The obtained models extrapolate well to non-linear dynamics, and reproduce back-propagating APs (bAPs), Ca2+ spikes (Larkum et al., 1999), and NMDA spikes (Schiller et al., 2000; Major et al., 2008) with few compartments. Additionally, we investigate whether a dendritic tree with given afferent spatial connectivity motifs can be simplified by ablating subtrees or branches and grouping synapses in the next proximal compartment. We find that effective weight-rescale factors for synapses can be computed if temporal conductance fluctuations stay below a limit that depends on the difference in input resistance between the ablated branch and the next proximal dendrite. Under application of these factors, voltage responses in the simplified tree are preserved. Finally, we demonstrate that our approach can fit reduced models to experimental recordings, without the need to reconstruct full morphologies. We have created a Python toolbox (NEural Analysis Toolbox – NEAT – (copy archived at swh:1:rev:1cb15f36aa0a764105348541d046c85ef38e21ee) that implements this method together with an extension to NEURON to simulate the obtained models.


A systematic simplification of complex neuron morphologies

Our simplification strategy fits compartments with voltage-gated channels to a reduced set of dendritic locations of interest (Figure 1A, left-middle). The locations, together with the morphology, provide a tree structure for the reduced model (Figure 1A, middle-right) where each node corresponds to a traditional compartment and each edge to a coupling conductance. The fit requires that the reduced model responds similarly to perturbative current steps as the full neuron, at the set of chosen locations (Figure 1B). We consider perturbations around any spatially uniform holding potential vh. Experimentally, vh may be reached by injecting a constant current ih, and thus vh is the effective equilibrium potential under this current injection. In models, we do not need to know ih; we simply assume that vh is the equilibrium potential around which to linearize the model. Note that if ih=0, vh is the resting membrane potential.

Figure 1 with 1 supplement see all
Flexible and accurate reduction methodology.

(A) For any set of locations on a given morphology (left, here an L2/3 pyramidal cell [Branco and Häusser, 2010]), a reduced compartmental model can be derived (middle), with an associated schematic representation (right). (B) Steps of our approach: (1) choice of locations at which the reduced model should reproduce the full model’s voltage, (2) coupling, leak and channel conductances are fitted to resistance matrices derived from the full model at different holding potentials, and (3) capacitances are fitted to mimic the largest eigenmode of the full model. (C) The resistance matrix of the passive full model (top) restricted to the five locations in A is approximated accurately by the inverse of the conductance matrix of the passive reduced model (bottom). Labels correspond to locations in A. (D) Example components of the quasi-active resistance matrix of the full model, equipped with a Na+-channel, as a function of the holding potential vh. Red lines show the four holding potentials at which our methodology evaluates the resistance matrix. Singularities correspond to holding potentials where the linearization is invalid and should be avoided in the fit. Labels correspond to locations in A. (E) Temporal shape of exemplar input impedance kernels of the full model (gray) and their reduced counterparts (blue, dashed). (F) Same as in E, but for transfer impedance kernels.

Suppose we have chosen N locations. Let δ𝐢 be a vector describing perturbative input currents to each of those N locations. The linearized voltage response of the full neuron, at those locations, is given in the Fourier domain for an input of a particular frequency ω by

(1) δ𝐯(ω)=Zvh(ω)δ𝐢(ω),

with Zvh the quasi-active (Koch and Poggio, 1985) N×N impedance matrix. In experiments, Zvh is extracted from δ𝐯 measured in response to δ𝐢. In models, Zvh is algorithmically computed based on the full morphology, the parameters of the passive currents, and the dynamics of the voltage-gated ion channels (Koch and Poggio, 1985). For our reduction method, we consider perturbative current steps, so that after an equilibration period, we obtain a steady-state voltage response δ𝐯. Note that this corresponds to evaluating (Equation 1) at zero frequency, and that the corresponding matrix Zvh(ω=0) contains the input and transfer resistances. Throughout the manuscript, we maintain the notation that z without argument signifies z(ω=0), representing an individual input or transfer resistance, and that Z signifies Z(ω=0), representing the resistance matrix.

Linearizing the reduced compartmental model around a holding potential (see Materials and methods — The conductance matrix) yields an N×N conductance matrix Gvh, with on its diagonal

(2) gLi+n𝒩igC,in+dig¯dild(vh),

where gLi is the unknown leak conductance of the i-th compartment, 𝒩i the set of nearest neighbor compartments, and gC,in the unknown coupling conductance between compartment i and compartment n𝒩i. The second sum runs over the set i of all linearized currents of ion channels present in compartment i, where g¯di is the unknown maximal conductance of ion channel d and ld(vh) is a factor that follows from linearizing the dynamics of channel d around the holding potential vh (see Materials and methods – Quasi-active channels). Thus, the unknowns gLi, gC,in for n𝒩i, and g¯di for di are to be fitted for each compartment, and the factors ld(vh) are known and determined by the ion-channel dynamics. The ij’th off-diagonal element of Gvh is -gC,ij – the negative of the coupling conductance between compartments i and j – if i and j are nearest neighbor compartments on the reduced tree structure, and zero otherwise. The conductance matrix relates the reduced model’s voltage response to the perturbative input current steps

(3) Gvhδ𝐯=δ𝐢.

The full neuron and reduced model will behave similarly for all possible perturbative input steps δ𝐢 if their responses δ𝐯 match. From (Equation 1) and (Equation 3), it follows that Zvh should be the inverse of Gvh. Consequently, we require that multiplying the known Zvh (measured or calculated) by the parametric Gvh yields the identity

(4) ZvhGvhI.

From (Equation 2), it can be seen that (Equation 4) is linear in the parameters that have to be fitted (leak, coupling, and maximal ion-channel conductances of the reduced compartments). By consequence, (Equation 4) can be cast into a least-squares problem and solved accurately (Figure 1C, Figure 1—figure supplement 1A).

Since the linearized ion channel activation ld(vh) depends on the holding potential, Zvh changes with vh (Figure 1D). The fit must disentangle the changes in Zvh induced by the various channels. In models, we first block all voltage-gated ion channels in the full and reduced models and fit Gpas to Zpas according to (Equation 4), and thus obtain leak and coupling conductances for each compartment. Then, we unblock one ion channel at a time and decompose the conductance matrix into Gvh=Gpas+Gvh,chan, with Gvh,chan a diagonal matrix containing the conductances of the unblocked channel, linearized around vh, at each compartment. Thus Gvh,chan depends linearly on the unknown maximal ion channel conductance parameters. With Zvh and Gpas known, we optimize these maximal conductances, so that the left-hand side of

(5) ZvhGvh,chanIZvhGpas.

matches the right-hand side, and that it does so for multiple holding potentials (we chose vh=-75,-55,-35, and 15 mV, Figure 1D, Figure 1—figure supplement 1D). Thus, we obtain an overdetermined system of equations in the unknown maximal channel conductances at each compartment and compute its solution in the least mean squares sense.

Having fixed the steady-state behavior of the reduced model, we are left with the capacitance of each compartment of the reduced model to match the temporal dynamics of the full model as closely as possible. In our reductions, the capacitance of each compartment is thus a parameter to be fitted. To do so, we found it sufficient to consider passive membrane dynamics. Blocking all voltage-gated ion channels in the reduced model (Gvh,chan=0), temporal voltage responses follow

(6) diag(𝐜)𝐯˙=-Gpas𝐯+𝐢,

with 𝐜 a vector containing the unknown capacitance of each compartment. We require that voltage-decay back to rest, as described by (Equation 6), matches voltage decay in the full neuron with all active channels blocked, at the N chosen locations. According to linear dynamical systems theory (Strogatz, 2000), this decay can be decomposed as a sum of exponentially decaying eigenmodes, each with an associated eigenvalue αk and eigenvector ϕk. The former gives the exponential time scale τk=1/αk of the decay and the latter the spatial profile of the mode. The most prominent of these modes has the largest time scale τ0=1/α0, and primarily models the voltage decay back to rest through trans-membrane currents (Holmes et al., 1992). In experiments, we extract this mode by fitting an exponential to the voltage decay. In full models, we compute the eigenvalue α0 and eigenvector ϕ0 – restricted to the N chosen locations – with the separation-of-variables method (Major et al., 1993). In the reduced model, eigenmodes are found as the eigenvalues and eigenvectors of the matrix S=-diag(𝐜)-1Gpas. To fit the capacitances, we require that S has α0 as eigenvalue and ϕ0 as corresponding eigenvector:

(7) Sϕ0=α0ϕ0.

This system of equations is linear in the reciprocals of the capacitances entering in S and can thus be solved efficiently.

To accurately reproduce the spatio-temporal voltage responses of the full model, the reduced model must approximate the ‘impedance kernels’ z(t) (Wybo et al., 2015), the inverse Fourier transforms of the frequency-dependent elements of the impedance matrix. Despite only fitting the largest time scale 1/α0, we accurately reproduce kernels at all times (Figure 1E,F, Figure 1—figure supplement 1B).

Unless a large number of closely spaced compartment sites is retained on the morphology, our reductions replace the true axial currents flowing along the dendritic tree, which are shaped by the distribution of ion channels in the dendritic membrane, with a strong abstraction: a single coupling conductance between compartments that may be far apart. Hence, it is not guaranteed that the resting voltage of reduced model will be equal to the resting voltage of the full model, at the N chosen locations. We ensure that this will be the case by fitting the leak reversals for each compartment. The product of leak conductance and reversal represents a constant current for each compartment, and thus does not influence model dynamics beyond fixing it’s resting voltage. To fit the leak reversals, we evaluate all ion-channel currents in the reduced model at the full model’s resting voltage and require that temporal voltage derivatives are zero. The obtained equations are linear in the reduced model’s leak reversals, and hence can be solved efficiently. Note that our methodology does not guarantee that the obtained reversals are within physiological limits. Thus, they should purely be seen as fit parameters that determine the resting membrane potential.

Reduced models match the voltage response of their full counterparts

We demonstrated the reduction on two computations that require accurate spatio-temporal interaction. We reproduced sequence discrimination in an L2/3 pyramidal cell model (Branco and Häusser, 2010), where a neuron responds more strongly to a centripetal sequence of inputs than to a centrifugal one, by only retaining compartments on a single branch (Figure 2A). We also reproduced input-order detection (Torben-Nielsen and Stiefel, 2009), where a neuron responds more strongly when one input arrives before the other, by only retaining compartments at the soma and the two input sites (Figure 2B).

Voltage-match between full and reduced models for spatio-temporal dendritic computations.

(A) Reduction of full model (right) to a single branch (middle) reproduces sequence discrimination (right), full model in gray, and reduced model colored for different time-steps between inputs, centripetal (dashed), and centrifugal (dash-dotted). (B) Full model (left) and three-compartment reduction (middle, bottom) discriminate temporal order of inputs, where the response to inputs (middle, top) ordered 12 is stronger than 21. Voltage responses (right) in full model (gray) and reduced model (colored). (C, D) Reductions of resp. L5 pyramidal and Purkinje cells with active ion channels at the soma, and excitatory (AMPA+NMDA) and inhibitory (GABA) synapses at dendritic compartment sites. From left to right: full model with compartment sites (soma S and a selected dendrite site D are labeled), reduced model, somatic voltage with zoom on a single AP (full model in gray and reduced model in purple, input spikes at the bottom), dendritic voltage at D (full model in gray and reduced model in green, input spikes at the bottom), and the relative root mean square voltage errors at each compartment (root mean squared error normalized by the standard deviation RRMSE(v)=Avg(vfull-vreduced)2/σvfull). Spike coincidence factor Γ is also shown.

We then tested our approach with non-linear currents concentrated at discrete points (hot-spots) along the morphology. In two biophysical models – an L5 neocortical pyramid (Hay et al., 2011; Figure 2C) and a Purkinje cell (Chen et al., 2013; Miyasho et al., 2001; Figure 2D) – we removed active dendritic channels, but added their total conductance at rest to the passive leak (we term this the 'passified’ model), while retaining all active channels at the soma. The soma, with its AP-generating channels, naturally forms such a non-linear hot-spot. We implemented dendritic hot-spots by clustering sets of excitatory (α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid [AMPA] + NMDA) and inhibitory (γ-aminobutyric acid [GABA]) synapses at randomly selected points on the morphology. Both the reduced pyramidal cell model and the reduced Purkinje cell model accurately reproduced the somatic and dendritic membrane voltage traces of their full counterparts (Figure 2C,D, middle panels). Furthermore, 97% and 89% of APs coincided within a 6 ms window, respectively (Jolivet et al., 2008; Figure 2C,D, right panels).

Models with passive dendritic membranes but non-linear synapse clusters can already reproduce much of the canonical computational repertoire of dendrites. Clusters of AMPA+NMDA synapses at the distal tips of basal dendrites can make neurons function as two-layer networks (Mel, 1993; Poirazi et al., 2003b; Poirazi et al., 2003a). Such reduced models also exhibit the same robustness to input noise as the full models (Poleg-Polsky, 2019), and can include modulation of the co-operativity between AMPA+NMDA clusters through compartments with shunting conductances (Wybo et al., 2019). Furthermore, because our reduced models reproduce dendritic input impedance properties, long time-scale NMDA currents activate strongly in distal compartments as high distal input impedance helps to overcome their voltage-dependent Mg2+ block (Jahr and Stevens, 1990a; Jahr and Stevens, 1990b), whereas short time-scale AMPA currents dominate in more proximal compartments. Thus passive models with non-linear synapse clusters also capture the shift in computational strategy from a proximal temporal code to a distal rate code (Branco and Häusser, 2011).

What level of spatial complexity reproduces characteristic responses generated by somatic and dendritic ion channels? These responses arise through the interplay of integrative properties of the morphology – captured in the impedance kernels z(t) – and channel dynamics. If a full model (here the L2/3 pyramid) is reduced to a single somatic compartment (Figure 3A), the impedance kernel has only one time scale – the membrane time-scale (Figure 3B, top – purple). The kernel of the full model contains additional shorter time scales (Figure 3B, top – gray), leading to a faster AP onset than in the reduced model (Figure 3B, bottom). This effect is thus a fundamental limitation of the single-compartment reduction. Extending the compartmental model with four nearby dendritic sites (Figure 3A) adds fast components to the somatic impedance kernel (Figure 3B, top – green) that model the spread of charge to the other compartments. These components increase AP amplitude and decrease AP delay (Figure 3C), leading to a better match with the full model.

Figure 3 with 1 supplement see all
Dendritic computations with active channels are captured by our reduced models.

(A–C) Effect of compartment distribution on AP dynamics in reduced models. (A) One- and five-compartment reductions of the L2/3 pyramid, equipped with somatic and dendritic ion channels. (B) Differences in short time-scale behavior in somatic input impedance kernels (top) between full model (gray) and one- (purple) and five-compartment (green) reductions result in different AP delays (bottom). (C) AP amplitude (top) and AP delay (bottom) for the three models (colors as in B). (D–F) Effect of compartment distribution on AP back-propagation in basal dendrites. (D) Four- and two-compartment reductions of a basal branch. (E) APs at soma (top) and most distal compartment site (bottom) for four models (full in gray, four compartments in blue, two compartments in orange, and full but with a passive dendrite in green). (F) Amplitude (top) and delay (bottom) for bAPs at different distances from soma (if compartment is present in model), averaged over all basal branches longer than 150 μm (error bars indicate standard deviation). (G–I) AP back-propagation in the apical dendrite of the L5 pyramid. (G) Reductions of the apical dendrite with increasing inter-compartment spacing. (H) Voltage waveform at soma (top, full in gray, Δx=100 μm in blue, Δx=200 μm in orange) and two dendritic sites (middle, bottom). (I) Waveform amplitude as a function of distance to soma for various inter-compartment spacings. (J, K) Ca2+-spike-mediated coincidence detection. (J) Reduction of the L5-pyramid’s apical dendrite to 11 compartments. (K) Response to a somatic current pulse (top, left), a dendritic synaptic current waveform (top, right), and the coincident arrival of both inputs (bottom).

Dendritic ion channels support the back-propagation of APs (Stuart et al., 1997). For each basal branch of at least 150 μm in the L2/3 pyramid, we derived reduced models with three compartments (at 50, 100, and 150 μm from the soma) and with one compartment. We compared bAP amplitudes and delays relative to the somatic AP (Figure 3E) and found that even models with a single distal compartment support active back-propagation (Figure 3E,F). In the apical dendrite of the L5 pyramid (Figure 3G), we found that a distance step of 100 μm between compartments was required to support bAPs to the same degree as in the full model (Figure 3H,I).

Finally, we considered the pairing of a somatic current pulse with a dendritic post-synaptic potential waveform. We found that an 11-compartment reduction (Figure 3J) – with compartments spaced at 100 μm through the apical trunk to support bAPs until the Ca2+ hot-zone – reproduced the Ca2+-mediated coincidence detection mechanism (Larkum et al., 1999; Pérez-Garci et al., 2013).

Conditions under which afferent spatial connectivity motifs can be simplified

Our method allows us to reduce morphological complexity by removing any branch or subtree without affecting the integrative properties of other branches or subtrees. Up until now, we have only considered reductions where the to-be removed branch or subtree does not receive direct synaptic input. However, if this branch or subtree does receive synaptic input, the morphological reduction also simplifies the afferent spatial connectivity motifs targeting it, as synapses have to be grouped at the nearest proximal intact compartment. We assess whether the computational repertoire of the neuron remains intact under such a simplification, or which elements of the repertoire might be lost. We study two proxies for this repertoire: the voltage difference between full and reduced models at the intact compartments – in order to assess whether the output of local computations in the ablated branch or subtree can be recovered, and the difference between full and reduced models in NMDA-spike threshold at the nearest proximal intact compartment (quantified as the smallest number of AMPA+NMDA synapses with a weight of 1 nS that need to be activated simultaneously to elicit an NMDA spike) – in order to assess whether local integrative properties at the intact compartments can also be retained under the shift of synapses. We allow ourselves the liberty of rescaling the weights of the shifted synapses by a temporally constant factor β (Figure 4A). Which spatial synapse distributions admit simplification in this sense, and under which input conditions?

Simplification of afferent spatial connectivity motifs.

(A) Removal of a branch with a synapse (red triangle) is considered possible if the correct voltages at the compartment sites (here, dendritic – green square and somatic – purple square) can be obtained by shifting the synapse to the compartment site and rescaling its weight with a fixed factor β for the all input conditions under consideration. (B) Comparison between weight-rescale factors for current-based (left) and conductance-based (right) input. Voltage trace for dendritic compartment without rescaling (gray) and with the current-based scale factor βcurr compensating for attenuation (blue). Bottom right panel shows voltage trace with conductance-based scale factor βcond. (C) Current- and conductance-based scale factors in the green dendritic branch, for a shift of a synapse at a given distance from the soma to the dendritic compartment site (green square). (D) Spatial peak voltage without (left, blue) and with (right, orange) application of βcond. (E) For a cluster of AMPA+NMDA synapses in isolation, scale factors for physiological constants can be obtained (see Materials and methods — Synaptic weight-rescale factors) that reproduce the correct voltage waveform. Colors as in D. (F) Maximal amplitude of NMDA-spike waveform upon activation of increasing numbers of synapses – NMDA-spike thresholds indicated with vertical lines. Colors as in D. (G–H) Removing a whole subtree and shifting multiple synapses (red triangles) to the next proximal compartment site (green square). NMDA-spike generation (gray voltage trace) at the compartment site through burst inputs to local AMPA+NMDA synapses (green inputs), with AMPA (blue) and GABA (red) background inputs spread throughout the subtree. Reductions shown without rescaling (blue, dash-dotted), with the analytical single-site rescaling rule (red, full) and the numerical multi-site rule (orange, dashed). (I) Error in NMDA-spike threshold for the three cases in H. (J) Dependence of the error in NMDA-spike threshold on the factor Δzg, with Δz the average input resistance difference between synapse sites and the compartment and g the average synaptic conductance.

For a single, current-based synapse, shifted from its original location on the ablated branch to the next proximal compartment, we analytically compute the weight-rescale factor and find βcurr=zcs/zcc, with zcs the transfer resistance from synapse site s to compartment site c and zcc the input resistance at c. Since zcs<zcc (Koch, 1998), the weight-rescale factor weakens the shifted synapse, so that it matches the attenuation of its distal counterpart. Nevertheless, when s is located more distally on the same branch or subtree as c, the transfer resistance zcs is often close to the input resistance zcc (Figure 1—figure supplement 1E). Thus, βcurr is often close to one, and hence negligible (Figure 4B,C). For a conductance-based synapse, rescaling weights by βcurr does not accurately fit the full model (Figure 4B). We decompose the synaptic conductance g(t) as g+δg(t), with g the temporal average and δg(t) the fluctuations around g. We then analytically compute (see Materials and methods — Synaptic weight-rescale factors) that the weight-rescale factor for a conductance-based synapse

(8) βcond=11+Δzg

recovers the full model’s voltage if

(9) Δzδg(t)1+Δzg,

where Δz=zss-zcc is the difference between input resistance zss at the original synapse site s and input resistance zcc at the nearest proximal compartment site c. Multiplying the synaptic weights by βcond weakens the synapse so that it matches the shunt effect of its distal counterpart, and so that the voltage at all compartments more proximal than c and in their respective side branches is reproduced (Figure 4D).

If s and c are close, Δz is small and (Equation 9) is satisfied for plausible conductance fluctuations. Thus, shifting the synapse to the next proximal compartment keeps the voltage response intact. If the separation between s and c increases, Δz increases accordingly, and the magnitude of tolerated conductance fluctuations shrinks (Equation 9). Thus, for a given magnitude of δg(t), it is possible that a thick branch can be removed, but not a thin branch, as the input resistance increase in the former is much more gradual than in the latter. Nevertheless, (Equation 9) suggests that even for large Δz the full model’s voltage can be recovered by downscaling synaptic weight by βcond, but only if δg(t)g. A tonic level of AMPA and/or GABA activation, as in the high-conductance state (Destexhe et al., 2003), can implement this input condition. For NMDA synapses, this input condition does not hold, as NMDA spikes arise through a strong, transient conductance increase. However, when a cluster of such synapses – to be moved to a proximal site – is considered in isolation, a workaround can be found by applying constant rescale factors not just to synaptic weight, but also to threshold and width of the Mg2+ block, as well as to the reversal (see Materials and methods — Synaptic weight-rescale factors). These factors recover NMDA-spike shape (Figure 4E) and threshold (Figure 4F) when the cluster is moved to from s to c.

Our analytical weight-rescale factor βcond treated reductions where a single site s was moved to c. We consider now the ablation of a whole subtree and the shift of all its synapses to c (Figure 4G). In our simulations, we distributed 100 AMPA and 100 GABA synapses on the to be ablated subtree, and activated them with a Poisson process with fixed rate. We studied how the reduction influenced the local integrative properties at the intact compartment c. As a proxy for those integrative properties, we checked the shape of the waveform elicited by activating AMPA+NMDA synapses at c in a short burst (Figure 4H), and we quantified the threshold for NMDA-spike generation. We found that neither shape nor threshold were preserved by the analytical weight-rescale factors βcond, derived for the case where only a single input site is present on the subtree. We numerically extend the derivation of the weight-rescale factors to multi-site inputs (see Materials and methods – Synaptic weight-rescale factors). These extended weight-rescale factors βcondext reproduced local voltage (Figure 4H) and NMDA-spike threshold (Figure 4I). We also investigate whether βcondext can still be conceptualized as depending on Δzavggavg – the product between (1) a difference in average input resistance between synapses on the to be ablated subtree and the compartment site and (2) the average conductance load exerted by those synapses. To do so, we activated the 100 AMPA and 100 GABA synapses on the to be ablated subtree at a wide range of input conditions (see Materials and methods – Simulation-specific parameters), and measured the error in threshold for NMDA-spike generation at c. We found that this error depended strongly on Δzavggavg and remained small for Δzavggavg1 (Figure 4J).

In conclusion, we can now assess whether given spatial synapse motifs on branches or subtrees admit simplification by shifting the synapses to the nearest proximal compartment site and ablating the branch or subtree. We find that reduction in this sense is possible for any motif if fluctuations are small, so that δg(t)g, by applying temporally constant weight-rescale factors. Otherwise, the synapses must be located sufficiently closely to the compartment site, so that Δz is small and that Δzδg remains below 1+Δzg for all relevant activation levels.

Active and passive reduced dendrites under synaptic bombardment

We next investigated how many compartments are needed to reproduce somatic and dendritic responses under a synaptic bombardment and quantified the contribution of active dendritic Na+, Ca2+, and K+ channels. We mimicked the in vivo state in the basal dendrites of the L2/3 pyramid by distributing 200 AMPA and 200 GABA synapses receiving Poisson inputs and 20 clusters of AMPA+NMDA synapses receiving bursts of inputs (Figure 5A). We sampled 10 different sets of meta-parameters governing synaptic activation (Poisson and burst rates and synaptic weights, see Materials and methods — Simulation-specific parameters), resulting in output spike rates between 0 and 8 Hz (Figure 5B). We derived reduced models, once based on the full model with all active channels (‘active reduced dendrite’ in Figure 5C) and once based on the passified full model (‘passive reduced dendrite’ in Figure 5C). We distributed increasing numbers of compartment sites on the basal branches and measured somatic and dendritic voltage responses for the active full model and for the reduced models with active and passive dendrites (Figure 5C, Figure 5—figure supplement 1 for dendritic traces). Spike coincidence (Figure 5D), subthreshold somatic voltage error (Figure 5E), and dendritic voltage error (Figure 5F) showed a significant improvement when compartment numbers increased from 0 (point-neuron – top row in Figure 5C) to 2 per 100 μm. In the presence of active channels, the error decreases until up to three compartments per 100 μm (Figure 5E,F). Spike coincidence factors were consistently ∼10% higher than in the passive case (Figure 5D).

Figure 5 with 1 supplement see all
Reductions with active and passive dendritic compartments under in vivo like conditions.

(A) AMPA and GABA synapses, and AMPA+NMDA synapse clusters are spread randomly throughout the basal dendrites of the L2/3 pyramid (inset). AMPA and GABA synapses receive Poisson inputs and AMPA+NMDA synapses receive bursting input (right). (B) Output spike rates of the full model for 10 different input configurations. (C) Reductions with increasing numbers of compartment sites along the basal dendrites (left). Somatic voltage traces and spike times for the full model (gray), for a reduction with passive dendrites (middle, blue) and a reduction with active dendrites (right, orange). (D–F) Spike coincidence factors (D), relative somatic voltage errors (E) and relative dendritic voltage errors (F) for reductions with increasing numbers of compartments, and with passive (blue) and active dendrites (orange), averaged over all 10 input configurations. Relative voltage errors are computed as in Figure 2C,D.

Fitting reduced models directly to experimental data

Traditionally, creating a morphological neuron model involves a reconstruction of the morphology and an evolutionary fit of the electrical parameters (Hay et al., 2011; Almog and Korngreen, 2014; Van Geit et al., 2016). Our method allows skipping these resource-intensive steps (Figure 6A). We adapt our method to a common experimental paradigm where hyper- and depolarizing current steps are injected under applications of various ion-channel blockers (Marti Mengual et al., 2020). We extract resistance matrices and holding potentials from voltage step amplitudes (Figure 6B), and the time scale of the largest eigenmode from the average voltage decay back to rest (Figure 6C). If the modeling goal is a reduced model, the traditional reconstruction and subsequent reduction both introduce errors (Figure 6D). Thus, we find that while our reduction of the full model is faithful within an expected accuracy margin (Figure 6E), fitting a reduced model directly to the experimental traces is more accurate (Figure 6F).

Fitting reduced models directly to experimental data.

(A) Comparison between the traditional neuron model creation paradigm (morphological reconstruction and evolutionary fit, possibly followed by a reduction), and the proposed direct experiment to reduced dendrite model paradigm. (B) Resistance matrices and holding potentials are extracted from the response amplitude to hyper- and depolarizing current step inputs, here measured once under application of an h-channel antagonist and once under control conditions. (C) Time scale of largest eigenmode is extracted from average decay back to the resting membrane potential. (D) Combined root mean square voltage error RMSE(v)=Avg(vfull-vreduced)2 of full model (gray) – fitted to current step data through an evolutionary algorithm, of the reduction of the full model (red), and of the direct fit of the reduced model to the data (brown). (E) Experimentally recorded traces (gray), traces from the full reconstruction (black), and its reduction (green and purple, dashed). (F) Experimentally recorded traces (gray) and traces from the directly fitted reduced model (green and purple, dashed).


We have presented a flexible yet accurate method to construct reduced neuron models from experimental data and morphological reconstructions. The method consists of linear parameter fits at holding potentials that are informative for the full non-linear dynamics. First, we derive leak and coupling conductances from the passified version of the full model. Second, we fit the maximal conductances of the active ion-channels, so that the reduced model, at various holding potentials, responds similarly to input perturbations as the full model. Third, we fit the capacitances to reproduce spatio-temporal integrative properties. Finally, we ensure that resting membrane potentials in the full and reduced models match. The resulting method can be adapted to a wide variety of use cases.

A fundamental use case is deriving reduced models that retain specific elements of the dendritic computational repertoire, to explore how those dendritic computations may improve network function. To that purpose, one can apply our method to compartments placed at input sites needed for the computation of interest. For instance, if a computation requires independent computational subunits (Poirazi et al., 2003b; Kastellakis et al., 2016), one would place compartments at the dendritic tips. If one wants to modulate co-operativity of excitatory inputs at distal sites, one would place compartments with shunting conductances at the bifurcations between the distal sites (Wybo et al., 2019). Detecting input sequences (Branco and Häusser, 2010) or implementing a dendritic integration gradient (Branco and Häusser, 2011) would require placing two or more compartments along the branches of interest. For many of these computations, a passive morphology with non-linear (NMDA) synapses suffices. Furthermore, our method can be combined with abstract models of AP generation (Pozzorini et al., 2013; Pozzorini et al., 2015).

A second use case is constructing models directly from electro-physiological recordings. Our method avoids the labor-intensive detour of reconstructing the morphology and optimizing model parameters with an evolutionary algorithm (Marti Mengual et al., 2020). In combination with advances in voltage-sensitive dye imaging, our method may thus form the basis of a high-throughput experiment-to-model paradigm.

A third use case is elucidating the effective complexity of a given dendritic tree. We assessed whether afferent spatial input motifs could be simplified by removing a branch or subtree and grouping the affected synapses at the next proximal compartment. We found that such simplification is possible when the input resistance of the affected synapse sites is close to the input resistance of the compartment site. If this condition does not hold, the simplification can still proceed if conductance fluctuations are small compared to the average conductance magnitude, but not otherwise. What is the simplest model for a given dendritic tree that captures its full computational repertoire? We find that spike coincidence and subthreshold voltage metrics reach satisfying accuracy at, but do not improve much beyond ∼3 compartments per 100 μm. Further research is required to tease apart mere numerical errors from potentially missed computational features.

The process of simplifying complex dynamical systems, such as morphological neuron models, is inherently approximate. Due to the vast variety of morphologies, voltage-gated ion-channel configurations and input patterns, it is beyond the scope of any single paper to give an exhaustive overview of possible inaccuracies. Rather, we propose three lines of inquiry if a reduction fails to reproduce a particular response pattern. First, we suggest checking the discrepancy between the impedance kernels z(t) of the full and reduced models. Often, response mismatches arise from slight inaccuracies in how ion-channel dynamics are integrated by the passive system (Figure 3A–C). Second, we suggest assessing to what degree afferent spatial input motifs, present in the full model, are simplified, and whether that simplification still admits a sufficient accuracy in reproducing local voltage-dependent responses, such as NMDA spikes (Figure 4). Finally, we suggest checking whether the distribution of compartments on the morphology admits ion channel dynamics that are sufficiently rich to reproduce the response characteristic in question. To illustrate these lines of inquiry, we have further simplified the L5 pyramid to five compartments, by omitting the compartments in the apical trunk (Figure 3—figure supplement 1A vs Figure 3J). This reduction results a change in AP shape, and in a failure to elicit AP bursts through the pairing of a somatic and a dendritic input, even though a Ca2+ spike can still be elicited (Figure 3—figure supplement 1B). Likely, the change in AP shape can be traced back to the difference in the somatic input impedance kernels (Figure 3—figure supplement 1C), while the failure to elicit an AP burst may be due to insufficient Na+-driven cross-talk between soma and apical tuft.

Computational tools have played a key role in neuroscience. Take NEURON (Carnevale and Hines, 2004), which accelerated our understanding of morphological neurons, or NEST (Gewaltig and Diesmann, 2007), which enabled simulation of large-scale spiking point-neuron networks. We have implemented a Python toolbox that automatizes the simplification process (NEural Analysis Toolkit – NEAT –; copy archived at swh:1:rev:1cb15f36aa0a764105348541d046c85ef38e21ee). The toolbox reads any morphology in the standard ‘.swc’ format (Ascoli, 2006) and returns the parameters of the reduced models, while also providing tools to export the models to NEURON (Carnevale and Hines, 2004).

Our method and toolbox fill a void in between the extremes of modeling large-scale networks with abstract models and modeling single cells in all their detail. By enabling the efficient systematic derivation of simplified neurons, amenable to simulation at the network level, our work bridges the gap between two branches of neuroscience that historically have remained separate.

Materials and methods


Request a detailed protocol

Three exemplar cell models were used for the analysis: a cortical L2/3 pyramidal cell (Branco and Häusser, 2010; Figure 2D), a cortical L5 pyramidal cell (Hay et al., 2011; Figure 2E), and a cerebellar Purkinje cell (Figure 2F). The latter morphology was retrieved from the repository (Ascoli, 2006) and the two others from the ModelDB repository (McDougal et al., 2017).

Physiological parameters

Request a detailed protocol

In our passive models, physiological parameters were set according to Major et al., 2008: the equilibrium potential was -75mV, the membrane conductance 100μS/cm2, the capacitance 0.8μF/cm2, and the intracellular resistivity 100Ωcm. For the active models, we took ion channels and parameters according to Branco and Häusser, 2010 for the L2/3 pyramid, Hay et al., 2011 for the L5 pyramid, and Miyasho et al., 2001 for the Purkinje cell.

AMPA and GABA synaptic input currents were implemented as the product of a conductance profile, here a double exponential (Rotter and Diesmann, 1999), with a driving force:

(10) isyn=g(er-v).

AMPA rise resp. decay times were τr=0.2ms, τd=3ms and AMPA reversal potential was e=0mV. For GABA, we used τr=0.2ms, τd=10ms, and e=-80mV. NMDA currents (Jahr and Stevens, 1990a) had the form:

(11) isyn=gσ(v)(e-v)

with rise resp. decay time τr=0.2ms, τd=43ms, and e=0mV, while σ(v), modeling the channel’s magnesium block, had the form (Behabadi and Mel, 2014):

(12) σ(v)=11+0.3e-0.1v.

The ‘conductance’ of an AMPA or GABA synapse signifies the maximum value of its conductance window. For an AMPA+NMDA synapse, the conductance is the maximal value of the AMPA conductance window, and the conductance of the NMDA component is determined by multiplying the AMPA conductance value with an NMDA ratio RNMDA, set to be either 2 or 3.

Biophysical models

Request a detailed protocol

We used the NEURON simulator (Carnevale and Hines, 2004) to implement biophysical and reduced models. For the biophysical models, the distance step was set according the lambda rule (Carnevale and Hines, 2004) or smaller.

Quasi-active channels

Request a detailed protocol

A voltage-dependent ion channel described by the Hodgkin–Huxley formalism can in general be written as follows:

(13) i=g¯f(y1,,yK)(v-e),y˙k=gk(yk,v) for k=1,,K,

where g¯ is the channel’s maximal conductance, e its reversal potential, y1,,yK its state variables, f() a function that depends on the channel type (e.g. for a typical sodium channel f(m,h)=m3h), and gk(yk,v) (k=1,,K) the functions governing state-variable activation (gk(yk,v) can usually be written as (yk(v)-yk)/τk(v) with yk(v) the state-variable’s activation and τk(v) its time scale). To obtain the channel’s quasi-active approximation (Mauro et al., 1970; Koch and Poggio, 1985) around a holding potential vh and a state variable expansion point y10,,yK0, we linearize (Equation 13):

(14) ilin=g¯[k=1Kfyk(y10,,yK0)(vhe)(ykyk0)+f(y10,,yK0)(vvh)],y˙k=gkyk(yk0,vh)(ykyk0)+gkv(yk0,vh)(vvh)for k=1,,K.

To obtain the zero-frequency contribution of (Equation 14) to the resistance and conductance matrices, we set y˙k=0, solve the linearized state-variable equations for (yk-yk0), and substitute the result in ilin:

(15) ilin=g¯[-k=1Kfykgkvgkyk(vh-e)+f(y10,,yK0)](v-vh).

Introducing a shorthand for the factor in square brackets

(16) l(vh)=-k=1Kfykgkvgkyk(vh-e)+f(y10,,yK0)

results in the linearized ion-channel current as described in (Equation 2).

The impedance matrix

Request a detailed protocol

The voltage response δvx(t) at a location x along the neuron to an input current perturbation δix(t) at location x can be computed as the convolution of an impedance kernel with δix(t):

(17) δvx(t)=zxx(t)δix(t).

The impedance kernel itself can be computed in the frequency domain from the quasi-active cable equation using Koch’s algorithm (Koch and Poggio, 1985). We may assume any a priori arbitrary set of holding potentials and ion-channel state variables to compute the quasi-active expansion (with Koch’s algorithm, the only constraint is that their distribution is uniform for each cylindrical segment). For our purpose, spatially uniform holding potentials and state-variable expansion points suffice. For a steady-state current δix, (Equation 17) simplifies to:

(18) δvx=zxxδix,

where we term zxx=0dtzxx(t) the ‘resistance’ (zxx is also known as the input resistance if x=x or the transfer resistance if xx). For passive membranes, impedance kernels can also be computed as a sum of exponentials (Holmes et al., 1992; Major et al., 1993):

(19) zxx(t)=l=0ϕl(x)ϕl(x)e-tτl,

where we adopt the ordering τ0τ1τ2 Essentially, this infinite sum can be thought of as the generalization of the eigenmodes for linear dynamical systems (Strogatz, 2000) to partial differential equations.

When a current perturbation is applied at multiple sites along the neuron (we write δ𝐢(t)=(δi1(t),,δin(t))), we obtain the voltage response δ𝐯(t)=(δv1(t),,δvn(t)) at those sites from:

(20) δ𝐯(t)=Z(t)δ𝐢(t),

with Z(t) the matrix of impedance kernels (the ij’th elements of this matrix is the impedance kernel between sites i and j). In the steady-state case, we call Z the resistance matrix.

Compartmental models

Request a detailed protocol

The voltage vi in a compartment i, connected to a set 𝒩i of nearest neighbor compartments, and with a set of ion channels i, is given by

(21) civ˙i+gLi(vi-eLi)+n𝒩igC,in(vi-vn)+diidi=ii,

where ci is the capacitance of compartment i, gLi resp. eLi its leak conductance resp. reversal, gCin its coupling conductance to the neighboring compartment n, idi the ion channel current (Equation 13) of channel d in compartment i and ii an arbitrary input current to compartment i.

The conductance matrix

Request a detailed protocol

To obtain the conductance matrix of the compartmental model for steady-state inputs, (Equation 21) is linearized around a certain holding potential vh and we assume that v˙i=0. After absorbing all constant terms into ii, we obtain

(22) gLiδvi+n𝒩igC,in(δvi-δvn)+dig¯dild(vh)δvi=ii,

where δvi=vi-vh, where g¯dild(vh) is the linearized ion-channel current (Equation 15), and where we used the notation (Equation 16). Summarizing the voltage responses for each compartment in a vector δ𝐯=(δv1,,δvN), (Equation 22) for each compartment i=1,,N can be recast as a matrix equation, yielding (Equation 3).

Simplification method details

Request a detailed protocol

For any given set of M sites on the morphology, we construct a simplified compartmental model whose connection structure follows a tree graph defined by the original morphology (Figure 1A). We do not allow triplet connections (e.g. sites 1–2, sites 2–3, and sites 3–1 all mutually connected) in our reduced models. Hence, to obtain accurate results, it is necessary to extend the original set of M sites with the B bifurcation points that lie in between (Figure 1—figure supplement 1C). With these N=M+B sites, we thus define a tree graph that provides the scaffold for our fit (and our reduced model).

The fitting process proceeds in four steps: (1) fit the passive leak and coupling conductances, (2) fit the capacitances, (3) fit the maximal conductances of the ion channels, and (4) fit the reversal potentials to obtain the same resting membrane potentials as in the biophysical model. The order of steps 2 and 3 is interchangeable, but not of the other steps. We describe these steps in detail below:

  1. If the biophysical model contains active conductances, we compute their opening probabilities at rest and add them to the leak. Otherwise, we simply take the leak as is. We then compute Z for the N compartment sites. In the passive case, Z is independent of the holding potential. In (Equation 21), we substitute vi˙=0, ic=0. The terms gLieLi are a constant contribution, and can thus be absorbed in ii. We obtain: 

    (23) gLivi+n𝒩igC,in(vi-vn)=ii,i=1,,n

    which, when recast in matrix form, yields (Equation 3). (Equation 4) then yields a system of N2 equations, linear in gLi and gC,in (note that this system always has 2N-1 unknowns). We recast this system of equations in a form A𝐠=𝐛 and solve it in the least-squares sense for 𝐠.

  2. To fit the capacitances, we set ic=0 and again absorb gLieLi in ii in (Equation 21) to obtain:

    (24) v˙i=1ci(-gLivi+n𝒩igC,in(vn-vi))+iici,i=1,,n.

    In matrix form, we obtain (Equation 6) with S=diag(𝐜)-1G. We require that this matrix has the same smallest eigenvalue (corresponding to the largest time scale τ0) α0=1/τ0 and eigenvector ϕ0=(ϕ0(x1),,ϕ0(xN)) as the biophysical model, as defined by (Equation 19). Hence, we obtain the system of (Equation 7) that can be solved for 𝐜.

  3. We compute the maximal conductances of each ion-channel type separately. Thus, next to setting vi˙=0, we set ic=0 for all but one of the channels – call the non-zero channel d – and replace id with its quasi-active, zero-frequency approximation (Equation 15). We choose a spatially uniform holding potential vh and expansion point (y10,,yK0) for ion channel d. We absorb constant terms in ii and obtain (Equation 21) in terms of δvi=(vi-vh):

    (25) gLiδvi+n𝒩igC,in(δvi-δvn)+g¯dild(vh)δvi=ii.

    We again recast this equation in matrix form, but distinguish the known passive component Gpas from the unknown diagonal matrix Gvh,chan d containing as its i’th diagonal element the channel term g¯dild(vh):

    (26) (Gpas+Gvh,chan d)δ𝐯=𝐢.

    From the biophysical model, we compute Z by setting all ion channel conductances except channel d to zero, to obtain the fit (Equation 5).

    By consequence, we have a system of N2 equations, linear in the N unknown maximal conductances 𝐠¯d=(g¯d1,,g¯dN), that can be recast in a form A𝐠¯d=𝐛, where A and 𝐛 depend on the holding potential and expansion points. On the one hand, we aim to compute these linearizations for a sufficiently large range of holding potentials and expansion points. On the other hand, the fit should be restricted to domains of the ion-channel phase space where the channel resides during normal input integration, so as to avoid over-fitting on domains of the phase space that are never reached. We choose four holding potentials around which to linearize: vh=-75,-55,-35 and 15 mV. For channels with a single state variable y, we choose y0=y(vh) for these four holding potentials. For channels with two state variables, we computed y10=y1(vh) and y20=y2(vh) for the four holding potentials, and choose sixteen expansion points as all possible combinations of these two state variables. We then weigh the matrices A and vectors b with the inverse of the open probability f() for their respective expansion point, concatenate the matrices A and vectors 𝐛 for each of these expansion points, and obtain the system Aext𝐠¯d=𝐛ext of EN2 equations (with E the number of expansion points) and N unknowns. We solve this system in the least-squares sense for 𝐠¯d and repeat this procedure for each voltage-dependent ion channel in the biophysical model.

  4. Finally, we fit eLi (i=1,,N) to reproduce the resting membrane potential 𝐯eq=(veq1,,veqN) of the biophysical model evaluated at the compartment sites. We substitute this potential in (Equation 21), set v˙i=0 and ii=0 for i=1,,N and obtain: 

    (27) gLi(veqi-eLi)+n𝒩igC,in(veqi-veqn)+c𝒞iic(veqi)=0,i=1,,n.

    This is a system of N equations, linear in the N unknowns eLi (i=1,,N), and can thus be solved with standard algebraic techniques.

Synaptic weight-rescale factors

Request a detailed protocol

To be able to analytically compute the weight-rescale factor when a synapse with temporal conductance g(t) is moved from its original site s to a compartment site c on a dendritic tree, we use the steady-state approximation (Equation 18). Because we implicitly convolve the rescaled synaptic current with the temporal input impedance kernel zcc(t) at the compartment site when we simulate the reduced model, the weight-rescale factor still yields accurate temporal voltage responses (Figure 4B,E,H). In the original configuration, we have

(28) vcvs=zcczcszsczssicg(t)(evs),

with ic an arbitrary current at the compartment site and g resp. e the synaptic conductance resp. reversal. Eliminating vs from this yields

(29) vc=(zcc+(zcczss-zcszsc)g(t))ic+zcsg(t)e1+zssg(t).

In the reduced configuration we have:

(30) vc=zcc(ic+βg(t)(e-vc)).

From requiring that vc=vc, β is found as follows:

(31) β=zcszcczcsic-e[1+(zss-zcs2zcc)g(t)]zccic-[1+(zss-zcs)g(t)]e.

When zcczcs, which is true in basal dendrites and a reasonable approximation in many apical dendrites if c is on the direct path from s to the soma, (Equation 31) reduces to

(32) β=11+(zss-zcc)g(t).

The weight-rescale factor hence depends on time. Decomposing the time-varying g(t) into g+δg(t), with g the temporal average and δg(t) the fluctuations around g, we find (Equation 8) and (Equation 9) from requiring that the denominator of (Equation 32) be constant in time:

(33) β=11+Δzg+Δzδg(t)=1(1+Δzg)(1+Δzδg(t)1+Δzg)11+ΔzgifΔzδg(t)1+Δzg1

where Δz=zss-zcc.

In the case of an NMDA synapse, the conductance depends sigmoidally on the local voltage. Using scale factor (Equation 8), we obtain for the rescaled synaptic current at the compartment site:

(34) is=gσ(vs;vT,Δv)1+(zss-zcc)gσ(vs;vT,Δv)(e-vc).

This current still depends on the voltage at the synaptic site vs in the sigmoid. A current at the compartment site that causes a voltage vc there would have caused a voltage

(35) vs=veq+zsszcc(vc-veq)

if it was placed at the synapse site. Hence, to retain the same activation level of the NMDA synapse, we substitute (Equation 35) in the sigmoid. This amounts to a change in threshold and width of the sigmoid:

(36) σ(vs;vT,Δv)σ(veq+zsszcc(vc-veq);vT,Δv)=σ(vc;zcczss(vT+veq)-veq,zcczssΔv).

We will denote this modified sigmoid as σ(v). We have not yet succeeded in reformulating (Equation 34) by only rescaling its parameters with constant rescale factors, since 1/(1+(zss-zcc)gσ(vc)) still depends on time through its dependence on vc. To obtain constant rescale factors, we substitute the (Equation 34) in (Equation 17):

(37) vc-veq=zcc(t)[gσ(vc)1+(zss-zcc)gσ(vc)(e-vc)].

Here, the current in square brackets is the synaptic current. The convolution to obtain vc is implemented implicitly by the compartmental model. We then approximately eliminate the denominator: 

(38) vc-veqzcc(t)[gσ(vc)(e-vc)-zss-zcczccgσ(vc)(vc-veq)]

to obtain for the synaptic current: 

(39) is=zsszccgσ(vc)(zcczss(e-veq)+veq-vc).

We thus find that the synaptic weight is rescaled by a constant factor zsszcc and the reversal potential is shifted to zcczss(e-veq)+veq. Note however that this reduction relies on two key assumptions: that there are no other inputs to the dendritic tree and that there are no fluctuations in vc generated by input currents not at site s, as this could lead to spurious activation or suppression of σ(vc).

Suppose now we shift M conductance synapses to N compartments. Let 𝐯=(v1,,vN,vN+1,,vN+M) be the vector with K=N+M components containing the voltages at the N compartment sites and at the M synapse sites. Similarly 𝐠=(0,,0,g1,,gM) is a vector of K components containing N zeros and the M synaptic conductances, while 𝐞=(0,,0,e1,,eM) contains the synaptic reversals. In matrix form, (Equation 18) for this system becomes:

(40) 𝐯=Zdiag(𝐠)(𝐞-𝐯),

and its solution:

(41) 𝐯=Γ𝐞,

with Γ a matrix given by:

(42) Γ=(I+Zdiag(g))1Zdiag(g).

In the reduced setting, we assign each synapse to one compartment site (a reasonable choice could be the closest site). We introduce an N×M compartment assignment matrix C where an element cnm is 1 if synapse m is assigned to compartment n and zero otherwise. We further introduce a vector of reduced compartment voltages 𝐯=(v1,,vN), a vector of rescaled synaptic conductances 𝐠β=(β1g1,,βMgM) and a vector of synaptic reversals 𝐞=(e1,,eM). In the reduced model, voltages are obtained from:

(43) 𝐯=ZCdiag(𝐠β)(𝐞-CT𝐯).

We then require that 𝐯N=𝐯, with 𝐯N a vector containing the first N components of 𝐯. We denote by ΓNM the matrix containing the first N rows and last M columns of Γ, and note that we can write (Equation 41) as 𝐯N=ΓNM𝐞. Substituting this in (Equation 43) yields:

(44) ΓNM𝐞=ZCdiag(𝐠β)(𝐞-CTΓNM𝐞).

We then fit the parameters β1,,βM (here absorbed in 𝐠β) in the least-square sense from:

(45) ΓNM=ZCdiag(gβ)(ICTΓNM),

which again is a linear fit.

Experimental recordings

Request a detailed protocol

Coronal slices (300 μm thick) containing the anterior cingulate cortex (ACC) were prepared from 10 to 12 week old C57BL/6 mice using a vibratome on a block angled at 15 degree to the horizontal in ice-cold oxygenated artificial cerebral spinal fluid (ACSF) and then maintained in the same solution at 37°C for 15–120 min. Normal ACSF contained (in mM) NaCl, 125; NaHCO3, 25; KCl, 2.5; NaH2PO4, 1.25; MgCl2, 1; glucose, 25; CaCl2, 2; pH 7.4. Individual neurons were visualized with a Nikon Eclipse E600FN fit with a combination of oblique infrared illumination optics and epifluorescence, the switch between optical configurations was software-triggered (Sieber et al., 2013). Pyramidal neurons were selected on the clearly visible, proximal apical dendrite. This selection criterion resulted in a homogeneous population of pyramidal neurons based on their firing properties and shape of the AP (i.e. all cells possessed a prominent after-hyperpolarization and a significant sag ratio at the soma). Dual somatic and dendritic whole-cell patch-clamp recordings were performed from identified L5 pyramidal neurons in the rostroventral ACC (1.1–1.4 mm below the pial surface, 1.1–0.2 mm rostral to the Bregma) using two Dagan BVC-700. During the experiments, the external recording solution (normal ACSF) was supplemented with 0.5 mM CNQX and 0.5 mM AP-5 to block excitatory glutamatergic synaptic transmission. Experiments were performed at physiological temperatures between 34–37°C. Whole-cell recording pipettes (somatic, 4 to 8 MΩ; dendritic, 12 to 32 MΩ), were pulled from borosilicate glass. The internal pipette solution consisted of (in mM) potassium gluconate, 135; KCl, 7; Hepes, 10; Na2-phosphocreatine, 10; Mg-ATP, 4; GTP, 0.3; 0.2% biocytin; pH 7.2 (with KOH); 291–293 mosmol l−1. For somatic recordings, 10–20 μm Alexa 594 was added to the intracellular solution: first, the soma was patched (whole-cell configuration by negative pressure); after 5 min of intracellular perfusion, the fluorescent signal allowed for the clear identification of the apical dendritic tree, then the dendritic region of interest was patched with a smaller pipette. Compensation was performed in current clamp mode by recovering the fast, initial square voltage response to a hyperpolarizing current injection (−100 pA, 50 ms). First the pipette capacitance was compensated to the level that the voltage response showed an immediate voltage drop due to the series resistance of the pipette that was adjusted subsequently. Compensation of dendritic series resistance yielded values between 30 and 60 MΩ for pipettes with a resistance between 17 and 21 MΩ. On average, series resistance was 2.3 times larger than the pipette resistance. Series resistance of both somatic and dendritic recording electrodes was monitored frequently, and experiments were terminated when proper compensation was not possible anymore (i.e. reached values of more than four times the pipette resistance). All cells were filled with biocytin, and PFA-fixed slices were developed with the avidin–biotin-peroxidase method for Neurolucida reconstructions (Egger et al., 2008). Data analysis was performed using Igor software (Wavemetrics) and Excel (Microsoft).

Simulation-specific parameters

Parameters Figure 2

Request a detailed protocol

For the sequence detection (Figure 1f), the synapse model is as in Branco and Häusser, 2010. The maximal conductance of the AMPA component is g¯AMPA=0.5 nS and its conductance window shape given by an alpha function (Rotter and Diesmann, 1999) with a time scale of 2 ms. The NMDA component is given by a kinetic model (Destexhe et al., 1998) with g¯NMDA=8 nS and external magnesium concentration of 1 mM, and receives an exponentially decaying (τ=.5 ms) neurotransmitter concentration with amplitude of 5 mM. For the input-order detection (Figure 1G), AMPA synapses 1 and 2 use the standard parameters and have respective maximal conductances of 10 and 5 nS. For the simulation with the L5 pyramidal cell, excitatory synapses have AMPA+NMDA components with RNMDA=2 and g¯=3 nS. For inhibitory synapses, g¯=2 nS. In the Purkinje cell, excitatory synapses only have AMPA components with g¯=10 nS, while for inhibitory synapses, g¯=5 nS. For the L2/3 pyramid, L5 pyramid, resp. Purkinje cell, excitatory Poisson rates are 2, five resp. 6 Hz and inhibitory firing rates are 4, one resp. 2 Hz. Simulation were run for 10,000 ms.

Parameters Figure 3

Request a detailed protocol

APs are evoked with a DC current pulse. For panels A–C, pulse amplitude is iamp=0.5 nA and pulse duration is tdur=5 ms. For panels D–F, we have iamp=3 nA and tdur=1 ms. For panels G–I, iamp=1.5 nA and tdur=5 ms. For panels J–K, the somatic current pulse had iamp=1.9 nA and tdur=5 ms, while the dendritic current injection had a double exponential waveform, with τr=0.5 ms and τd=5 ms and amplitude iamp=0.5 nA. Onset of the somatic current pulse precedes onset of the dendritic current injection by 5 ms.

Parameters Figure 4

Request a detailed protocol

In panel B, we inject Ornstein–Uhlenbeck (OU) processes for current (mean μ=0.08 nA and standard deviation σ=0.025 nA) and conductance (μ=0.005 μS, σ=0.0025 μS). Reversal was 0 mV for the excitatory conductance and −80 mV for the inhibitory conductance. All OU processes had a time scale of 30 ms. In panels D–F, we use AMPA+NMDA synapses with g¯=0.5 nS and RNMDA=3. Burst firing is mimicked by drawing synapse activation times from a Gaussian distribution with as mean the burst time and a standard deviation of 2 ms.

In panels G–J, AMPA+NMDA synapses at the compartment site (green square) have g¯=1 nS and RNMDA=2. To determine the NMDA-spike threshold, we activate between 0 and 200 synapses in a burst. The burst is again modeled by drawing spike times from a Gaussian distribution with a 2 ms standard deviation. We then multiply the amplitude of the resulting waveform with its half-width, and average this quantity over five trials for each number of activated synapses. The number of synapses at which this quantity increases most is taken to be the NMDA-spike threshold.

We aim to test the synaptic weight-rescale factors for the AMPA and GABA synapses βcondext under a wide range of input conditions. To do so, we conduct simulations with different total time-averaged conductance loads gavg, exerted by the AMPA and a GABA synapse on the subtree that is to be reduced. These synapses are activated with spike times drawn from a homogeneous Poisson process. We also also change the size of the conductance fluctuations, by conducting simulations with different firing rates for the homogeneous Poisson processes (since the total conductance load for a simulation is fixed, a small firing rate means that individual synaptic weights have to be increased to reach the desired total conductance load, hence resulting in larger conductance fluctuations). We furthermore change the number of locations on the subtree nloc over which the conductance load is distributed. Note that by consequence, the conductance load at each location is gavg/nloc. We also adapt these locations to a specified average difference in input resistance Δzavg between the compartment site and the location. Finally, we change the average voltage level achieved by the combined AMPA and GABA input by changing the ‘nudging potential’ en (Urbanczik and Senn, 2014)


where gAMPA resp. gGABA are the time-average conductances of individual AMPA resp. GABA synapses, and gAMPA=0 mV resp. gGABA=-80 mV their respective reversal potentials.

We thus draw five meta-parameters for each simulation: gavg, ravg, nloc, Δzavg and en. The time-averaged conductance at each location is given by 

(47) gavg/nloc=gAMPA+gGABA.

Together with (Equation 46), this fixes the time-averaged conductance of each individual AMPA synapse and GABA synapse. Since the temporal conductance profile of a synapse (e.g. the AMPA synapse) following the arrival of a single input AP is modeled as the product of a weight factor wAMPA and a unitary conductance window guAMPA(t), the time-averaged conductance of that synapses activated by a homogeneous Poisson process input is 

(48) gAMPA=wAMPAravg0dtguAMPA(t).

Thus, given ravg and gAMPA, the synaptic weight wAMPA can be extracted (and similarly for the GABA synapse).

To perform simulations across a wide range of input conditions, we draw 200 samples of these five meta-parameters from the intervals: en[-80 mV,-50 mV], gavg[0.01 nS,300 nS], ravg[1 Hz,100 Hz], nloc[1,20], and Δzavg[0 MΩ,1023.2 MΩ] following the Latin hypercube (LH) method (Press et al., 2007). gavg and ravg are sampled on a log scale, whereas for all other parameters, we use a linear scale.

Parameters Figure 5

Request a detailed protocol

To obtain compartment sites, we divide the longest branch in each basal subtree in n parts of equal length, thus giving us n distances from the soma, where we increase n from 0 (point-neuron) to 10. In each subtree, we distribute compartments at all sites at these distances, and also add all bifurcation sites in between compartment sites. In this way, we obtain 11 reductions that we quantify according to ‘no. of compartments per 100 μm of dendrite’. We implement 400 ‘background’ synapses; 200 AMPA and 200 GABA synapses with an average conductance of gAMPA resp. gGABA=4gAMPA and firing rate rAMPA=ravg resp. rGABA=ravg. We implement 20 synapse clusters, consisting of AMPA+NMDA synapses with maximal conductance g¯A+N and RNMDA=2. These clusters are activated with a Poissonian burst rate rburst, the number of spikes per burst is drawn from a Poisson distribution with parameter navg, and spike times are drawn for each burst time according to a Gaussian distribution with standard deviation of 5 ms. For each of the parameters not given previously, 10 LH samples were drawn from gAMPA[0.4nS,0.8nS], ravg[1.5Hz,3.0Hz], g¯A+N[0.5nS,1.5nS], navg[10,30], and rburst[0.25Hz,0.6Hz], resulting in 10 different output spike rates shown in B. For each parameter set, a simulation was run for 10,000 ms.

Parameters Figure 6

Request a detailed protocol

Hyper- resp. depolarizing current steps have iamp=-.3 nA resp. iamp=.1 nA and tdur=500 ms. The full model was optimized with an evolutionary algorithm using the BluePyOpt library (Van Geit et al., 2016), where we ran 100 iterations with and offspring size of 100. Goodness of fit was evaluated in a multi-objective manner as the root mean square error of the resting voltage (average voltage 100 ms before each current step), the final step voltage amplitudes after sag (average voltage during the last 100 ms of the DC current injection), and the voltage root mean square error of the full trace. We optimized the specific membrane capacitance and the conductance densities for passive leak and h, Kir and Km channels. The membrane currents followed an exponential distribution g(x)=g0ex/dx, with x the distance from the soma, and as parameters g0 – the conductance at the soma – and dx – the length constant of the distribution.

Data and software availability

Request a detailed protocol

NEAT (NEural Analysis Toolbox), our open-source Python toolbox to obtain reduced models, is available on  (copy archived at swh:1:rev:1cb15f36aa0a764105348541d046c85ef38e21ee).

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files.

The following previously published data sets were used
    1. Hay E
    2. Hill S
    3. Schürmann F
    4. Markram H
    5. Segev I
    (2011) ModelDB
    ID 139653. L5b PC model constrained for BAC firing and perisomatic current step firing.
    1. Branco T
    2. Clark BA
    3. Häusser M
    (2010) ModelDB
    ID 140828. Dendritic Discrimination of Temporal Input Sequences.
    1. Chen XR
    2. Heck N
    3. Lohof AM
    4. Rochefort C
    5. Morel M
    6. Rosine W
    7. Doulazmi M
    8. Marty S
    9. Cannaya V
    10. Avci HX
    11. Mariani J
    12. Rondi-Reig L
    13. Vodjdani G
    14. Sherrard RM
    ID NMO_100072. Mature Purkinje Cells Require the Retinoic Acid-Related Orphan Receptor-α (RORα) to Maintain Climbing Fiber Mono-Innervation and Other Adult Characteristics.


  1. Book
    1. Destexhe A
    2. Mainen ZF
    3. Sejnowski TJ
    Kinetic models of synaptic transmission
    In: Koch C, Segev I, editors. Methods in Neuronal Modeling. MIT Press. pp. 1–26.
  2. Book
    1. DiLorenzo PM
    2. Victor JD
    Can we predict every spike?
    In: DiLorenzo P. M, editors. Spike Timing. CRC Press. pp. 70–72.
  3. Book
    1. Press WH
    2. Teukolsky SA
    3. Vetterling WT
    4. Flannery BP
    Numerical Recipes 3rd Edition: The Art of Scientific Computing
    Cambridge University Press.
  4. Book
    1. Strogatz SH
    Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering
    Westview Press.

Article and author information

Author details

  1. Willem AM Wybo

    Department of Physiology, University of Bern, Bern, Switzerland
    Conceptualization, Resources, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1385-4980
  2. Jakob Jordan

    Department of Physiology, University of Bern, Bern, Switzerland
    Resources, Software, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3438-5001
  3. Benjamin Ellenberger

    Department of Physiology, University of Bern, Bern, Switzerland
    Resources, Software
    Competing interests
    No competing interests declared
  4. Ulisses Marti Mengual

    Department of Physiology, University of Bern, Bern, Switzerland
    Resources, Data curation
    Competing interests
    No competing interests declared
  5. Thomas Nevian

    Department of Physiology, University of Bern, Bern, Switzerland
    Supervision, Funding acquisition, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9804-608X
  6. Walter Senn

    Department of Physiology, University of Bern, Bern, Switzerland
    Conceptualization, Supervision, Funding acquisition, Visualization, Writing - original draft, Project administration, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3622-0497


H2020 European Research Council (720270)

  • Walter Senn

Swiss National Science Foundation (180316)

  • Walter Senn

Swiss National Science Foundation (159872)

  • Thomas Nevian

H2020 European Research Council (785907)

  • Walter Senn

H2020 European Research Council (945539)

  • Walter Senn

H2020 European Research Council (682905)

  • Thomas Nevian

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


This work was supported by the Swiss National Science Foundation (Grant 159872 to TN; Grant 180316 to WS – with F Helmchen), the European Union’s Horizon 2020 Framework Programme (Grant 720270, 785907 and 945539 to WS – Human Brain Project SGA 1–3) and the European Research Council (Grant 682905 to TN). We thank all our lab colleagues for helpful discussions and critical comments on the figures.

Version history

  1. Received: July 10, 2020
  2. Accepted: January 4, 2021
  3. Version of Record published: January 26, 2021 (version 1)


© 2021, Wybo et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 1,623
  • 241
  • 21

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Willem AM Wybo
  2. Jakob Jordan
  3. Benjamin Ellenberger
  4. Ulisses Marti Mengual
  5. Thomas Nevian
  6. Walter Senn
Data-driven reduction of dendritic morphologies with preserved dendro-somatic responses
eLife 10:e60936.

Share this article

Further reading

    1. Neuroscience
    Elissavet Chartampila, Karim S Elayouby ... Helen E Scharfman
    Research Article

    Maternal choline supplementation (MCS) improves cognition in Alzheimer’s disease (AD) models. However, the effects of MCS on neuronal hyperexcitability in AD are unknown. We investigated the effects of MCS in a well-established mouse model of AD with hyperexcitability, the Tg2576 mouse. The most common type of hyperexcitability in Tg2576 mice are generalized EEG spikes (interictal spikes [IIS]). IIS also are common in other mouse models and occur in AD patients. In mouse models, hyperexcitability is also reflected by elevated expression of the transcription factor ∆FosB in the granule cells (GCs) of the dentate gyrus (DG), which are the principal cell type. Therefore, we studied ΔFosB expression in GCs. We also studied the neuronal marker NeuN within hilar neurons of the DG because reduced NeuN protein expression is a sign of oxidative stress or other pathology. This is potentially important because hilar neurons regulate GC excitability. Tg2576 breeding pairs received a diet with a relatively low, intermediate, or high concentration of choline. After weaning, all mice received the intermediate diet. In offspring of mice fed the high choline diet, IIS frequency declined, GC ∆FosB expression was reduced, and hilar NeuN expression was restored. Using the novel object location task, spatial memory improved. In contrast, offspring exposed to the relatively low choline diet had several adverse effects, such as increased mortality. They had the weakest hilar NeuN immunoreactivity and greatest GC ΔFosB protein expression. However, their IIS frequency was low, which was surprising. The results provide new evidence that a diet high in choline in early life can improve outcomes in a mouse model of AD, and relatively low choline can have mixed effects. This is the first study showing that dietary choline can regulate hyperexcitability, hilar neurons, ΔFosB, and spatial memory in an animal model of AD.

    1. Neuroscience
    Guozheng Feng, Yiwen Wang ... Ni Shu
    Research Article

    Brain structural circuitry shapes a richly patterned functional synchronization, supporting for complex cognitive and behavioural abilities. However, how coupling of structural connectome (SC) and functional connectome (FC) develops and its relationships with cognitive functions and transcriptomic architecture remain unclear. We used multimodal magnetic resonance imaging data from 439 participants aged 5.7–21.9 years to predict functional connectivity by incorporating intracortical and extracortical structural connectivity, characterizing SC–FC coupling. Our findings revealed that SC–FC coupling was strongest in the visual and somatomotor networks, consistent with evolutionary expansion, myelin content, and functional principal gradient. As development progressed, SC–FC coupling exhibited heterogeneous alterations dominated by an increase in cortical regions, broadly distributed across the somatomotor, frontoparietal, dorsal attention, and default mode networks. Moreover, we discovered that SC–FC coupling significantly predicted individual variability in general intelligence, mainly influencing frontoparietal and default mode networks. Finally, our results demonstrated that the heterogeneous development of SC–FC coupling is positively associated with genes in oligodendrocyte-related pathways and negatively associated with astrocyte-related genes. This study offers insight into the maturational principles of SC–FC coupling in typical development.