How neuronal morphology impacts the synchronisation state of neuronal networks

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.


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][4][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 actionpotential 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][11][12][13][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][16][17][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-ornone 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][20][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][23][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.

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 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 G in .The passive input conductance is given by the constant level of external input current I ext required to change the voltage from its steady state value (when all active conductances are zero) by an amount Δv In the DS model, we apply I ext to the soma and use the voltage change Δv at the soma for G in .In each case G in is the sum of its somatic (G σ ) and dendritic (G δ ) contributions.For a dendrite of effective length When ℓ � 1, we can treat the dendrite as semi-infinite and the passive input conductance of the DS model becomes For the S model G in = G L , the total leak conductance, allowing for a straightforward comparison between the two morphologies.The total capacitance of the S model is denoted by C m .
While it is possible to adjust the dendritic parameters (ρ, λ, L) to match G in for the S and DS models, this is not possible for the passive input admittance Y in (the reciprocal of the input impedance, Y in � Z À 1 in ).Denoting the angular frequency of the input signal as ω, for the S model Y in is simply given as a first-order filter while for the semi-infinite and finite DS models with somatic capacitance C σ we have where γ(ω) is defined as Here setting C m = C σ or choosing any other frequency-independent value of C m will lead to Y F;1 in and Y S in 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, Y in can be fully parametrised by G in and τ δ .From the form of Y in , we see that τ δ is a key parameter for comparing the DS and S models.When τ δ = 0, γ = 1 at all frequencies and hence Y F in and Y 1 in become equal to Y S in in this limit for G in = G L .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 Z in .Fig 2A shows that for any given G in , |Z in | is higher in the S model than the semi-infinite DS model at all non-zero frequencies.In Fig 2B, we see that |Z in | 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 G in constant while decreasing ℓ, then |Z in | decreases at all non-zero frequencies.
In our analyses of the DS model, we use G in and τ δ to compare dendritic arbours rather than the underlying electrophysiological parameters.G in 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 G in without affecting τ δ .Changing τ δ , on the other hand, requires changing the per-area membrane properties of the arbour.While τ δ can be modified without changing G in 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, G in 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 Y in 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 G in and decrease τ δ [32][33][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 Y in 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-Buzsa ´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 (G in , τ δ ) for the semi-infinite dendrite, it naturally follows that we should choose (G in , τ δ ) as bifurcation parameters along with I ext 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 G in above 2 nS denote the conductance load added to the neuron by the dendrite.As in our analysis of the passive impedance properties, here I ext is always applied to the soma; analytical treatment of I ext applied at any dendritic location is described in S1 Text.In short, applying I ext along the dendrite systematically increases the value of I ext for all bifurcations but leaves the other bifurcation parameters, including G in , unchanged.Thus all transitions in the spiking onset type have the same values of G in and τ δ irrespective of the applied current location; only the values of I ext change by a known amount.
An equivalent approach to varying the input conductance G in and applied current I ext can be taken by using the relative dendritic input conductance, ρ = G δ /G σ , and the applied potential, μ ext = I ext /G σ , as bifurcation parameters, as detailed in S1 Text.We show the absolute quantities G in and I ext 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 (I ext , G in ) 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 G in = G L .There are two saddle-node (SN) bifurcations at lower and higher I ext , which we will hereafter refer to as "lower" (dashed) and "higher" (solid line) respectively.These two SN bifurcations converge at the cusp as G in increases.Three fixed points exist for I SN;low ext < I ext < I SN;high ext , 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 (I ext , G in ) in the S and DS models.
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 G in until it reaches the cusp at t BTC d ¼ 12:9 ms (hereafter denoted as the BTC point).Increasing τ δ beyond t BTC d moves the BT point onto the lower SN bifurcation curve (dashed) and G BT in 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 .When the Hopf bifurcation that emerges from the BT point is subcritical, it can be switched to supercritical by increasing G in .This criticality switch moves to lower G in as τ δ increases towards t BTC d .At all τ δ , there is a fold of Hopf bifurcations when G in 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 G in 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 G in in the DS model for t d < t BTC d .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][44][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 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 G SNL in also increases with τ δ when t d < t BTC d .Indeed, this is what we see in Fig 4, with not only G SNL in increasing with τ δ but also the conductance difference between the BT and SNL points decreasing before eventually converging at t BTC d .Thus both the SNIC !HOM switch and the G in range for homoclinic onset are affected by the dendritic time constant.For t d > t BTC d , 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 I FLC ext < I HOM ext , 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 I HOM ext , making it hard to distinguish from class II excitability; and (3) the difference between I FLC ext and I HOM ext is typically extremely small [46] and thus I HOM ext 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 in 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 G in .As a result, one would expect that changing τ δ changes the value of G SNL in .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 G in ) is varied, and that the big SNL bifurcation is reached first when the timescale is shortened (e.g.increasing G in or decreasing τ δ ), followed by the BT bifurcation.Taken together, this means that G SNIC in < G SNL in < G BT in for the two-dimensional conductance-based neuron model, and hence the switching between dynamic spiking types of SNIC !HOM !subcritical Hopf with increasing G in (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 G in .The switching between dynamical spiking types is qualitatively similar to increasing G L in the single-compartment model, however the G in values for the SNIC !HOM switch and subcritical Hopf onset increase with τ δ , whilst the value of G in 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][20][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 (G in , τ δ ).We chose the (G in , τ δ ) values to be in the neighbourhood of the SNL, BT, and cusp bifurcations, as shown by the coloured points in  not found, as predicted from the disappearance of the SNL bifurcation in the previous section.The G in 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  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] 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 (ψ*) < 0. Furthermore, when the synapses are instantaneous and current-based, H(ψ) is equal to twice the odd part of the PRC We simulate one pair of neurons with low dendritic load (G in = 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 (G in = 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.
Starting the initial phase difference at ψ 0 = 0. why these different phase-locked states are achieved: while H(ψ = 1/2) = 0 for both pairs with H 0 (ψ) < 0, the zero-crossing at H HOM (1/2) is much further from its other crossings than the zero-crossings of H SNIC are from 1/2.This gives the HOM case a larger basin of attraction.
Furthermore, H SNIC has lower amplitude across the phase range than H HOM , 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 ffi ffi ffi ffi N p R ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi This expression can be simplified by noting that ∑ i ψ i = 1 and hψi = 1/N: R ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 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].
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 G in is close to G BT in , then one must use more sophisticated approximations of the input impedance (for example [65][66][67]).Our method showing how to calculate G BT in 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][4][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][71][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-Buzsa ´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 t d < t BTC d .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 G in 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][33][34][74][75][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][79][80].Weakly active dendritic currents, such as I h and small-conductance Ca 2+ -activated potassium current I SK , 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][83][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][88][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.

Active soma and passive dendrite model
A single-compartment (S) conductance-based neuron model consists of a leak current, B voltage-activated currents I a,j and an externally applied current where each voltage-activated current depends on a set of activation/inactivation gating variables (hereafter termed "active variables"), a j ¼ ða j1 ; :::; a jb j Þ. Altogether this gives K active variables indexed a i 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 a i,1 (v) saturate as v !±1.Each active current has a reversal potential E j , a maximal conductance G j , and depends on the product of its active variables where p i is the gating exponent associated with active variable a i .Finally, the each active variable evolves in time with a voltage-dependent time constant τ i (v) da i dt ¼ a i;1 ðvÞ À a i t i ðvÞ ¼ f a;i ða 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 where c δ is the dendritic membrane capacitance per unit area, g δ is the leak membrane conductance per unit area, r a 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 f S , this means that the somatic voltage v σ evolves as where the last term represents the axial current flowing from the dendrite to the soma.The dendritic potential obeys the passive cable equation 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 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 (I ext , G in , τ δ ) 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 a * i ¼ a i;1 ðv s Þ at equilibrium.The cable equation of Eq 17 at equilibrium becomes the second-order ODE This has the following solution in terms of which in the semi-infinite limit reduces to Thus all the active variables and the dendritic potential for all x at equilibrium are given in terms of the somatic resting potential v * s .Defining the steady-state current as we can substitute in the dendritic voltage derivative at x = 0 to give in the semi-infinite case The somatic equilibrium potential v * s is obtained by numerically solving I 1 = 0. We note that the steady-state current equation (Eq 23) is the same as what we would find for the S model if G L = G σ (1 + ρ) as τ δ is not present.
Saddle-node (SN) bifurcation.At a saddle-node bifurcation, a saddle and node fixedpoint meet.Therefore we have a repeated root v SN s of I 1 = 0, which means at the saddle-node For codimension one bifurcations such as the saddle-node bifurcation, we choose I ext as the bifurcation parameter as we are interested in the onset current for spiking.Hence we can numerically solve Eq 24 to obtain v SN s as it does not contain I ext .I SN ext is obtained by rearranging I 1 = 0. Note that for the same G in , 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 For codimension-two bifurcations such as the cusp, we will use (I ext , G in ) as our bifurcation parameters.This is analogous to the prior use of (I ext , G L ) in point-neuron models [13].Since I 00 ext ðv s Þ does not depend on G in or I ext , we numerically solve Eq 25 to obtain v C s .G C in and I C ext 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 ðI C ext ; G C in Þ 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 � q, this means that we have From this pair of equations, we can derive two simultaneous nonlinear equations in terms of ðv H s ; o H Þ 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 v H s and then I H ext is obtained from from I 1 = 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 saddlenode 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 (I ext , G in ) 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 lJ ¼ 0; Jr ¼ 0; lr ¼ 0: ð27Þ These equations can be reduced to a single equation where for the semi-infinite dendrite α 0 = 1/2.Since Eq 28 has no dependence on our bifurcation parameters (I ext , G in ), we can numerically solve it to find the equilibrium potential v BT s .G BT in can then be obtained by rearranging the saddle-node condition (Eq 24) and I BT ext from the equilibrium condition I 1 = 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 (I ext , G L , C m ) as the bifurcation parameters, while in [6], Al-Darabsah and Campbell use the M-current conductance instead of C m .Here we use the bifurcation parameters (I ext , G in , τ δ ) as they are common to all conductancebased neuron models and it includes two dendritic parameters.Furthermore, τ δ behaves similarly to C m , 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 v BTC s , the saddle-node condition (Eq 24) to give G BTC in , and the BT condition (Eq 28) to yield t BTC d .

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.jlpackage [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 I ext 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 t k 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.

Fig 1 .
Fig 1.A comparison of the single-compartment (S) model with a dendrite-and-soma (DS) model with equivalent passive input conductance G in .In the DS model, the constant input current I ext is applied to the somatic compartment.In the S model, G in is equivalent to the lumped leak conductance G L , while for the DS model, G in is the sum of the passive somatic (G σ ) and dendritic (G δ ) contributions.When we increase G in in the DS model, G σ is kept constant while G δ is increased.https://doi.org/10.1371/journal.pcbi.1011874.g001

Fig 2 .
Fig 2. (A) When one fixes the input conductance G in between the S model and a semi-infinite DS model, the magnitude of the input impedance Z in will be higher for the S model for non-zero frequencies.τ δ = 10 ms.(B) Increasing τ δ of the DS model decreases |Z in | at all non-zero frequencies.G in = 8 nS.(C) Decreasing the effective dendritic length ℓ while maintaining the same G in also decreases |Z in | at all non-zero frequencies.τ δ = 10 ms.https://doi.org/10.1371/journal.pcbi.1011874.g002

Fig 3 .
Fig 3. (A) Two-parameter local bifurcation diagram in terms of (I ext , G in ) 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 G in .Increasing τ δ shifts the BT bifurcation and shrinks the Hopf bifurcation curve.(B) shows that initially increasing τ δ moves the BT point to higher G in until it reaches the cusp at the BTC point.Increasing τ δ beyond t BTC d moves the BT point to the lower SN bifurcation curve (dashed) and thus decreases G BT in .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 I SN,high , while (D) shows the Hopf bifurcation when we measure the external input current relative to the lower saddle-node value I SN,low .The saddle-node and cusp bifurcations are depicted in black because they are unaffected by changes to τ δ .https://doi.org/10.1371/journal.pcbi.1011874.g003

Fig 4 .
Fig 4. Bifurcation diagram of the Morris-Lecar DS model in terms of (τ δ , G in ) focussing on dynamical switches.Here we have taken I ext as the onset current for every value of (τ δ , G in ).Points indicate values of (τ δ , G in ) 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 G SNL in and G BT in until the BT and SNL bifurcations meet the cusp at the BTC point at t BTC d .The difference between G BT in and G SNL in decreases as τ δ increases, meaning that the range of G in for which HOM onset exists becomes smaller.For t d > t BTC d , the SNL bifurcation no longer exists and G BT in decreases.Furthermore, the Hopf onset for G in > G BT in switches from subcritical to supercritical.G in for subcritical Hopf onset increases with τ δ until t BTC d , while G in for supercritical Hopf onset decreases with τ δ throughout the whole range.https://doi.org/10.1371/journal.pcbi.1011874.g004

Fig 4 .
At G in = 3 nS, the Morris-Lecar neuron has SNIC spiking onset for all τ δ .Hence Fig 5A shows symmetric positive-valued at this input conductance.At G in = 4.7 nS (Fig 5B), the neuron is operating in the HOM regime for lower τ δ and thus we see asymmetric positivelyskewed 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 G in (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 G in 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.Thus the procession of PRCs of the DS model as G in is similar to the S model when t d < t BTC d .For t d > t BTC d , the asymmetric positive-valued PRC associated with HOM onset is

Fig 5 .
Fig 5. Phase-response curves (PRCs) of the Morris-Lecar DS model for various τ δ and G in .Values of G in and τ δ have been chosen to be around the dynamical switches in the Morris-Lecar DS model.For example at τ δ = 0, G BT in ¼ 4:77 nS and for all τ δ the cusp bifurcation lies at G C in ¼ 5:53 nS.When t d < t BTC d (t BTC d ¼ 12:9 ms), increasing G in 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 G in at which the SNIC !HOM transition occurs and also decreases the G in value of the HOM !Hopf transition.For t d > t BTC d , the PRC transitions straight from SNIC to Hopf-like.DS parameters used ℓ = 5. https://doi.org/10.1371/journal.pcbi.1011874.g005 with only the soma (representing the S model with default class I parameters), we increased G in 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 G in 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, G in is low and we have a symmetric PRC indicative of SNIC onset.Increasing G in first gives rise to an asymmetric HOMlike PRC (Fig 6B and 6C) before eventually yielding a PRC with substantial negative regions representing Hopf onset (Fig 6D).

Fig 6 .
Fig 6.The PRCs obtained by simulating detailed multicompartment Purnkinje cell neuron show that the results of the simplified DS morphology are applicable to complicated and realistic dendritic arbours.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.https://doi.org/10.1371/journal.pcbi.1011874.g006

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

Fig 7 .
Fig 7. A comparison of two pairs of neurons with excitatory coupling between, one in which the neurons have SNIC onset and the other in which the neurons have HOM onset.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.https://doi.org/10.1371/journal.pcbi.1011874.g007

Fig 8 .
Fig 8.A comparison of two all-to-all networks with excitatory coupling between 5 neurons, one in which all the neurons have SNIC onset and the other in which all the neurons have HOM onset.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.https://doi.org/10.1371/journal.pcbi.1011874.g008