https://orcid.org/0000-0003-3913-5650
Figures
Abstract
The biophysical properties of neurons not only affect how information is processed within cells, they can also impact the dynamical states of the network. Specifically, the cellular dynamics of action-potential generation have shown relevance for setting the (de)synchronisation state of the network. The dynamics of tonically spiking neurons typically fall into one of three qualitatively distinct types that arise from distinct mathematical bifurcations of voltage dynamics at the onset of spiking. Accordingly, changes in ion channel composition or even external factors, like temperature, have been demonstrated to switch network behaviour via changes in the spike onset bifurcation and hence its associated dynamical type. A thus far less addressed modulator of neuronal dynamics is cellular morphology. Based on simplified and anatomically realistic mathematical neuron models, we show here that the extent of dendritic arborisation has an influence on the neuronal dynamical spiking type and therefore on the (de)synchronisation state of the network. Specifically, larger dendritic trees prime neuronal dynamics for in-phase-synchronised or splayed-out activity in weakly coupled networks, in contrast to cells with otherwise identical properties yet smaller dendrites. Our biophysical insights hold for generic multicompartmental classes of spiking neuron models (from ball-and-stick-type to anatomically reconstructed models) and establish a connection between neuronal morphology and the susceptibility of neural tissue to synchronisation in health and disease.
Author summary
Cellular morphology varies widely across different cell types and brain areas. In this study, we provide a mechanistic link between neuronal morphology and the dynamics of electrical activity arising at the network level. Based on mathematical modelling, we demonstrate that modifications of the size of dendritic arbours alone suffice to switch the behaviour of otherwise identical networks from synchronised to asynchronous activity. Specifically, neurons with larger dendritic trees tend to produce more stable phase relations of spiking across neurons. Given the generality of the approach, we provide a morphology-based hypothesis that explains the differential sensitivity of tissue to epilepsy in different brain areas and assigns relevance to cellular morphology in healthy network computation.
Citation: Gowers RP, Schreiber S (2024) How neuronal morphology impacts the synchronisation state of neuronal networks. PLoS Comput Biol 20(3): e1011874. https://doi.org/10.1371/journal.pcbi.1011874
Editor: Lyle J. Graham, Université Paris Descartes, Centre National de la Recherche Scientifique, FRANCE
Received: June 18, 2023; Accepted: January 30, 2024; Published: March 4, 2024
Copyright: © 2024 Gowers, Schreiber. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All data and code are publically available: https://github.com/rpgowers/morph-excite-code.
Funding: This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no 864243 to SS) and from the Einstein Foundation Berlin (EP-2021-621 to SS). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Network states are instrumental for neural computation: they correlate with the behavioural state in healthy animals, such as in central pattern generator circuits for movement [1, 2], and are also often altered in neuropathologies [3–5]. Recent theoretical and experimental work highlights that it is not only the connectivity between neurons which plays a role in determining network behaviour, but that neuron-intrinsic properties also exert an influence [2, 6, 7]. Mechanistically, these influences arise not only from indirect effects on connectivity, such as plastic changes in synaptic transmission or modulation of plasticity rules, but also from direct effects on the very core of processing by a single neuron: the qualitative dynamics of action-potential generation that define a neuron’s excitability class as described by Hodgkin [8]. Along these lines, weakly coupled inhibitory neurons with class 1 cell-intrinsic excitability do not foster synchronous network states [9], while neurons arranged in the same network topology with homoclinic-type action-potential generation favour in-phase synchronisation [7, 10]. Among the parameters that alter the cellular excitability class, we find those that directly affect ion channel dynamics, including channel composition, extracellular ion concentration, and modulators such as temperature [6, 7, 10–14]. However, these are not all.
Here, we explore the implications of a neuronal property that has received comparatively less attention in the context of network dynamics, presumably due to its more inflexible, less variable nature: neuronal morphology. Although previous work has shown that even passive dendritic arbours can play an important role in processing inputs [15–18], their effect on the network state has not yet been explored. We demonstrate that differences in the extent of dendritic arborisation alone are sufficient to induce network behaviour that is either synchronised, with stable phase relationships between neuronal firing, or asynchronous. Differences in neuronal morphology can thus contribute to the differential susceptibility of neuronal tissue to synchronisation, which is likely to be of relevance also for pathological phenomena such as epileptiform activity or spreading depolarisation.
Our approach exploits the fact that the dynamics of regularly spiking neurons with all-or-none action potentials come in at least three qualitatively different flavours—hereafter referred to as dynamical spiking types—corresponding to the three different mathematical bifurcations that underlie the onset of spiking at threshold. Despite the large diversity of neuronal properties, including a zoo of ion channels that shape a cell’s conductance portfolio as well as extrinsic modulators such as ionic concentrations and temperature, regular spiking in neurons is initiated by either a saddle-node on an invariant cycle (SNIC) bifurcation, a subcritical Hopf bifurcation (in conjunction with a fold of limit cycles bifurcation), or a saddle homoclinic orbit (HOM) bifurcation. The different spike onset dynamics of these dynamical spiking types have been shown to influence the temporal relationships of spikes across neurons in weakly connected networks [19–21]. Indeed, it is the combination of synaptic and cellular voltage dynamics that determines the state of the network, with the influence of cellular voltage dynamics being particularly pronounced in fast synaptic transmission. In this context, passive properties such as neuronal morphology also have an influence on the effective dynamics of spike generation. Larger dendritic branches induce a larger leak, which—as we demonstrate in this study—is sufficient to alter the dynamical spiking type, and consequently, the network synchronisation state.
Given the considerable heterogeneity between dendritic arbours—even in neurons of the same class and layer [22–24]—deciphering their influence on properties of neural computation and network states is a worthy endeavour. In this study, we methodically demonstrate that the morphology of dendritic arbours can, via their effect on a cell’s input impedance, be tuned to promote (de)synchronisation of network activity. Using a neuron model of an active conductance-based soma attached to a passive dendrite, we first recapitulate how the input impedance differs between single-compartment and dendrite-and-soma models. We then show how to calculate the local bifurcation structure (and hence the dynamical spiking type) analytically when morphology is included. Of particular interest is the Bogdanov-Takens (BT) bifurcation, which acts as an organising centre for different spiking onset bifurcations and hence dynamical types [13, 25]. In the next step, we show that the change in spike onset bifurcation from the dendritic load is reflected in the spiking timing response near the onset bifurcation. We further verify our results in detailed, anatomically reconstructed neuron models with varying degrees of dendritic arborisation, demonstrating that our results from the simplified dendrite-and-soma model are quantitatively precise, generalise and are widely applicable. Finally, we demonstrate the strong effect dendritic arborisation can have on network synchronisation via network simulations of multicompartment neurons with different dendritic extents. This exemplifies the potential contribution of cellular morphology to the susceptibility of neuronal tissue to specific network states.
Results
Passive impedance properties
Network behaviour is influenced by cellular properties of the constituent neurons via the dynamical spiking type they grant to each neuron. This is because distinct dynamical spiking types result in distinct timing sensitivities and consequently synchronisation properties of the network. Bringing neuronal morphology into play in a network environment therefore consists in first understanding how morphology impacts the dynamical spiking type.
Given that varying the speed at which the neuron’s membrane potential changes has been shown to induce all known dynamical switches for regular spiking [10, 13], we start our investigation by analysing how dendrites affect temporal filtering of the neuron’s membrane potential. To differentiate the voltage filtering effects caused by a dendritic arbour versus those captured by an isopotential point neuron, we first calculate the passive impedance properties of the single-compartment (S) and dendrite-and-soma models (DS). While the passive input impedance of an S model neuron is a first-order low-pass filter, spatial neuron models yield qualitatively different impedance profiles that depend on the dendritic properties [26, 27].
As we show in Fig 1, the DS model has an active soma attached to a passive dendrite which represents the equivalent cable of a branched dendritic arbour [28, 29]. If the active conductances in the soma have the same valued parameters and equations as the S model, then the two models differ only in their passive properties. The passive properties of the dendrite can be parametrised in terms of its length L, electrotonic length constant λ, passive time constant τδ, and dendritic dominance factor ρ [29] (the ratio between the characteristic dendritic conductance and the somatic leak conductance).
In the S model, Gin is equivalent to the lumped leak conductance GL, while for the DS model, Gin is the sum of the passive somatic (Gσ) and dendritic (Gδ) contributions. When we increase Gin in the DS model, Gσ is kept constant while Gδ is increased.
In this study, we compare different degrees of neuronal arborisation as captured by the dendritic contribution to the passive input conductance (Gδ). Such a comparison can be thought of either as a means of comparing arbours from different neurons, or as the effect of dynamically changing the dendritic properties of an existing neuron. A neuron with more dendrites radiating from the soma will have a thicker equivalent cable and thus contribute more to the total passive input conductance Gin. The passive input conductance is given by the constant level of external input current Iext required to change the voltage from its steady state value (when all active conductances are zero) by an amount Δv (1)
In the DS model, we apply Iext to the soma and use the voltage change Δv at the soma for Gin. In each case Gin is the sum of its somatic (Gσ) and dendritic (Gδ) contributions. For a dendrite of effective length ℓ = L/λ
(2)
When ℓ ≫ 1, we can treat the dendrite as semi-infinite and the passive input conductance of the DS model becomes
(3)
For the S model Gin = GL, the total leak conductance, allowing for a straightforward comparison between the two morphologies. The total capacitance of the S model is denoted by Cm.
While it is possible to adjust the dendritic parameters (ρ, λ, L) to match Gin for the S and DS models, this is not possible for the passive input admittance Yin (the reciprocal of the input impedance, ). Denoting the angular frequency of the input signal as ω, for the S model Yin is simply given as a first-order filter
(4)
while for the semi-infinite and finite DS models with somatic capacitance Cσ we have
(5)
(6)
where γ(ω) is defined as
(7)
Here setting Cm = Cσ or choosing any other frequency-independent value of Cm will lead to and
differing for almost the whole frequency range.
This mismatch in the passive input admittance encapsulates the difference between the S and DS models, and we will explore the implications that it has on the neuronal dynamical spiking type. For the semi-infinite DS model, Yin can be fully parametrised by Gin and τδ. From the form of Yin, we see that τδ is a key parameter for comparing the DS and S models. When τδ = 0, γ = 1 at all frequencies and hence and
become equal to
in this limit for Gin = GL. Thus the voltage dynamics of vσ in the DS model become identical to the dynamics of v in the S model when τδ = 0.
We illustrate the differences between the passive input admittances of the S and DS models via its more commonly measured reciprocal, the input impedance Zin. Fig 2A shows that for any given Gin, |Zin| is higher in the S model than the semi-infinite DS model at all non-zero frequencies. In Fig 2B, we see that |Zin| decreases at all non-zero frequencies when we increase τδ in the semi-infinite DS model. In the finite DS model, Fig 2C shows that if we hold Gin constant while decreasing ℓ, then |Zin| decreases at all non-zero frequencies.
(A) When one fixes the input conductance Gin between the S model and a semi-infinite DS model, the magnitude of the input impedance Zin will be higher for the S model for non-zero frequencies. τδ = 10 ms. (B) Increasing τδ of the DS model decreases |Zin| at all non-zero frequencies. Gin = 8 nS. (C) Decreasing the effective dendritic length ℓ while maintaining the same Gin also decreases |Zin| at all non-zero frequencies. τδ = 10 ms.
In our analyses of the DS model, we use Gin and τδ to compare dendritic arbours rather than the underlying electrophysiological parameters. Gin and τδ are defined in terms of Eqs 3 and 15. We chose these parameters because they can be readily measured experimentally by looking at the transient and long-term response to a step-current input. However, differences in these parameters have several biophysical interpretations. Firstly, increasing the dendritic diameter d increases Gin without affecting τδ. Changing τδ, on the other hand, requires changing the per-area membrane properties of the arbour. While τδ can be modified without changing Gin by increasing the dendritic membrane capacitance per unit area cδ, there are conflicting findings on how much neuronal membrane capacitance per area varies [30, 31]. Hence differences in τδ are more likely to arise from differences in the dendritic per-area conductance gδ. In either case, Gin can change from both the dendritic size and the conductive membrane properties, while τδ can only change from the per-area membrane properties.
While the two biophysical changes discussed so far are typically set once the neuron has developed and will not vary over short time scales, it is important to recognise that Yin can be changed dynamically. One example of this is that increased synaptic bombardment distributed across the dendrite will increase the per-area leak conductance of the dendrite gδ. This in turn will increase both Gin and decrease τδ [32–34]. Increasing d or decreasing gδ also increases the length constant λ (Eq 15). For short dendrites, this is an important consideration, however for long dendrites with inputs applied to the soma, λ does not appear in Yin or in any of the calculations of the local bifurcations as we show in the Methods section.
Bifurcation structure
We now make the DS model active by giving the soma the dynamics of the Morris-Lecar model [35] while keeping the dendrite passive. Dendritic arbours in general will contain various active channels [36, 37], yet we limit ourselves here to passive dendrites throughout this article for reasons of both mathematical tractability and to focus on the effects of morphological extent. The effects of active dendritic arbours are detailed further in the Discussion, and some simulations with an active dendrite are provided in Supporting information: S1 Text.
The Morris-Lecar model was chosen because it has a single time-dependent gating variable, making it one of the simpler conductance-based neuron models, and also because it has been extensively studied for its ability to change the dynamical spiking type upon parameter variation [12, 13, 19, 38, 39]. We chose the default class I excitability parameter set of the Morris-Lecar model with Gσ = 2 nS, where details of the model’s dynamics and parameters are stated in Table A in S1 Text. Other higher-dimensional conductance-based neuron models with class I excitability, such as the Wang-Buzsáki model [40], could have been chosen and are amenable to the analysis presented here and would yield similar results. The analyses of these models in other studies [6, 7, 10, 13], along with our mathematical derivations relying on the general bifurcation structure in these models (detailed in S1 Text), allow us to claim that the results we obtain from the lower-dimensional Morris-Lecar model will be transferable to higher-dimensional neuron models.
Given that the DS model differs from the S model in terms of its passive input impedance, and that the input impedance can be fully described in terms of (Gin, τδ) for the semi-infinite dendrite, it naturally follows that we should choose (Gin, τδ) as bifurcation parameters along with Iext in assessing the effect of morphology on the dynamical spiking type. One can interpret variations in the bifurcation parameters as either varying the passive properties of an existing dendritic arbour or alternatively as a means of comparing dendritic arbours belonging to different neurons.
Since the passive leak conductance in the soma Gσ = 2 nS, values of Gin above 2 nS denote the conductance load added to the neuron by the dendrite. As in our analysis of the passive impedance properties, here Iext is always applied to the soma; analytical treatment of Iext applied at any dendritic location is described in S1 Text. In short, applying Iext along the dendrite systematically increases the value of Iext for all bifurcations but leaves the other bifurcation parameters, including Gin, unchanged. Thus all transitions in the spiking onset type have the same values of Gin and τδ irrespective of the applied current location; only the values of Iext change by a known amount.
An equivalent approach to varying the input conductance Gin and applied current Iext can be taken by using the relative dendritic input conductance, ρ = Gδ/Gσ, and the applied potential, μext = Iext/Gσ, as bifurcation parameters, as detailed in S1 Text. We show the absolute quantities Gin and Iext here for intuitive clarity, but note that it is the relative size of the dendritic input conductance and external input current in comparison to the somatic leak conductance Gσ that determines the shape of the bifurcation diagram.
In Fig 3, we show several two-dimensional bifurcation diagrams in terms of (Iext, Gin) plotted for various values of the third bifurcation parameter τδ. As discussed earlier in for the analysis of the input impedance, the case of τδ = 0 is equivalent to the S model where Gin = GL. There are two saddle-node (SN) bifurcations at lower and higher Iext, which we will hereafter refer to as “lower” (dashed) and “higher” (solid line) respectively. These two SN bifurcations converge at the cusp as Gin increases. Three fixed points exist for , one stable and two unstable. When spiking onset is caused by a SNIC bifurcation, this occurs on the higher SN bifurcation. These bifurcations do not vary with τδ, and hence occur at exactly the same values of (Iext, Gin) in the S and DS models.
(A) Two-parameter local bifurcation diagram in terms of (Iext, Gin) of the Morris-Lecar DS model for various dendritic time constants τδ. At τδ = 0, the DS model is equivalent to an S model with a leak conductance of Gin. Increasing τδ shifts the BT bifurcation and shrinks the Hopf bifurcation curve. (B) shows that initially increasing τδ moves the BT point to higher Gin until it reaches the cusp at the BTC point. Increasing τδ beyond moves the BT point to the lower SN bifurcation curve (dashed) and thus decreases
. The Hopf bifurcation emerging from the BT point also switches from subcritical to supercritical when τδ passes the BTC point. (C) shows the Hopf bifurcation emerging more clearly from a BT point when we measure the external input current relative to the higher saddle-node value ISN,high, while (D) shows the Hopf bifurcation when we measure the external input current relative to the lower saddle-node value ISN,low. The saddle-node and cusp bifurcations are depicted in black because they are unaffected by changes to τδ.
At a Bogdanov-Takens (BT) bifurcation, a saddle-node, Hopf and homoclinic bifurcation converge to a single point [41]. Since bifurcations which can produce spiking onset meet at this bifurcation, we can often find different spiking onset types by varying parameters near the BT point. Given the multiple conditions required for the bifurcation, two bifurcation parameters are necessary to locate the BT point. For τδ = 0, the BT point is located on the higher SN bifurcation (solid line) and has a subcritical Hopf bifurcation emerging from it. This Hopf bifurcation permits class II excitability. Thus the BT point has often been used heuristically as separating class I and class II excitability. Increasing τδ from zero initially shifts the BT bifurcation to higher Gin until it reaches the cusp at ms (hereafter denoted as the BTC point). Increasing τδ beyond
moves the BT point onto the lower SN bifurcation curve (dashed) and
now decreases as τδ increases. As the BT point passes the cusp as τδ increases, the criticality of the emerging Hopf bifurcation switches from subcritical to supercritical at the BTC point. The transfer of the BT point from the higher to lower SN bifurcation with increasing τδ is detailed more clearly in the inset and Fig 4.
Here we have taken Iext as the onset current for every value of (τδ, Gin). Points indicate values of (τδ, Gin) we later use for the spike timing response. At the saddle-node loop (SNL) bifurcation, the dynamical spiking type switches from SNIC to HOM, while Hopf onset becomes possible at the BT bifurcation. Increasing τδ above zero increases both and
until the BT and SNL bifurcations meet the cusp at the BTC point at
. The difference between
and
decreases as τδ increases, meaning that the range of Gin for which HOM onset exists becomes smaller. For
, the SNL bifurcation no longer exists and
decreases. Furthermore, the Hopf onset for
switches from subcritical to supercritical. Gin for subcritical Hopf onset increases with τδ until
, while Gin for supercritical Hopf onset decreases with τδ throughout the whole range.
When the Hopf bifurcation that emerges from the BT point is subcritical, it can be switched to supercritical by increasing Gin. This criticality switch moves to lower Gin as τδ increases towards . At all τδ, there is a fold of Hopf bifurcations when Gin becomes sufficiently high, and spiking is no longer possible at any applied current. This is shown as the maxima of the Hopf bifurcation curves in Fig 3A. The value of Gin for the fold of Hopf bifurcations decreases as τδ increases. The changes to the BT and Hopf bifurcations with τδ mean that we would expect the transition between class I and class II excitability to occur at higher Gin in the DS model for
.
As the HOM bifurcation that emerges from the BT point involves an unstable limit cycle for the ML model, this means that it is not responsible for HOM onset and we do not show it here (see [42] for mathematical details). Instead, to identify the input conductance for the SNIC to HOM switch, we must look at the saddle-node-loop (SNL) bifurcation [10, 43–45]. At the SNL bifurcation, an SN and HOM bifurcation merge as the homoclinic orbit created by the HOM bifurcation goes through the saddle-node caused by the SN bifurcation. This means that we can look for the switch from SNIC to HOM by looking along the SN bifurcation curve.
In Fig 4, we show the bifurcation types responsible for dynamical switching as a function (τδ, Gin). For all τδ, , with SNIC onset for
. The limit cycle formed by the HOM bifurcation in this case contains all three fixed points of the system (big-HOM). Since the SNL bifurcation converges to the BTC point, it therefore makes sense that
also increases with τδ when
. Indeed, this is what we see in Fig 4, with not only
increasing with τδ but also the conductance difference between the BT and SNL points decreasing before eventually converging at
. Thus both the SNIC → HOM switch and the Gin range for homoclinic onset are affected by the dendritic time constant. For
, we do not find an SNL bifurcation for the range of τδ we test here (up to τδ = 20 ms).
Strictly speaking there can exist another global bifurcation in the region for which a fold of limit cycles (FLC) appears with a bifurcation current
, thus making the FLC the onset bifurcation [13, 42]. While this bifurcation can determine the switch between class I and class II excitability, we neglect to show it for the following reasons: (1) the spike timing perturbation response for FLC onset in this region is extremely similar to HOM onset; (2) though HOM onset in theory has class I excitability, its f–I curve is extremely steep at
, making it hard to distinguish from class II excitability; and (3) the difference between
and
is typically extremely small [46] and thus
gives an accurate approximation of the onset current in this region. Thus due to simplicity and the high degree of functional similarity, we refer to the onset type for the whole of
as being HOM.
We see that some bifurcations are affected by τδ while others are not. Both the SN and cusp bifurcations are calculated from the existence of fixed points alone, and so the timescale of the dynamics (such as from τδ) will not affect them. On the other hand, the Hopf and BT bifurcations involve the stability switch of fixed points, for which knowledge of the timescale of the dynamics is necessary. For the SNL bifurcation, one can extend the reasoning in [10] to show that changing the timescale of the dynamics will break any existing homoclinic orbits at a given Gin. As a result, one would expect that changing τδ changes the value of .
The switching in dynamical spiking type from SNIC to HOM to Hopf with increasing dendritic load is a pattern that is expected to be found in all conductance-based neuron models which start from a SNIC onset soma and have the BT point on the higher SN bifurcation branch. This follows from the analysis of two-dimensional conductanced-based models by Hesse et al [10], in which it was shown that a SNIC bifurcation is always enclosed by two SNL bifurcations when a timescale parameter (such as τδ or Gin) is varied, and that the big SNL bifurcation is reached first when the timescale is shortened (e.g. increasing Gin or decreasing τδ), followed by the BT bifurcation. Taken together, this means that for the two-dimensional conductance-based neuron model, and hence the switching between dynamic spiking types of SNIC → HOM → subcritical Hopf with increasing Gin (and equivalently, with decreasing τδ). A related explanation using the properties of the BTC bifurcation is also given by Kirst [13, 46]. Since this explanation is valid for conductance-based neuron models with an arbitrary number of dimensions, we anticipate that the argument by Hesse et al [10] could be suitably extended to apply to all conductance-base neuron models.
Our bifurcation analysis indicates that switching between dynamical spiking types can be induced by increasing the dendritic arborisation as parametrised by Gin. The switching between dynamical spiking types is qualitatively similar to increasing GL in the single-compartment model, however the Gin values for the SNIC → HOM switch and subcritical Hopf onset increase with τδ, whilst the value of Gin for supercritical Hopf onset decreases with τδ.
Phase-response curves
The dynamics of a neuronal spike affects how the neuron responds to external perturbations, which can come from chemical synapses, gap junction coupling, local field potentials, or externally applied currents. This has been found in both experimental [47, 48] and in modelling studies [19, 49]. The change in spike time caused by a perturbation to a tonically spiking neuron is described by the phase-response curve (PRC). In many cases, such as when neurons are weakly coupled or subject to weakly correlated inputs, one can use an individual neuron’s PRC to infer synchronisation conditions and the overall network state [19–21]. To see how the different dynamical spiking types affect the neuron’s spiking response, we calculated the PRCs for the S and DS models for a range of (Gin, τδ). We chose the (Gin, τδ) values to be in the neighbourhood of the SNL, BT, and cusp bifurcations, as shown by the coloured points in Fig 4.
At Gin = 3 nS, the Morris-Lecar neuron has SNIC spiking onset for all τδ. Hence Fig 5A shows symmetric positive-valued at this input conductance. At Gin = 4.7 nS (Fig 5B), the neuron is operating in the HOM regime for lower τδ and thus we see asymmetric positively-skewed PRCs associated with HOM onset [10, 50, 51], while for higher τδ the SNL bifurcation has not yet been reached and we still have symmetric SNIC PRCs. Further increasing Gin (Fig 5C) causes the neuron with τδ = 10 ms to pass its SNL point, inducing an asymmetric HOM PRC, while the PRC for τδ = 20 ms has negative regions after passing its BT bifurcation and adopting supercritical Hopf onset. Finally when Gin is above the cusp value (Fig 5D), all the neurons have Hopf onset with negative regions in their PRCs. However, the neuron with τδ > 20 ms has a far more sinusoidal PRC with a greater negative region due to its onset being via a supercritical rather than subcritical Hopf bifurcation.
For example at τδ = 0, nS and for all τδ the cusp bifurcation lies at
nS. When
(
ms), increasing Gin switches the onset PRC shape first from a symmetric SNIC PRC (A) to an asymmetric HOM PRC (B-C), and later to a Hopf-like PRC (D). Increasing τδ increases the value of Gin at which the SNIC → HOM transition occurs and also decreases the Gin value of the HOM → Hopf transition. For
, the PRC transitions straight from SNIC to Hopf-like. DS parameters used ℓ = 5.
Thus the procession of PRCs of the DS model as Gin is similar to the S model when . For
, the asymmetric positive-valued PRC associated with HOM onset is not found, as predicted from the disappearance of the SNL bifurcation in the previous section. The Gin bifurcation values found earlier are thus useful predictors of the neuron’s spiking response in the dendrite-and-soma model.
Full morphology test
To demonstrate that switches in the dynamical spiking type from increased input conductance are applicable to more complex and realistic neuronal morphologies, we calculated the PRCs from simulations of reconstructed dendritic arbours. In this case, we used the reconstructed dendritic arbour of a real Purkinje cell from NeuroMorpho.Org [52, 53]. We kept the somatic properties the same as the DS model investigated earlier and set all the dendritic compartments to have passive dynamics with τδ = 2.5 ms.
Starting with only the soma (representing the S model with default class I parameters), we increased Gin by adding dendritic compartments to reconstruct the full dendritic arbour in stages, as shown in the top row of Fig 6. At each stage of arborisation, we found the somatic onset current and measured the PRC. If we choose the Gin of the full dendritic arbour to be what would be in the subcritical Hopf regime of the DS model, then we can find different PRCs representing their respective dynamical spiking types by tuning the extent of dendritic arborisation. For instance, at the small arborisation in Fig 6A, Gin is low and we have a symmetric PRC indicative of SNIC onset. Increasing Gin first gives rise to an asymmetric HOM-like PRC (Fig 6B and 6C) before eventually yielding a PRC with substantial negative regions representing Hopf onset (Fig 6D).
Here we increase the input conductance of the multicompartment model by “growing” the dendritic arbour. We can see this switches the PRC shape from SNIC (A) to HOM (B-C) to Hopf (D) as in the DS model earlier. PRCs obtained from the equivalent DS model closely agree with the full morphology at each of the morphological stages.
We then compared the detailed-morphology PRCs with PRCs obtained from an equivalent DS model. For this DS model, the dendritic length L and length constant λ were extracted for each stage of the Purkinje cell by reducing the arbour to an equivalent cable [28, 54]. Due to the different number of compartments between the detailed arbours and their DS equivalent models, we normalised the PRCs when comparing them by setting each PRC peak value to one. Fig 6 shows that the relative PRCs of the equivalent DS model in each stage agree very closely with the PRCs obtained from the full morphology, with small deviations of PRC peak location in the HOM region being the only noticeable difference.
While we have not attempted to accurately capture the complex active dynamics of Purkinje cells (whether in the soma or in the dendrites), we have used its highly branched and complex dendritic arbour to demonstrate that our modelling of a single equivalent dendrite and soma is adequate to approximately capture the affect on neuronal dynamical spiking type of realistic dendritic arbours. Our analysis is thus not limited to the equivalent dendrite morphology in the DS model, but is applicable to realistic morphologies with highly bifurcated arbours with tapering dendrites.
Network simulations
We now illustrate how the dendritic conductance load can affect the network synchronisation states of a small population of these neurons via a switch in dynamical spiking type. To show most clearly how the PRC affects the network state, we first examine pairs of coupled neurons. For a pair of weakly coupled, spiking neurons, a phase-reduction of the model leads to the phase difference ψ evolving in time as [9, 19, 55, 56]
(8)
where Δω is the difference in the uncoupled neuronal spiking frequencies between the two neurons and H(ψ) is the coupling function. The phase locked states ψ* occur when dψ/dt = 0 and for neurons with identical uncoupled spiking frequencies (Δω = 0), this means H(ψ) = 0. Phase-locked states are stable when the gradient of the coupling function at ψ* is negative H′(ψ*) < 0. Furthermore, when the synapses are instantaneous and current-based, H(ψ) is equal to twice the odd part of the PRC
(9)
We simulate one pair of neurons with low dendritic load (Gin = 3.45 nS), placing them in the region of SNIC onset, and another pair of neurons with higher dendritic load in the region of HOM onset (Gin = 4.53 nS). Each pair of coupled neurons was connected to the other with excitatory, instantaneous current-based synapses. The synapses were connected to the soma of each neuron with zero transmission delay. The synaptic strengths were chosen such that the maximum phase advance from a single synaptic input as predicted by the PRC is ∼0.1, as shown in Fig 7C. The firing rate of each uncoupled neuron in both networks was set at 1 Hz. Thus the two pairs differ in their dynamical spiking type and not the spiking frequency or the effective synaptic strength.
Both pairs had the same initial phase relation of ψ0 = 0.3 with the final network states shown in (A-B). (A) The SNIC pair almost achieves a synchronous network state, though this is weakly stable due to the near-symmetry of the PRC. (B) The HOM network robustly achieves an anti-phase state in which the neurons are phase-locked to fire half a cycle apart. (C) The PRCs of each neuron in the SNIC and HOM networks. (D) Coupling functions of the SNIC and HOM pairs. Phase-locked states in a pair of neurons exist where the coupling function is zero. In panels (A-B), time is measured in terms of the number of interspike-intervals (ISIs) of the network spiking rate.
Starting the initial phase difference at ψ0 = 0.3, Fig 7A shows that the SNIC pair converges to a phase locked state close to in-phase synchrony, while the HOM pair in Fig 7B converges to clear anti-phase synchronisation. The coupling functions of the two pairs in Fig 7D reveal why these different phase-locked states are achieved: while H(ψ = 1/2) = 0 for both pairs with H′(ψ) < 0, the zero-crossing at HHOM(1/2) is much further from its other crossings than the zero-crossings of HSNIC are from 1/2. This gives the HOM case a larger basin of attraction. Furthermore, HSNIC has lower amplitude across the phase range than HHOM, meaning that any phase-locked states from the SNIC pair will be achieved more slowly.
We next expand this analysis by looking at networks each with 5 neurons with all-to-all connectivity and maintaining the excitatory, instantaneous current-based synapses as in the neuronal pairs. As before, all neurons in a given network have identical properties. For a group of neurons where each member has a similar spiking frequency, we measure the synchronicity of the network by the standard deviation of the phase differences ψi between the N neurons scaled by
(10)
This expression can be simplified by noting that ∑i ψi = 1 and 〈ψ〉 = 1/N:
(11)
For N neurons all synchronised in-phase, R = 1, while in the splayed-out state with ψi = 1/N for all i, R = 0.
We see in Fig 8A that the phase relations between neurons approaches a synchronous state but Fig 8C shows that this network does not settle to stable in-phase synchrony. As in the pairwise coupling case, the lack of a strongly attractive stable network spiking state for this network arises from the near-perfect symmetry of the SNIC PRC [9].
Both networks had the same initial phase relations with the final network states shown in (A-B). (A) The SNIC network almost achieves a synchronous network state, though this is weakly stable due to the near-symmetry of the PRC. (B) The HOM network robustly achieves a splay state in which the neurons are phase-locked in which neuron i has a constant phase difference of 1/5 with neuron i + 1. (C) The synchrony measure over time shows that the SNIC network gradually and non-monotonically approaches the synchronous state while the HOM network converges to the splay state far more rapidly. In panels (A-C), time is measured in terms of the number of interspike-intervals (ISIs) of the network spiking rate.
Fig 8B shows that the HOM network converges to a splayed-out network spiking state which is approached rapidly as shown in Fig 8C. This splay state is made stable by the fact that the HOM PRC is asymmetric with a peak at a phase less than 1/N = 0.2 [57].
Thus we have demonstrated how changes to the neuronal dynamics conferred by a more extensive neuronal morphology can affect the network behaviour. While these network states will be altered further by synaptic dynamics, transmission delays, and heterogeneities in neuronal properties, in the weak-coupling regime the neuron’s PRC will remain an essential component in determining the stable network states. Inclusion of synaptic dynamics, for example, involves integrating the time course of the synaptic conductance with the PRC [19, 58]. Heterogeneity in neuronal properties complicates the network state as the PRC and uncoupled firing rate of each neuron can then differ. Whether the splayed state (or any other synchronous state) is preserved in this case, can in principle be assessed by considering the difference between all the firing rates [7, 9], with greater differences typically leading to a reduction in phase-locked states.
Finally expanding the size of the neuronal network to higher N (while tuning the synaptic strength such that we remain in the weakly coupled regime) would mean that the full splay state with phase differences 1/N is no longer stable for the HOM network [57], but subsets of the network in phase-locked neurons in splay or anti-phase synchrony could exist. In addition, the anti-phase state shown in Fig 7B, in which some neurons are locked 1/2 out-of-phase from the rest, will remain stable irrespective of the size of the network if the gradient of the PRC at θ = 0 and θ = 1/2 is negative [56, 59]. The basin of attraction for the anti-phase state, along with all other phase-locked states, will tend to decrease as N increases however, making it less robust to perturbations.
Discussion
In this article we have shown that not only can passive dendrites be included in the analysis of conductance-based neuron models, but also that the addition of a dendrite switches the dynamical spiking type of the neuron. In particular, we find that the reduction in input impedance caused by the dendrite can induce the SNIC → HOM dynamical switch, changing the onset PRC from symmetric to asymmetric. This dynamical spiking switch not only affects how information is processed within the individual neurons, it also implies that alterations in dendritic arborisation allow the network dynamics to achieve stable (de)synchronisation states. Specifically, our network simulations show that neurons with greater dendritic load (i.e. more extensive dendritic arbours) can achieve splay states for excitatory coupling due to their HOM onset. In contrast, these splay states could not be achieved by neurons with lower dendritic load (i.e. smaller dendritic arbours) which had SNIC onset. Different network (de)synchronisation states can thus be achieved by tuning the passive dendritic properties common to all neurons. Moreover, our findings establish a direct relationship between the susceptibility of neural tissue to synchronous network states and the morphology of the cells involved, helping us to better understand principles of neural design as well as the effect of deviations thereof in neuropathologies.
Dendritic modelling studies of conductance-based models have been performed previously. These include investigation of the effect of dendritic load on the firing frequency [60] and burst firing [15, 61], the effect of dendritic perturbations on spike timing [62, 63], and the effect of dendritic coupling between spike-generating zones in the same neuron [64]. However, the work presented here shows via bifurcation analysis how morphology alters the dynamical spiking type and ultimately the (de)synchronisation state of the network.
The differential effect on the local bifurcations and PRCs when comparing the single-compartment and dendrite-and-soma models has utility for reducing the number of compartments in neuronal models for larger network simulations. If the input conductance is in the SNIC onset regime and far from the BT point, then a single-compartment model with leak conductance equal to the input conductance of the morphology is appropriate. However, if Gin is close to , then one must use more sophisticated approximations of the input impedance (for example [65–67]). Our method showing how to calculate
in spatially extended models thus informs one when these more sophisticated approximations are required.
Our results also prompt the hypothesis that neuronal morphology may contribute to the differential susceptibility of brain areas to pathological network states like epileptiform activity or spreading depolarisation [5, 68]. This could occur in two related ways. Firstly, pathological spiking patterns have been produced in modelling studies of single cells by changing the spiking dynamics [68, 69]. Because we have shown that the dendritic impedance alters spiking dynamics, changing the dendritic arborisation could therefore move the dynamical state of the neuron to a pathological region (via a switch in its neuronal dynamical spiking type) where network synchrony is changed. Second, pathological behaviour is often associated with either surplus or insufficient synchronous activity [3–5]. Since we have shown that the dynamical spiking switches resulting from added dendritic leak affect (de)synchronisation states, it follows that dendritic arborisation may help move the degree of network synchrony to or from a pathological state.
Presently there is much research interest in how dendrites contribute to neuronal computation. This has largely focused on how either the nonlinear active dendritic channels [70–72] or nonlinear dendritic synapses in conjunction with passive dendritic compartments [73] affect the voltage or firing rate at the soma (i.e. voltage or rate coding). In demonstrating how the passive dendritic contribution changes the dynamical spiking type generated by the spiking compartment, we have added insights on one additional mechanism how dendrites can affect the temporal encoding of neuronal networks.
Methodologically, bifurcation analysis allowed us to calculate important local bifurcations of the system from a model consisting of an active soma attached to a spatially continuous passive dendrite. This approach is not restricted to the Morris-Lecar model examined in this article, but to any conductance-based neuron model with independent voltage-activated ion channels (e.g. the Wang-Buzsáki model [40]). In addition, the reasoning in [10, 42, 46] for related parameters that change the input conductance and the neuronal timescale implies that the bifurcation structure found is general across the whole class of conductance-based neuron models. The calculation of the BT and BTC bifurcations enabled us to predict how the dendrite affects the dynamical type, as these bifurcations act as organising centres for different dynamical types and different switches between dynamical types respectively.
Specifically, using the external input current, passive input conductance and dendritic time constant as bifurcation parameters, we found that the dendrite differentially affects the local bifurcations of the system. The saddle-node and cusp bifurcations are unaffected by the dendritic time constant of the system, meaning that all dendrite-and-soma models have the same saddle-node bifurcation locations as their corresponding single-compartment model with equivalent leak conductance. In contrast, the BT bifurcation moves in an “anticlockwise” direction about the cusp as the dendritic time constant increases, with the emerging Hopf bifurcation switching from subcritical to supercritcal at the BTC bifurcation. This change in the BT bifurcation with the dendritic time constant demonstrated that the effect of passive dendritic load on the dynamical spiking type cannot be fully replicated with a single-compartment model with an increased leak conductance.
On a mathematical note, the “anticlockwise” shift in the BT bifurcation with increasing τδ meant that the switch between class I and class II excitability occurs at higher input conductances for . Furthermore, it means that the switch between SNIC and HOM onset at the SNL bifurcation occurs at higher input conductances, and that the input conductance region of HOM onset is smaller. Examining the PRCs confirmed the predicted spiking onset types from the bifurcation analysis: the SNIC → HOM switch is shifted to higher Gin when τδ is increased in the Morris-Lecar model. When τδ is above the BTC value, the HOM PRC region was eliminated.
Interestingly, the temporal sensitivity of neurons, as captured by the PRCs obtained from the morphologically detailed Purkinje cell reconstruction and its simplified DS model demonstrated the quantitative validity of our reduction approach in the earlier part of the analysis: PRCs of reconstructed neurons and those from the reduction of the dendritic arbour to a single equivalent cylinder were in excellent agreement. This demonstrates that the dynamical spiking type of a morphologically detailed neuron with passive dendrites can be predicted by knowing just its equivalent dendrite reduction and its active spike-generating currents. Furthermore, the equivalent dendrite reduction implies that an arbour that has more dendrites branching off the soma and thicker dendrites will have a greater impact on the neuron’s dynamical spiking type.
We note that this work also establishes a framework from which other influences of the dendritic arbour on a neuron’s spiking dynamics can be explored. As a first example, it allows us to analyse the effect of inputs applied to the dendrites. For a steady external current applied at an arbitrary location, our approach to calculating the local bifurcation structure can still be applied, as detailed in S1 Text.
Meanwhile for transient dendritic perturbations, the PRCs from dendritic input have been simulated and observed experimentally in previous studies [17, 62, 63], where it is found that the amplitude of the PRC decreases with the distance of the perturbation from the soma. On the other hand, the PRC shape is unchanged for class I excitability PRCs and is shifted for class II excitability PRCs. Our work is thus complementary to studies on PRCs measured from dendritic perturbations, as it explores the question of how morphological size and membrane properties affect the PRC. In principle, our approach can be combined with the research in [17, 62, 63] which focuses on how synaptic position affects the PRC.
Furthermore, background synaptic activity targeting the dendritic arbour also affects integrative properties of the neuron [32–34, 74–76]. This can be accommodated in our analysis via changes to the dendritic leak conductance and time constant.
The effect of some dendritic active currents can also, in principle, be included in our modelling approach. If the currents are linearisable under the quasi-active approximation [77], then the effect of the active currents can included in the filtering properties of the dendrite [64, 78–80]. Weakly active dendritic currents, such as Ih and small-conductance Ca2+-activated potassium current ISK, have been found to change the PRCs from being positive-valued everywhere to having negative regions [63, 81], and thus may change the network synchronisation state. If the active currents are strong enough to elicit dendritic spikes [37], then this falls outside of the framework described in our work, firstly because there are now multiple spiking compartments, and second because the interaction between strong dendritic nonlinearities (e.g. spikes, plateau potentials) and axosomatic action potentials can induce bursting behaviour outside of the regular spiking regime [82–84]. However, some aspects of our work may still be used if one considers the strongly active dendritic region to be another oscillator coupled dendritically to the axosomatic compartment [64, 85]. In S1 Text, we extract simulated PRCs from an active dendrite-and-soma model and show that changes in morphology affect the PRC apart the case where the dendrite and soma have identically valued active properties (Fig A in S1 Text). This is because increasing the relative size of the dendrite in comparison to the soma shifts the overall dynamics of the system to be more like that of the dendrite in isolation. Furthermore, it is physiologically unexpected that both the channel types and densities would be identical throughout the neuron, with, for example, Na+ channels being typically found in much higher densities in the AIS than the dendrites [86].
Lastly, there has been much interest in how the geometry of the axonal initial segment (AIS) affects spike threshold and spike shape [87–89] and encoding of time-varying input [16, 90]. Extending our framework to include features of the axon, such as the spatial separation between the soma and the spike-initiation point on the AIS, would allow exploration of how the geometry of the system affects the dynamical spiking type. Changes in the dynamical type would in turn give further insight into how axonal geometry affects the neuron’s temporal encoding of input.
In summary, we demonstrate that neuronal morphology can significantly influence the state of neuronal networks. Moreover, our results describe a method by which one can directly assess the impact of a dendritic arbour on neuronal excitability. Our approach is flexible in allowing for any conductance-based neuron model and external input currents applied to any location. Thus, the proposed approach enables further exploration of how dendritic arbour impacts the tuning of temporal encoding via changes to the network (de)synchronisation state. Our work highlights neuronal morphology as a contributor to neural function via changes to network synchrony, making it therefore likely that morphology is relevant for the differential sensitivity of neuronal tissue to synchronisation in health and disease.
Model and methods
Active soma and passive dendrite model
A single-compartment (S) conductance-based neuron model consists of a leak current, B voltage-activated currents Ia,j and an externally applied current Iext (12)
where each voltage-activated current depends on a set of activation/inactivation gating variables (hereafter termed “active variables”),
. Altogether this gives K active variables indexed ai from i = 1, …, K. As in prior research regarding conductance-based models [6, 13, 51], we assume that the active variables are independent of each other and that their steady state distributions ai,∞(v) saturate as v → ±∞. Each active current has a reversal potential Ej, a maximal conductance Gj, and depends on the product of its active variables
(13)
where pi is the gating exponent associated with active variable ai. Finally, the each active variable evolves in time with a voltage-dependent time constant τi(v)
(14)
Using Rall’s principle of an equivalent cylinder, a passive dendritic arbour emanating from the soma can be simplified to a single equivalent dendrite [28, 29]. We therefore modify the above model to include the dendritic morphology by attaching a passive equivalent dendrite to a conductance-based active soma. The dendrite is spatially continuous, with the coordinate x denoting the distance away from the soma and vδ(x) denoting the dendritic potential at location x.
The dendrite is parametrised by its electrotonic length constant λ, its passive time constant τδ and the dendritic dominance factor ρ (the ratio between the characteristic dendritic conductance and the somatic leak conductance) [29]. These are each defined in terms of electrophysiological parameters
(15)
where cδ is the dendritic membrane capacitance per unit area, gδ is the leak membrane conductance per unit area, ra is the axial resistivity, d is the dendritic diameter and Gσ is the leak conductance of the soma.
Recalling that we denoted dv/dt in Eq 12 as fS, this means that the somatic voltage vσ evolves as
(16)
where the last term represents the axial current flowing from the dendrite to the soma.
The dendritic potential obeys the passive cable equation
(17)
where for simplicity we have assumed that the leak reversal potential is the same in the dendrite as the soma. The dendritic potential is subject to continuity of potential at x = 0 so vδ(x = 0) = vσ, and a sealed-end boundary condition at x = L
(18)
Defining the electrotonically normalised length as ℓ = L/λ, when ℓ ≫ 1, the distal dendritic end is too far away to be influenced by somatic activity, and we can simplify the model by making the dendrite semi-infinite in extent.
Calculation of local bifurcations
Here we will describe how the local bifurcations of the DS system can be calculated for any conductance-based soma model with (Iext, Gin, τδ) as the bifurcation parameters. Although we primarily focus on the semi-infinite dendrite in this article, the method we outline here can be adapted for the finite dendrite, as given in S1 Text. The equations we derive for each bifurcation are applicable to a spatially continuous dendrite rather than for a specific finite number of dendritic compartments.
Fixed points.
Since local bifurcations are defined by the properties of a fixed point, we first outline how to calculate fixed points of the DS system. The fixed points are defined as the values (a*, v*) for which all time derivatives of the system (Eqs 14, 16 and 17) are equal to zero. From Eq 14, each of the active variables satisfies at equilibrium. The cable equation of Eq 17 at equilibrium becomes the second-order ODE
(19)
This has the following solution in terms of vσ (20)
which in the semi-infinite limit reduces to
(21)
Thus all the active variables and the dendritic potential for all x at equilibrium are given in terms of the somatic resting potential . Defining the steady-state current as
(22)
we can substitute in the dendritic voltage derivative at x = 0 to give in the semi-infinite case
(23)
The somatic equilibrium potential is obtained by numerically solving I∞ = 0. We note that the steady-state current equation (Eq 23) is the same as what we would find for the S model if GL = Gσ(1 + ρ) as τδ is not present.
Saddle-node (SN) bifurcation.
At a saddle-node bifurcation, a saddle and node fixed-point meet. Therefore we have a repeated root of I∞ = 0, which means at the saddle-node
(24)
For codimension one bifurcations such as the saddle-node bifurcation, we choose Iext as the bifurcation parameter as we are interested in the onset current for spiking. Hence we can numerically solve Eq 24 to obtain as it does not contain Iext.
is obtained by rearranging I∞ = 0. Note that for the same Gin, Eq 24 is also identical for both S and DS models.
Cusp bifurcation.
At the cusp, two saddle-node bifurcations meet. Thus, the cusp bifurcation satisfies the condition
(25)
For codimension-two bifurcations such as the cusp, we will use (Iext, Gin) as our bifurcation parameters. This is analogous to the prior use of (Iext, GL) in point-neuron models [13]. Since does not depend on Gin or Iext, we numerically solve Eq 25 to obtain
.
and
are obtained from rearranging Eqs 23 and 24 respectively. As with the saddle-node bifurcation, the cusp bifurcation condition does not depend on τδ, thus the cusp bifurcation values
will be the same for the S and DS models.
Hopf bifurcation.
At a Hopf bifurcation, a fixed-point changes stability and a limit cycle appears. The criticality of the Hopf bifurcation determines the stability of the limit cycle involved; at a subcritical Hopf bifurcation we have the transition: stable FP + unstable LC → unstable FP, whilst at a supercritical Hopf bifurcation we have: stable FP → unstable FP + stable LC.
The Jacobian evaluated at a Hopf bifurcation has two purely imaginary conjugate eigenvalues ±iωH. Denoting the corresponding right-eigenvector as q and its complex conjugate as , this means that we have
(26)
From this pair of equations, we can derive two simultaneous nonlinear equations in terms of for the discretised system. Then we take the continuum limit to obtain two nonlinear equations for the continuous dendrite for which further details are found in S1 Text, Eqs BH and BI. These equations can be obtained to find
and then
is obtained from from I∞ = 0.
Unlike all the previous bifurcations listed, the Hopf bifurcation depends on τδ. Therefore, we should expect the bifurcation values of the Hopf bifurcation to differ between the S and DS models. The criticality of the Hopf bifurcation can be calculated using the approach outlined in [41], as detailed in S1 Text, Eq (BO).
Bogdanov-takens (BT) bifurcation.
At the BT point, a saddle-node, Hopf and saddle-node homoclinic orbit (HOM) bifurcations meet. Since the global HOM bifurcation can also determine the spiking onset type of the neuron [10, 13, 25, 45], the BT point acts as an organising centre in two parameter dimensions for different spiking onset types. Varying the bifurcation parameters (Iext, Gin) in the vicinity of the BT point can therefore induce spiking onsets associated with the three codimension one bifurcations that meet here.
The Jacobian of the system has two zero eigenvalues at the BT point. This means that J will have orthogonal left, l, and right, r, eigenvectors satisfying
(27)
These equations can be reduced to a single equation
(28)
where for the semi-infinite dendrite α0 = 1/2. Since Eq 28 has no dependence on our bifurcation parameters (Iext, Gin), we can numerically solve it to find the equilibrium potential
.
can then be obtained by rearranging the saddle-node condition (Eq 24) and
from the equilibrium condition I∞ = 0. We again see that the BT bifurcation depends on τδ and thus we should expect the location of the BT bifurcation to change with τδ.
Bogdanov-takens-cusp (BTC) bifurcation.
Finally, we outline how to find the BTC point. The BTC point is of particular interest because bifurcations associated with spike onset transitions coincide at this point. At a spiking onset transition, the spiking changes from one type to another. We should distinguish this from the lower codimension BT point by stating that the BT point organises bifurcations responsible for spiking onsets, while the BTC point organises bifurcations responsible for transitions of spiking onsets. The bifurcations that meet at the BTC point include the Bautin (where a Hopf bifurcation changes stability), neutral saddle-node (where a homoclinic bifurcation and a fold of limit cycles meet), saddle-node loop (where saddle-node and homoclinic bifurcations meet), and BT bifurcations [13, 25, 41].
Since the BTC point acts as an organising centre for spiking onset transitions, deviations in the bifurcation parameters around it will lead to variations in the spiking onset structure in two parameter dimensions. The BTC point also has the advantage of being easier to calculate than the global bifurcations which coalesce at it, as it is calculable from the properties of fixed-points.
The BTC point has codimension 3, meaning that its location must be expressed in terms of three bifurcation parameters. In the work of Kirst et al [13], an equation for the BTC point for point-neuron conductance-based models is given with (Iext, GL, Cm) as the bifurcation parameters, while in [6], Al-Darabsah and Campbell use the M-current conductance instead of Cm. Here we use the bifurcation parameters (Iext, Gin, τδ) as they are common to all conductance-based neuron models and it includes two dendritic parameters. Furthermore, τδ behaves similarly to Cm, in that increasing it slows the response of the dendrite to somatic activity.
As the BTC point is where a cusp and BT bifurcation coincide, one can use the cusp condition (Eq 25) to obtain the equilibrium voltage , the saddle-node condition (Eq 24) to give
, and the BT condition (Eq 28) to yield
.
Simulation details
To simulate the DS models, a cable of length L = 1000 μm was discretised into M = 50 evenly spaced spatial compartments with step size Δx = 20 μm. The soma occupied a separate compartment at x = 0. All simulations utilised the DifferentialEquations.jl package [91] in the Julia programming language [92], from which the Tsitouras 5/4 Runge-Kutta method was used.
For the SNL bifurcation and PRCs, we required an estimate of the onset current at which stable spiking begins. Here we used the value of Iext which produced regular spiking at the lowest possible frequency above 1 Hz. This means that for class I onset the spiking frequency was approximately 1 Hz, while for class II onset this was typically greater.
At this onset current, PRCs were obtained by increasing the somatic potential by an amount Δv at a single time tk corresponding to a phase θk. To ensure that the phase shift is approximately in the linear regime, Δv was chosen for each neuron such that the PRC maximum was never greater than 0.1. This procedure was repeated for input phases corresponding to θk = 0.01k − 0.005, k = 1, …, 100.
Code used to run all simulations, along with that used to calculate the local and global bifurcations, is available in the “morph-excite-code” GitHub repository we created [93].
Supporting information
S1 Text. Document containing the equations and parameters of the Morris-Lecar model (Table A); an alternative, relative conductance formulation of the model; active dendrite simulations (Fig A); mathematical derivations of how to find the local bifurcations; procedures for finding the global bifurcations; and the approach and parameters used to calculate the PRCs for reconstructed neuronal morphologies (Table B).
https://doi.org/10.1371/journal.pcbi.1011874.s001
(PDF)
Acknowledgments
We would like to thank Jan-Hendrik Schleimer and Paul Pfeiffer for insightful discussions about excitability switching.
References
- 1. Marder E, Bucher D. Central pattern generators and the control of rhythmic movements. Current biology. 2001;11(23):R986–96. pmid:11728329
- 2. Hürkey S, Niemeyer N, Schleimer JH, Ryglewski S, Schreiber S, Duch C. Gap junctions desynchronize a neural circuit to stabilize insect flight. Nature. 2023:1–8. pmid:37225999
- 3. Bergman H, Deuschl G. Pathophysiology of Parkinson’s disease: From clinical neurology to basic neuroscience and back. Movement Disorders. 2002 Mar;17(S3):S28–40. pmid:11948753
- 4. Uhlhaas PJ, Singer W. Abnormal neural oscillations and synchrony in schizophrenia. Nature Reviews Neuroscience. 2010;11(2):100–13. pmid:20087360
- 5. Jiruska P, De Curtis M, Jefferys JG, Schevon CA, Schiff SJ, Schindler K. Synchronization and desynchronization in epilepsy: controversies and hypotheses. The Journal of physiology. 2013;591(4):787–97. pmid:23184516
- 6. Al-Darabsah I, Campbell SA. M-current induced Bogdanov–Takens bifurcation and switching of neuron excitability class. The Journal of Mathematical Neuroscience. 2021;11(1):1–26. pmid:33587210
- 7. Hesse J, Schleimer JH, Maier N, Schmitz D, Schreiber S. Temperature elevations can induce switches to homoclinic action potentials that alter neural encoding and synchronization. Nature Communications. 2022;13(1):1–12. pmid:35803913
- 8. Hodgkin AL. The local electric changes associated with repetitive action in a non-medullated axon. The Journal of physiology. 1948;107(2):165. pmid:16991796
- 9.
Izhikevich EM. Synchronization. In: Dynamical systems in neuroscience. MIT press; 2007. p. 443–505.
- 10. Hesse J, Schleimer JH, Schreiber S. Qualitative changes in phase-response curve and synchronization at the saddle-node-loop bifurcation. Physical Review E. 2017;95(5):052203. pmid:28618541
- 11. Stiefel KM, Gutkin BS, Sejnowski TJ. Cholinergic neuromodulation changes phase response curve shape and type in cortical pyramidal neurons. PloS one. 2008;3(12):e3947. pmid:19079601
- 12. Prescott SA, De Koninck Y, Sejnowski TJ. Biophysical basis for three distinct dynamical mechanisms of action potential initiation. PLoS computational biology. 2008;4(10):e1000198. pmid:18846205
- 13. Kirst C, Ammer J, Felmy F, Herz A, Stemmler M. GABA regulates resonance and spike rate encoding via a universal mechanism that underlies the modulation of action potential generation. bioRxiv. 2017:206581.
- 14. Contreras SA, Schleimer JH, Gulledge AT, Schreiber S. Activity-mediated accumulation of potassium induces a switch in firing pattern and neuronal excitability type. PLoS Computational Biology. 2021;17(5):e1008510. pmid:34043638
- 15. Mainen ZF, Sejnowski TJ. Influence of dendritic structure on firing pattern in model neocortical neurons. Nature. 1996;382(6589):363–6. pmid:8684467
- 16. Eyal G, Mansvelder HD, de Kock CP, Segev I. Dendrites impact the encoding capabilities of the axon. Journal of Neuroscience. 2014;34(24):8063–71. pmid:24920612
- 17. Tiroshi L, Goldberg JA. Population dynamics and entrainment of basal ganglia pacemakers are shaped by their dendritic arbors. PLOS Computational Biology. 2019;15(2):e1006782. pmid:30730886
- 18. Beniaguev D, Shapira S, Segev I, London M. Multiple Synaptic Contacts combined with Dendritic Filtering enhance Spatio-Temporal Pattern Recognition capabilities of Single Neurons. bioRxiv. 2022.
- 19. Ermentrout B. Type I membranes, phase resetting curves, and synchrony. Neural computation. 1996;8(5):979–1001. pmid:8697231
- 20. Teramae Jn, Tanaka D. Robustness of the noise-induced phase synchronization in a general class of limit cycle oscillators. Physical review letters. 2004;93(20):204103. pmid:15600929
- 21. Ko TW, Ermentrout GB. Phase-response curves of coupled oscillators. Physical Review E. 2009;79(1):016211. pmid:19257126
- 22. Larkman A, Mason A. Correlations between morphology and electrophysiology of pyramidal neurons in slices of rat visual cortex. I. Establishment of cell classes. Journal of Neuroscience. 1990;10(5):1407–14. pmid:2332787
- 23. Mohan H, Verhoog MB, Doreswamy KK, Eyal G, Aardse R, Lodder BN, et al. Dendritic and axonal architecture of individual pyramidal neurons across layers of adult human neocortex. Cerebral Cortex. 2015;25(12):4839–53. pmid:26318661
- 24. Van Aerde KI, Feldmeyer D. Morphological and physiological characterization of pyramidal neuron subtypes in rat medial prefrontal cortex. Cerebral cortex. 2015;25(3):788–805. pmid:24108807
- 25. Pereira U, Coullet P, Tirapegui E. The Bogdanov–Takens normal form: a minimal model for single neuron dynamics. Entropy. 2015;17(12):7859–74.
- 26. Rall W. Membrane potential transients and membrane time constant of motoneurons. Experimental neurology. 1960;2(5):503–32. pmid:13739270
- 27.
Tuckwell HC. Introduction to theoretical neurobiology: linear cable theory and dendritic structure. vol. 1. Cambridge University Press; 1988.
- 28. Rall W, Rinzel J. Branch input resistance and steady attenuation for input to one branch of a dendritic neuron model. Biophysical journal. 1973;13(7):648–88. pmid:4715583
- 29.
Rall W. Core conductor theory and cable properties of neurons. In: Comprehensive Physiology. Wiley Online Library; 2011. p. 39–97.
- 30. Eyal G, Verhoog MB, Testa-Silva G, Deitcher Y, Lodder JC, Benavides-Piccione R, et al. Unique membrane properties and enhanced signal processing in human neocortical neurons. Elife. 2016;5:e16553. pmid:27710767
- 31. Beaulieu-Laroche L, Toloza EH, Van der Goes MS, Lafourcade M, Barnagian D, Williams ZM, et al. Enhanced dendritic compartmentalization in human cortical neurons. Cell. 2018;175(3):643–51. pmid:30340039
- 32. Bindman L, Meyer T, Prince C. Comparison of the electrical properties of neocortical neurones in slices in vitro and in the anaesthetized rat. Experimental brain research. 1988;69(3):489–96. pmid:3371433
- 33. Paré D, Shink E, Gaudreau H, Destexhe A, Lang EJ. Impact of spontaneous synaptic activity on the resting properties of cat neocortical pyramidal neurons in vivo. Journal of neurophysiology. 1998;79(3):1450–60. pmid:9497424
- 34. Destexhe A, Paré D. Impact of network activity on the integrative properties of neocortical pyramidal neurons in vivo. Journal of neurophysiology. 1999;81(4):1531–47. pmid:10200189
- 35. Morris C, Lecar H. Voltage oscillations in the barnacle giant muscle fiber. Biophysical journal. 1981;35(1):193–213. pmid:7260316
- 36. Major G, Larkum ME, Schiller J. Active properties of neocortical pyramidal neuron dendrites. Annual review of neuroscience. 2013;36:1–24. pmid:23841837
- 37. Larkum ME, Wu J, Duverdin SA, Gidon A. The guide to dendritic spikes of the mammalian cortex in vitro and in vivo. Neuroscience. 2022. pmid:35182699
- 38. Tsumoto K, Kitajima H, Yoshinaga T, Aihara K, Kawakami H. Bifurcations in Morris–Lecar neuron model. Neurocomputing. 2006;69(4-6):293–316.
- 39. Liu C, Liu X, Liu S. Bifurcation analysis of a Morris–Lecar neuron model. Biological cybernetics. 2014;108(1):75–84. pmid:24435761
- 40. Wang XJ, Buzsáki G. Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. Journal of neuroscience. 1996;16(20):6402–13. pmid:8815919
- 41.
Kuznetsov YA. Elements of applied bifurcation theory. vol. 112. Springer Science & Business Media; 2013. p. 293–392.
- 42. De Maesschalck P, Wechselberger M. Neural excitability and singular bifurcations. The Journal of Mathematical Neuroscience (JMN). 2015;5(1):1–32. pmid:26246435
- 43. Schecter S. The saddle-node separatrix-loop bifurcation. SIAM journal on mathematical analysis. 1987;18(4):1142–56.
- 44. Guckenheimer J, Labouriau J. Bifurcation of the Hodgkin and Huxley equations: a new twist. Bulletin of Mathematical Biology. 1993;55(5):937–52.
- 45.
Izhikevich EM. Bifurcations. MIT press; 2007.
- 46.
Kirst C. Synchronization, Neuronal Excitability, and Information Flow in Networks of Neuronal Oscillators. Georg-August-Universität Göttingen; 2011.
- 47. Phoka E, Cuntz H, Roth A, Häusser M. A new approach for determining phase response curves reveals that Purkinje cells can act as perfect integrators. PLoS computational biology. 2010;6(4):e1000768. pmid:20442875
- 48. Wang S, Musharoff MM, Canavier CC, Gasparini S. Hippocampal CA1 pyramidal neurons exhibit type 1 phase-response curves and type 1 excitability. Journal of neurophysiology. 2013;109(11):2757–66. pmid:23468392
- 49. Hansel D, Mato G, Meunier C. Synchrony in excitatory neural networks. Neural computation. 1995;7(2):307–37. pmid:8974733
- 50. Brown E, Moehlis J, Holmes P. On the phase reduction and response dynamics of neural oscillator populations. Neural computation. 2004;16(4):673–715. pmid:15025826
- 51. Schleimer JH, Schreiber S. Phase-response curves of ion channel gating kinetics. Mathematical Methods in the Applied Sciences. 2018;41(18):8844–58.
- 52. Vetter P, Roth A, Häusser M. Propagation of action potentials in dendrites depends on dendritic morphology. Journal of neurophysiology. 2001;85(2):926–37. pmid:11160523
- 53. Ascoli GA, Donohue DE, Halavi M. NeuroMorpho. Org: a central resource for neuronal morphologies. Journal of Neuroscience. 2007;27(35):9247–51. pmid:17728438
- 54. Amsalem O, Eyal G, Rogozinski N, Gevaert M, Kumbhar P, Schürmann F, et al. An efficient analytical reduction of detailed nonlinear neuron models. Nature communications. 2020;11(1):1–13. pmid:31941884
- 55.
Kuramoto Y. Chemical oscillations, waves and turbulence. Springer; 1984.
- 56.
Ermentrout GB, Terman DH. Mathematical Foundations of Neuroscience. vol. 35. Springer Science & Business Media; 2010.
- 57. Ashwin P, Bick C, Burylko O. Identical phase oscillator networks: Bifurcations, symmetry and reversibility for generalized coupling. Frontiers in Applied Mathematics and Statistics. 2016;2:7.
- 58. Van Vreeswijk C, Abbott L, Ermentrout GB. When inhibition not excitation synchronizes neural firing. Journal of computational neuroscience. 1994;1(4):313–21. pmid:8792237
- 59. Ermentrout GB. Stable periodic solutions to discrete and continuum arrays of weakly coupled nonlinear oscillators. SIAM Journal on Applied Mathematics. 1992;52(6):1665–87.
- 60. Schwemmer MA, Lewis TJ. Effects of dendritic load on the firing frequency of oscillating neurons. Physical Review E. 2011;83(3):031906. pmid:21517524
- 61. Psarrou M, Stefanou SS, Papoutsi A, Tzilivaki A, Cutsuridis V, Poirazi P. A simulation study on the effects of dendritic morphology on layer V prefrontal pyramidal cell firing behavior. Frontiers in cellular neuroscience. 2014;8:287. pmid:25278837
- 62. Crook SM, Ermentrout GB, Bower JM. Dendritic and synaptic effects in systems of coupled cortical oscillators. Journal of computational neuroscience. 1998;5(3):315–29. pmid:9663554
- 63. Goldberg JA, Deister CA, Wilson CJ. Response Properties and Synchronization of Rhythmically Firing Dendritic Neurons. Journal of Neurophysiology. 2007 Jan;97(1):208–19. pmid:16956986
- 64. Remme MW, Lengyel M, Gutkin BS. The role of ongoing dendritic oscillations in single-neuron dynamics. PLoS Comput Biol. 2009;5(9):e1000493. pmid:19730677
- 65. Major G, Evans JD, Jack J. Solutions for transients in arbitrarily branching cables: I. Voltage recording with a somatic shunt. Biophysical journal. 1993;65(1):423–49. pmid:8369447
- 66. Wybo WA, Stiefel KM, Torben-Nielsen B. The Green’s function formalism as a bridge between single-and multi-compartmental modeling. Biological cybernetics. 2013;107(6):685–94. pmid:24037222
- 67. Wybo WA, Jordan J, Ellenberger B, Mengual UM, Nevian T, Senn W. Data-driven reduction of dendritic morphologies with preserved dendro-somatic responses. Elife. 2021;10. pmid:33494860
- 68. Wei Y, Ullah G, Schiff SJ. Unification of neuronal spikes, seizures, and spreading depression. Journal of Neuroscience. 2014;34(35):11733–43. pmid:25164668
- 69. Depannemaecker D, Ivanov A, Lillo D, Spek L, Bernard C, Jirsa V. A unified physiological framework of transitions between seizures, sustained ictal activity and depolarization block at the single neuron level. Journal of computational neuroscience. 2022;50(1):33–49. pmid:35031915
- 70. Behabadi BF, Polsky A, Jadi M, Schiller J, Mel BW. Location-dependent excitatory synaptic interactions in pyramidal neuron dendrites. PLoS computational biology. 2012;8(7):e1002599. pmid:22829759
- 71. Tzilivaki A, Kastellakis G, Poirazi P. Challenging the point neuron dogma: FS basket cells as 2-stage nonlinear integrators. Nature communications. 2019;10(1):1–14. pmid:31413258
- 72. Gidon A, Zolnik TA, Fidzinski P, Bolduan F, Papoutsi A, Poirazi P, et al. Dendritic action potentials and computation in human layer 2/3 cortical neurons. Science. 2020;367(6473):83–7. pmid:31896716
- 73. Stöckel A, Eliasmith C. Computational properties of multi-compartment LIF neurons with passive dendrites. Neuromorphic Computing and Engineering. 2022.
- 74. Holmes WR, Woody CD. Effects of uniform and non-uniform synaptic ‘activation-distributions’ on the cable properties of modeled cortical pyramidal neurons. Brain research. 1989;505(1):12–22. pmid:2611664
- 75. Häusser M, Clark BA. Tonic synaptic inhibition modulates neuronal output pattern and spatiotemporal synaptic integration. Neuron. 1997;19(3):665–78. pmid:9331356
- 76. Rudolph M, Destexhe A. A Fast-Conducting, Stochastic Integrative Mode for Neocortical Neurons InVivo. Journal of Neuroscience. 2003;23(6):2466–76. pmid:12657707
- 77. Koch C. Cable theory in neurons with active, linearized membranes. Biological Cybernetics. 1984;50(1):15–33. pmid:6324889
- 78. Coombes S, Timofeeva Y, Svensson CM, Lord GJ, Josić K, Cox SJ, et al. Branching dendrites with resonant membrane: a “sum-over-trips” approach. Biological Cybernetics. 2007;97(2):137–49. pmid:17534649
- 79. Remme MW, Rinzel J. Role of active dendritic conductances in subthreshold input integration. Journal of computational neuroscience. 2011;31(1):13–30. pmid:21127955
- 80. Zhuchkova E, Remme MW, Schreiber S. Somatic versus dendritic resonance: differential filtering of inputs through non-uniform distributions of active conductances. PLoS One. 2013;8(11):e78908. pmid:24223864
- 81. Schultheiss NW, Edgerton JR, Jaeger D. Phase response curve analysis of a full morphological globus pallidus neuron model reveals distinct perisomatic and dendritic modes of synaptic integration. Journal of Neuroscience. 2010;30(7):2767–82. pmid:20164360
- 82. Larkum ME, Zhu JJ, Sakmann B. A new cellular mechanism for coupling inputs arriving at different cortical layers. Nature. 1999;398(6725):338–41. pmid:10192334
- 83. Williams SR, Stuart GJ. Mechanisms and consequences of action potential burst firing in rat neocortical pyramidal neurons. The Journal of physiology. 1999;521(Pt 2):467. pmid:10581316
- 84. Grienberger C, Chen X, Konnerth A. NMDA receptor-dependent multidendrite Ca2+ spikes required for hippocampal burst firing in vivo. Neuron. 2014;81(6):1274–81. pmid:24560703
- 85. Yi G, Wang J, Wei X, Deng B. Action potential initiation in a two-compartment model of pyramidal neuron mediated by dendritic Ca2+ spike. Scientific Reports. 2017;7(1):1–16.
- 86. Kole MH, Ilschner SU, Kampa BM, Williams SR, Ruben PC, Stuart GJ. Action potential generation requires a high sodium channel density in the axon initial segment. Nature neuroscience. 2008;11(2):178–86. pmid:18204443
- 87. Kole MHP, Stuart GJ. Signal Processing in the Axon Initial Segment. Neuron. 2012;73(2):235–47. pmid:22284179
- 88. Hamada MS, Goethals S, Vries SID, Brette R, Kole MHP. Covariation of axon initial segment location and dendritic tree normalizes the somatic action potential. PNAS. 2016;113(51). pmid:27930291
- 89. Goethals S, Brette R. Theoretical relation between axon initial segment geometry and excitability. Elife. 2020;9:e53432. pmid:32223890
- 90. Zhang C, Hofmann D, Neef A, Wolf F. Ultrafast population coding and axo-somatic compartmentalization. PLoS computational biology. 2022;18(1):e1009775. pmid:35041645
- 91. Rackauckas C, Nie Q. Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia. Journal of Open Research Software. 2017;5(1).
- 92. Bezanson J, Edelman A, Karpinski S, Shah VB. Julia: A fresh approach to numerical computing. SIAM review. 2017;59(1):65–98.
- 93.
Gowers R. morph-excite-code GitHub repository; 2024. Available from: https://github.com/rpgowers/morph-excite-code.