Introduction

Characterizing the relationship between the dynamics of individual ion channels and the neuronal activity in vivo is crucial for mechanistic understanding of biological computation, to develop more realistic circuit models of brain activity, and to identify potential therapeutic targets.

Ion channels have been extensively studied both in isolation (Hille, 2001) and also in their interactions with other channels to explain the emergence of diverse neuronal activity patterns observed (Hodgkin and Huxley, 1952; Koch, 1999). A dominant view emerged that neurons form non-linear dynamical systems, with state variables corresponding to membrane potential, gate openings, and second messenger concentrations (Izhikevich, 2007). Studying the impact of particular ion channels on the neuronal response in such a dynamical system under in vivo conditions is hindered by two major obstacles: First, synaptic inputs to neurons arrive in complex temporal patterns often engaging a combination of ion channels that interact with each other in potentially highly nonlinear, state-dependent manner. Measuring channel-contributions under in vitro conditions (Losonczy and Magee, 2006) or by pharmacological manipulations (Palmer et al., 2014) may not engage the ion channels in the same way as during the high-conductance state in vivo, (Destexhe et al., 2003; Ujfalussy et al., 2018). In response to these challenges, the currentscape technique has recently been developed to provide an intuitive and simple way to visualize the high-dimensional current dynamics during naturalistic neural activity in one-compartmental biophysical neuron models (Alonso and Marder, 2019).

Second, synaptic inputs target the spatially extended dendritic tree of the neurons (Fig. 1B) where electrical and biochemical compartmentalization renders state variables local (Stuart and Häusser, 2001; London and Häusser, 2005; Branco and Häusser, 2010). With local state variables, the recruitment of ion channels can vary greatly between different dendritic domains (Häusser and Mel, 2003; Poirazi et al., 2003), which makes it especially difficult to understand or visualise current-interactions in neurons with a large dendritic tree under in vivo conditions. There are three simple ways to adapt the currentscape technique to neurons with realistic morphology: First, one could visualize the current dynamics in each compartment separately (Guet-McCreight and Skinner, 2020; Linaro et al., 2022). This strategy could work for few compartmental models but is clearly not feasible beyond a critical size and does not provide an intuitive description of the current propagation in neurons. Second, it is possible to sum up all membrane currents throughout the dendritic tree (Fig. 1A-D). However, this approach can not reveal the contribution of currents to the output of the cell, as even large-amplitude distal dendritic currents often fail to propagate to the soma (white arrowhead in Fig. 1D). Third, one could focus on a single compartment and study the input currents locally, for example at the somatic action potential generation site (Linaro et al., 2022). However, most input currents can be currents flowing axially as a consequence of the activation of synaptic or intrinsic currents in other compartments (Fig. 1D). Here, we develop a computational technique, the extended currentscape, which is able to trace the axial currents back to their origin: the location and type of the membrane currents that generated them.

Fig. 1:

To demonstrate the ability of the extended currentscape method to provide an intuitive visualisation of input integration in neuron models with complex morphology under in vivo-like conditions, we apply it to study the generation of complex spike bursts (CSB) in hippocampal CA1 pyramidal neurons (PCs). In PCs CSBs are associated with Ca²⁺-spikes in the apical dendrites (Larkum et al., 1999; Grienberger et al., 2014; Magó et al., 2021) and are believed to play a critical role in learning (Payeur et al., 2019) and memory formation (Bittner et al., 2015, 2017; Grienberger and Magee, 2022). However, the precise input conditions necessary to trigger Ca²⁺-spikes in these cells remained controversial: On the one hand, Ca²⁺-spikes can be readily initiated by distal dendritic, but not by somatic current injection (Golding et al., 1999; Takahashi and Magee, 2009), the timing of the Ca²⁺-spikes coincides with inputs to the most distal branches during theta oscillation (Bittner et al., 2015) and blocking entorhinal inputs decreases the probability of CSB firing (Bittner et al., 2015). On the other hand, CSBs can be triggered by somatic current injection in most PCs at random spatial locations (Bittner et al., 2015, 2017) and distal tuft dendrites are variably recruited during putative CSB events (O’Hare et al., 2025) in vivo.

In this paper, we first describe the extended currentscape technique that generalizes the standard currentscape visualization method to neurons with spatially extended dendritic trees. Next, we show that the extended currentscape accurately and intuitively captures the origin and type of dendritic events driving place cell-like activity in a biophysical model CA1 PC. We further show that the membrane potential throughout the entire dendritic arbor has low-dimensional dynamics dominated by global activity even when the cell is driven by a few dendritic branches receiving strong input. Next, we analyze CSBs in the model using the extended currentscape method and show that they can be started with highly variable initiation dynamics. By contrasting CSBs and isolated single spikes we reveal that although distal inputs facilitate the generation of CSBs, there is no need for exceptionally strong tuft inputs nor Ca²⁺ hotspots for their initiation.

Results

The extended currentscape method

We start by describing the extended currentscape visualization method which aims to capture how membrane currents in a distal reference (e.g., dendritic) compartment influence the activity of a target (e.g., soma) compartment. Our method relies on the equivalent circuit model of neurons (Koch and Segev, 2000) and requires measuring the membrane and axial currents throughout the dendritic tree of a neuron (in every node of the circuit) to partition the axial currents flowing into the target compartment by the underlying membrane currents. Since this kind of data is not available in real neurons, we used a multicompartmental biophysical model (Jarsky et al., 2005; Ujfalussy and Makara, 2020) to develop and test our method. With the advent of novel single-neuron voltage imaging techniques (Brooks et al., 2024; Park et al., 2025), our method could potentially also be adapted to data from real neurons in the near future.

The intuition behind our method is that distal currents can directly influence the membrane potential of the target compartment only if there is a continuous flow of axial current from reference to target. In order to characterise the relationship between the axial currents and the membrane currents we start with Kirchhoff’s current law, stating that the sum of all currents flowing into or out of a node in an electrical circuit must be zero. Applying it to a small section of a neuronal process, we have:

Here is the capacitive current, I^mem denotes the membrane currents, including all synaptic and intrinsic currents, I ^ax denotes the axial currents and I_inj is the current injected into the cell through a stimulation electrode. Eq. (1) states that if the axial current between the reference and the target compartment is negative (i.e., it is an outward current for the reference compartment as in Fig. 2A) then the sum of the inward currents must be larger than the sum of the remaining outward currents in the reference compartment. In this case, we postulate that this excess inward current is responsible for the outward flowing axial current, and we partition the axial current in the proportion of the inward currents in the reference compartment (Fig. 2A). In contrast, when the direction of the axial current between the reference and the target reverses, we partition it proportionally to the outward currents in the reference compartment (Fig. 2B). Note that the sign of the membrane currents can also reverse when the membrane potential crosses its reversal potential. In our simulations this often happened with the leak current, which had a reversal potential E_leak = −66 mV, falling between the reversal potential of excitatory (E_NMDA = 0 mV) and inhibitory inputs E_{GABA, slow} = −80 mV. This way, the same current type can contribute to inward or outward axial currents at different time points or in different compartments.

Fig. 2:

In a multicompartmental model, the partitioning of the axial currents can be applied recursively starting from the most distal compartment and proceeding towards the target iteratively once all axial currents arriving at the intermediate targets from the more distal branches had been partitioned (Fig. 2C-D). In the next step, the partitioned axial currents are added to the membrane currents in the intermediate target (Fig. 2E) before partitioning the axial current flowing towards the next target (Fig. 2F-G).

The topology of the model is represented by an acyclic graph, where each node is a compartment, the directed edges represent the flow of the axial current and the root node of the graph is the target (Fig. 2H). The nodes where the direction of the axial current reverses are collision points (Fig. 2I). Since the propagation of the axial current is blocked at the collision points, the membrane currents of the more distal nodes do not have a direct influence on the activity at the target. Therefore, to make the algorithm computationally more efficient, in every time step we truncate the graph at the collision points and apply the partitioning recursion only to the truncated graph (Fig. 2J-K).

After partitioning the axial currents of the target compartment we can use the standard currentscape plot (Alonso and Marder, 2019) to visualize the contribution of distal membrane currents to the membrane potential dynamics of the target. The entire process must be repeated when a new target compartment (e.g. a dendritic branch) is selected. Throughout the paper, we will use two different variants of the partitioning algorithm: we will either partition axial currents by the type of membrane currents (e.g., current flowing through Ca²⁺ or Na⁺ channels) or by the dendritic region of the reference compartments (i.e., basal, oblique or tuft dendrites). However, alternative variants of the partitioning algorithm can also be applied (see Methods). Next, we will apply the extended currentscape technique to analyze input integration and the mechanism of burst firing in hippocampal place cells.

Currentscape analysis of dendritic integration CSBs in the CA1 pyramidal neuron model

To study the synaptic input conditions that lead to burst firing under in vivo-like conditions we extended a previous CA1 PC model (Jarsky et al., 2005; Ujfalussy and Makara, 2020) with a calcium and a potassium channel. Experimental data indicates that R-type Ca²⁺-channels are mainly responsible for both burst firing and for dendritic plateau potentials in this cell type (Magee and Carruth, 1999; Metz et al., 2005; Takahashi and Magee, 2009). We used a novel Ca²⁺-channel model that displayed similar activation and inactivation kinetics to the Ca²⁺ currents recorded in the apical dendrites in CA1 PCs (Fig. S1A; Magee and Johnston 1995) and distributed it uniformly across all apical dendrites (trunk, obliques, and tuft) of the cell (Magee and Johnston, 1995; Poirazi et al., 2003).

Various potassium channels contribute to the repolarization after Ca²⁺-spikes, including Ca²⁺-activated potassium channels (Golding et al., 1999; King et al., 2015), but the precise biophysical characterization of these channels is currently lacking. To model their overall effect, we used a high voltage activated potassium current with relatively slow kinetics (Fig. S1B) in the apical dendrites.

Equipped with these channels, our model was able to fire a short burst of 3-4 action potentials with decreasing amplitude riding on a depolarizing wave upon dendritic current injection (Fig. S1C). These somatic bursts were accompanied by dendritic Ca²⁺-spikes, so we refer to them as complex spike bursts or CSBs. In agreement with the experimental data (Golding et al., 1999) the model had a lower threshold for initiating Ca²⁺-spikes in the apical trunk than in the soma (Fig. S1D). Under in vivo-like synaptic input conditions (see below and Methods), dendritic Ca²⁺-spikes could also be evoked by somatic current injection (Fig. S1E), as in Bittner et al. (2015).

To illustrate the ability of the extended currentscape method to provide a compact and intuitive summary of dendritic events in the model, we first tested it under simple, spatially and temporally restricted input conditions. We first used the model without dendritic Ca²⁺-channels and stimulated an increasing number of synapses in a single oblique branch Fig. 3A-C). We found that the somatic response became superlinear beyond a branch-specific threshold, with a sigmoid superlinearity (Fig. 3A-B), reminiscent of experimental findings with 2-photon glutamate uncaging (Losonczy and Magee, 2006). Partitioning the somatic currents by the current type (Fig. 3C, middle) indicated that before stimulation, the inward currents were mainly leak currents (dark green), with some contribution from the somatic and axonal Na⁺-currents (red). Partitioning the somatic currents by the region of their origin (Fig. 3C, bottom) showed, that current flew from the axon and basal dendrites towards the soma and further towards the apical dendrites, which was more hyperpolarized than the soma due to the larger density of K⁺ channels in the apical trunk. Upon stimulation, the inward currents show first a fast, AMPA-mediated and later a slower NMDA-mediated component (yellow and olive in Fig. 3C, middle), both originating from the stratum radiatum (levander in Fig. 3C, bottom). At n=20 inputs a Na⁺-spike that remains local to the stimulated branch is responsible for the spikelet appearing in the somatic response (pink in Fig. 3C, right). During these events, the inward current arrived transiently from the apical dendrites (see shades of levander in Fig. 3C, bottom), and left the soma towards the basal dendrites (pink in Fig. 3C). At the end of the events the NMDA contribution disappeared from the soma as the axial current reversed to its original direction, flowing from basal towards apical dendrites.

Fig. 3:

Next we repeated the same stimulation protocol after including the Ca²⁺-channels to the apical dendrites of the neuron. The addition of the Ca²⁺-channels rendered the shape of the superlinearity more step-like, as in single-photon uncaging or synaptic stimulation experiments (Fig. 3D-E; Wei et al. (2001); Ariav et al. (2003); Cai et al. (2004)). Under these conditions the Ca²⁺-spikes did not lead to somatic action potential firing but remained localised to the stimulated branch. The currentscape analysis revealed that although Ca²⁺-channels first activated at near-threshold responses (at n=10 inputs, orange in Fig. 3F), they were greatly boosted by local Na⁺-spikes (at n=15 or 20 inputs; Fig. 3F). During these events NMDA-mediated- and Ca²⁺-currents contributied similarly to the somatic depolarization (Fig. 3). Taken together, these simulations confirmed that extended currentscapes provide a rich form of visualization of the sequence of dendritic events leading to somatic responses under relatively simple input conditions. Next, we used this model to investigate the synaptic input patterns leading to CSB firing during in vivo-like inputs.

CSBs in the CA1 pyramidal neuron model

To create in vivo-like input condition we simulated the activity of 2000 excitatory and 200 inhibitory presynaptic neurons during the traversal of a 2 m long linear track as described previously (Fig. 4A; (Ujfalussy and Makara, 2020; Kim et al., 2023)). Briefly, each excitatory neuron had a ∼ 20 cm long place field and neurons showed theta modulation and displayed phase precession. Place fields were distributed uniformly across the linear track. Inhibitory neurons showed theta modulation but were spatially untuned. The majority of the 2000 excitatory synapses were placed randomly on the dendritic tree (avoiding the soma and the proximal apical trunk; light green in Fig. 1B), while the remaining 240 inputs active in the middle of the track were organized into 12 functional synaptic clusters where presynaptic neurons with similar place fields (dark green circles in Fig. 4A) innervated neighbouring dendritic locations (Fig. 1B) and had larger synaptic conductance (see Methods). These synaptic clusters provided strong drive to the cell in the middle of the track by activating a few dendritic branches (Ujfalussy and Makara, 2020). Half of the inhibitory cells targeted the soma and the proximal apical trunk, the other half was randomly distributed throughout the dendritic tree (blue circles in Fig. 1B).

Fig. 4:

First, we analyzed the somatic and dendritic membrane potential response of the model. The neuron responded to synaptic stimulation with a train of action potentials in the middle of the simulated track consisting of a variable number of single spikes and CSBs (Fig. 4B). During CSBs the distal apical dendrites showed a membrane potential response similar to the Ca²⁺-spikes observed in vitro (Fig. 4B, bottom trace). In a set of 16 simulations with variable cluster placement and synaptic input patterns, the membrane potential of the entire neuron was highly correlated (Fig. 4C), with a single factor explaining ∼70% of the total variance across the dendritic branches (Fig. 4D). The weights associated with this first factor uniformly covered the entire dendritic tree (Fig. 4E, see also Fig. S2) indicating the dominance of the global modulation of the dendritic membrane potential mainly due to back-propagating action potentials (bAPs, Spruston et al. 1995). In contrast, the weights associated with the second principal component increased consistently along the apical dendrite, demonstrating a greater level of independence from the soma in distal tuft branches (Fig. 4E).

When we inspected the dendritic membrane potential during CSBs we found that in our model CSBs are not preceded by large, global Ca²⁺-spikes in the tuft. Instead, dendritic Ca²⁺-spikes typically started asynchronously (O’Hare et al., 2025), often after a local Na⁺-spike or a back-propagating action potential (Fig. 4F, Park et al. 2025). Although in distal dendrites we could occasionally observe local Na⁺ or Ca²⁺-spikes not coupled to somatic activity (Fig. 5), local dendritic membrane potential peaks closely followed somatic APs in the majority of cases (>90%; Fig. 4G). This indicates that dendritic events tend to follow somatic APs even when the neuron is mainly driven by strong, functionally clustered synaptic inputs to a few dendritic branches (Ujfalussy and Makara, 2020).

Fig. 5:

Currentscape analysis of place field dynamics

To better understand how individual dendritic events influenced the somatic response of the cell, we calculated the extended currentscape of the simulated neuron in the time window around its activity in the place field (Fig. 5). Partitioning the somatic axial currents by current type revealed that before entering the place field, the cell was driven mainly by synaptic currents (Fig. 5Ba-c). Voltage dependent Na⁺ or Ca²⁺ channels only activated occasionally and briefly, and their effect often remained localized, not being able to efficiently propagate to the soma (Fig. 5Ba-c). When the effect of dendritic Na⁺-spikes reached the soma, they either appeared as small spikelets (Fig. 5Ba) or as fast rise in the somatic membrane potential (Fig. 5Ca). Partitioning the somatic axial currents by their origin indicated that the soma was driven mainly by the basal dendrites (Fig. 5Bd). Performing similar partitioning of the axial currents in the tuft by their type (Fig. 5Be) and by their origin (Fig. 5Bf) revealed that the tuft region was driven by local synaptic currents and was largely decoupled from the rest of the cell.

Within the place field, somatic action potentials were evoked directly by synaptic inputs (first AP in Fig. 5C) or triggered by Na⁺-spikes propagating from nearby basal dendritic branches (second AP in Fig. 5C). At the beginning of the action potential, a brief current originating from the axon can be observed (Fig. 5Cd), indicating that APs are initiated in the axon. In our model, active backpropagation of the APs is limited to the apical trunk, since the Na⁺-channels in the higher order dendritic branches have a higher voltage threshold (Ujfalussy and Makara, 2020). Indeed, back-propagating APs do not recruit local, dendritic-type Na⁺-channels in the tuft (Fig. 5Ce). In contrast, Ca²⁺-currents are recruited in the apical dendrites and contribute to the spike after depolarization (Fig. 5Ca-f).

The model neuron also fired complex spike bursts (CSBs) within its place field (Fig. 5D). When we magnified the current dynamics around the second CSB, we found that it was preceded by strong Ca²⁺-currents and dendritic Na⁺-spikes that eventually reached the soma and elicited a somatic AP (Fig. 5Da-d). In this case, the cell was driven by a mixture of synaptic and intrinsic inputs from the apical dendrites already before the CSB (Fig. 5Dd). The first AP propagated back to the dendrites and amplified the dendritic Ca²⁺- and Na⁺-currents (Fig. 5Db-d). These intrinsic currents, originating from the apical dendrite, became the dominant inward drive in the soma and drove the cell to fire further APs (Fig. 5Db-d) until the slow potassium current in the dendrites terminated the burst (Fig. 5De-f). Taken together, the extended currentscape plots provide compact and intuitive visualization of the complex current dynamics underlying neuronal activity during spatially and temporally structured synaptic inputs. Next, we studied the conditions leading to CSB firing in our model.

Input conditions for complex spike burst generation

To identify the input conditions leading to complex spike burst firing in our model, we focused on 41 CSB events from 16 simulations with different cluster location and input patterns, and contrasted them with 58 isolated action potentials (iAP; no other APs within 30 ms). First, we focussed on the membrane potential dynamics in long (>60 µm) terminal branches in the basal, oblique or tuft region. We found that during CSBs the membrane potential was highly variable across events, but was typically substantially depolarized for ∼40 ms in both tuft and oblique dendrites (Fig. 6A), indicating the presence of dendritic Ca²⁺-spikes in these branches. The large depolarizations associated with local Ca²⁺-spikes were absent in basal dendrites or during iAPs (Fig. 6A-B).

Fig. 6:

Interestingly, the depolarization started earlier in the tuft than in the oblique branches, where it reached its maximum level only after the second AP in the burst. During CSBs most tuft branches participated in Ca²⁺-spike firing (Fig. 6C), whereas the prevalence of dendritic Ca²⁺-spikes remained low during isolated spikes (Fig. 6D). In the tuft branches, but not in obliques, both the mean membrane potential and the probability of Ca²⁺-spikes diverged between CSBs and iAPs already ∼10 ms before the start of the event (Fig. 6B,D). Note that the synapse density, the ion channel mechanisms and the input statistics are identical for tuft and oblique branches, suggesting that the morphology of the cell could be the primary factor underlying the increased excitability of the tuft for Ca²⁺-spikes.

To more directly test the involvement of tuft dendrites in CSB generation, we compared the location of active synaptic clusters during CSBs and iAPs. We found that the average cluster distance from the soma was significantly larger during CSBs than iAPs (Fig. 6E). These observations confirm earlier results suggesting the special importance of tuft branches in controlling CSBs in CA1 neurons(Takahashi and Magee, 2009; Bittner et al., 2015; Grienberger and Magee, 2022; Park et al., 2025).

To further analyze the synaptic and dendritic events leading to CSB firing, we turned to our extended currentscape method and calculated the average currentscapes for both CSBs and iAPs (Fig. 7A-B). The somatic currentscapes by current type and current origin revealed the strong Ca²⁺- and Na⁺-currents driving burst firing in the soma during CSBs and being responsible for the relative depolarization after iAPs. However, to our surprise, there was almost no difference between CSBs and iAPs until the end of the first spike neither in the magnitude, nor in the origin or the composition of the somatic currents.

Fig. 7:

Since we observed a large variability in the membrane potential, we speculated that averaging across many events could potentially conceal important differences between iAPs and CSBs. Therefore, we calculated the total magnitude of each current type in the soma for each event in the 8 ms period before the first AP of iAPs or CSBs (Fig. 7C). We found a large variability in the currents underlying both iAPs and CSBs with the distributions entirely overlapping between the two types of events. We got similar results when analysing the origin of the axial currents in the soma (Fig. 7D). This indicates that it is impossible to distinguish CSBs and iAPs before the first AP from the somatic state of the neuron, and suggests that the soma can not control burst firing.

Next, we analyzed currents in the apical dendrites from the distal trunk region and from the apical tuft. Although we found some small, but significant alterations in the current types in both the trunk (not shown) and the tuft (Fig. 7E), the largest differences were in the origin of the currents driving these regions (Fig. 7F): before CSBs the tuft did not receive inward (excitatory) currents from proximal regions but was driven solely by local excitatory intrinsic and synaptic currents. In contrast, outward currents from the tuft propagated towards the trunk and obliques significantly more during CSBs than during iAPs. However, we still observed a substantial variability between the events indicating that both iAPs and CSBs could be highly diverse. We were wondering whether this diversity could be captured by a few dominant factors along which iAPs and CSBs were segregated.

In order to test this idea, we applied factor analysis to the dataset containing the magnitude of 42 current types in the 3 different dendritic domains (soma, trunk, tuft; see Methods) before the first spike of the event. The factor analysis revealed that the data indeed was low-dimensional: The first two factors together explained almost 80% of the total variance with additional factors providing negligible contribution (Fig. 8A). When we projected the data into the space defined by the first two factors, the datapoints corresponding to CSBs and iAPs tended to occupy different regions of the state space (Fig. 8B) indicating that these factors capture characteristic differences between these event types.

Fig. 8:

Next, we checked how different current types contributed to these first two factors. We found that the outward currents (GABA and K⁺) in the tuft region had the largest positive loadings, while the distal inward currents (NMDA and Ca²⁺ currents) had the largest negative loadings to the first factor (Fig. 8C; note that inward currents are negative). Thus, strong excitatory or inhibitory tuft inputs both increase the first factor. We thus interpreted the first factor as capturing the strength of the inputs to the tuft region. Similarly, the second factor had the strongest positive loadings on inhibitory, and negative weights on excitatory currents in the trunk region (Fig. 8C). Interestingly, proximal NMDA and Ca²⁺-current had opposite effects (Fig. 8C): somatic Ca²⁺-current, originating from the apical dendrites, increased both factors, while somatic NMDA currents, indicating dominant excitation from basal branches, decreased them. We interpret the second factor as capturing the strength of inputs to the oblique/trunk region of the cell.

This analysis revealed that most of the variability between iAPs is confined to a single dimension with low tuft input and a varying degree of oblique or basal excitation (Fig. 8B,D-E). Moreover, iAPs occur only exceptionally at high tuft input levels (Fig. 8F). In contrast, CSBs occupy a greater fraction of the state space with varying tuft and trunk input, but generally avoiding the lowest levels of tuft excitation (Fig. 8B, G-I). Note that occasionally highly similar currents preceded iAPs and CSBs and which of these occurred could depend on later synaptic inputs or subtle differences in the state of intrinsic conductances (Fig. 8B, D-I; see also Fig. S3).

To further elucidate the role of the inputs targeting different dendritic domains in evoking spikes or CSBs, we selectively decreased the input rates at basal, oblique or tuft branches, simulating optogenetic inhibition experiments (Bittner et al., 2015). Specifically, we randomly selected 300 input synapses targeting either of these domains, decreased their input rates by 75%, and measured the changes in postsynaptic spiking and CSB firing. We found that decreasing the tuft input had the strongest effect in reducing both the CSBs and the spike rate of the cell (Fig. S4A). This effect was not explained by the more clustered location of the inhibited synapses in the tuft (Fig. S4A), further emphasizing the crucial role of the tuft in controlling the excitability of the cell. In agreement with the experimental results (Bittner et al., 2015), the tuft inhibition caused a larger reduction in the CSBs rates (≈ 20% of the control) than in the spike rates (≈ 35% of the control; Fig. S4). Nevertheless, inhibition of basal or oblique synapses also decreased the number of CSBs, which confirms that complex spike bursts are not exclusively triggered by distal inputs in our model.

Taken together, these analyzes demonstrated that, in agreement with previous experimental results (Takahashi and Magee, 2009; Bittner et al., 2015; Grienberger and Magee, 2022; Park et al., 2025) strong distal inputs facilitate the generation of CSBs. However, in our model there was no need for exceptionally high input synchrony at the tuft synapses neither for Ca²⁺-hotspots in the dendrites to evoke CSBs. Instead, CSBs could be evoked under a large variety of synaptic and intrinsic current conditions by input to apical dendrites even at relatively weaker tuft excitation.

Discussion

Summary

In this paper, we first introduced the extended currentscape technique to quantify and visualize the effect of dendritic events on the somatic response of a biophysical model neuron. Next, we illustrated the extended currentscape method in a standard experiment probing dendritic integration in vitro. Our method could compactly and intuitively capture the sequence of dendritic events leading to a superlinear somatic response. After validating the model using spatio-temporally restricted stimuli, we turned towards in vivo-like inputs and examined dendritic activity during simulated place field dynamics. We found that although correlated activity dominated the dendritic membrane potential dynamics of the model, the extended currentscape method could identify various dendritic events underlying somatic spikelets, action potentials, and complex spike bursts. Next, we contrasted the somatic and dendritic current dynamics preceding iAPs and CSBs. Although we found a large diversity of current contributions before the two event types, the distribution of the currents were surprisingly similar before CSBs and iAPs. Finally, we applied factor analysis to the somatic and dendritic currents and showed that the large variability is dominated by two factors, corresponding to inputs to the proxial and distal apical dendrites. We found that CSBs and iAPs could occur at variable oblique and basal input levels, but strong distal inputs facilitated the generation of CSBs.

The extended currentscape

Our analysis is based on the extended currentscape method to identify the type and origin of membrane currents that drive the activity of any compartment in a neuron with spatially extended morphology. Our technique relies on simulating the neuron using a standard biophysical modeling software (Hines and Carnevale, 1997) and measuring the membrane currents and the membrane potential in each model segment at every timestep. We used Kirchhoff’s current law to partition the axial currents between neighbouring pairs of segments proportionally to the composition of the incoming or the outgoing membrane currents, depending on the direction of the axial current. We applied this partitioning algorithm recursively starting from the tip of the dendrites (leaf nodes) and proceeding towards the target segment.

In this paper, we partitioned the axial currents either by the channel type or by the dendritic region of the associated membrane currents. Partitioning the axial currents by the type of the membrane currents generating them provides a compact and intuitive summary of the sequence of events mediated by different intrinsic and synaptic conductances during various neuronal activity patterns, illustrated here for the case of theta activity, spiking and bursting in place cells. In contrast, partitioning the axial currents by their dendritic origin reveals the relative contribution of the different dendritic domains to driving the neural response at any given time. Here we distinguished only a few major regions (axon, soma, basal, oblique and tuft dendrites) within the cell, but different divisions, with more fine grain compartmentalization or combination with partitioning by current type, would be an interesting future application of the method.

Previous approaches dissecting the contribution of individual ion channel types to the somatic response of the model relied on targeted perturbation of the available biophysical mechanisms both in real neurons (Wei et al., 2001; Ariav et al., 2003; Losonczy and Magee, 2006; Smith et al., 2013; Palmer et al., 2014) and in biophysical models (Gasparini et al., 2004; Gómez González et al., 2011; Hay et al., 2011; Behabadi et al., 2012; Grienberger et al., 2014; Goetz et al., 2021; Kim et al., 2023). However, interpreting the result of these perturbations is hampered by the potential interactions between various biophysical mechanisms, especially under in vivo conditions (Smith et al., 2013; Palmer et al., 2014). For example, Losonczy and Magee (2006) found that blocking dendritic Na⁺-spikes inhibited both the fast and the slow components of the superlinear somatic response. They posited that the supralinear input-output function critically depends on the facilitation of NMDA spikes by local dendritic Na⁺spikes. Conversely, Smith et al. (2013) found the blocking NMDA receptors greatly reduced the number of local dendritic fast spikes in layer 2/3 PCs, suggesting a reciprocal interaction between NMDA– and Na⁺–spikes. Here we propose an alternative approach that does not require perturbations but calculates the contribution of ion channel types directly from the recorded axial and membrane currents, albeit currently only in biophysical models. Importantly, our method works equally well under a wide variety of input conditions, including stimulation localised to a single dendritic branch or distributed throughout the entire dendritic tree.

Biophysical model / Implications for burst firing

Our biophysical model CA1 neuron is a modified version of the model of Jarsky et al. (2005) originally developed to capture the generation and propagation of Na⁺ spikes along the apical trunk and later modified to account for the integration of synaptic inputs via Na⁺-, and NMDA-spikes in basal and oblique dendrites (Losonczy and Magee, 2006; Ujfalussy and Makara, 2020). As Ca²⁺-spikes in CA1 PCs are primarily mediated by R-type channels (Takahashi and Magee, 2009), we added Ca²⁺channels using kinetic scheme taken from (Magee and Johnston, 1995) to model dendritic Ca²⁺spikes. Our model was able to generate complex spike bursts, with associated dendritic plateau potentials upon current injection or synaptic stimulation into the dendrites. However, this model is not yet able to capture all, potentially important, properties of Ca²⁺-spikes in CA1 neurons, including their failure to generate Ca²⁺-spike upon strong somatic current injection in most cells, the delayed initiation of the Ca²⁺spikes, the inhibitory effects of the Na⁺-APs on Ca²⁺ spike initiation (Golding et al., 1999) and the generation of somatic plateau potentials (Bittner et al., 2015). Currently we are not aware of any CA1 PC model that would be able to reproduce these experimental observations (see the references in Sáray et al. 2021) indicating the need for further research revealing the biophysical mechanisms of dendritic Ca²⁺spikes.

Recognizing the incompleteness of our biophysical model, we believe that our simulations still provide several important insights into the generation of Ca²⁺-spikes under in vivo-like conditions. First, it was unclear whether increased Ca²⁺-channel density in the tuft region or special perforant path input patterns were necessary for evoking burst firing (Bittner et al., 2015; Grienberger and Magee, 2022). Our simulations demonstrated that distal tuft branches are especially well suited to promote burst firing even with uniform Ca²⁺channel density and similar input statistics along the apical subtree. This observation suggests that it is the morphology of the cell and the presence of strong inhibitory conductance load on the apical trunk that makes tuft branches more excitable than oblique dendrites.

Second, we observed a high variability between the activation of individual tuft and oblique branches during CSB firing. This observation is reminiscent of Ca²⁺ imaging data from the distal dendrites showing a variable recruitment of tuft dendrites during place field formation and subsequent traversals (O’Hare et al., 2025). Our simulations clarified that tuft Ca²⁺-plateaus are not necessarily associated with uniform tuft activation and the high variability between branches can be caused by diversity in the synaptic inputs.

Third, we found a large heterogeneity among CSBs in the input currents preceding the events: while some CSBs were elicited at unusually high tuft currents, many were triggered at lower or intermediate level of distal inputs. This suggests that the presence of distal inputs facilitates rather than controls the emergence of CSBs. This finding explains why it is possible to elicit plateau potentials in vivousing somatic current injection during BTSP induction (Bittner et al., 2015) where the same stimulus usually fails to evoke Ca²⁺-spikes in vitro (Golding et al., 1999) and makes a direct prediction that reduction of CA3 inputs to the oblique or basal dendrites would also decrease CSB rate, BTSP induction and the associated reward zone over-representation similarly to the effect of inhibiting the more distal EC inputs (Fig. S4; Bittner et al. 2015; Grienberger and Magee 2022; Fan et al. 2023). The facilitatory action of distal tuft inputs on CSBs is reminiscent of the conjunctive bursting mode, where bursts are generated by a synergistic interaction between different input streams (Larkum et al., 1999; Naud and Sprekeler, 2018) and is consistent with the conclusions of the recent study postulating that dendritic plateaus are initiated within the distal regions of stratum radiatum by strong inputs to both distal tuft and radial oblique dendrites (O’Hare et al., 2025).

Although we could frequently observe complex spike bursts in our simulations, we did not observe large amplitude (>20 mV) prolonged (> 100 ms), plateau-like depolarization events in the soma with substantially reduced AP amplitude (Bittner et al., 2015) during naturalistic synaptic inputs. However, a response more similar to plateau potentials could be evoked in our model by direct somatic current injection (Fig. S1E). We thus speculate that the generation of plateau potentials might require strong perisomatic excitatory currents. Our preliminary simulations suggest that Ca²⁺-channels added to the basal dendrites can provide this additional excitation. The biophysical mechanism and the natural synaptic input conditions that lead to CSB versus plateau potentials in CA1 pyramidal neurons is a promising subject for future research.

Limitations of the study

In order to make the partitioning algorithm self-consistent, we had to also include the capacitive current among the membrane currents in Eq. (1): It is possible that in a given segment the only inward current may be the capacitive current but the axial current is still outward. The capacitive current acts as a delay line in the membrane equation, representing the indirect effect of currents that were active earlier and have charged the membrane capacitance. In practice, we found that the capacitive current may sometimes have a large contribution to the membrane dynamics (Fig. 5Ce) and can thus mask the delayed contribution of various membrane currents. A potential extension of our method would be to also partition the capacitive current, but it is not clear how this could be achieved self-consistently and is beyond the scope of our current paper.

Our partitioning algorithm identifies only the direct contribution of the membrane currents to membrane potential changes in the target compartments, and ignores all indirect effects. This property follows from the presumption that there must be a direct chain of links between cause and effect (Pearl, 2009), in particular, there must be a continuous flow of axial current between the membrane current in a distal dendrite and the change it triggers in the target. Thus, our method is blind to indirect effects otherwise known to be present in complex dendritic trees, such as off-the-path inhibition (Gidon and Segev, 2012). Revealing the contribution of such indirect causes under complex synaptic stimulation would require causal perturbation methods that have been proposed in the context of synaptic effects (Bicknell and Häusser, 2021).

Testing the model: electrophysiology or voltage imaging

Experimental testing of the partitioning algorithm presents considerable challenges. Accurate measurements of the contribution of membrane currents to the somatic activity of a neuron require simultaneous measurement of all membrane currents and potentials, which is not feasible with current experimental methods.

Although measuring individual current contributions is not feasible in real neurons, a recent study performed whole-cell recordings in visual cortical neurons in vivo and used the systematic change of the input resistance with depolarization to estimate the contribution of intrinsic and synaptic currents to neuronal responses (Li et al., 2020). They found that during visual stimulation the intrinsic and synaptic conductances have comparable contribution to the subthreshold membrane potential changes of the cell, with intrinsic channels amplifying the synaptic response. Similar analysis on in vivo hippocampal recordings could test whether synaptic and intrinsic conductances amplify or counteract each other in vivo (Hoffman et al., 1997; Bittner et al., 2015; Epsztein et al., 2011; Cohen et al., 2017). However, in the absence of dendritic recordings, this study was unable to identify the spatial origin of synaptic and intrinsic changes within the dendritic tree.

More direct test of current propagation in real neurons would require measuring the membrane and axial currents at multiple spatial locations. As a recent step towards achieving this Meszéna et al. (2023) combined somatic patch-clamp recordings and multichannel extracellular recordings to reconstruct the spatiotemporal distribution of the membrane currents and the membrane potential of a single neuron during AP generation. Alternatively, one could use in vivo voltage imaging to monitor the membrane potential of dendrites and the soma simultaneously (Abdelfattah et al., 2023; Park et al., 2025). From the local membrane potential one can calculate the axial currents after estimating the intracellular resistivity and the dendritic cable diameters. Our method of partitioning the axial current by its origin within the dendritic tree can be applied directly to this kind of data. Therefore, such data could be used to directly test both the behaviour of the biophysical model under in vivo-like input conditions (Fig. 4F-G) and the large diversity of the origin of the input currents before single spikes and CSBs in CA1 neurons (Fig. 8).

Methods

Biophysical models

All simulations were performed with the NEURON simulation environment (version 8.2) embedded in Python (version 3.9). Code for simulating the biophysical model, preprocessing, axial currentscape calculation and visualization are available at https://github.com/bencefogel/currentscapes_invitro_demo and https://github.com/bencefogel/currentscapes_invivo_demo.

CA1 neuron model

We used a modified version of a previously developed CA1 PC model (Jarsky et al., 2005; Ujfalussy and Makara, 2020), maintaining the same ion channel distributions and densities for Na⁺, K⁺ and proximal and distal type K⁺ channels as described in Ujfalussy and Makara (2020) (Fig. S5).

To capture dendritic Ca²⁺ spikes and somatic CSBs, we added an R-type Ca²⁺ channel to the model. The Ca²⁺ channel kinetics were based on the steady-state activation and inactivation curves fitted by Magee and Johnston (1995), with a slight modification to allow for a larger window current (Fig. S1A). The Ca²⁺ channel was expressed in the apical dendritic tree with a uniform distribution and an ion channel density of 0.006 S/cm².

Furthermore, we added a high-voltage activated, slow K⁺ channel that was used to simulate the combined effect of multiple K⁺ channels including Ca²⁺-activated K⁺ channels, contributing to the repolarization after the Ca²⁺ spikes. The slow K⁺ channel was expressed in the apical dendritic tree with a uniform distribution and an ion channel density of 0.001 S/cm².

The CA1 cell received inputs from 2000 excitatory and 200 inhibitory presynaptic neurons. We simulate the activity of excitatory and inhibitory neurons during the movement of a mouse on a 2-m long circular track with constant velocity of 20 cm/s. Each presynaptic excitatory neuron had a single idealised place field and represented the CA3 inputs received by the postsynaptic CA1 neuron. Excitatory synapses were placed randomly with a uniform distribution on the entire dendritic tree, except 240 inputs, active in the middle of the track, that were selected to form functional synaptic clusters. There were a total of 12 clusters, and each cluster had 20 synapses. Synaptic clusters were placed on the middle of terminal dendritic branches longer than 60 µm with 1 µm distance between the synapses. Inhibitory presynaptic inputs were modulated by theta oscillation but were not spatially tuned. These synapses were divided into two groups with 80 synapses targeting the soma and the apical trunk and the remaining synapses distributed randomly along the entire dendritic tree.

The model included AMPA and NMDA excitation and slow and fast GABA inhibition. Synaptic parameters were kept identical to those described in Ujfalussy and Makara (2020) (see Fig. S6 for kinetic parameters). To induce reliable firing within the neuron’s place field, we increased the AMPA (from g_max = 0.6 nS to g_max = 1 nS) and NMDA (from g_max = 0.8 nS to g_max = 1.2 nS) conductance associated with the clustered synaptic inputs. Inhibitory synapses had a maximal conductance of g_max = 0.2 nS.

We simulated place cell activity for 10 seconds as the hypothetical animal completed a single lap on the track. A total of 16 simulations were run, each with different random configurations of synapse placement and presynaptic input patterns. During each simulation, we recorded the membrane potential (mV), intrinsic (mA/cm²) and synaptic currents (nA) of all segments of the neuron (N_sec = 161 sections and N_seg = 1529 segments in total, including the 2N_sec internal segments; see below). This was the raw dataset we used to calculate the currentscapes. We also saved the model connectivity structure between segments together with the axial resistance values (MΩ) to calculate axial currents offline. Finally, we measured the area (µm²) of each segment, which is needed during the preprocessing of membrane currents.

Preprocessing

The goal of preprocessing was to convert the raw dataset, saved by NEURON into two tables that we can use to calculate the partitioning of the axial currents. The first table contains the membrane currents for each segment and for each time step. The second table contains the axial currents between the segment pairs.

During the simulation, state variables were computed using NEURON’s built-in multi order variable time-step integration method. For subsequent preprocessing and currentscape calculation, the recorded output vectors were downsampled to 5 kHz.

Preprocessing membrane and synaptic currents

In the NEURON software, neurons are modelled as a series of sections, which represent unbranched lengths of continuous cable, such as dendrites or axon. These sections are connected together to form an acyclic graph, according to the morphology of the neuron. Each section is divided into multiple segments of equal length.

NEURON uses a normalized distance to express locations along a section, where 0 represents the start (closest to the parent section) and 1 represents the distal end. To compute membrane dynamics, NEURON calculates membrane potentials and membrane currents at discrete positions, known as nodes, located at the center of the segments. The nodes are evenly spaced internal points within a section. The number of segments within a section was chosen to ensure accurate description of signal propagation within the cable with reasonable computational efficiency. As a rule of thumb, we found that axial current partitioning is accurate if the magnitude of the membrane currents within the segment is small compared to the magnitude of the axial current, which was on average ≈ 1/17 for non-terminal branches in our case.

All segments that belong to the target section were merged for further analysis (that is, the membrane currents of the segments were summed and no axial current was calculated within the target section). The connections between the segments of the target section and its child branches were reassigned to the segment of the merged target.

In addition to the internal nodes, there are external nodes located at the 0 and 1 ends of each section. These external nodes are only used to connect segments, but no membrane mechanisms (synapses or ion channels) are associated with these external nodes (see below).

Each synapse is represented by a conductance-based mechanism that injects current into the postsynaptic compartment when triggered via a NetCon object. Since each NetCon delivers spikes independently, and multiple synapses can target the same segment, synaptic currents (in nA) were summed per segment and per synaptic type (AMPA, NMDA, GABA_fast, GABA_slow) to accurately capture the total synaptic input.

To ensure consistency in subsequent calculations, intrinsic currents were converted from units of mA/cm² to nA. This was done by multiplying the recorded current by the corresponding segment area (in µm²) and applying a scaling factor of 10⁻². Finally, intrinsic and synaptic currents were combined into a single data structure, where each row corresponds to a unique segment and current type pair, and each column is a simulation time point.

Preprocessing axial currents

To compute axial currents we first created a list with the adjacent segment pairs between which axial currents will be calculated. Connections between sections (branches) are implemented in NEURON through a pair of external nodes: one at the distal end of the parent section and another at the proximal end of the child section. Although these nodes formally belong to different sections (parent and child), they have identical membrane potential as they represent the same point within the cable, so we merged them to remove the duplicates from our primary dataset.

The axial current between a parent and child segment is computed using Ohm’s law:

Here V_parent and V_child denote the membrane potentials of the parent and child segments, respectively, and R_i denotes the axial resistance between them.

Extended Currentscape Calculation

The extended currentscape calculation requires the following inputs: the preprocessed membrane and axial currents across all segments, a specified target segment, and a defined partitioning strategy, such as current type-specific or region-specific partitioning.

The first step of calculating the extended currentscape was to separate positive and negative membrane currents. Certain currents, such as capacitive currents, can be both inward and outward at the same time in different compartments. Without separation, such positive and negative currents would partially cancel each other, reducing the contribution of the current type during the extended currentscape calculation. This effect is even more pronounced when partitioning the axial currents by the origin (see below): Strong positive and negative currents from two different basal branches can cancel each other, largely concealing the basal contribution to the somatic activity. By treating the positive and negative membrane currents independently, we ensured that the outward and inward membrane current components of the target’s axial current are calculated separately.

The core algorithm underlying the extended currentscape calculation is the axial current partitioning algorithm. By default, this algorithm decomposes the axial current of a target compartment into underlying membrane current components. However, alternative variants of the partitioning algorithm can be applied if the membrane current data structure is reindexed and recalculated. In region-specific partitioning, membrane currents are reassigned based on their location within different subcellular regions (for example, axon or dendrite). These different approaches can also be combined in a flexible way to measure e.g., the somatic influence of NMDA channels in a particular dendritic branch. In this study, we demonstrate two such variants: Membrane current type-specific partitioning and a combined partitioning where we distinguished 5 subcellular regions (axon, soma, basal dendrites, trunk+oblique dendrites and tuft dendrites) and 2 current types (synaptic and intrinsic currents). The intrinsic category incorporated all non-synaptic currents, including capacitive and passive currents.

Following preprocessing, including reindexing and recalculation of membrane currents when required by the region-specific partitioning, a directed graph was constructed at each time step. The nodes of the graph represent individual neuronal segments and the directed edges indicate the direction of axial currents flowing between them.

To determine the partitioning order, we performed a depth-first search (DFS) on the directed graph starting from the target segment along each dendritic branch separately. The search proceeds through the directed graph until it encounters collision points, where axial currents reverse direction. At these points, current propagation is blocked, and distal nodes beyond the collision points are excluded from further calculations, as they do not directly influence the target compartment. We create two separate subgraphs from the target segment: one for the inwardly flowing and one for the outward axial currents. The two subgraphs identified by this search define the structures upon which axial current partitioning is performed. Partitioning is calculated independently for the two subgraphs, with the inward (outward) axial graph partitioned according to the inward (outward) membrane currents, respectively.

Partitioning is carried out recursively, beginning from the leaf nodes of the subgraph and progressing toward the target compartment. At each time step, the axial current is partitioned into inward (negative) or outward (positive) membrane current components depending on the direction of the axial current flow. Specifically, for a given time point t, the positive or negative membrane currents I ^m,ref of the reference compartment are first normalised by their sum. The resulting fractions are multiplied by the magnitude of the axial current I ^ax flowing between the reference and target compartment. The contribution of membrane current of type i to the axial current can be expressed as:

The partitioned axial current components I ^part(t) are then added to the corresponding membrane currents of the parent section according to:

where is magnitude of the membrane current type i measured in the parent segment. This recursive procedure is repeated for each pair of nodes along the subgraph until the target compartment is reached.

The final output of the extended currentscape calculation consists of two data structures, each including the partitioned axial currents of the target compartment: one containing the inward membrane currents and the other containing the outward membrane currents in each time step.

Data Analysis

CSBs and isolated spikes

CSBs were detected by first smoothing the somatic membrane potential with a 20 ms Gaussian kernel and then detecting the crossing of a -55 mV threshold of the smoothed subthreshold Vm (Fig. 4B). The start of the CSB was defined as the time of the first spike after crossing this threshold. Isolated spikes were defined as action potentials with a distance of at least 30 ms from the closest AP.

Ca²⁺-spike detection

Dendritic Ca²⁺-spikes (Fig. 6C-D) were detected by smoothing the dendritic membrane potential with a 4 ms Gaussian kernel and then detecting the crossing of a -35 mV threshold of the smoothed subthreshold Vm (Bittner et al., 2015).

PCA of dendritic voltage

We recorded the membrane potential in the middle segment of all dendritic and somatic compartments of the model and calculated the covariance matrix of the z-scored dendritic voltages in each run separately. Fig. 4C shows the average of these 16 covariance matrices. PCA was then performed by calculating the eigenvalues and eigenvectors of the covariance matrices separately (Fig. 4D-E).

Dendritic events

Dendritic Vm peaks were detected by the find_peaks function from the signal modul of the scipy python package with parameters height=-30 mV and prominence=30 mV. Fig. 4F shows the histogram calculated from the time difference between dendritic events and the closest somatic action potential.

Dendritic domains

In Fig. 6A-B we show the average membrane potential in different dendritic domains during CSBs and isolated APs. We included all terminal branches with L > 60 µm in this analysis (basals: N_basal = 30, average distance from the soma: D_basal = 135 µm; obliques: N_oblique = 19, D_oblique = 313 µm; N_tuft = 9, D_tuft = 729 µm;).

Average currentscapes

To calculate the average currentscapes in Fig. 6 we first calculated the percentage of the different current shares of each event and then averaged the percentages across the events.

Factor analysis

Factor analysis in Fig. 8 was performed using the partitioned somatic, trunk, and tuft current types before the first spike of the CSB or iAP events. Specifically, we averaged the magnitude of the incoming current types in the interval t = {−10, −2} before each event. We performed the averaging separately for the positive and negative currents for all types with nonzero variance (e.g., NMDA currents are always inward, so their contribution to outward currents has zero mean and variance). For this analysis, we used the raw current magnitudes instead of the percentages. We performed factor analysis using the FactorAnalysis function using sklearn’s decomposition module.

Simulating optogenetic inhibition

We simulated targeted optogenetic inhibition of different dendritic domains (Fig. S4) by randomly selecting a similar number (E[n] = 300) of presynaptic inputs and decreasing their firing rate to 25% of their baseline value. In different experiments, we selected inputs targeting either basal, oblique, or tuft dendrites.

Supplementary figures

Fig. S1:

Fig. S2:

Fig. S3:

Fig. S4:

Fig. S5:

Fig. S6:

Acknowledgements

We thank Judit K Makara, Szabolcs Káli and János Brunner and all members of the Neuronal Signalling and the Biological Computation Research Groups for discussions and for their comments on the manuscript. B.B.U. was supported by the National Research, Development and Innovation Office of Hungary (FK-125324) and by the Hungarian Research Network.