Numerical simulations were performed with the NEURON 7.3 simulator [33; RRID: nif-0000-00081], using its backward Euler integration method and 25 μs time steps. Fixed time step integration was used throughout. Preliminary, feasibility testing with a variable time step method (CVODE in NEURON, [33]), even with a high error tolerance, did not produce significant speeding of simulations (data not shown).

Reduction algorithm

Our reduction algorithm is a simple extension to an existing reduction algorithm [34,35]. Like the original formulation, it collapses the elaborate dendritic arbour into a smaller number of compartments whilst conserving axial resistance (Ra ). However, unlike prior formulations it collapses the dendritic tree to just a single compartment. Furthermore, it has an additional step whereby the length of the dendrite compartment is increased, with a compensatory change in radius to keep the compartment’s volume constant. This added step is required in order for the reduced Purkinje neuron model to replicate the behavior of the full Purkinje neuron model and express the trimodal activity pattern. The optimal dendrite length is found by manual tuning.

Recently an alternative extension to the algorithm of [34,35] has been published, in which the collapsing function incorporates a Strahler analysis of the dendritic tree [36,37]. However, this method typically produces reduced models of the cerebellar Purkinje cell that have tens of compartments and run, on average, ~20 times faster than the full model [37]. Our reduction method produces a faithful 2 compartment model that runs >400 times faster.

Our algorithm merges successive dendritic branches into an equivalent cylinder, preserving the axial resistance of the original branches. The radius (R) of the equivalent cylinder is given by:

$$ R\kern0.5em =\kern0.5em \sqrt{\underset{i}{\Sigma \kern0.5em {r}_i^2}} $$ (1)

where r i are the radii of the collapsed branches for the equivalent cylinder. The length (l) of the equivalent cylinder is taken as an average of the lengths of the collapsed branches (l i ) for that cylinder, weighted by their respective radii (r i ):

$$ l\kern0.5em =\kern0.5em \frac{\underset{i}{\Sigma}{l}_i{r}_i}{\underset{i}{\Sigma}{r}_i} $$ (2)

The volume of the cylinder (V) is:

$$ V\kern0.5em =\kern0.5em \pi \kern0.5em \cdot \kern0.5em {R}^2\kern0.5em \cdot \kern0.5em l $$ (3)

We can then change the length of the cylinder (l) to a desired value and can compensate for this change, to keep the cylinder volume (V) constant, by adjusting the radius (R) value accordingly:

$$ R\kern0.5em =\kern0.5em \sqrt{\frac{V}{\pi \kern0.5em \cdot \kern0.5em l}} $$ (4)

This final manipulation is to confront the fact that the collapsing function does not conserve dendritic length. It gives the equivalent cylinder a short Length of 120.9 μm (and Diameter of 6.74 μm). This is a problem because Purkinje cell functioning requires a degree of uncoupling (distance) between somatic and dendritic events. With this length of dendritic compartment, the model cannot replicate the Purkinje cell’s trimodal firing pattern (data not shown). To address this, we set the equivalent cylinder with a Length of 529.29 μm and adjusted its Radius, to 1.61 μm, in order to keep the volume constant (Eq. 4). This optimal length was found by manual tuning and was with a specific axial resistivity (Ra) of 35.4 Ωcm (the default NEURON setting). If Ra is different, then the optimal dendrite compartment dimensions are different. For example, with a Ra value of 250 Ωcm, the dimensions are instead: Length = 240 μm; Radius = 4.78 μm.

We have previously published a 41 compartment and a 5 compartment model of the cerebellar Purkinje neuron [11]. These models can replicate the Purkinje cell’s trimodal pattern of firing. However, the model of this current manuscript is superior to these because a model with just 2 compartments is the more elegant solution; not least because it runs faster. In addition, we describe the 2 compartment model here in immense detail to maximize its utility to other researchers. This level of detail is missing in [11].

With our collapsing method, the membrane surface area of the reduced model is less than the original because, although axial resistance is conserved, surface area (and hence the membrane resistance and capacitance) is not. This is compensated for by introducing in the equivalent cylinder, a dendritic correction factor (C d ), which rescales the current/pump/exchanger density values (g i ) and membrane capacitance (C m ) in the dendrites such that:

$$ {g}_i^{\hbox{'}}\kern0.5em =\kern0.5em {C}_d{g}_i $$ (5)

$$ {C}_i^{\hbox{'}}\kern0.5em =\kern0.5em {C}_d{C}_m $$ (6)

The dendritic correction factor C d is the ratio of the total surface area of the dendritic segments to their equivalent cylinder (42,310/6,874 = 6.16). The model’s somatic compartment is not introduced into the collapsing algorithm and so in the reduced model it has the same dimensions as in the full model, and the C d correction factor is not applied to any of its parameters.

Following the methodology of [35], our 2 compartment Purkinje cell model has C d applied to the model parameter, depth, which sets the depth of the sub-membrane shell that Ca2+ diffuses in within the model dendrite.

$$ dept{h}^{\hbox{'}}\kern0.5em =\kern0.5em {C}_d depth $$ (7)

The Q parameter is a setting in the model’s description of extracellular K+ dynamics and in the full, 1089 compartment model it is 0.06. This parameter isn’t set or constrained by any feature of the reduction algorithm used and is not constrained by any experimental literature. It is a free parameter. We manually tuned this parameter to confer a best fit between reduced and full model output, ultimately assigning it a value of 0.0119. This is distinct from the value of 0.37 obtained by a C d manipulation (0.06 * C d = 0.37, where C d is 6.16) but this is of no consequence, as there is no directive in the literature for applying C d to such a parameter. If a different value of Ra is assigned, different dendrite dimensions are to be used (as aforementioned), and the Q parameter needs to be re-tuned. For example, if Ra = 250 Ωcm, Q = 0.0103 gives a best fit.

In the 1089 compartment model of Forrest et al. [11], spines were accounted for (implicitly) by setting the 1003 spiny dendrite compartments with a higher specific membrane capacitance (C m ) value than the 84 smooth dendrite compartments (C m = 0.8 μF/cm2 for the smooth dendrites, C m = 1.5 μF/cm2 for the spiny dendrites). Our 2 compartment model makes no such account of spines; for its dendritic compartment: C m = 0.8 μF/cm2 (before the C d manipulation).

So, to summarize the reduction methodology: the full morphology is collapsed into just a soma compartment and a single dendrite compartment. This is done using the algorithm/NEURON code of [35], but modified to collapse the dendrites to just one compartment instead of three (this modified code is included with our model entry in ModelDB [12]). The length of this single dendrite compartment is then manually tuned to confer the best fit between reduced and full model output. For an inputted change in dendrite length, coding is used to automatically update other linked variables such that the user only has to tune this single parameter. For a change in dendrite length, the dendrite radius is automatically, accordingly updated to keep the dendrite volume constant. Furthermore, the surface area of the updated dendrite compartment is measured and a new, appropriate C d value is automatically calculated and applied. This coding form is in the NEURON code of the model, available in ModelDB [12]. So, to re-iterate, only a single variable needs be adjusted to tune the reduced model; such that it can best reproduce the behaviour of the full model. We adjusted it manually but with it being just a single parameter, it is likely tractable to automatic tuning methods. Thus described is the principal, crucial step to tuning. However, as aforementioned, after the optimal dendrite length has been found, an even closer fit between reduced and model output can then be reached by going on to tune the Q model parameter. This latter step is specific to this cerebellar Purkinje neuron model but we anticipate that the prior steps are more general and could be used to produce 2 compartment models for other neuron types.

Model equations

C m is the membrane capacitance (0.8 μF/cm2), I is the current, V is the membrane potential in mV, t is time, T is temperature (36° C), Ra is the specific axial resistivity (35.4 Ωcm and g max is the maximal conductance (“current density”).

The model currents have equations and kinetic parameters as described in their source literature with the exception of g max values, which have been modified. g max values, for the different currents, are shown in Table 1. These values are drawn from the detailed Purkinje model of [11] and were found by manual tuning. Below, as we describe each current, we also state the publication that it is drawn from.

Table 1 Maximal current conductances (g max ; mS/cm 2 ) for the Full, 1089 Compartment model and the Simplified, Two Compartment model Full size table

m, h and z are Hodgkin-Huxley “particles”/gates [38]; for example, for the m Hodgkin-Huxley gate:

$$ \frac{dm}{dt}=\frac{m_{\infty }-m}{\tau_m} $$ (8)

The voltage (and/or intracellular calcium) dependence of a Hodgkin-Huxley (H-H) current [38] can be expressed by stating, for each H-H gate (e.g. for the m gate), either [m ∞ , m ∞ ] OR [α m , β m ]. These entities are voltage (and/or intracellular calcium) dependent. The latter set can give the former set through the relations:

$$ {m}_{\infty }={\alpha}_m/\left({\alpha}_m+{\beta}_m\right) $$ (9)

$$ {\tau}_m=1/\left({\alpha}_m+{\beta}_m\right) $$ (10)

Soma

The soma is a cylinder (length = 22 μm, diameter = 22 μm). E K is the reversal potential for K+ (initiated at -88 mV), E Na is the reversal potential for Na+ (initiated at +70 mV), E Ca is the reversal potential for Ca2+ (initiated by the NEURON default value; +132 mV), E L is the reversal potential for the Leak current (-70 mV), E h is the reversal potential for the hyperpolarisation activated cation current (-30 mV), Intracellular Ca2+ concentration is initiated at the NEURON default of 5*10−5 mM; Extracellular Ca2+ concentration is initiated at the NEURON default of 2 mM. The somatic membrane voltage (V) is initiated at the NEURON default of -65 mV.

The soma has highly TEA sensitive (I K_fast ), moderately TEA sensitive (I K_mid ) and TEA insensitive (I K_slow ) voltage-gated K+ currents, a BK voltage-and-Ca2+-gated K+ current (I BK ), a resurgent Na+ current (I Na-R ), a P-type Ca2+ current (I CaP ), a hyperpolarization activated cation current (I H ), a leak current (I L ), a SK Ca2+-gated K+ current (I SK ) and an intracellular Ca2+ dynamics abstraction. The soma has two Na+/K+ pump descriptions (Is pump and is pump ) and a Na+/Ca2+ exchanger mechanism Is ex .

$$ {C}_m\cdot \frac{dV}{dt}=-\left({I}_{K\_ fast}+{I}_{K\_mid}+{I}_{K\_ slow}+{I}_{BK}+{I}_{CaP}+{I}_H+{I}_L+{I}_{NaR}+{I}_{SK}+{i}_{pump}^s+{I}_{pump}^s+{I}_{ex}^s+{I}_{transfer\_DS}\right) $$ (11)

Dendrite – Soma electrotonic current [39]

$$ {I}_{transfer\_DS}=\frac{\left({V}_D-{V}_S\right)}{R_{DS}} $$ (12)

$$ {R}_{DS}=\frac{Ra\cdot \left({L}_S/2\right)}{\pi \cdot {r_S}^2}+\frac{Ra\cdot \left({L}_D/2\right)}{\pi \cdot {r_D}^2} $$ (13)

V D is the membrane voltage at the centre of the dendrite compartment, V S is the membrane voltage at the centre of the soma compartment and R DS is the axial Resistance between the two. Ra is the specific axial Resistivity, L S and r S are the length and radius of the soma respectively; L D and r D are the length and radius of the dendrite respectively.

Highly TEA sensitive K+ current [40]

$$ {I}_{K\_ fast}={g}_{\max}\cdot {m}^3\cdot h\cdot \left(V-{E}_K\right) $$ (14)

$$ {m}_{\infty }=\frac{1}{ \exp \left(-\frac{V--24}{15.4}\right)} $$ (15)

$$ {\tau}_m=\left\{\begin{array}{l}0.000103+0.0149\ast \exp \left(0.035\ast V\right)......................\left[V<-35 mV\right]\\ {}0.000129+1/\left[ \exp \left(\frac{V+100.7}{12.9}\right)+ \exp \left(\frac{V-56}{-23.1}\right)\right]......\left[V\ge -35 mV\right]\end{array}\right. $$ (16)

$$ {h}_{\infty }=0.31+\frac{1-0.31}{ \exp \left(-\frac{V--5.8}{-11.2}\right)} $$ (17)

$$ {\tau}_h=\left\{\begin{array}{l}1.22*{10}^{-5}+0.012\ast \exp \left[-{\left(\frac{V+56.3}{49.6}\right)}^2\right]......................\left[V\le 0 mV\right]\\ {}0.0012+0.0023* \exp \left(-0.141*V\right)..............................\left[V>0 mV\right]\end{array}\right. $$ (18)

Moderately TEA sensitive K+ current [40]

$$ {I}_{K\_mid}={g}_{\max}\cdot {m}^4\cdot \left(V-{E}_K\right) $$ (19)

$$ {m}_{\infty }=\frac{1}{ \exp \left(-\frac{V--24}{20.4}\right)} $$ (20)

$$ {\tau}_m=\left\{\begin{array}{l}0.000688+1/\left[ \exp \left(\frac{V+64.2}{6.5}\right)+ \exp \left(\frac{V-141.5}{-34.8}\right)\right]..............\left[V<-20 mV\right]\\ {}0.00016+0.0008* \exp \left(-0.0267*V\right)..............................\left[V\ge -20 mV\right]\end{array}\right. $$ (21)

TEA insensitive K+ current [40]

$$ {I}_{K\_ slow}={g}_{\max}\cdot {m}^4\cdot \left(V-{E}_K\right) $$ (22)

$$ {m}_{\infty }=\frac{1}{ \exp \left(-\frac{V--16.5}{18.4}\right)} $$ (23)

$$ {\tau}_m=0.000796+1/\left[ \exp \left(\frac{V+73.2}{11.7}\right)+ \exp \left(\frac{V-306.7}{-74.2}\right)\right] $$ (24)

P-type Ca2+ current [40]

$$ {I}_{CaP}={g}_{\max }*m*ghk $$ (25)

Goldman-Hodgkin-Katz (ghk) equation:

$$ ghk=\left(4*{P}_{C{a}^{2+}}\right)*\frac{V\cdot {F}^2}{R\cdot T}*\frac{{\left[C{a}^{2+}\right]}_i-{\left[C{a}^{2+}\right]}_o* \exp \left(\frac{-2\cdot F\cdot V}{R\cdot T}\right)}{1- \exp \left(\frac{-2\cdot F\cdot V}{R\cdot T}\right)} $$ (26)

P Ca 2+ is 5*10−5cm/sec, [Ca2+] i = 100 nM, [Ca2+] o = 2 mM, T = 295 K, F is the Faraday constant and R is the gas constant. [Ca2+] i and [Ca2+] o are fixed constants, as seen by this equation – it does not access the changing value of [Ca2+] i as set by the intracellular Ca2+ equations (given later).

$$ {m}_{\infty }=\frac{1}{ \exp \left(-\frac{V--19}{5.5}\right)} $$ (27)

$$ {\tau}_m=\left\{\begin{array}{l}0.000264+0.128* \exp \left(0.103*V\right).....................\left[V\le -50 mV\right]\\ {}0.000191+0.00376* \exp \left[-{\left(\frac{V+11.9}{27.8}\right)}^2\right]...........\left[V>-50 mV\right]\end{array}\right. $$ (28)

Hyperpolarisation activated cation current [40]

$$ {I}_H={g}_{\max}\cdot m\cdot \left(V-{E}_h\right) $$ (29)

$$ {m}_{\infty }=\frac{1}{ \exp \left(-\frac{V--90.1}{-9.9}\right)} $$ (30)

$$ {\tau}_m=0.19+0.72* \exp \left[-{\left(\frac{V+81.5}{11.9}\right)}^2\right] $$ (31)

BK type K+ current [40]

$$ {I}_{BK}={g}_{\max}\cdot {m}^3\cdot {z}^2\cdot h\cdot \left(V-{E}_K\right) $$ (32)

$$ {m}_{\infty }=\frac{1}{ \exp \left(-\frac{V--28.9}{6.2}\right)} $$ (33)

$$ {h}_{\infty }=0.085+\frac{1-0.085}{ \exp \left(-\frac{V--32}{-5.8}\right)} $$ (34)

$$ {\tau}_m=0.000505+1/\left[ \exp \left(\frac{V+86.4}{10.1}\right)+ \exp \left(\frac{V-33.3}{-10}\right)\right] $$ (35)

$$ {\tau}_h=0.0019+1/\left[ \exp \left(\frac{V+48.5}{5.2}\right)+ \exp \left(\frac{V-54.2}{-12.9}\right)\right] $$ (36)

$$ {z}_{\infty }=\frac{1}{1+\frac{0.001}{\left[C{a}^{2+}\right]}} $$ (37)

$$ {\tau}_z=1 $$ (38)

Leak current [40]

$$ {I}_L={g}_{\max }*\left(V-{E}_L\right) $$ (39)

Intracellular Ca2+ concentration [40]

[Ca2+] is calculated for the intracellular space within 100 nm of the membrane. [Ca2+] changes as I Ca 2+ (negative by convention; inward currents are negative) brings Ca2+ into this space and as Ca2+ leaves by diffusion to the bulk cytoplasm. The diffusion rate constant,β, is set to 1/msec.

$$ \frac{d\left[C{a}^{2+}\right]}{dt}=\beta *\left[C{a}^{2+}\right] $$ (40)

[Ca2+] at time step, t:

$$ {\left[C{a}^{2+}\right]}_t={\left[C{a}^{2+}\right]}_{t-1}+\Delta t*\left(\frac{-(100)*{I}_{C{a}^{2+}}}{\left(2\cdot F\right)*\left( depth\cdot Area\right)}-\beta *{\left[C{a}^{2+}\right]}_{t-1}\right) $$ (41)

F is the Faraday constant, depth = 0.1 μm and membrane surface Area = 1,521μm 2. [Ca2+] was constrained to not fall below 100 nM by coding of the form:

$$ if\left(\left[C{a}^{2+}\right]<100\right)\left\{\left[C{a}^{2+}\right]=100\right\} $$ (42)

Resurgent Na+ current [40,41]

$$ {I}_{NaR}=\kern0.5em {g}_{\max}\kern0.5em *\kern0.5em O\kern0.5em *\kern0.5em \left(V\kern0.5em \hbox{-} \kern0.5em {E}_{Na}\right)\ O\kern0.5em \mathrm{is}\ \mathrm{the}\ \mathrm{occupancy}\ \mathrm{of}\ \mathrm{the}\ \mathrm{Open}\ \mathrm{state}. $$ (43)

This current is described by a Markov scheme, shown in Figure 1. The rate constants, labelled in Figure 1, are (ms−1):

Figure 1 The Resurgent Na+ current is described by a Markov scheme [40,41]. [C1 to C5] denote sequential Closed states; O denotes the Open state. [I1 to I6] denote Inactivated states. OB denotes the state entered by a second mechanism of inactivation, which is hypothesized to be equivalent to Open Channel Block. The rate constants between states are given in Eq. [44], Eq. [45], Eq. [46], Eq. [47], Eq. [48] and Eq. [49]. Full size image

$$ \alpha =150* \exp \left(\frac{V}{20}\right) $$ (44)

$$ \beta =3* \exp \left(\frac{2\cdot V}{20}\right) $$ (45)

$$ \gamma =150;\kern0.5em \delta =40;\kern0.5em Con=0.005;\kern0.5em Coff=0.5;\kern0.5em Oon=0.75;\kern0.5em Ooff=0.005 $$

$$ a={\left(\frac{Oon}{Con}\right)}^{1/4} $$ (46)

$$ b={\left(\frac{Ooff}{Coff}\right)}^{1/4} $$ (47)

$$ \varepsilon =1.75 $$ (48)

$$ \zeta =0.03* \exp \left(\frac{2\cdot V}{25}\right) $$ (49)

SK type K+ current [42]

$$ {I}_{SK}={g}_{\max }*z*\left(V-{E}_K\right) $$ (50)

$$ z=\frac{1}{1+{\left(0.00019/\left[C{a}^{2+}\right]\right)}^4} $$ (51)

Na+/K+ pump [11]

The somatic Na+/K+ pump (density = ds pump , 1 mA/cm2) transports 3 Na+ out (is pump_Na ) for every 2 K+ in (is pump_K ). It has a fixed voltage (V) dependency and an exponential relation to intracellular Na+ concentration ([Na+] i ). The affinity constant for Na+, K Na = 40 mM.

$$ \left\{\begin{array}{l}{\mathrm{i}}_{\mathrm{pump}}^{\mathrm{s}} = {\mathrm{d}}_{\mathrm{pump}}^{\mathrm{s}}\left(\mathrm{V} + 75\right)\ /\ \left[\left(\mathrm{V} + 80\right)\left(1+ \exp \left(\frac{{\mathrm{K}}_{\mathrm{Na}}\hbox{-} {\left[{\mathrm{Na}}^{+}\right]}_{\mathrm{i}}}{1}\right)\right)\right]\\ {}{\mathrm{i}}_{\mathrm{pump}\_\mathrm{N}\mathrm{a}}^{\mathrm{s}} = 3{\mathrm{i}}_{\mathrm{pump}}^{\mathrm{s}}\\ {}{\mathrm{i}}_{\mathrm{pump}\_\mathrm{K}}^{\mathrm{s}} = \hbox{-} 2\ {\mathrm{i}}_{\mathrm{pump}}^{\mathrm{s}}\end{array}\right. $$ (52)

Na+/Ca2+ exchanger current and an electrically counterbalancing Na+/K+ pump current [11]

The soma has a simple Na+/Ca2+ exchanger mechanism (Eq. 53) and a simple Na+/K+ pump mechanism (Eq. 54). The Na+/Ca2+ exchanger current (Is ex_net ) is net depolarizing (-1), inwardly passing 3 singly positive Na+ ions (3*[+1]) for the extrusion of every doubly positive Ca2+ ion (1*[+2]). By contrast, the Na+/K+ pump current (Is pump_net ) is net hyperpolarizing (+1) in its transport of 3 Na+ out (3*[+1]) for every 2 K+ in (2*[+1]). The exchanger density is gs ex ; the pump density is gs pump .

$$ \left\{\begin{array}{l}{\mathrm{I}}_{\mathrm{ex}\_\mathrm{N}\mathrm{a}}^s=-3\cdot \left[+1\right]\cdot {g}_{ex}^s\\ {}{\mathrm{I}}_{\mathrm{ex}\_\mathrm{C}\mathrm{a}}^{\mathrm{s}}=1\cdot \left[+2\right]\cdot {g}_{ex}^s\\ {}{I}_{ex\_net}^s=\left({I}_{ex\_Ca}^s-{I}_{ex\_Na}^s\right)\Rightarrow {g}_{ex}^s\cdot \left[-1\right]\Rightarrow -{g}_{ex}^s\end{array}\right. $$ (53)

$$ \left\{\begin{array}{l}{\mathrm{I}}_{\mathrm{pump}\_\mathrm{N}\mathrm{a}}^s=3\cdot \left[+1\right]\cdot {g}_{pump}^s\\ {}{\mathrm{I}}_{\mathrm{pump}\_\mathrm{K}}^{\mathrm{s}}=-2\cdot \left[+1\right]\cdot {g}_{pump}^s\\ {}{I}_{pump\_net}^s=\left({I}_{pump\_Na}^s-{I}_{pump\_K}^s\right)\Rightarrow {g}_{pump}^s\cdot \left[+1\right]\Rightarrow +{g}_{pump}^s\end{array}\right. $$ (54)

The Na+/K+ pump current of Eq. 54 largely, but incompletely, counterbalances the Na+/Ca2+ exchanger current of Eq. 53 – there is a slight mismatch [gs ex = 0.511 mA/cm2, gs pump = 0.5 mA/cm2] (Eq. 55) which permits a small net influx of Na+ ions.

$$ \left[{\mathrm{g}}_{\mathrm{ex}}^{\mathrm{s}} = {\mathrm{g}}_{\mathrm{pump}}^{\mathrm{s}} + 0.011\right]\kern0.5em \Rightarrow \kern0.5em \left[\hbox{-} {\mathrm{I}}_{\mathrm{ex}\_\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{s}}\kern0.5em \approx \kern0.5em +{\mathrm{I}}_{\mathrm{pump}\_\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{s}}\right] $$ (55)

Intracellular Na+ concentration [11]

[Na+] i is initiated at 10 mM and then changes in time,

$$ \frac{\partial {\left[N{a}^{+}\right]}_i}{\partial t}=\frac{I_{Na\_net}}{\left[d\cdot F\cdot (10000)\right]/4}+\frac{D{\partial}^2{\left[N{a}^{+}\right]}_i}{{\left(\partial x\right)}^2} $$ (56)

$$ {\left({I}_{Na\_net}\right)}_{\left[t\right]}={\left({I}_{Na\_ in}-{I}_{Na\_ out}\right)}_{\left[t-\tau \right]},\tau = 5s $$ (57)

$$ {I}_{Na\_ in}={I}_{Na-R}+{I}_{ex\_Na}^s $$ (58)

$$ {I}_{Na\_ out}\kern0.5em =\kern0.5em {i^s}_{pump\_Na}\kern0.5em +\kern0.5em {I_{pump}^s}_{\_Na} $$ (59)

F is the Faraday constant, d is the somatic diameter and D is the diffusivity constant (0.6 μm2/ms). The second term on the Right Hand Side (RHS) of Eq. 56 accounts for longitudinal diffusion of Na+ out of the soma compartment, along the longitudinal distance (x). The effects of this term are fairly negligible and it could be dropped to quicken simulation speeds. I Na_net is the difference between Na+ current flowing into the soma (I Na_in ; through I Na-R and I s ex_Na ) and Na+ current pumped out of the soma by the Na+/K+ pump (I Na_out ), lagged by parameter τ = 5 s. Intracellular Na+ stimulates the Na+/K+ pump and this lag τ accounts for the duration of sodium’s diffusion from channels to pumps. is pump_Na and Is pump_Na are Na+ currents produced by the pumping action of the Na+/K+ pump at the soma, set by Eq. 52 and Eq. 54 respectively. “Catch coding” is applied:

$$ \mathrm{if}\ \left({\left[\mathrm{N}{\mathrm{a}}^{+}\right]}_{\mathrm{i}} < 10\ \right)\ \left[{\left[\mathrm{N}{\mathrm{a}}^{+}\right]}_{\mathrm{i}} = 10\ \right] $$ (60)

$$ \mathrm{if}\ \left({\mathrm{E}}_{\mathrm{Na}} < 70\ \right)\ \left[{\mathrm{E}}_{\mathrm{Na}} = 70\right] $$ (61)

Dendrite

The dendrite is a cylinder (length = 529.29 μm, diameter = 3.22 μm). E K is the reversal potential for K+ (initiated by the NEURON default value; -77 mV), E Na is the reversal potential for Na+ (initiated by the NEURON default value; +50 mV), E Ca is the reversal potential for Ca2+ (initiated by the NEURON default value; +132 mV), E L is the reversal potential for the Leak current (-80 mV), E h is the reversal potential for the hyperpolarisation activated cation current (-32.9 mV). Intracellular Ca2+ concentration is initiated at 4*10−5 mM; Extracellular Ca2+ concentration is initiated at 2.4 mM. Extracellular Na+ concentration is set by the NEURON default value; +140 mM.

All the mechanisms in the dendrite are distinct from those in the soma. The dendrite has hyperpolarization activated cation current (I H ); T-type (I CaT ), Class-E (I CaE ) and P-type (I CaP ) voltage-gated Ca2+ currents; a leak current (I L ); A-type (I KA ), D-type (I KD ), M-type (I KM ), Delayed Rectifier (I DR ) and Kv1.2 (I Kv1.2 ) voltage-gated K+ currents; BK (I BK ) and K2 (I K2 ) type voltage-and-Ca2+-gated K+ currents and an intracellular Ca2+ dynamics abstraction. The dendrite has two Na+/K+ pump descriptions (Id pump and id pump ) and a Na+/Ca2+ exchanger mechanism Id ex .

$$ {C}_m\cdot \frac{dV}{dt}=-\left(\begin{array}{l}{I}_{CaT}+{I}_{CaE}+{I}_{CaP}+{I}_H+{I}_{Kv1.2}+{I}_{KA}+{I}_{KM}+{I}_{KD}+{I}_{DR}+{I}_{BK}+{I}_{K2}+{I}_L+{i}_{pump}^d\\ {}+{I}_{pump}^d+{I}_{ex}^d+{I}_{transfer\_SD}\end{array}\right) $$ (62)

Soma - dendrite electrotonic current [39]

$$ {I}_{transfer\_SD}=\frac{\left({V}_S-{V}_D\right)}{R_{SD}} $$ (63)

$$ {R}_{SD}=\frac{Ra\cdot \left({L}_S/2\right)}{\pi \cdot {r_S}^2}+\frac{Ra\cdot \left({L}_D/2\right)}{\pi \cdot {r_D}^2} $$ (64)

V S is the membrane voltage at the centre of the soma compartment, V D is the membrane voltage at the centre of the dendrite compartment and R SD is the axial Resistance between the two. Ra is the specific axial Resistivity, L S and r S are the length and radius of the soma respectively; L D and r D are the length and radius of the dendrite respectively.

T-type Ca2+ current [43]

$$ {I}_{CaT}={g}_{\max}\cdot m\cdot h\cdot \left(V-{E}_{Ca}\right)\kern1.75em ;\kern1.25em {\mathrm{E}}_{\mathrm{Ca}}\mathrm{is}\ \mathrm{f}\mathrm{ixed}\ \mathrm{at} + 135\ \mathrm{mV}\ \mathrm{f}\mathrm{o}\mathrm{r}\ \mathrm{this}\ \mathrm{current}. $$ (65)

$$ {\alpha}_m=\frac{2.6}{1+ \exp \left(\frac{V+21}{-8}\right)} $$ (66)

$$ {\beta}_m=\frac{0.18}{1+ \exp \left(\frac{V+40}{4}\right)} $$ (67)

$$ {\alpha}_h=\frac{0.0025}{1+ \exp \left(\frac{V+40}{8}\right)} $$ (68)

$$ {\beta}_h=\frac{0.19}{1+ \exp \left(\frac{V+50}{-10}\right)} $$ (69)

\( mt={3}^{\frac{T-37}{10}} \); T is temperature in degrees centigrade (36).

$$ {\tau}_m=\frac{1}{\left({\alpha}_m+{\beta}_m\right)\cdot mt} $$ (70)

$$ {\tau}_h=\frac{1}{\left({\alpha}_h+{\beta}_h\right)\cdot mt} $$ (71)

E-type Ca2+ current [43]

$$ {I}_{CaE}={g}_{\max}\cdot m\cdot h\cdot \left(V-{E}_{Ca}\right);{\mathrm{E}}_{\mathrm{Ca}}\mathrm{is}\ \mathrm{f}\mathrm{ixed}\ \mathrm{at} + 135\ \mathrm{mV}\ \mathrm{f}\mathrm{o}\mathrm{r}\ \mathrm{this}\ \mathrm{current}. $$ (72)

$$ {\alpha}_m=\frac{2.6}{1+ \exp \left[\left(V+7\right)/-8\right]} $$ (73)

$$ {\beta}_m=\frac{0.18}{1+ \exp \left[\left(V+26\right)/4\right]} $$ (74)

$$ {\alpha}_h=\frac{0.0025}{1+ \exp \left[\left(V+32\right)/8\right]} $$ (75)

$$ {\beta}_h=\frac{0.19}{1+ \exp \left[\left(V+42\right)/-10\right]} $$ (76)

\( mt={3}^{\frac{T-37}{10}} \); T is temperature in degrees centigrade (36).

$$ {m}_{\exp }=1- \exp \left(\frac{\left[-dt*mt\right]\cdot \left[{\alpha}_m+{\beta}_m\right]}{4}\right) $$ (77)

$$ {h}_{\exp }=1- \exp \left(\frac{\left[-dt*mt\right]\cdot \left[{\alpha}_h+{\beta}_h\right]}{10}\right) $$ (78)

P-type Ca2+ current [43]

$$ {I}_{CaP}={g}_{\max}\cdot m\cdot \left(V-{E}_{Ca}\right);{\mathrm{E}}_{\mathrm{Ca}}\mathrm{is}\ \mathrm{f}\mathrm{ixed}\ \mathrm{at} + 135\ \mathrm{mV}\ \mathrm{f}\mathrm{o}\mathrm{r}\ \mathrm{this}\ \mathrm{current}. $$ (79)

$$ {\alpha}_m=\frac{8.5}{1+ \exp \left(\left[V+-8\right]/-12.5\right)} $$ (80)

$$ {\beta}_m=\frac{35}{1+ \exp \left(\left[V+74\right]/14.5\right)} $$ (81)

\( mt={3}^{\frac{T-37}{10}} \); T is temperature in degrees centigrade (36).

$$ {\tau}_m=\frac{1}{\left({\alpha}_m+{\beta}_m\right)\cdot mt} $$ (82)

Hyperpolarisation activated cation current [44]

$$ {I}_h={g}_{\max}\cdot m\cdot \left(V-{E}_h\right);\ {\mathrm{E}}_{\mathrm{h}} = \hbox{-} 32.9\ \mathrm{mV} $$ (83)

$$ {\tau}_m=\frac{1}{ \exp \left(-17.9-0.116\cdot V\right)+ \exp \left(-1.84+0.09\cdot V\right)}+100 $$ (84)

$$ {m}_{\infty }=\frac{1}{1+ \exp \left[\left(V+84.1\right)/10.2\right]} $$ (85)

$$ {m}_{\exp }=1- \exp \left(\frac{-dt}{\tau_m}\right) $$ (86)

Kv1.2 K+ current [45]

$$ {I}_{Kv1.2}={g}_{\max}\cdot {m}^4\cdot \left(V-{E}_K\right) $$ (87)

$$ {\alpha}_m=0.12899* \exp \left(\frac{-\left(V+45\right)}{-33.90877}\right) $$ (88)

$$ {\beta}_m=0.12899* \exp \left(\frac{-\left(V+45\right)}{12.42101}\right) $$ (89)

\( mt={3}^{\frac{T-22}{10}} \); T is temperature in degrees centigrade (36).

$$ {\tau}_m=\frac{1}{\left({\alpha}_m+{\beta}_m\right)\cdot mt} $$ (90)

A-type K+ current [43]

$$ {I}_{KA}={g}_{\max}\cdot {m}^4\cdot h\cdot \left(V-{E}_K\right) $$ (91)

$$ {\alpha}_m=\frac{1.4}{1+ \exp \left(\left[V+27\right]/-12\right)} $$ (92)

$$ {\beta}_m=\frac{0.49}{1+ \exp \left(\right[V+30/4\Big)} $$ (93)

$$ {\alpha}_h=\frac{0.00175}{1+ \exp \left(\right[V+50/8\Big)} $$ (94)

$$ {\beta}_h=\frac{0.49}{1+ \exp \left(\right[V+13/-10\Big)} $$ (95)

\( mt={3}^{\frac{T-37}{10}} \); T is temperature in degrees centigrade (36).

$$ {\tau}_m=\frac{1}{\left({\alpha}_m+{\beta}_m\right)\cdot mt} $$ (96)

$$ {\tau}_h=\frac{1}{\left({\alpha}_h+{\beta}_h\right)\cdot mt} $$ (97)

M-type K+ current [43]

$$ {I}_{KM}={g}_{\max}\cdot m\cdot \left(V-{E}_K\right) $$ (98)

\( ft={2.3}^{\frac{T-36}{10}} \); T is temperature in degrees centigrade (36).

$$ {\tau}_m=\frac{1000/ ft}{3.3\cdot \left({e}^{+\left(V+35\right)/40}+{e}^{-\left(V+35\right)/20}\right)} $$ (99)

$$ {m}_{\infty }=\frac{1}{1+{e}^{-\left(V+35\right)/10}} $$ (100)

D-type K+ current [43]

$$ {I}_{KD}={g}_{\max}\cdot m\cdot h\cdot \left(V-{E}_K\right) $$ (101)

$$ {\alpha}_m=\frac{8.5}{1+ \exp \left(\left[V+17\right]/-12.5\right)} $$ (102)

$$ {\beta}_m=\frac{35}{1+ \exp \left(\left[V+99\right]/14.5\right)} $$ (103)

$$ {\alpha}_h=\frac{0.0015}{1+ \exp \left(\left[V+89\right]/8\right)} $$ (104)

$$ {\beta}_h=\frac{0.0055}{1+ \exp \left(\left[V+83\right]/-8\right)} $$ (105)

$$ {m}_{\infty }={\alpha}_m/\left({\alpha}_m+{\beta}_m\right) $$ (106)

$$ {h}_{\infty }={\alpha}_h/\left({\alpha}_h+{\beta}_h\right) $$ (107)

\( mt={3}^{\frac{T-37}{10}} \); T is temperature in degrees centigrade (36).

$$ {m}_{\exp }=1- \exp \left(\frac{\left[-dt*mt\right]\cdot \left[{\alpha}_m+{\beta}_m\right]}{10}\right) $$ (108)

$$ {h}_{\exp }=1- \exp \left(\left[-dt*mt\right]\cdot \left[{\alpha}_h+{\beta}_h\right]\cdot 1.6\right) $$ (109)

These equations are based upon the NEURON code of [43], as opposed to the equations in their associated journal paper; there is some minor discrepancy between the two. This code was sourced from ModelDB [12].

Delayed Rectifier type K+ current [43]

$$ {I}_{DR}={g}_{\max}\cdot {m}^4\cdot \left(V-{E}_K\right) $$ (110)

$$ {\alpha}_m=0.1\cdot vtrap $$ (111)

$$ catch= fabs\left(\frac{-\left(V+55\right)}{10}\right) $$ (112)

Where fabs(x) returns the absolute value of a floating point number; the absolute value of its argument (|x|).

$$ vtrap=\left\{\begin{array}{l}10\cdot \left(1-\frac{-\left(V+55\right)}{10}/2\right)....................\left[ catch<1{e}^{-6}\right]\\ {}\frac{-\left(V+55\right)}{ \exp \left(\frac{-\left(V+55\right)}{10}\right)-1}.......................\left[ catch\ge 1{e}^{-6}\right]\end{array}\right. $$ (113)

$$ {\beta}_m=0.125\cdot \exp \left(\frac{-\left(V+65\right)}{80}\right) $$ (114)

\( mt={3}^{\frac{T-37}{10}} \); T is temperature in degrees centigrade (36).

$$ {m}_{\exp }=1- \exp \left(-dt\cdot mt\cdot \left[{\alpha}_m+{\beta}_m\right]\right) $$ (115)

BK type K+ current [43]

$$ {I}_{BK}={g}_{\max}\cdot m\cdot {z}^2\cdot \left(V-{E}_K\right) $$ (116)

$$ {\alpha}_m=7.5 $$ (117)

$$ {\beta}_m=\frac{0.11}{ \exp \left(\left[V+-35\right]/14.9\right)} $$ (118)

$$ {m}_{\exp }=1- \exp \left(-dt\cdot \left[{\alpha}_m+{\beta}_m\right]\right) $$ (119)

$$ {\alpha}_z=1 $$ (120)

$$ {\beta}_z=\frac{400}{\left[C{a}^{2+}\right]*1000} $$ (121)

$$ {\tau}_z=10 $$ (122)

$$ {z}_{\exp }=1- \exp \left(\frac{-dt}{\tau_z}\right) $$ (123)

These equations are based upon the NEURON code of [43], as opposed to the equations in their associated journal paper; there is some minor discrepancy between the two. This code was sourced from ModelDB [12].

K2 type K+ current [43]

$$ {I}_{K2}={g}_{\max}\cdot m\cdot {z}^2\cdot \left(V-{E}_K\right) $$ (124)

$$ {\alpha}_m=25 $$ (125)

$$ {\beta}_m=\frac{0.075}{ \exp \left(\left[V+5\right]/10\right)} $$ (126)

$$ {m}_{\exp }=1- \exp \left(-dt\cdot \left[{\alpha}_m+{\beta}_m\right]\right) $$ (127)

$$ {\alpha}_z=1 $$ (128)

$$ {\beta}_z=\frac{20}{\left[C{a}^{2+}\right]*1000} $$ (129)

$$ {\tau}_z=10 $$ (130)

$$ {z}_{\exp }=1- \exp \left(\frac{-dt}{\tau_z}\right) $$ (132)

These equations are based upon the NEURON code of [43], as opposed to the equations in their associated journal paper; there is some minor discrepancy between the two. This code was sourced from ModelDB [12].

Leak current [43]

$$ {I}_L=\kern0.5em {g}_{\max }*\kern0.5em \left(V-{E}_L\right);{\mathrm{E}}_{\mathrm{L}}\mathrm{is}\ \hbox{-} 80\ \mathrm{mV}\ \mathrm{f}\mathrm{o}\mathrm{r}\ \mathrm{this}\ \mathrm{current}\ \mathrm{in}\ \mathrm{the}\ \mathrm{dendrite}. $$ (133)

Intracellular Ca2+ dynamics [43]

$$ \frac{d{\left[C{a}^{2+}\right]}_i}{dt}\kern0.5em =\kern0.5em chan\kern0.5em +\kern0.5em \left(\frac{-\kern0.5em kt\kern0.5em *\kern0.5em {\left[C{a}^{2+}\right]}_i}{{\left[C{a}^{2+}\right]}_i\kern0.5em +\kern0.5em kd}\right)\kern0.5em +\kern0.5em \left(\frac{y\kern0.5em -\kern0.5em {\left[C{a}^{2+}\right]}_i}{ta{u}_r}\right) $$ (134)

$$ chan=\left(\frac{-(10000)*{I}_{C{a}^{2+}}}{2*F* depth}\right) $$ (135)

$$ if\left( chan<0\right)\left\{ chan=0\right\} $$ (136)

where [Ca2+] i is the intracellular Ca2+ concentration in a supra-membrane shell of depth = 0.1 μm, F is the Faraday constant, \( {I}_{C{a}^{2+}} \) is the Ca2+ membrane current (negative by convention; inward currents are negative), kt = 4*10−5 mM/ms, kd = 4*10−5 mM, tau r = 2 ms and y = 4*10−5 mM.

Na+/K+ pump [11]

The dendritic Na+/K+ pump (density = dd pump , 0.001 mA/cm2) has a 3Na+:2 K+ stoichiometry, no voltage dependency and a hyperbolic relation to extracellular K+ concentration ([K+] o .

$$ \left\{\begin{array}{l}{\mathrm{i}}_{\mathrm{pump}}^d = {\mathrm{d}}_{\mathrm{pump}}^d/\ \left(1+\left(2.245/\ {\left[{\mathrm{K}}^{+}\right]}_{\mathrm{o}}\right)\right)\\ {}{\mathrm{i}}_{\mathrm{pump}\_\mathrm{N}\mathrm{a}}^d = 3{\mathrm{i}}_{\mathrm{pump}}^d\\ {}{\mathrm{i}}_{\mathrm{pump}\_\mathrm{K}}^d = \hbox{-} 2\ {\mathrm{i}}_{\mathrm{pump}}^d\end{array}\right. $$ (137)

Na+/Ca2+ exchanger current and an electrically counterbalancing Na+/K+ pump current [11]

The dendrite has a simple Na+/Ca2+ exchanger mechanism (Eq. 138) and a simple Na+/K+ pump mechanism (Eq. 139). The Na+/Ca2+ exchanger current is net depolarizing (-1), inwardly passing 3 singly positive Na+ ions (3*[+1]) for the extrusion of every doubly positive Ca2+ ion (1*[+2]). By contrast, the Na+/K+ pump current is net hyperpolarizing (+1) in its transport of 3 Na+ out (3*[+1]) for every 2 K+ in (2*[+1]).

$$ \left\{\begin{array}{l}{\mathrm{i}}_{\mathrm{ex}\_\mathrm{N}\mathrm{a}}^d = \hbox{-} 3\cdot \kern0.5em \left[+1\right]\cdot {g}_{ex}^d\\ {}{\mathrm{i}}_{\mathrm{ex}\_\mathrm{C}\mathrm{a}}^{\mathrm{d}} = 1\cdot \left[+2\right]\cdot {g}_{ex}^d\\ {}{\mathrm{i}}_{ex\_Net}^{\mathrm{d}} = \left({\mathrm{i}}_{ex\_Ca}^{\mathrm{d}}\kern0.5em -\kern0.5em {\mathrm{i}}_{ex\_Na}^{\mathrm{d}}\right)\kern0.5em \Rightarrow \kern0.5em {g}_{ex}^d\kern0.5em \cdot \left[-1\right]\kern0.5em \Rightarrow \kern0.5em -\kern0.5em {g}_{ex}^d\end{array}\right. $$ (138)

$$ \left\{\begin{array}{l}{\mathrm{I}}_{\mathrm{pump}\_\mathrm{N}\mathrm{a}}^d = 3\cdot \kern0.5em \left[+1\right]\cdot {g}_{pump}^d\\ {}{\mathrm{I}}_{\mathrm{pump}\_\mathrm{K}}^{\mathrm{d}} = 2\cdot \left[+1\right]\cdot {g}_{pump}^d\\ {}{I}_{pump\_net}^d = \left({I}_{pump\_Na}^d\kern0.5em -\kern0.5em {I}_{pump\_K}^d\right)\kern0.5em \Rightarrow \kern0.5em {g}_{pump}^d\kern0.5em \cdot \left[+1\right]\kern0.5em \Rightarrow \kern0.5em -\kern0.5em {g}_{pump}^d\end{array}\right. $$ (139)

gd ex and gd pump are Na+/Ca2+ exchanger and Na+/K+ pump membrane current densities (respectively) in the dendrite and their equality at 0.0021 mA/cm2 ensures an electrical counterbalance So, the model dendrite has a Na+/Ca2+ exchanger current and an electrically counterbalancing Na+/K+ pump current:

$$ \left[-{g}_{\mathrm{ex}}^{\mathrm{d}}\kern0.5em =\kern0.5em +\kern0.5em {g}_{\mathrm{pump}}^{\mathrm{d}}\right]\kern0.5em \Rightarrow \kern0.5em \left[-{I}_{\mathrm{ex}\_\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{d}}\kern0.5em =\kern0.5em +\kern0.5em {I}_{\mathrm{pump}\_\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{d}}\right] $$ (140)

Extracellular K+ dynamics [11]

Extracellular K+ concentration ([K+] o ) to the dendrite is initiated at 2 mM and then changes in time, t, according to the relationship:

$$ \frac{d{\left[{K}^{+}\right]}_0}{dt}=\frac{(1000)\cdot Q\cdot {I}_{K\_net}}{F\cdot wid} $$ (141)

$$ {I}_{K\_net}\kern0.5em =\kern0.5em {I}_{K\_ out}\kern0.5em -\kern0.5em {I}_{K\_ in} $$ (142)

$$ {I}_{K\_ out}={I}_{KD}+{I}_{KA}+{I}_{KM}+{I}_{DR}+{I}_{BK}+{I}_{K2}+{I}_{Kv1.2} $$ (143)

$$ {I}_{K\_ in}\kern0.5em =\kern0.5em {i}_{pump\_K}^d\kern0.5em +\kern0.5em {I}_{pump\_K}^d $$ (144)

where F is the Faraday constant, wid is the thickness of an extracellular region around the compartment that K+ accumulates in (70*10−3 μm), Q is a K+ accumulation factor (0.0119) and I K_net (Eq. 142) is the difference between K+ current flowing out of the compartment (I K_out ; through gated K+ conductances, Eq. 143) and K+ current pumped into the compartment (I K_in ; by the Na+/K+ pump, Eq. 144). id pump_K is set by Eq. 137; Id pump_K is set by Eq. 139. “Catch coding” is applied to [K+] o :

$$ \mathrm{if}\ \left({\left[{\mathrm{K}}^{+}\right]}_{\mathrm{o}}>3.03\right)\ \left[{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{o}} = 3.03\right] $$ (145)

$$ \mathrm{if}\ \left({\left[{\mathrm{K}}^{+}\right]}_{\mathrm{o}}<2\right)\ \left[{\left[{\mathrm{K}}^{+}\right]}_{\mathrm{o}} = 2\right] $$ (146)

Eq. 141 is based upon the NEURON code of [11], as opposed to the relevant equation in their associated journal paper; there is a minor discrepancy between the two (the paper does not have the dimensionality factor (1000) in its equation).