Bifurcation analysis of a two-compartment hippocampal pyramidal cell model

Journal of Computational Neuroscience, May 2016

The Pinsky-Rinzel model is a non-smooth 2-compartmental CA3 pyramidal cell model that has been used widely within the field of neuroscience. Here we propose a modified (smooth) system that captures the qualitative behaviour of the original model, while allowing the use of available, numerical continuation methods to perform full-system bifurcation and fast-slow analysis. We study the bifurcation structure of the full system as a function of the applied current and the maximal calcium conductance. We identify the bifurcations that shape the transitions between resting, bursting and spiking behaviours, and which lead to the disappearance of bursting when the calcium conductance is reduced. Insights gained from this analysis, are then used to firstly illustrate how the irregular spiking activity found between bursting and stable spiking states, can be influenced by phase differences in the calcium and dendritic voltage, which lead to corresponding changes in the calcium-sensitive potassium current. Furthermore, we use fast-slow analysis to investigate the mechanisms of bursting and show that bursting in the model is dependent on the intermediately slow variable, calcium, while the other slow variable, the activation gate of the afterhyperpolarisation current, does not contribute to setting the intraburst dynamics but participates in setting the interburst interval. Finally, we discuss how some of the described bifurcations affect spiking behaviour, during sharp-wave ripples, in a larger network of Pinsky-Rinzel cells.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs10827-016-0606-8.pdf

Bifurcation analysis of a two-compartment hippocampal pyramidal cell model

J Comput Neurosci Bifurcation analysis of a two-compartment hippocampal pyramidal cell model Laura A. Atherton 0 1 2 Luke Y. Prince 0 1 2 Krasimira Tsaneva-Atanasova 0 1 2 Action Editor: J. Rinzel 0 1 2 0 Department of Mathematics, College of Engineering, Mathematics and Physical Sciences, & EPSRC Centre for Predictive Modelling in Healthcare, University of Exeter , Exeter, Devon, EX4 4QF , UK 1 Physiology, Pharmacology & Neuroscience, University of Bristol , Bristol, England , UK 2 Engineering Mathematics, and Physiology, Pharmacology & Neuroscience, University of Bristol , Bristol, England , UK The Pinsky-Rinzel model is a non-smooth 2compartmental CA3 pyramidal cell model that has been used widely within the field of neuroscience. Here we propose a modified (smooth) system that captures the qualitative behaviour of the original model, while allowing the use of available, numerical continuation methods to perform full-system bifurcation and fast-slow analysis. We study the bifurcation structure of the full system as a function of the applied current and the maximal calcium conductance. We identify the bifurcations that shape the transitions between resting, bursting and spiking behaviours, and which lead to the disappearance of bursting when the calcium conductance is reduced. Insights gained from this analysis, are then used to firstly illustrate how the irregular spiking activity found between bursting and stable spiking states, can be influenced by phase differences in the calcium and dendritic voltage, which lead to corresponding changes in the calcium-sensitive potassium current. Furthermore, we use fast-slow analysis to investigate the mechanisms of bursting and show that bursting in the model is dependent on the intermediately slow variable, calcium, while the other slow variable, the activation gate of the afterhyperpolarisation current, does not contribute to setting the intraburst dynamics but participates in setting the interburst interval. Finally, we discuss how some of the described bifurcations affect spiking behaviour, during sharp-wave ripples, in a larger network of Pinsky-Rinzel cells. Dynamical system; Bifurcation analysis; Bursting and spiking; Numerical continuation; Parameter dependence 1 Introduction An extensive body of literature has implicated the hippocampal formation in spatial navigation and episodic memory (O’Keefe and Nadal 1978; Burgess et al. 2002; Morris et al. 1982; Riedel et al. 1999; Squire 1992) . Within the hippocampus, the recurrent connectivity within the autoassociative CA3 network on the one hand confers onto the region a strong computational capacity which can mediate behavioural processes such as pattern completion (McNaughton and Morris 1987; Rolls and Kesner 2006; Nakazawa et al. 2002) while on the other hand predisposes the region to epileptiform activity (Traub and Wong 1982; Miles and Wong 1983; Le Duigou et al. 2014) . At the cellular level, CA3 pyramidal cells intrinsically display a variety of firing patterns, ranging from single action potentials to complex bursts (Kandel and Spencer 1961; Traub et al. 1991; Mizuseki et al. 2012; Wong and Prince 1978; 1981; Spruston and McBain 2007) . Such bursting is important for place cell activity (Harvey et al. 2009; Epsztein et al. 2011; Bittner et al. 2015) , signal propagation and the induction of synaptic plasticity (Lisman 1997; Buchanan and Mellor 2010) . In an attempt to understand the behaviour of isolated cells and the coupled CA3 network, numerous computational models of CA3 pyramidal cells have been proposed, ranging from detailed, multi-compartmental models down to single compartmental models (Traub et al. 1991; Pinsky and Rinzel 1994; Migliore et al. 1995; Grobler et al. 1998; Xu and Clancy 2008; Nowacki et al. 2011) . The Pinsky-Rinzel model (Pinsky and Rinzel 1994) was originally formulated over 20 years ago as a two-compartment reduction of the 19-compartment Traub CA3 cell model, developed earlier (Traub et al. 1991) . Since then variants of Pinsky-Rinzel model have been used to investigate hippocampal sharp wave ripple oscillations (Taxidis et al. 2012) ; carbacholinduced gamma oscillations (Tiesinga et al. 2001) ; rate and temporal coding of place cells (Kamondi et al. 1998; Booth and Bose 2001) ; and the influence of dendritic morphology on firing patterns (Mainen and Sejnowski 1996), to name but a few examples. Whilst parameter regimes for tonic spiking and bursting activity were investigated in the original Pinsky-Rinzel paper (Pinsky and Rinzel 1994) , a detailed mathematical analysis of the dynamical regimes of the model has only ever been performed once (Hahn and Durand 2001) . This is likely because the Pinsky-Rinzel model is non-smooth, and so traditional techniques used to study dynamical systems cannot be employed. In this previous paper, bifurcation analysis (Strogatz 2001) was used to investigate the transitions between resting, bursting and spiking states as the size of the extracellular potassium concentration was increased (Hahn and Durand 2001) . However, a detailed description of how this analysis was performed, given the non-smooth nature of the system, was never provided. Due to the complexity of the model, a great deal of extra insight can be gained by analysing how some of the many other parameters shape the dynamical landscape of the model. This can then inform parameter choices and potentially explain dynamic behaviour in larger networks of Pinsky-Rinzel cells, such as Taxidis et al. (2012) , which are much more difficult to analyse. Therefore we recast the original model equations using fully continuous functions. This permitted the use of available numerical continuation methods to perform bifurcation analysis using three notable bifurcation parameters: the applied somatic and dendritic currents, and the maximum calcium conductance. These first two parameters were chosen in order to ask what the mathematical mechanisms were for the originally observed transitions between resting, bursting and spiking as the applied current was increased (Pinsky and Rinzel 1994) . This is of particular interest for Pinsky-Rinzel cell behaviour in a larger network, where excitatory and inhibitory inputs may impinge onto either the somatic or dendritic compartment, or both. Would the mechanisms behind the transitions be qualitatively similar from current applied to either compartment? Or would the mechanisms differ? Having understood these transitions, the third bifurcation parameter was chosen because lowering the maximum calcium conductance can switch behaviour from bursting to spiking (Pinsky and Rinzel 1994; Traub et al. 1991) , as has been used to better represent the firing properties of CA1 cells (Taxidis et al. 2012). How does this occur from dynamical systems point of view? Is it just the bursting that disappears, or do other behaviours also change as the calcium conductance is reduced, and if so how might calcium be involved in shaping some of these behaviours? As an example of how answers to the above questions could be used in informing understanding of network activity, we then discuss the implications for some of the bifurcations identified, in relation to a larger network of Pinsky-Rinzel cells, exhibiting sharp-wave ripple oscillations (Taxidis et al. 2012) . Finally, given the relationship between the maximum calcium conductance and bursting, fast slow analysis (Rinzel 1987) was used to further investigate the intraburst dynamics. Traditionally this has been used to isolate computationally the important variables, responsible for bursting, given the difficulty of teasing apart the system experimentally. When this technique was applied to two separate pyramidal cell models based on, but slightly different to, the Pinsky-Rinzel model; bursting behaviour was either of the square-wave/fold-homoclinic type and dependent on the activation variable of the slow potassium current (Kepecs and Wang 2000) or of the parabolic type and dependent on both the slow autocatalytic activation variable for a Ttype calcium channel and the activation variable for the slow calcium-dependent potassium current (Xu and Clancy 2008) . Given these differences, it is difficult to interpolate which mechanism of bursting might exist in the original model, for which fast-slow analysis has never been performed. Moreover, while in the original paper, the authors describe burst initiation as occurring when the slow activation variable, (q), for the potassium afterhyperpolarisation current fell below a threshold value (Pinsky and Rinzel 1994) ; subsequent studies have shown instead, using phaseplane analysis on a piecewise reduced system of the original model, that although q is important for controlling the interburst interval, it is not important for setting the burst initiation threshold (Booth and Bose 2001; Bose and Booth 2005) . Therefore we perform fast-slow analysis, using the modified continuous version of the Pinsky-Rinzel model we propose here, in an attempt to clarify the role of q and also investigate its interaction with the other slow variable, calcium (Ca), in controlling the burst dynamics. 2 Model The Pinsky-Rinzel model characterises a typical pyramidal cell as comprising a single axosomatic and a single dendritic compartment, where the somatic compartment contains a transient sodium INa, delayed rectifier potassium IKDR , and leak current. The dendritic compartment contains a persistent calcium ICa, calcium activated potassium IKCa , after-hyperpolarisation potassium current IKAHP , and leak current. The two compartments are coupled by a coupling current ISD or IDS and their membrane voltages summarised by the following differential equations: dVs Cm dt dVd Cm dt = −ILeak − INa − IKDR − = −ILeak − ICa − IKCa − IKAHP + IDS ISapp p + p ISD IDapp , + (1 − p) + (1 − p) that evolve in time (measured in milliseconds). All currents are conductance-based, using the HodgkinHuxley formalism (Hodgkin and Huxley 1952) of activation and inactivation gates m, h, n, s, c, q dependent on voltage or intracellular calcium that drive the current. Additionally, IKCA has a saturation function dependent on calcium χ (Ca). INa = gNa m2∞(Vs) h (Vs − VNa) IKDR = gKDR n (Vs − VK) ICa = gCa s2 (Vd − VCa) IKCa = gKCa c χ (Ca) (Vd − VK) IKAHP = gKAHP q (Vd − VK) ISD = −IDS = gC (Vd − Vs) ILeak = gL (V − VL) Maximal conductance parameters were taken (in μ S/cm2) as gNa = 30, gKDR = 15, gKCa = 15, gKAHP = 0.8, gCa = 10, gL = 0.1 and gC = 2.1, while reversal potentials were taken (in mV) as VNa = 60, VK = −75, VCa = 80, and VL = −60. The size of the axosomatic compartment as a proportion of the entire cell was given by p = 0.5 and that of the dendritic compartment as 1 − p. The activation and inactivation gates evolve as a function of their steady state activation x , and time constant τx curves, ∞ where U represents the membrane potentials Vs or Vd, or the intracellular calcium Ca. dx dt = x∞(U ) − x τx (U ) x∞ and τx are often expressed in terms of forward and backward rate functions α and β. αx (U ) x∞(U ) = αx (U ) + βx (U ) 1 τx (U ) = αx (U ) + βx (U ) . In the original formulation, the m, h, n, and s variables are driven solely by continuous rate functions, whereas the c, q, and χ are given as discontinuous rate functions, where H () is the Heaviside step function, 0.32(−46.9 − Vs) αm(Vs) = (exp((−46.9 − Vs)/4) − 1) 0.28(Vs + 19.9) βm(Vs) = (exp((Vs + 19.9)/5) − 1) 0.016(−24.9 − Vs) αn(Vs) = (exp((−24.9 − Vs)/5) − 1) βn(Vs) = 0.25 exp(−1 − 0.025Vs) αh(Vs) = 0.128 exp( (−43 − Vs) ) 18 4 βh(Vs) = (1 + exp((−20 − Vs)/5)) 1.6 αs (Vd) = (1 + exp(−0.072(Vd − 5))) 0.02(Vd + 8.9) βs (Vd) = (exp((Vd + 8.9)/5) − 1) (1 − H (Vd + 10)) exp((Vd + 50)/11 − (Vd + 53.5)/27) αc(Vd) = 18.975 +H (Vd + 10) (2 exp( (−53.5 − Vd) )) 27 βc(Vd) = (1 − H (Vd + 10))(2 exp( (−53.5 − Vd) ) − αc(Vd)) 27 To allow us to perform bifurcation analysis, we approximated the discontinuous functions by fitting continuous functions directly to the steady state and time activation curves. A comparison of the original and fitted curves are shown in Fig. 1, which shows a small relative error for each fitted curve, of the order 10−2 − 10−3. To confirm that the approximated system displayed the same behaviour as the original system, we reproduced the firing patterns and Fig. 1 Approximation of discontinuous function. Approximated functions (a, c, e) and corresponding error of approximation (b, d, f) for c∞, (blue) and τc (red) (A-B), q∞ (blue) and τq (red) (c–d), and χ (blue) (e). For the approximated functions, the original discontinuous functions are plotted in black. Relative error was used for the time constant functions (τc and τq ), where error was computed over a large range of values, while absolute error was used for c∞, q∞ and χ a c e b d f f/I curve of the original model (Fig. 2). The approximated functions are as follows, c∞(Vd) = (1/(1 + exp((−10.1 − Vd)/0.1016)))0.00925, τc(Vd) = 3.627 exp(0.03704 Vd), q∞(Ca) = (0.7894 exp(0.0002726 Ca)) −(0.7292 exp(−0.01672 Ca)), τq (Ca) = (657.9 exp(−0.02023 Ca)) +(301.8 exp(−0.002381 Ca)), χ (Ca) = 1.073 sin(0.003453 Ca + 0.08095) +0.08408 sin(0.01634 Ca − 2.34) +0.01811 sin(0.0348 Ca − 0.9918). Intracellular free calcium in the effective sub-membrane shell is given as dimensionless, since in the original model it was reasoned that the size of the effective shell was undetermined, so appropriate scaling could not be given. Nevertheless, this evolves according to d Ca = −0.13ICa − 0.075Ca. d t Bifurcation and continuation analysis was conducted in XPPAUT, a tool for simulating and analysing dynamical systems (Ermentrout 2002) . One- and two-parameter bifurcation diagrams were constructed using AUTO within XPPAUT. For fast-slow analysis, the differential equations describing Ca and q were independently removed from the full system of equations and instead Ca and q were treated as parameters. All Figures were constructed in MATLAB. The code for the model simulations is available on ModelDB and could be found here: https://senselab.med.yale. edu/ModelDB/ShowModel.cshtml?model=189088. 3 Results 3.1 Response to steady somatic applied current In the original paper, as the applied current to the somatic or dendritic compartment was increased, the authors saw various qualitative shifts in behaviour between resting, bursting and spiking (Pinsky and Rinzel 1994) . Understanding the dynamic mechanisms behind these transitions is important because, within a neural network, individual Pinsky-Rinzel cells do not operate in isolation but rather there are various factors which might affect the net current impinging onto the somatic and dendritic compartments, most notably excitatory and inhibitory synaptic input and also the presence of neuromodulatory factors. Therefore we begin by investigating the effect of somatically or dendritically applied current on the dynamical behaviour of the cell using bifurcation analysis. The bifurcation diagram using ISapp as the initial bifurcation parameter is formed of an S-shaped curve of steady states and a curve of periodic orbits (Fig. 3a). The lower branch of the S-curve consists of stable nodes (see Fig. 3di, ISapp = −1) which correspond to the hyperpolarised resting state of the cell. The system is driven to these nodes predominantly by the potassium-based ILeak, in accordance with experimental work (Bernstein 1902; Hodgkin and Horowicz 1959) . Stability is lost via a saddle-node bifurcation (SN1) at ISapp = 0.02651, which represents the rheobase for ISapp and gives rise to a branch of saddles, which forms the middle and part of the upper branches of the S-curve. This branch of saddles turns around at another saddle-node bifurcation (SN2) at ISapp = −81.57 before regaining stability at the supercritical Hopf bifurcation (HB) for ISapp = 23.69. For increasing ISapp there remains (IDapp = −0.5, gNMDA = 1.25, gC = 1.425). b f/I curves in the original model (black) or approximated model (red) for the isolated soma or dendrite (upper panels) and the 2-compartment model with standard gC = 2.1 (middle panels), infinite gC (lower left panel) and gC = 1.85 (lower right panel) for comparison with the Traub model (Traub et al. 1991) in blue. Circles indicate burst frequency, triangles are spike frequency and squares are mixed burst-spike frequency a branch of stable nodes, corresponding to a depolarised resting state of around -30mV (see Fig. 3di, ISapp = 25). A branch of stable periodic solutions emerges from the HB with low amplitude and high frequency (see Fig. 3dvi, ISapp = 22). These solutions are only apparent in a small range of ISapp since stability is lost for decreasing ISapp between two torus bifurcations (TR1 at ISapp = 21.14 and TR2 at ISapp = 15.87). Within this range, an undulating spiking pattern is observed with increasing and then decreasing spike amplitude (see Fig. 3dv, ISapp = 17). For Fig. 3 Bifurcation analysis for applied somatic current. a Bifurcation diagram with ISapp as the bifurcation parameter. Stable nodes (black line), saddles (red dashed line), stable periodic orbit maxima and minima (green dots), unstable periodic orbit maxima and minima (blue dots), bifurcation points (magenta dots). SN1 saddle node bifurcation 1, SN2 saddle node bifurcation 2, HB supercritical hopf bifurcation, TR1 torus bifurcation 1, TR2 torus bifurcation 2, PD period doubling bifurcation, HC homoclinic bifurcation. b Period of limit cycle solutions. c Phase portrait in Vs and h at the homoclinic bifurcation point. SN1 (magenta), saddle at homoclinic bifurcation (cyan dot). D) Representative traces of somatic voltage at different values of ISapp. Dashed line represents −75 mV for reference a b ISapp below that at TR2, a second branch of stable periodic solutions emerges, representing regular spiking activity (see Fig. 3div, ISapp = 3) with greater spike amplitude but lower frequency than the first branch of stable spiking solutions (Fig. 3b), as described in the original model (Pinsky and Rinzel 1994) . Such spiking is dominated by the dynamics of INa and IK in the somatic compartment, as well described by seminal work by Hodgkin and Huxley (1952) . The aperiodic and spike doubling behaviour for lower ISapp (see Fig. 3diii , ISapp = 2) is brought about via the period doubling (PD) bifurcation which occurs at ISapp = 2.288 and leads to a final branch of unstable periodic solutions. For further decreasing ISapp, aperiodic and spike doubling behaviour transitions to very low frequency (VLF) bursting as in the original model (see Fig. 3dii). The currents involved in shaping these latter behaviours are discussed in more detail later in this paper. Finally, the periodic solutions disappear via a homoclinic bifurcation (HC) at I = −12.35 at which point the period is infinitely large (Fig. 3b). A homoclinic rather than saddle-node on infinite circle (SNIC) bifurcation was confirmed by showing that a very high period (1.27 × 107) limit cycle joins together with a saddle for this value of ISapp, rather than coalescing with the saddle at the SN1 bifurcation (Fig. 3c). It is important to note that for the numerical continuation results presented in Figs. 3a and 4a, we consider a very large range of the bifurcation parameter values, namely ISapp/Dapp ∈ [−500, 500], in order to locate the Hopf bifurcations giving rise to stable spiking solutions. The large range of the bifurcation parameter increases the likelihood of finding turning points for the curve of steady states. Interpreting the entirety of this range as meaningful should be cautioned. It is extremely unlikely that the currents required to reach the solutions between the Torus bifurcations and between the Torus and Hopf bifurcations in the model would ever be physiologically relevant, given that the unitary amplitude of a single CA3-CA3 monosynaptic EPSC is of the order of 20–30 pA (Debanne et al. 1995) . Therefore these solutions are not considered any further. 3.2 Response to steady dendritic applied current The bifurcation diagram, using IDapp as the bifurcation parameter (Fig. 4a), is qualitatively very similar to that for ISapp, although there are some differences mainly in the curve of steady states. The bottom branch of this curve consists of stable nodes, representing the hyperpolarised resting state (see Fig. 4di, IDapp = −1). Stability is again lost via a Fig. 4 Bifurcation analysis for applied dendritic current. a Bifurcation diagram with IDapp as the bifurcation parameter. Stable nodes (black line), saddles (red dashed line), stable periodic orbit maxima and minima (green dots), unstable periodic orbit maxima and minima (blue dots), bifurcation points (magenta dots). SN1 saddle node bifurcation 1, SN2 saddle node bifurcation 2, HB supercritical hopf bifurcation, SN3 saddle node bifurcation 3, TR1 torus bifurcation 1, TR2 torus bifurcation 2, HC homoclinic bifurcation. b Period of limit cycle solutions. c Phase portrait in Vs and h at the homoclinic bifurcation point. SN1 represented by magenta dot and saddle at homoclinic bifurcation represented by cyan dot. d Example traces of somatic voltage at different values of IDapp. Dashed line represents −75 mV for reference a b c saddle-node bifurcation (SN1) at IDapp = 0.02728, leading to a branch of saddles which again turns around at another saddle-node bifurcation (SN2) at IDapp = −83.33, before regaining stability again via a supercritical Hopf bifurcation (HB) at IDapp = 99.78. Notably, IDapp at the HB is much greater than that for ISapp, and thus there is a larger range of IDapp for which there are periodic solutions, despite the rheobase being similar. The branch of stable nodes arising from the HB, which represent the depolarised resting state (see Fig. 4di, IDapp = 120) this time loses stability at a third saddle-node bifurcation (SN3) at IDapp = 127.6, giving rise to a final branch of saddles. The branch of stable periodic spiking solutions emanating from the HB again grows in amplitude and period with decreasing IDapp (see Fig. 4dvi, IDapp = 70) until stability is lost between the two torus bifurcations (TR1 at IDapp = 15.59 and TR2 at IDapp = 28.75). Within this unstable branch, there is again a form of waxing and waning spike amplitude (see Fig. 4dv, IDapp = 25). For decreasing IDapp below that of TR2, a branch of stable, regular periodic spiking solutions appears (see Fig. 4div, IDapp = 12) before stability is lost again for IDapp < 9.127 at a period doubling bifurcation (PD). The final unstable branch of periodics for low IDapp disappears at IDapp = −3.486 again via a homoclinic bifurcation (HC). Within this unstable branch, the system transitions for decreasing IDapp from doublet activity for IDapp between ≈ 7-9, to a mixture of single spikes and doublets for IDapp of ≈ 6, to then aperiodic activity with a varying number of spikes per burst between IDapp of ≈ 2.6 and 5 (see Fig. 4diii, IDapp = 4), to finally periodic VLF bursting between HC and IDapp of ≈ 2.5 (see Fig. 4dii, IDapp = 0.3). Again the HC was confirmed by showing a very high period limit cycle solution (1.77 × 107) for this value of IDapp (Fig. 4b), which we use as an approximation of the homoclinic orbit, and with which the SN1 point did not coalesce (Fig. 4c). It is not entirely surprising that the bifurcation diagrams presented in Figs. 3a and 4a are qualitatively very similar, given that the two compartments are coupled quite closely. Nevertheless, there is clearly a larger range of IDapp within which the system is between the hyperpolarised resting and stable spiking states. This suggests that preferentially dendritic depolarisation could engage a wider repertoire of active behaviours and that dendritic currents could be in a more privileged position to shape these behaviours. 3.3 Dependence of model cell behaviour on gCa Indeed lowering the maximal conductance through the dendritically-located, persistent calcium channel (gCa) has been shown to switch bursting to spiking activity for a set value of applied current (Pinsky and Rinzel 1994; Traub et al. 1991) . A change in gCa from 10 to 7 is sufficient to switch the behaviour of the Pinksy-Rinzel cell from intrinsic bursting to tonic firing with frequency accommodation to better represent the firing properties of CA1 cells (Taxidis et al. 2012). We therefore investigated how gCa affects the bifurcation structure of the system, in an attempt to understand this change mathematically. We initially re-computed the bifurcation diagrams in ISapp and IDapp for the CA1 value of gCa = 7 (Fig. 5). Using ISapp as the bifurcation parameter, the bifurcation diagram (Fig. 5a) is extremely similar to that of the CA3 cell with one noticeable difference; instead of the periodic solutions disappearing for low ISapp via a HC bifurcation, they disappear via a saddle-node on an invariant circle (SNIC) bifurcation. The saddle-node bifurcation of equilibrium solutions corresponding to this value of ISapp is that of SN1. This was confirmed by depicting in the Vs − h phase plane that the SN1 equilibrium bifurcation lies on the limit cycle at the SNIC bifurcation (Fig. 5c) and that this limit cycle used to approximate the SNIC orbit has high period (5.65 × 105) (Fig. 5b). In addition, whereas a PD bifurcation was associated with the loss of stability of periodic solutions in the CA3 cell, no such transition was detected for the CA1 cell. Using IDapp as the bifurcation parameter, there are two noticeable differences between the CA3 and CA1 cell. Firstly while the upper branch of unstable steady states did not bifurcate back to a stable solution in the range of −500 ≤ IDapp ≥ 500 for the CA3 cell, in the CA1 cell a b d e c f Fig. 5 Bifurcation analysis of the CA1 cell for applied somatic and dendritic current. NB. For these diagrams gCa = 7 a Projection of the bifurcation diagram with ISapp as the bifurcation parameter. Stable nodes (black line), saddles (red dashed line), stable periodic orbit maxima and minima (green dots), unstable periodic orbit maxima and minima (blue dots), bifurcation points (magenta dots). SN1 saddle node bifurcation 1 (ISapp = 0.0557), SN2 saddle node bifurcation 2 (ISapp = −81.11), HB supercritical hopf bifurcation (ISapp = 24.01), TR1 torus bifurcation 1 (ISapp = 18.73), TR2 torus bifurcation 2 (ISapp = 17.37), SNIC saddle-node on invariant circle bifurcation (ISapp = 0.0556). b Period of limit cycle solutions for ISapp as the bifurcation parameter. c Phase portrait in Vs and h at the SNIC bifurcation for ISapp. SN1 is represented by a magenta dot. d Projection of the bifurcation diagram with IDapp as the bifurcation parameter. SN1 at ISapp = 0.05745, SN2 at IDapp = −83.33, HB at IDapp = 141.0, SN3 saddle node bifurcation 3 at IDapp = 288.3, SN4 saddle node bifurcation 4 at IDapp = −175.2 and SNIC at IDapp = 0.0574. e Period of limit cycle solutions for IDapp as the bifurcation parameter. f Phase portrait in Vs and h at the SNIC bifurcation for IDapp. SN1 is represented by a magenta dot a fourth saddle-node bifurcation was detected for IDapp = −175.2, where the steady state solutions turn around again, giving rise to a second depolarised steady state. This means that in the range of IDapp between SN1 and SN3, the system is bistable (either periodic spiking and depolarised steady states between SN1 and HB, or two depolarised steady states between HB and SN3, depending on the initial parameters). It is noted, however, that a similar behaviour might be observed in the CA3 cell if the bifurcation analysis was extended beyond −500 ≤ IDapp ≤ 500 range. Secondly, no torus bifurcations were detected for the periodic spiking solutions in the CA1 cell, whereas a pair of torus bifurcations were present in the CA3 cell. Whilst this may be of mathematical interest, as discussed above, the range of IDapp for which this difference manifests suggests it has no physiological role. As for ISapp, the periodic solution disappears for low IDapp via a SNIC bifurcation instead of a HC bifurcation, as was the case for the CA3 cell. Therefore for all ISapp/IDapp between the SNIC at ISapp = 0.0556, IDapp = 0.0574 and the HB at ISapp = 24.01, IDapp = 141.0 (excluding the regions between the pair of torus bifurcations for ISapp), the CA1 cell limit cycle solutions will be periodic spiking with decreasing amplitude and period as ISapp/IDapp increase, with importantly no regimes of bursting or aperiodic and spike doubling behaviours. The disappearance of these regimes is likely directly attribute to the presence of a codimension-2 bifurcation in the (I , gCa) parameter space (Fig. 6a-b), which causes the periodic solution to switch from vanishing via a HC to vanishing via a SNIC, as gCa is decreased. The dependence of this regime change on gCa also implies that Ca plays an important role in shaping both irregular spiking and bursting behaviours. To examine this further, we initially explored how switches between singlet activity and doublet activity during the irregular spiking regime may be shaped by Ca, for example when ISapp = 2. As shown in Fig. 7, during the transition between repetitive singlet activity and doublet activity, the Vd and Ca peaks phase advance but to different extents, leading to a reduction in the (Ca−Vd) phase difference and an increase in the peak IKCa which drives the membrane potential to a more hyperpolarised level. This train of events to shift singlet activity to doublet activity seems plausible given changes in gKCa , and hence IKCa , have previously been shown to shift a cell’s behaviour from spiking to bursting (Tsaneva-Atanasova et al. 2007; Nowacki et al. 2010; Szalai et al. 2011; Iosub et al. 2015) . 3.4 Mechanisms of bursting Finally we study how Ca and the Ca-dependent q dynamics shape the bursting behaviour in the model. In the original paper, the VLF bursting dynamics of the CA3 cell are described to be largely determined by q which, when passes below a threshold, was thought to trigger the initial somatic spike of the burst. As such, bursting was not thought to persist when q was replaced by its mean value (Pinsky and Rinzel 1994) . Since then, the idea that q acts as the threshold for burst initiation has been called into question from phaseplane analysis on a piecewise reduced model (Booth and Bose 2001; Bose and Booth 2005) . We therefore performed fast-slow analysis to probe the mathematical mechanism behind the bursting dynamics further. Since both Ca and q decay on slower timescales (13.33ms and 100 − 1000ms respectively) than the other gating variables (< 6ms within the effective ranges of Vs and Vd see Fig. 8a) (Pinsky and Rinzel 1994)) and therefore the ratio of timescales between fast and slow variables is << 1, this analysis is appropriate Fig. 6 Projection of the 2 parameter continuation in ISapp (a) or IDapp (b) and gCa . Solid black line - SN1, dotted black line - SN2, dashed black line SN3, green line - HB, blue line HC, magenta dot - SNIC bifurcation for CA1 cell. Inset zoomed in at the codimension-2 SNIC bifurcation where the HC and SN1 curves meet. To construct these diagrams, the bifurcations found in I were then continued in the gCa parameter space such that a horizontal slice through these diagrams (e.g. gCa = 10) would be the same as the 1-parameter bifurcation diagram in I for this value of gCa (i.e. that of Figs. 3–4) a b Fig. 7 Role of Ca in irregular spiking activity. a Example trace of a irregular spiking activity whereby a train of singlets transition into doublet and triplet activity. Vd (red), Ca (blue) and IKCa (magenta). The oscillation used to calculate phase is shown above (black). Demarked atop of the oscillation are the early singlet Vd (red x) and Ca (blue x) peak times. The oscillation period was set by the mean inter-spike interval of these early singlet Vd peaks, and was shifted to minimise the distance between the oscillation peak and the early singlet Vd peaks. Notice how the last singlet Vd (red o) and Ca (blue o) peaks in the singlet train are phase-advanced in the oscillation. b Polar plot of the Vd and Ca peak phases taken from 7 cases where there are at least 5 consecutive singlets before transition to doublet activity. 0 degrees is the peak of the oscillation and 180 is the trough. The radial axis (0−1) shows the normalised Vd or Ca amplitude (relative to maximum range). The last singlet Vd and Ca peaks (o) are consistently phase-advanced compared to the early singlet peaks (x). c Average Ca − Vd phase difference against the average peak IKCa . There is a fairly consistent Ca − Vd phase difference and IKCa for the early singlets (x). However for the last singlet in the train (o), the Ca −Vd phase difference is reduced and this corresponds to a higher IKCa to investigate how both the Ca and q dynamics influence bursting. The fast-slow method was originally pioneered by Rinzel (1987) , and involves separating the full system into a fast subsystem and a slow subsystem, where the slow variables (in this case Ca and q) can be used as parameters for bifurcation analysis of the fast subsystem. We analysed the VLF bursting dynamics for both applied somatic current and applied dendritic current. ISapp or IDapp were set to 0.3 with the other equal to 0 (values that generates VLF bursting in the full system). The same qualitative bifurcation structure was observed in both cases. Therefore we present only the bifurcation diagrams for ISapp, whilst describing both bifurcation values for the fast subsystem, where ISapp = 0.3/IDapp = 0.3. We begin our fast-slow analysis using q as the bifurcation parameter. For steady applied current, the bifurcation diagram of the fast subsystem (Fig. 8b) is formed of a curve of steady states and a curve of bursting periodic orbits. Since the system continues to burst while q is fixed, the intraburst dynamics are independent of q, and this variable (q) cannot act as a threshold for burst initiation, confirming results in (Booth and Bose 2001; Bose and Booth 2005) . To elaborate, the curve of steady states has a lower branch of stable nodes, representing the hyperpolarised resting state, which loses stability for decreasing q via a saddle node on invariant circle (SNIC) bifurcation at q = 0.1136/q = 0.1119 (see Fig. 8d) and gives rise to an upper branch of saddles. From the SNIC, a stable Fig. 8 Bifurcation structure in q for CA3 bursting. a Gating variable time constants within the effective range of Vs (i), for τn (black) and τh (red); Vd (ii) for τs (blue) and τc (green); and Ca (iii) for τq (magenta). b Bifurcation diagram of the fast subsystem with q as the bifurcation parameter. Inset zoomed in around the values of interest for q. Stable nodes (black line), saddles (dotted red line), unstable bursting orbit maxima and minima (blue lines), stable bursting orbit (green lines), bifurcation points (magenta dots), and burst trajecory in the q − Vs plane (cyan). SNIC saddle node on invariant cicle bifurcation; PD12 period doubling bifurcations 1-2; SNP saddle node bifurcation of periodic solutions; HC homoclinic bifurcation. c Time course of the burst for Vs (black) and q (blue). A dotted red line representing the value of q at the SNIC is overlayed. Phase portrait in Vs and h at the SNIC bifurcation (magenta dot) (d) and HC bifurcation (e) - note the lack of overlap between the SNIC (magenta dot) and the saddle at the HC (cyan dot). f Period of bursting solutions a b c d f e branch of bursting periodic solutions also emerges with initially high period (Fig. 8f). As q decreases, the amplitude of the burst remains fairly constant while the period decreases. Stability is lost at q = −0.05923/q = −0.06627 via a period doubling bifurcation (PD1) then regained briefly via a second period doubling bifurcation (PD2 at q = −0.4779/q = −0.4844) before being lost again via a saddle node of periodics (SNP) bifurcation at q = −0.5024/ q = −0.528. At this point, the bursting periodic orbit branch turns around and increases in period for increasing q (Fig. 8f) until it terminates via a homoclinic bifurcation (HC at q = 0.2524/q = 0.2653, see Fig. 8e). Following the burst trajectory (shown in cyan in Fig. 8b), the phase point tightly follows the branch of stable nodes for decreasing q and then jumps to the maxima of the bursting trajectory once past the SNIC. In the full system, for this value of applied current, this transition takes several hundred ms (Fig. 8c), again confirming that the threshold for burst initiation is not described by the q dynamics. Midway through the burst, q then begins to increase which pushes the system back through the SNIC. After this point, the phase point transiently oscillates before returning to the branch of stable nodes. Therefore since the transition below the SNIC in q permits the occurrence of the burst, the time between q increasing above the SNIC and then decreasing below the SNIC strongly controls at least part of the interburst interval; a finding consistent with previous work (Bose and Booth 2005) . Given that q could not fully describe the intraburst dynamics, we then turned our attention to the second slow variable, Ca and performed fast-slow analysis using Ca as the bifurcation parameter. For steady applied current, the bifurcation diagram in Ca (Fig. 9a) has two long branches of stable nodes, an upper branch representing a depolarised steady state for values of Ca < 127.5/Ca < 127.6 and a lower branch representing a hyperpolarised steady state for Ca > 4.263/Ca > 4.117. Stability in the lower branch is lost for decreasing Ca via a saddle node bifurcation (SN1) at Ca = 4.263/Ca = 4.117. The resulting unstable branch of saddles turns around and continues for larger Ca values. Stability in the upper branch is also lost via a saddle node bifurcation (SN2) at Ca = 127.5/Ca = 127.6. The branch turns around briefly and is formed of a series of saddles until stability is regained momentarily between another couple of saddle node bifurcations (SN3 at Ca = 112.5/Ca = 112.6 and SN4 at Ca = 127.2/Ca = 127.4). After SN4, the upper branch consists of a second series of saddles which turn around at a fifth saddle-node bifurcation (SN5 at Ca = 62.76/Ca = 63.73) before continuing on to larger values of Ca. Between SN4 and SN5, there is a subcritical Hopf bifurcation (HB) at Ca = 112.7/Ca = 113.9 which gives rise to the unstable periodic solutions. For decreasing Ca, the amplitude and period of the solution increases (Fig. 9c), before turning around at the saddle node of periodics (SNP a at Ca = 11.21/Ca = 11.92) and terminating via a homoclinic bifurcation (HC at Ca = 14.58/Ca = 15.23) (Fig. 9d). Following the burst trajectory, it can be seen that the phase point, loosely follows the lower branch of stable nodes for decreasing Ca and overshoots SN5 before jumping to the upper branch of stable nodes. Again the attraction is weak, because the Ca dynamics are only intermediately slow, and so the phase point transiently oscillates around this upper branch. Even when the phase point is past the HB, there are smaller oscillations around the unstable upper branch before returning to the attraction of the lower branch of steady states. VLF bursting is therefore governed by the Ca dynamics and is of the fold/subHopf type, adoping the nomenclature in Izhikevich (2000) , which has been reported also for biophysical models of pituitary somatotrophs and lactotrophs (Stern et al. 2008; Tsaneva-Atanasova et al. 2007; Tabak et al. 2007) . 4 Discussion 4.1 Mechanisms of spiking The discontinuous functions in the original Pinsky-Rinzel model were modified to facilitate numerical bifurcation analysis on the full system (Fig. 1), whilst maintaining the qualitative behaviour of the original model (Fig. 2). We show that spiking behaviour emerges in the model from b c d Fig. 9 Bifurcation structure in Ca for CA3 bursting. a Bifurcation diagram of the fast subsystem with Ca as the bifurcation parameter. Inset - zoomed in around the hopf bifurcation. Stable nodes (black line), saddles (dotted red line), unstable periodic orbit maxima and minima (blue line), bifurcation points (magenta dots), and burst trajectory in the Ca − Vs plane (cyan). SN1-5 saddle node bifurcations 1-5; HB subcritical hopf bifurcation; SNP saddle node bifurcation of periodic solutions; HC homoclinic bifurcation. b Time course of the burst for Vs (black) and Ca (blue). C) Period of limit cycle solutions. d Phase portrait in Vs and h at the homoclinic bifurcation point. Note the lack of overlap between SN1 (magenta dot) and the saddle at the HC (cyan dot) a supercritical Hopf bifurcation and that the rheobase is formed from a saddle-node bifurcation (Figs. 3 and 4). This mechanism of spiking is therefore similar to the bifurcation structure, described earlier, for a single compartmental CA3 model (Grobler et al. 1998) . As gCa is decreased, a codimension-two SNIC bifurcation switches the periodic spiking solutions from terminating via a homoclinic bifurcation to terminating via a SNIC bifurcation (Fig. 6). This eliminates the unstable periodic regime, where bursting and irregular spiking solutions exist, and thus explains why a reduction in gCa leads to unperturbed spiking for increasing IDapp and a larger parameter range of ISapp (Fig. 5), as previously observed (Pinsky and Rinzel 1994; Traub et al. 1991) . The disappearance of irregular spiking as gCa decreases, is suggestive that Ca, as well as shaping bursting behaviour, is also important in shaping spiking in the model. Indeed we show that a reduction in the (Ca-Vd) phase difference, and corresponding increase in IKCa is likely to drive activity from singlet spiking to doublet/triplet spiking in this irregular spiking regime (Fig. 7). As an example of how understanding where PinskyRinzel cells sit dynamically is important for understanding spiking in larger networks, we turn out attention to a combined CA3 and CA1 hippocampal network model, of Pinsky-Rinzel cells that exhibits sharp wave ripple (SWR) oscillations (Taxidis et al. 2012) . SWRs are transient network events that occur during periods of awake immobility and slow wave sleep (Buzsaki et al. 1983; Buzsaki 1986; O’Neill et al. 2010) , and are known to be important for spatial learning and memory (Girardeau et al. 2009; EgoStengel and Wilson 2010; Jadhav et al. 2012) . During sharp wave ripples, hippocampal place cells that fire in discrete locations of an environment (O’Keefe and Dostrovsky 1971) , reactivate such that co-active place cell activity, representing a discrete environmental location, or sequences of place cell activity, representing a trajectory through an environment, are replayed (Wilson and McNaughton 1994; Nadasdy et al. 1999; Lee and Wilson 2002; Foster and Wilson 2006; Csicsvari et al. 2007; Diba and Buzsaki 2007; Davidson et al. 2009; Karlsson and Frank 2009) . However the mechanisms that select which place cells are active in SWRs is still unclear, although the intrinsic excitability of the cells and the local synaptic connectivity structure of the network within which they reside are both likely to play a role (for review see Atherton et al. (2015)). In the Taxidis model (Taxidis et al. 2012) , CA1 cells are driven to spike by CA3 cell activity. From a dynamical systems perspective, the occurrence of a single CA1 spike could be influenced by CA3 in two ways. Firstly, if the CA1 cell was in a spiking limit cycle regime, excitatory input from CA3 could phase shift CA1 cell spike times. Secondly, if the CA1 cell was in a hyperpolarised resting state, excitatory input from CA3 could push the CA1 cell trajectory through the phase space in a manner that resembles a spike, but which does not actually push the CA1 cell trajectory onto a stable limit cycle i.e. causes a transient. The bifurcation analysis, conducted here, facilitates the differentiation between these quite different mechanisms of CA1 spiking, namely periodic attractor driven verses transient oscillations. The analysis shows that for the parameter values of ISapp and gCa used in the model (0 and 7 respectively), the CA1 cells are in a hyperpolarised resting regime below the SNIC. Therefore at least in this model, it can be concluded that CA1 cell spiking in SWRs is driven by transient oscillations, induced by presynaptic CA3 cell activity, and not by a periodic attractor driven excitability. 4.2 Mechanisms of bursting Experimental findings have shown that the intrinsic bursting of CA3 cells is driven by several processes: the intital fast spike is driven by activity through sodium channels; the ensuing after depolarising potential (ADP) subsequently drives the burst; and the burst is then terminated by an after hyperpolarising potential (AHP) (Wong and Prince 1978; 1981) . Although the ionic bases for the ADP and AHP that occur during bursting in CA3 pyramidal cells are not fully understood, they are both calcium driven (Wong and Prince 1978; 1981; Schwartzkroin and Stafstrom 1980) and are likely mediated by T/R-type calcium channels and calciumactivated potassium channels respectively, as has been found in some cases for CA1 cell bursting (Magee and Carruth 1999; Metz et al. 2005; Alger and Nicoll 1980) , although see (Azouz et al. 1996). Kv7/KCNQ/M-channels also play a role in terminating the ADP (Vervaeke et al. 2006; Brown and Randall 2009) . Previous studies using fast-slow analysis of CA3 pyramidal cell models have proposed different mechanisms as underlying bursting activity (Kepecs and Wang 2000; Xu and Clancy 2008) . In a reduced version of the Pinsky-Rinzel model with similar somatic currents but only a persistent sodium current and a slow potassium current in the dendritic compartment; bursting activity was found to be of the square-wave/fold-homoclinic type and dependent on the activation variable of the slow potassium current (Kepecs and Wang 2000). In contrast, using fast-slow analysis on a single-compartmental CA3 model, incorporating various additional currents to the Pinsky-Rinzel model, bursting was found to be of the parabolic type and dependent on two variables: the slow autocatalytic activation variable for a T-type calcium channel and the activation variable for the slow calcium-dependent potassium current (Xu and Clancy 2008) . It remained to be established which, if either, of these mechanisms is true for the original Pinsky-Rinzel model, where the two candidate slow variables for controlling bursting are calcium and the activation gate, q, of the dendritic after hyperpolarisation current. Here, we identify a mechanism intermediate between the two cases, whereby bursting is of the fold/sub-Hopf type and relies on just one variable, calcium, which is dependent on the dendritic calcium current, ICa, and in fact controls both the calciumdependent potassium current, IKCa , and the activation gate, q, of the dendritic after hyperpolarisation current, IKAHP . The rise of Ca through the subcritical Hopf bifurcation mediates the slow termination of the burst (Fig. 9), thereby biophysically explaining previous results on a piecewise reduced system, showing disruption of the ping-pong mechanism and termination of the burst as ICa was activated (Booth and Bose 2001; Bose and Booth 2005) . The results, here, are also consistent with the blockade of calcium channels perturbing CA3 bursting Wong and Prince (1978, 1981) which could not be explained by bursting dependent on activation of a slow potassium current alone (Kepecs and Wang 2000). The original Pinsky-Rinzel paper proposed that q sets the burst initiation threshold (Pinsky and Rinzel 1994) . Subsequent work, using the piecewise reduced system, showed this threshold was independent of q and instead set by Vd (Booth and Bose 2001; Bose and Booth 2005) . We clarify this point here by showing, using fast-slow analysis, that while q plays no role in setting the threshold for the initial somatic spike in the burst or indeed the intraburst dynamics, it does set part of the interburst interval (Fig. 8). The biophysical implications of this is that the dendritic afterhyperpolarisation current plays an important role in setting the interburst interval. Nevertheless, it it important to note that the dependence of bursting on just one slow variable is likely due to the bistability of the fast subsystem in the Pinsky-Rinzel model which is driven by the window INa (Golomb et al. 2006) . This derives from the overlap in the steady-state activation and inactivation curves for the INa current, where the inactivation curve is shifted to more depolarised values (> 20mV ) compared to experimental results (Sah et al. 1988) . When this window current was reduced, bursting in a different CA3 model was instead dependent on two slow variables and was of the parabolic type (Xu and Clancy 2008) . 5 Conclusion By proposing a modified smooth system of differential equations that can be used to further analyse the full PinskyRinzel model, we hope that future studies using PinskyRinzel cells as part of a network, or as a means to analyse firing patterns will have a tool with which to systematically identify parameter regimes of interest to their research question. Characterising the bifurcation structure for applied current, calcium conductance and applying fast-slow analysis to investigate the bursting regime in the systems already provides insight into the parameter dependence of the model dynamics. Acknowledgments LAA is supported by the Engineering and Physical Sciences Research Council (EPSRC) and Eli Lilly & Company; LYP is supported by the Wellcome Trust; and KT-A is supported by grant EP/N014391/1 of the EPSRC. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Compliance with ethical standards Conflict of interest The authors declare that they have no conflict of interest. Alger , B.E. , & Nicoll , R.A. ( 1980 ). Epileptiform burst afterhyperpolarization - calcium-dependent potassium potential in hippocampal CA1-pyramidal cells . Science , 210 ( 4474 ), 1122 - 1124 . Atherton , L.A. , Dupret , D. , & Mellor , J.R. ( 2015 ). Memory trace replay: the shaping of memory consolidation by neuromodulation . Trends in Neurosciences , 38 ( 9 ), 560 - 570 . Azouz , R. , Jensen , M.S. , & Yaari , Y. ( 1996 ). Ionic basis of spike afterdepolarization and burst generation in adult rat hippocampal CA1 pyramidal cells . Journal of Physiology-London , 492 ( 1 ), 211 - 223 . Bernstein , J. ( 1902 ). Tests on the thermodynamics of bioelectric currents . Pflu¨gers Archiv , 92 ( 10 - 12 ), 521 - 562 . Bittner , K.C. , Grienberger , C. , Vaidya , S.P. , Milstein , A.D. , Macklin , J.J. , Suh , J. , Tonegawa , S. , & Magee , J.C. ( 2015 ). Conjunctive input processing drives feature selectivity in hippocampal CA1 neurons . Nature Neuroscience , 18 ( 8 ), 1133 . Booth , V. , & Bose , A. ( 2001 ). Neural mechanisms for generating rate and temporal codes in model CA3 pyramidal cells . Journal of Neurophysiology , 85 ( 6 ), 2432 - 2445 . Bose , A. , & Booth , V. ( 2005 ). Bursting in 2-compartment neurons: a case study of the pinsky-rinzel model . Bursting: The Genesis of Rhythm in the Nervous System , 123 - 144 . Brown , J.T. , & Randall , A.D. ( 2009 ). Activity-dependent depression of the spike after-depolarization generates long-lasting intrinsic plasticity in hippocampal CA3 pyramidal neurons . Journal of Physiology-London , 587 ( 6 ), 1265 - 1281 . Buchanan , K.A. , & Mellor , J.R. ( 2010 ). The activity requirements for spike timing-dependent plasticity in the hippocampus . Frontiers in synaptic neuroscience , 2 , 11 - 11 . Burgess , N. , Maguire , E.A. , & O'Keefe , J. ( 2002 ). The human hippocampus and spatial and episodic memory . Neuron , 35 ( 4 ), 625 - 641 . Buzsaki , G. ( 1986 ). Hippocampal sharp waves - their origin and significance Brain Research 398 ( 2 ). Buzsaki , G. , Leung , L.W. , & Vanderwolf , C.H. ( 1983 ). Cellular bases of hippocampal eeg in the behaving rat . Brain research , 287 ( 2 ), 139 - 71 . Csicsvari , J. , O 'Neill , J. , Allen , K. , & Senior , T. ( 2007 ). Place-selective firing contributes to the reverse-order reactivation of CA1 pyramidal cells during sharp waves in open-field exploration . European Journal of Neuroscience , 26 ( 3 ), 704 - 716 . Davidson , T.J. , Kloosterman , F. , & Wilson, M.A. ( 2009 ). Hippocampal replay of extended experience Neuron 63 ( 4 ). Debanne , D. , Guerineau , N.C. , Gahwiler , B.H. , & Thompson , S.M. ( 1995 ). Physiology and pharmacology of unitary synaptic connections between pairs of cells in areas ca3 and ca1 of rat hippocampal slice cultures . Journal of Neurophysiology , 73 ( 3 ), 1282 - 1294 . Diba , K. , & Buzsaki , G. ( 2007 ). Forward and reverse hippocampal place-cell sequences during ripples . Nature Neuroscience , 10 ( 10 ), 1241 - 1242 . Ego-Stengel , V. , & Wilson, M.A. ( 2010 ). Disruption of rippleassociated hippocampal activity during rest impairs spatial learning in the rat . Hippocampus , 20 ( 1 ), 1 - 10 . Epsztein , J. , Brecht , M. , & Lee , A.K. ( 2011 ). Intracellular determinants of hippocampal CA1 place and silent cell activity in a novel environment . Neuron , 70 ( 1 ), 109 - 120 . Ermentrout , B. ( 2002 ). Simulating, analyzing, and animating dynamical systems: a guide to XPPAUT forresearchers and students . Society for Industrial and Applied Mathematics Philadelphia. Foster , D.J. , & Wilson, M.A. ( 2006 ). Reverse replay of behavioural sequences in hippocampal place cells during the awake state Nature 440 ( 7084 ). Girardeau , G. , Benchenane , K. , Wiener , S.I. , Buzsaki , G. , & Zugaro , M.B. ( 2009 ). Selective suppression of hippocampal ripples impairs spatial memory . Nature Neuroscience , 12 ( 10 ), 1222 - 1223 . Golomb , D. , Yue , C. , & Yaari , Y. ( 2006 ). Contribution of persistent na+ current and M-type K+ current to somatic bursting in CA1 pyramidal cells: Combined experimental and modeling study . Journal of Neurophysiology , 96 ( 4 ), 1912 - 1926 . Grobler , T. , Barna , G. , & Erdi , P. ( 1998 ). Statistical model of the hippocampal CA3 region - i. the single-cell module: bursting model of the pyramidal cell . Biological Cybernetics , 79 ( 4 ), 301 - 308 . Hahn , P.J. , & Durand , D.M. ( 2001 ). Bistability dynamics in simulations of neural activity in high-extracellular-potassium conditions . Journal of Computational Neuroscience , 11 ( 1 ), 5 - 18 . Harvey , C.D. , Collman , F. , Dombeck , D.A. , & Tank , D.W. ( 2009 ). Intracellular dynamics of hippocampal place cells during virtual navigation . Nature , 461 ( 7266 ), 941 - U196 . Hodgkin , A.L. , & Horowicz , P. ( 1959 ). The influence of potassium and chloride ions on the membrane potential of single muscle fibres . Journal of Physiology-London , 148 ( 1 ), 127 - 160 . Hodgkin , A.L. , & Huxley , A.F. ( 1952 ). A quantitative description of membrane current and its application to conduction and excitation in nerve . Journal of Physiology-London , 117 ( 4 ), 500 - 544 . Iosub , R. , Avitabile , D. , Grant , L. , Tsaneva-Atanasova , K. , & Kennedy , H.J. ( 2015 ). Calcium-induced calcium release during action potential firing in developing inner hair cells . Biophysical Journal , 108 ( 5 ), 1003 - 1012 . Izhikevich , E.M. ( 2000 ). Neural excitability, spiking and bursting . International Journal of Bifurcation and Chaos , 10 ( 6 ), 1171 - 1266 . Jadhav , S.P. , Kemere , C. , German , P.W. , & Frank , L.M. ( 2012 ). Awake hippocampal sharp-wave ripples support spatial memory . Science , 336 ( 6087 ), 1454 - 1458 . Kamondi , A. , Acsady , L. , Wang , X.J. , & Buzsaki , G. ( 1998 ). Theta oscillations in somata and dendrites of hippocampal pyramidal cells in vivo: Activity-dependent phase-precession of action potentials . Hippocampus , 8 ( 3 ), 244 - 261 . Kandel , E.R. , & Spencer , W.A. ( 1961 ). Electrophysiology of hippocampal neurons .2. after-potentials and repetitive firing . Journal of Neurophysiology , 24 ( 3 ), 243 . Karlsson , M.P. , & Frank , L.M. ( 2009 ). Awake replay of remote experiences in the hippocampus . Nature Neuroscience , 12 ( 7 ), 913 - 32 . Kepecs , A. , & Wang , X.J. ( 2000 ). Analysis of complex bursting in cortical pyramidal neuron models . Neurocomputing , 32 , 181 - 187 . Le Duigou , C. , Simonnet , J. , Telenczuk , M.T. , Fricker , D. , & Miles , R. ( 2014 ). Recurrent synapses and circuits in the CA3 region of the hippocampus: an associative network . Frontiers in Cellular Neuroscience , 7 . Lee , A.K. , & Wilson, M.A. ( 2002 ). Memory of sequential experience in the hippocampus during slow wave sleep Neuron 36(6 ). Lisman , J.E. ( 1997 ). Bursts as a unit of neural information: Making unreliable synapses reliable . Trends in Neurosciences , 20 ( 1 ), 38 - 43 . Magee , J.C. , & Carruth , M. ( 1999 ). Dendritic voltage-gated ion channels regulate the action potential firing mode of hippocampal CA1 pyramidal neurons . Journal of Neurophysiology , 82 ( 4 ), 1895 - 1901 . Mainen , Z.F. , & Sejnowski, T.J. ( 1996 ). Influence of dendritic structure on firing pattern in model neocortical neurons . Nature , 382 ( 6589 ), 363 - 366 . McNaughton , B.L. , & Morris , R.G.M. ( 1987 ). Hippocampal synaptic enhancement and information-storage within a distributed memory system . Trends in Neurosciences , 10 ( 10 ), 408 - 415 . Metz , A.E. , Jarsky , T. , Martina , M. , & Spruston , N. ( 2005 ). R-type calcium channels contribute to afterdepolarization and bursting in hippocampal CA1 pyramidal neurons . Journal of Neuroscience , 25 ( 24 ), 5763 - 5773 . Migliore , M. , Cook , E.P. , Jaffe , D.B. , Turner , D.A. , & Johnston , D. ( 1995 ). Computer-simulations of morphologically reconstructed CA3 hippocampal-neurons . Journal of Neurophysiology , 73 ( 3 ), 1157 - 1168 . Miles , R. , & Wong , R.K.S. ( 1983 ). Single neurons can initiate synchronized population discharge in the hippocampus . Nature , 306 ( 5941 ), 371 - 373 . Mizuseki , K. , Royer , S. , Diba , K. , & Buzsaki , G. ( 2012 ). Activity dynamics and behavioral correlates of CA3 and CA1 hippocampal pyramidal neurons . Hippocampus , 22 ( 8 ), 1659 - 1680 . Morris , R.G.M. , Garrud , P. , Rawlins , J.N.P. , & Okeefe , J. ( 1982 ). Place navigation impaired in rats with hippocampal-lesions Nature 297 ( 5868 ). Nadasdy , Z. , Hirase , H. , Czurko , A. , Csicsvari , J. , & Buzsaki , G. ( 1999 ). Replay and time compression of recurring spike sequences in the hippocampus . Journal of Neuroscience , 19 ( 21 ), 9497 - 9507 . Nakazawa , K. , Quirk , M.C. , Chitwood , R.A. , Watanabe , M. , Yeckel , M.F. , Sun , L.D. , Kato , A. , Carr , C.A. , Johnston , D. , Wilson, M.A. , & Tonegawa , S. ( 2002 ). Requirement for hippocampal CA3 nmda receptors in associative memory recall . Science , 297 ( 5579 ), 211 - 218 . Nowacki , J. , Mazlan , S. , Osinga , H.M. , & Tsaneva-Atanasova , K. ( 2010 ). The role of large-conductance calcium-activated (bk) channels in shaping bursting oscillations of a somatotroph cell model . Physica D: Nonlinear Phenomena , 239 ( 9 ), 485 - 493 . Nowacki , J. , Osinga , H.M. , Brown , J.T., Randall , A.D. , & TsanevaAtanasova , K. ( 2011 ). A unified model of CA1/3 pyramidal cells: an investigation into excitability . Progress in Biophysics & Molecular Biology , 105 ( 1-2 ), 34 - 48 . O'Keefe , J. , & Dostrovsky , J. ( 1971 ). Hippocampus as a spatial map - preliminary evidence from unit activity in freely-moving rat . Brain Research , 34 ( 1 ), 171 - 175 . O'Keefe , J. , & Nadal , L. ( 1978 ). The hippocampus as a cognitive map . Oxford: Clarendon Press. O'Neill , J. , Pleydell-Bouverie , B. , Dupret , D. , & Csicsvari , J. ( 2010 ). Play it again: reactivation of waking experience and memory . Trends in Neurosciences , 33 ( 5 ), 220 - 229 . Pinsky , P.F. , & Rinzel , J. ( 1994 ). Intrinsic and network rhythmogenesis in a reduced traub model for CA3 neurons . Journal of computational neuroscience , 1 ( 1-2 ), 39 - 60 . Riedel , G. , Micheau , J. , Lam , A.G.M. , Roloff , E.V. , Martin , S.J. , Bridge , H., de Hoz, L. , Poeschel , B. , McCulloch , J. , & Morris , R.G.M. ( 1999 ). Reversible neural inactivation reveals hippocampal participation in several memory processes . Nature Neuroscience , 2 ( 10 ), 898 - 905 . Rinzel , J. ( 1987 ). A formal classification of bursting mechanisms in excitable systems , In Teramoto, E., & Yamaguti , M. (Eds.) Lecture Notes in Biomathematics, Vol. 71 . Mathematical Topics in Population Biology, Morphogenesis and Neurosciences; International Symposium, Kyoto, Japan, November 10-15 , 1985 . Ix+348p. Springer-Verlag: New York, New York, USA; Berlin, West Germany. Illus. Paper pp 267 - 281 . Rolls , E.T. , & Kesner , R.P. ( 2006 ). A computational theory of hippocampal function, and empirical tests of the theory . Progress in Neurobiology, 79 ( 1 ), 1 - 48 . Sah , P. , Gibb , A.J. , & Gage , P.W. ( 1988 ). The sodium current underlying action-potentials in Guinea-pig hippocampal CA1 neurons . Journal of General Physiology , 91 ( 3 ), 373 - 398 . Schwartzkroin , P.A. , & Stafstrom , C.E. ( 1980 ). Effects of egta on the calcium-activated afterhyperpolarization in hippocampal CA3- pyramidal cells . Science , 210 ( 4474 ), 1125 - 1126 . Spruston , N. , & McBain , C. ( 2007 ). Structural and functional properties of hippocampal neurons , pp. 133 - 201 . New York: Oxford University Press. Squire , L.R. ( 1992 ). Memory and the hippocampus - a synthesis from findings with rats, monkeys, and humans . Psychological Review , 99 ( 2 ), 195 - 231 . Stern , J.V. , Osinga , H.M. , LeBeau , A. , & Sherman , A. ( 2008 ). Resetting behavior in a model of bursting in secretory pituitary cells: Distinguishing plateaus from pseudo-plateaus . Bulletin of Mathematical Biology , 70 ( 1 ), 68 - 88 . Strogatz , S.H. ( 2001 ). Nonlinear Dynamics and Chaos: With applications to Physics, Biology, Chemistry and Engineering Perseus Books. Szalai , R. , Tsaneva-Atanasova , K. , Homer , M.E. , Champneys , A.R. , Kennedy , H.J. , & Cooper , N.P. ( 2011 ). Nonlinear models of development, amplification and compression in the mamMalian cochlea . Philosophical Transactions of the Royal Society aMathematical Physical and Engineering Sciences , 369 ( 1954 ), 4183 - 4204 . Tabak , J. , Toporikova , N. , Freeman , M.E. , & Bertram , R. ( 2007 ). Low dose of dopamine may stimulate prolactin secretion by increasing fast potassium currents . Journal of Computational Neuroscience , 22 ( 2 ), 211 - 222 . Taxidis , J. , Coombes , S. , Mason , R. , & Owen , M.R. ( 2012 ). Modeling sharp wave-ripple complexes through a CA3-CA1 network model with chemical synapses . Hippocampus , 22 ( 5 ), 995 - 1017 . Tiesinga , P.H.E. , Fellous , J.M. , Jose, J.V. , & Sejnowski, T.J. ( 2001 ). Computational model of carbachol-induced delta, theta, and gamma oscillations in the hippocampus . Hippocampus , 11 ( 3 ), 251 - 274 . Traub , R.D. , & Wong , R.K.S. ( 1982 ). Cellular mechanism of neuronal synchronization in epilepsy . Science , 216 ( 4547 ), 745 - 747 . Traub , R.D. , Wong , R.K.S. , Miles , R. , & Michelson , H. ( 1991 ). A model of a CA3 hippocampal pyramidal neuron incorporating voltage-clamp data on intrinsic conductances . Journal of Neurophysiology , 66 ( 2 ), 635 - 650 . Tsaneva-Atanasova , K. , Sherman , A., van Goor , F. , & Stojilkovic , S.S. ( 2007 ). Mechanism of spontaneous and receptor-controlled electrical activity in pituitary somatotrophs: Experiments and theory . Journal of Neurophysiology , 98 ( 1 ), 131 - 144 . Vervaeke , K. , Gu , N. , Agdestein , C. , Hu , H. , & Storm , J.F. ( 2006 ). Kv7/kcnq/m-channels in rat glutamatergic hippocampal axons and their role in regulation of excitability and transmitter release . Journal of Physiology-London , 576 ( 1 ), 235 - 256 . Wilson , M.A. , & McNaughton , B.L. ( 1994 ). Reactivation of hippocampal ensemble memories during sleep . Science , 265 ( 5172 ), 676 - 679 . Wong , R.K.S. , & Prince , D.A. ( 1978 ). Participation of calcium spikes during intrinsic burst firing in hippocampal neurons . Brain Research , 159 ( 2 ), 385 - 390 . Wong , R.K.S. , & Prince , D.A. ( 1981 ). After-potential generation in hippocampal pyramidal cells . Journal of Neurophysiology , 45 ( 1 ), 86 - 97 . Xu , J. , & Clancy , C.E. ( 2008 ). Ionic mechanisms of endogenous bursting in CA3 hippocampal pyramidal neurons: A model study Plos One 3(4).


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10827-016-0606-8.pdf

Laura A. Atherton, Luke Y. Prince, Krasimira Tsaneva-Atanasova. Bifurcation analysis of a two-compartment hippocampal pyramidal cell model, Journal of Computational Neuroscience, 2016, 91-106, DOI: 10.1007/s10827-016-0606-8