J Comput Neurosci DOI 10.1007/s10827-015-0561-9

Action potential initiation in a multi-compartmental model with cooperatively gating Na channels in the axon initial segment ¨ 1 · Min Huang2 · Fred Wolf3 Pinar Oz

Received: 10 September 2014 / Revised: 11 March 2015 / Accepted: 13 April 2015 © Springer Science+Business Media New York 2015

Abstract Somatic action potentials (AP) of cortical pyramidal neurons have characteristically high onset-rapidness. The onset of the AP waveform is an indirect measure for the ability of a neuron to respond to temporally fastchanging stimuli. Theoretical studies on the pyramidal neuron response usually involves a canonical Hodgkin-Huxley (HH) type ion channel gating model, which assumes statistically independent gating of each individual channel. However, cooperative activity of ion channels are observed for various cell types, meaning that the activity (e.g. opening) of one channel triggers the activity (e.g. opening) of a certain fraction of its neighbors and hence, these groups of channels behave as a unit. In this study, we describe a multi-compartmental conductance-based model with cooperatively gating voltage-gated Na channels in the axon initial segment. Our model successfully reproduced the somatic sharp AP onsets of cortical pyramidal neurons. The onset latencies from the initiation site to the soma and

Action Editor: Alain Destexhe

¨  Pinar Oz [email protected] 1

Department of Molecular Biology ¨ udar University, Istanbul, Turkey and Genetics, Usk¨

2

State Key Laboratory of Cognitive Neuroscience and Learning, Beijing Normal University, Beijing, China

3

Theoretical Neurophysics Group, Max Planck Institute for Self-Organization and Dynamics, G¨ottingen, Germany

the conduction velocities were also in agreement with the previous experimental studies. Keywords Action potential · Onset rapidness · Cooperativity hypothesis · Na channel · Axon initial segment · Multi-compartmental neuron model

1 Introduction Action potential (AP) generation and conduction in neurons form the basis for neuronal signalling and interneuronal communication. The collective behavior of individual neurons in a neural circuit, e.g. synchronization, strongly depends on the AP-generating dynamics of a neuron (Ermentrout et al. 2001), while initiation and propagation of APs depend on the activity of ion channels scattered on the neuronal membrane. Therefore, understanding the structure-function relationship of these ionic channels is crucial to understand the AP-generating mechanisms. AP waveforms recorded at the somata of cortical pyramidal neurons exhibit distinct sharp AP onsets. This phenomenon was explained by the back-propagation of spikes (McCormick et al. 2007; Shu et al. 2007) from the axon initial segment, which is the AP initiation site, to the soma (Colbert and Johnston 1996; Colbert and Pan 2002; H¨ausser et al. 1995; Kole et al. 2008; Kole and Stuart 2008; Palmer and Stuart 2006; Stuart and Schiller 1997). Several multicompartmental models based on canonical Hodgkin-Huxley type models (Hodgkin and Huxley 1952) with high sodium channel density in the AIS reproduced this phenomenon (Mainen et al. 1995). However, the channel densities used for these models were usually higher than the suggested densities in AIS (Kole et al. 2008; Fleidervish et al. 2010; Colbert and Pan 2002) and the models didn’t reproduce the

J Comput Neurosci

high cut-off frequencies observed in the cortical pyramidal neurons in vivo or in vitro (K¨ondgen et al. 2008; Naundorf et al. 2006). As reported by a recent theoretical study (Wei and Wolf 2011), smaller timescales at the AP onset help to reproduce the high cut-off frequencies of cortical pyramidal neurons. The findings imply that onset rapidness can be used as an indirect measure of cut-off frequencies. Therefore, we focus on the onset rapidness of AP in our study. The sharp somatic AP onsets can be described as the abrupt increase in the somatic membrane potential due to reaching a threshold interval for the Na channels at the initiation site (Brette 2013). Another way to explain the sharp AP onsets is the cooperativity hypothesis, which uses that description and further suggests the cooperative activity of Nav channels in AIS. In this hypothesis, the activity (e.g. opening) of one channel triggers the activity (e.g. opening) of a certain fraction of its neighbors and hence, these groups of channels behave as a unit (Naundorf et al. 2006; Huang et al. 2012). Although there is no direct evidence supporting the hypothesis up to date, there are many reports of similar type of activity for various kinds of ion channels. Considering the special molecular structure of the AIS cytoskeleton, localized distribution of Nav1.6 and Nav1.2 channels with higher densities compared to other neuronal regions, the possibility of such cooperative activity still holds for future investigation. As a third way, it was suggested that the geometrical discontunity between soma and AIS can account for the observations of sharp AP onsets at the soma, rather than back-propagation of spikes or the cooperativity of Na channels at AIS. This model design utilizes a compartmental structure and theoretical framework, which is very close to what is presented in this study. However, instead of the cooperativity, this model utilizes the distribution of Nav1.2 and Nav1.6 relatively denser at the proximal and distal AIS respectively. It was suggested that the difference in the activation kinetics of these two channels and their localizations, together with the compartmentalization of soma and AIS can result in the sharp AP onsets at the soma (Brette 2013). In this study, we describe a multi-compartmental conductance-based model with cooperative Na channel kinetics at AIS. Our results show that this model successfully reproduces the kink of the somatic AP onsets in cortical pyramidal neurons. 1.1 Cooperatively gating ion channels The distance between the ion channels is usually far enough to exclude a possibility of interaction. However, some ionic channels tend to form clusters on cellular membrane. For different cell types from Archea to mammalian neurons,

inter-channel coupling was shown to exist in these clusters (Dekker and Yellen 2006; Dixon et al. 2012; Gabelli et al. 2014; Huang 2006; Marx et al. 1998; Marx et al. 2001; Molina et al. 2006; Saito et al. 1988; Sonkusare et al. 2014; Undrovinas and Makielski 1992). Coupled gating dynamics of neighboring ion channels can be observed in a variety of channels - from the one with the simplest structure to much more complex ones. For instance, a K+ channel from the soil bacteria Streptomyces lividans, KcsA, was shown to exhibit unexpectedly complex opening and closing patterns (Molina et al. 2006), although it has the simplest K+ channel structure known up to date (Doyle et al. 1998). The patterns of KcsA channel opening suggest positive coupling of two to five channels and this result was also supported by confocal microscopy and FRET measurements (Molina et al. 2006). Inter-channel coupling was also shown for yet another K+ channel : hyperpolarization-activated cyclic nucleotide-gated (HCN) channel (Dekker and Yellen 2006). The study on HEK 293 cells showed that the single channel recordings in response to a voltage step resulted in variable delays with no channel openings, followed by synchronous activity which appeared highly correlated in time. This coupled activity was also evident from the analysis of stochastic fluctuations in microscopic currents. The authors suggested the allosteric communication between channels as the simplest hypothesis that explains the results (Dekker and Yellen 2006). Another example of cooperative gating was reported for ryanodine receptors. In a study that utilized electron microscopy to characterize the structure and kinetics of ryanodine receptors of the endoplasmic and sarcoplasmic reticulum (Saito et al. 1988), the ryanodine receptors were shown to form clusters. As mentioned previously, the clustering of membrane proteins may lead to functional interprotein coupling (Huang 2006). Therefore, such interaction could also be expected for ryanodine receptors. This prediction was confirmed by following studies (Marx et al. 1998; Marx et al. 2001). The results pointed out that ≈ 10% of ryanodine receptors act as functionally coupled. It was also intriguing that the electron micrographs of purified RyR1 homotetramers, that were obtained in the previous study (Saito et al. 1988), showed that ≈ 10% of the structures are physically connected to form contacting pairs. A recent study indicated that L-type CaV 1.2 channels in mouse ventricular myocytes also exhibit oligomerization via their C termini (Dixon et al. 2012). This oligomerization resulted in the amplification of Ca2+ influx into the cells. The authors proposed that the binding at C-tails of two neighbouring CaV 1.2 channels makes them allosteric activators of adjoining channels. This binding was suggested

J Comput Neurosci

to increase the cooperativity between neighboring channels, and therefore, observed multi-channel openings of these channels result from the coordinated opening of adjoined channels and coincident acivation of independently gating channels. The cooperative gating of NaV channels was shown on cardiomyocytes from rat and rabbit (Undrovinas and Makielski 1992). The study indicated that lysophosphatidylcholine (LPC)-treated NaV channels exhibit synchronized opening patterns that can only be explained by inter-channel coupling. Mechanisms such as changes in the membrane fluidity or in the lipid protein interaction, the activation of mechanisms, that are frequently inactive in non-LPC-treated cells (e.g. increase in cAMP due to increasing levels of LPC at heart tissue) and the disruption of cytoskeleton due to LPC-induced ischemia, which leads to further changes in the membrane structure were suggested as possible explanations for LPC-modified cooperative gating of NaV channels. In addition to the cardiomyocyte NaV channels, a recent structural study (Gabelli et al. 2014) showed that two NaV 1.5 channels can form an asymmetrical homodimer via the interactions of calmodulins that are structurally linked to the carboxy terminals of NaV 1.5 channels, forming a complex. Although biological evidences about the cooperative activity of ionic channels are provided, details of the underlying molecular mechanisms for such phenomena are still unknown. To have the level of synchronous activity in a population of ionic channels as presented in the previously mentioned studies, it is clear that the gating state changes of one channel must be transmitted to the neighboring channels. For this type of allosteric communication between ionic channels, several mechanisms are possible candidates : 1. A protein-protein interface between channels, i.e. due to dimer, trimer or tetramer formations (Neumcke and St¨ampfli 1983; Almers and Stirling 1984), 2. Reconciliation by an adaptor protein, 3. Allosteric changes in a linking cytoskeleton (Undrovinas and Makielski 1992; Post et al. 1985; Grubb and Burrone 2010), 4. Interactions between subdomains of channel proteins, e.g. that are the binding sites of cytosolic secondary messengers like cAMP, Ca2+ or PIP2 (Dekker and Yellen 2006). Other than allosteric communication, the coupling through the changes in the local electric field is yet another possibility. Furthermore, the voltage-dependent or independent elastic changes of the lipid membrane may also lead to conformational changes in the ionic channels, that

may contribute to inter-channel coupling (Ursell et al. 2007; Silberberg and Magleby 1997). The specific molecular structure of axon initial segment, which is the initiation site of action potentials for neurons in general, an allosteric communication through linking cytoskeletal structures might coordinate cooperative NaV channel gating in this region. 1.2 A neuron model with cooperatively gating Na channels Voltage-gated Na+ channels play the major role in the onset and the upstroke of the action potentials. Therefore, the cooperative gating of these channels would result in significant changes at the action potential initiation timescale. It was shown that the onset rapidness could be an indirect measure of cut-off frequency (Fourcaud-Trocme et al. 2003). Furthermore, the recent in vitro studies provided information about high cut-off frequencies of real neurons (K¨ondgen et al. 2008; Boucsein et al. 2009) and a recent theoretical work reported that smaller timescales at the AP onset help to reproduce such high cut-off frequencies (Wei and Wolf 2011). Combining these results, it can be concluded that very rapid onsets, i.e. in cortical neurons (Naundorf et al. 2005a, 2006) can be associated high cut-off frequencies in physiological conditions. It is also reported that cooperative gating kinetics, when implemented in classical Hodgkin-Huxley formulations, can indeed reproduce two most important effects of the cortical neuron dynamics: rapid AP onsets and threshold variability (Huang et al. 2012). Therefore, it is alluring to test this hypothesis on a more realistic neuron model, and i.e. implementing the dynamics in AIS. To test the cooperative gating hypothesis, a conductancebased multi-compartmental neuron model was designed such that Hodgkin-Huxley type conductance algorithms are the skeleton of the model and the cooperative Na+ channel gating kinetics are implemented only at the AIS. In addition to analyzing the effects of this implementation, the multi-compartmental cooperative axon initial segment (CAIS) model can be used to investigate some important aspects of real neurons, such as : • • •

The action potential initiation dynamics and the site of initiation, The variations in the Na+ channel density distribution in different parts of a neuron, The action potential propagation and the conduction velocity.

The Na+ channel density in the site of action potential initiation is still a hot debate. The suggestions of very high Na+ channel densities in AIS compared to soma was based

J Comput Neurosci

on the expectation that, when similar ion gating properties are assumed, a greater density of sodium channels in AIS should lead to a significantly larger Na+ current in this site . This expectation is challenged with the patch clamp studies, which revealed that there is only a small difference (3-10 fold) in Na+ current between soma, proximal dendrites and AIS (Colbert and Pan 2002; Kole et al. 2008). These findings are incompatible with what should be expected from large density differences (e.g. 50-1000 fold of soma Kole et al. 2008 Mainen et al. 1995). Therefore, in our simulations, we assumed a Na+ channel density that is only 3-10 fold of soma. Implementation of the cooperative gating of Na+ channels in the axon initial segment requires a model for the coupled activity of channels. The basis of the model in this paper was previously established by Huang et al. (2012) to study the effects of cooperativity on a single-compartment model. First of all, it is assumed that a NaV channel is coupled to K neighboring channels. Each member of this coupled neighborhood increases the opening probability of its neighbors. Second, a coupling constant in mV, J ,is introduced as a measure of coupling strength. J is defined by the shift along the activation curve that would increase the opening probability of a single, isolated channel the same amount as a coupled channel. The coupling strength depends on both the coupling constant J and the number of coupled neighbors K. For simplicity, the coupling strength will be referred as KJ. Third, it is assumed that only a fraction (percentage, p) of the Na+ channels cooperatively gate in the given membrane surface area. The activation kinetics of the cooperative NaV channels was explained in detail in the previous paper (Huang et al. 2012). Our study extends the findings presented in the previous paper by implementing the framework for cooperative activity in each of the AIS segments in our multi-compartmental model. A high value of KJ would imply strong coupling between neighboring channels. Depending on the cooperativity percentage, the channels that are coupled in the unit area increase and this indirectly affects the coupling strength, as the number of the channels involved is also included in its estimation. Therefore, if there are 100 channels in the unit area, a coupling strength of 1000 mV with 10 % cooperativity would mean that the inter-channel coupling of two channels would create an effect in the gating of the channel as much as a 100 mV increase in the membrane potential would create. Depending on their functional roles, many neurons are intrinsically very noisy, therefore, exhibit spontaneous firing even in the absence of an external stimulus. Also, the outgoing spike train was considered as a random process even in the presence of a stimulus. Both the stimulus and

intrinsic noise trigger the spike trains in a cooperative fashion (Wiesenfeld and Moss 1995). Therefore, in a realistic simulation, the input current is presented with a noisy component in addition to the deterministic part.

2 Materials and methods The model was implemented using C++ in Linux environment. The implicit backward Euler integration method was used with a 10 μs time step. To maintain a constant ν0 between trials for specific experimental designs, a current adapter was also included such that μI of the injected current was altered automatically until the desired ν0 is obtained. For the analysis of the action potential onset rapidness , V˙ (tthreshold ) = 10V /s is chosen to define the threshold of an action potential. Then, the slope of the phase plot (V vs V˙ ) at V˙ (tthreshold ) is measured to obtain . Data analysis were performed in Matlab with user-defined protocols. The source code for the model is also available online at ModelDB website. 2.1 Geometry and passive properties of the model The multi-compartmental cooperative axon initial segment model consists of successive cylindirical compartments,that vary in size and electrical properties to reflect the characteristics of a particular part of the model. The basic parts can be listed as –

– –

–

a thin cable-like extension with homogeneous thickness and homogeneous electrical properties to represent the dendritic assembly, an octagon-shaped ball with homogeneous electrical properties to represent soma, a thin cable-like extension with heterogeneous thickness and heterogeneous electrical properties to represent the axon, a small octagon-shaped bleb to end the axon.

The axon itself can be divided into further subgroups as the axon hillock (AH), the axon initial site (AIS), the myelinated compartments and the nodes of Ranvier (NR). The details of the model are given in Fig. 1 and Table 1. The membrane capacitance ( Cm ) was 0.1 nF / mm2 . The membrane resistance (rm ) was 1 M.mm2 , longitudinal membrane resistance, rL , was either 3 or 5 M.mm and the membrane time constant was estimated as 10 ms. A summary of the passive parameters for each part is given in Table 1.

J Comput Neurosci

step by instantaneous activation. This model exhibits type-I excitability and it has been employed in various theoretical studies to explore the dynamical response properties of cortical neurons (Fourcaud-Trocme et al. 2003; Wang and Buzs`aki 1996; Fourcaud-Trocme and Brunel 2005; Naundorf et al. 2005b). The channel gating kinetics for a gating particle z is modelled by the equation

τz

dz = z∞ − z(t) dt

(1)

where τz is the time constant and z∞ is the steady-state opening probability (or the limiting value). Each gating particle is assumed to switch between open and close states with an opening rate, α(V ), and a closing rate, β(V ), respectively. Using these rates, the dynamics of the time constant and the steady-state probability is formulated as

τz =

 , α(V ) + β(V )

z∞ =

Fig. 1 The geometry of Multi-Compartmental Cooperative AIS (CAIS) model

2.2 Channel gating kinetics I : Statistically independent gating CAIS model utilized a modified version of the conductancebased Wang-Buzsaki (WB) model (Wang and Buzs`aki 1996), that is a Hodgkin-Huxley type model and therefore, an essential assumption of the model is that the activity of individual channels is independent from its neighboring channels. In WB model, Na+ channels respond to a voltage

(2)

α(V ) , α(V ) + β(V )

(3)

where  is a constant.To maintain a maximum activation time constant of 50 μs, which is suggested by the sodium current measurements on cortical neurons (Baranauskas and Martina 2006),  was set to 0.1. The opening probability, m(t), and the closing probability, h(t), of NaV constructs the general formulations for the NaV conductance as gNa = m(t)3 h(t)g¯ Na

(4)

where g¯ Na is the unit conductance for the channel. The activation of Na+ channel is instantaneous, therefore m is substituted by the steady-state function, m∞ . In our model,

Table 1 The parameters of multi-compartmental cooperative AIS (CAIS) model. Please also see Fig. 1 and text PART

Ncomp

lcomp (μm)

a∗in (μm)

a∗out (μm)

Cm (nF/μm2 )

Dendrite 1 Dendrite 2 Soma Axon Hillock Axon Initial Segment Myelinated Segment Node of Ranvier Bleb

6 10 20 n 25 - n 5 1 5

50 10 2 2 2 10 2 2

0.5 0.5 1 1.5 1.5- ax 0.5 0.5 0.5

0.5 1 1.5 1.5 - ax 0.5 0.5 0.5 0.5

10−5 10−5 10−5 10−5 10−5 2*10−6 10−5 10−5

J Comput Neurosci

K+ channel conductance is also included with the opening probabiliyt n(t) and the unit conductance g¯ K

and coop dhcoop

gK = n(t) g¯ K . 4

(5)

The opening and closing rates are dependent only on local membrane potential, V, as given in Table 2.

τh

coop

τh

dt

coop

= h∞

− hcoop (t),

(10)

(V ) = τh (V + Vshif t ),

(11)

2.3 Channel gating kinetics II : Cooperative gating coop

The cooperative gating model assumes that a single channel is coupled to K neighbor channels, such that the opening of one channel increases the opening probability of its neighbors. It might also be assumed that this cooperative activity will have similar effects with a voltage shift in the membrane potential. Therefore, a unit coupling strength in mV, J , can be defined such that it represents a voltage shift in the membrane potential that results in the same amount of increase in the opening probability of a single channel. Then, cooperativity can be implemented in the classical model, by giving the voltage shift Vshif t = KJ (mcoop )3 hcoop

(6)

and adding it in the canonical equations, coop dmcoop

τm

coop

τm

dt

coop

= m∞

− mcoop (t),

(7)

(V ) = τm (V + Vshif t ),

(8)

coop

m∞ (V ) = m(V + Vshif t ),

(9)

h∞ (V ) = h(V + Vshif t ).

(12)

The opening of individual channels are modeled as a Markov process; and the eqs. (7) and (10) represent the mean field approximation of cooperative channel gating among a coupled population (Naundorf et al. 2006). It is also assumed that only some fraction p of the channel population in a certain area couple. Therefore, the cooperative gNa can be rewritten as   gNa = (1 − p)m3 h + pm3coop hcoop g¯ Na

(13)

and p represents the percentage of coupling channels in the population. Cooperative channel gating was implemented with the assumption that only the Nav channels of AIS were act cooperatively. Therefore, the cooperativity parameters (cooperativity percentage p and coupling strength KJ ) were set to zero in the remaining compartments. The changes in mcoop and m during an AP can be seen in Fig. 10. Somatic g¯ Na , gsoma , was equal to 300 pS/μm2 and remained constant. The conductance in the the dendrites (gdend ), axon hillock (gAH ), myelinated parts (gmye ) and terminal bleb (gbleb ) also remained constant and were set as shown in Table 3. The conductances in the other parts, namely in AIS (gAI S ) and nodes of Ranvier (gNR ) were varied from 1-10-fold gsoma .

Table 2 The voltage dependence of Na+ and K+ type channels in Wang-Buzsaki model (Wang and Buzs`aki 1996) Table 3 The channel density distribution for CAIS Channel Type

Gating Particle

V-dependence

Na+ channel

m

αm =

0.1(V +35) 1−exp[−0.1(V +35)]

βm = 4exp h

K+ channel

n



−(V +60) 18

αh = 0.35exp





−(V +58) 20

βh =

5 1+exp[−0.1(V +28)]

αn =

0.05(V +34) 1−exp[−0.1(V +34)]

βn = 0.625exp





−(V +44) 80



PART

gN a (pS/μm2 )

gN a gK gLeak (gsoma ) (pS/μm2 ) (pS/μm2 )

Dendrite Soma Axon Hillock Axon Initial Segment Myelinated Segment Node of Ranvier Bleb

100, 300 300 300, 3000 300 − 3000 0 300 − 3000 300

0.3 ,1 1 1, 10 1 − 10 0 1 −10 1

150 150 150 150 150 150 150

1 1 1 1 0.2 1 1

A dash between two values means that the value of parameter was chosen in this interval. A comma between two values means that the parameter took either the first or the second value

J Comput Neurosci

2.4 Current injection In the brain, neurons receive information under a bombardment of excitatory and inhibitory synaptic inputs. The activity of numerous synapses in the cortex creates a noisy background for the information processing of a single neuron. This background noise can be characterized by Gaussian statistics and by an autocorrelation function exponentially decaying with the time constant τI . Fast synaptic currents such as AMPA- and GABAA -mediated currents have a τI in the range of 5-20 ms. Therefore, this range is used in our simulations. The noisy current (Inoisy ) was generated as a realization of an Ornstein-Uhlenbeck stochastic process with zero mean and variance σI2 to mimic the synaptic noisy input as following : ˙ = −Inoisy + κζt . τI Inoisy

(14)

where ζt is a random variable drawn at every time step from a Gaussian distribution with a zero mean, called white noise; κ = X/ t and ζt = ξt / t. Integrating the deterministic part gives us 

t Inoisy (t + t) = Inoisy (t)exp − τI

 + Xξt

3-5 K.mm to observe its effect on the action potential propagation in the model. Setting intracellular resistivity(rL ) to 5 K.mm and the sodium channel density at axon initial segment (gAI S ) and Node of Ranvier (gNR ) as 10-fold of the density at soma (gsoma ), the initiation site of the model was shown to be approximately 45 μm away from the soma, which is inside the borders of the axon initial segment (AIS) in our model. The action potentials initiated at the AIS reached to the soma by retrograde propagation in 200 μs (Fig. 2). This value of delay is in aggrement with the previous reports (Yu et al. 2008; Palmer and Stuart 2006). When gAI S and gNR were both 3-fold of gsoma , the latency to soma was decreased to 100 μs. However, with the implementation of the cooperative sodium channel activity at axon initial segment, with a cooperativity fraction of 10% and a coupling strength of 400 mV, the latency was again elevated to 200 μs level (Fig. 3). The orthodromic and antidromic propagation of action potentials were analyzed using the onset latencies ( tonset ) or peak latencies ( tpeak ) of the first action potentials that were initiated at AIS. The latencies are measured as the difference between the onset/peak time at AIS and the

(15)

and X can be set as  X = σI



2 t 1 − exp τI

 (16)

3 Results 3.1 Action potential initiation site in CAIS model depends on the ratio of gAI S /gsoma and gAI S /gNR The experimental findings up to date support the view that AIS is the initiation site in most of the myelinated neurons of the central nervous system. The experimental and theoretical studies emphasize the importance of ionic, i.e. Na+ , channel density distribution and gating kinetics in determining the location of the action potential initiation in a neuron (Colbert and Pan 2002; Mainen et al. 1995). Therefore, Na+ channel density and the parameters that affect and alter the channel gating kinetics in a compartment were altered systematically to investigate their impact on the site of initiation. Furthermore, the cooperativity parameters, such as the percentage of Na+ channels that cooperatively gate, p, and the coupling strength of these channels, KJ , were also altered to understand their effect on the action potential initiation. In addition to active electrical parameters, the intracellular resistivity was also altered between

Fig. 2 Action potential initiation and propagation for gAI S = gN R = 10*gsoma (non-cooperative model). Upper: The contour plot for the propagation of an action potential on CAIS model. Lower: The onset latency-distance plot, that was obtained from the same simulation. The somatic onset latency is around 200 μs, which is in agreement with previous studies (Palmer and Stuart 2006; Kole et al. 2008; Hu et al. 2009). A constant current was injected at the soma to produce a steady-state firing rate ν0 = 10 Hz. gAI S = gN R = 10*gsoma = 3000 pS/μm2 , rL = 5 K.mm. The model was non-cooperative, and dendritic conductances were set as same with soma

J Comput Neurosci Table 4 Axonal conduction velocities with respect to channel densities in AIS and NR

Fig. 3 Action potential initiation and propagation for gAI S = gN R = 3*gsoma (cooperative and non-cooperative model). The dots represent the numerical results. These results were fit using a linear function of latency and distance (solid blue line) and the conduction velocity of each case was estimated from the slopes of the fitting functions (as given next to the solid line). The conduction velocities for both cases were the same.Red : cooperative (p = 10 %, KJ = 400 mV) and rL = 5 K.mm; black : non-cooperative and rL = 5 K.mm. A noisy current with τI = 5 ms and an automatically adjusting μI was injected to maintain ν0 = 10 Hz. σI was set to give σV ∼ = 5 mV

onset/peak time at other compartments. The conduction velocities were estimated from fitting the slope of onset latency versus the distance that the action potential travels on both directions (Figs. 3 and 4). Given same relative Na+ channel densities in AIS, NR and soma (gAI S :gNR :gsoma ) and same rL , the axonal conduction velocity was 0.53 m/s for both cooperative and non-cooperative cases (Fig. 3). When the intracellular resistivity rL was 3 K.mm instead of 5 K.mm, the conduction velocity increased to 0.9 m/s (Fig. 4). A summary for the effects of different channel density values in AIS, NR and soma on the axonal and dendritic conduction velocities is given in Tables 4 and 5.

Channel Density

gN R = gsoma

gN R = 3*gsoma

gAI S = gsoma gAI S = 3*gsoma

0.48 m/s 1 m/s

0.48 m/s 0.9 m/s

The axonal conduction velocities were estimated with the slope of increase in the onset latency between x1 = 102 μm and x2 =362 μm for the combinations of channel densities in AIS and NR. The model was non-cooperative. A noisy current with τI = 5 ms and an automatically adjusting μI was injected to maintain ν0 = 10 Hz. σI was set to give σV ∼ = 5 mV

The results indicated that the only factor that affected the site of initiation is the channel density distribution (Fig. 4). A trivial prediction for the effect of rL on the conduction velocity would be that as the rL increases, the conduction velocity should decrease. This prediction was approved by the simulation results. For the non-cooperative model, the conduction velocity for rL = 3 K.mm was 0.9 m/s and for rL = 5 K.mm was 0.5 m/s. The cooperativity parameters (p and KJ ) were ineffective on the conduction velocity regardless of the intracellular resistivity. However, the coupling strength had a clear impact on the somatic onset latency. 3.2 Onset Rapidness of Axonal and Somatic Action Potentials The AP waveforms (Fig. 5) of CAIS model clearly indicate the difference in onset rapidness () between canonical and cooperative models : In cooperative model the onset grows piecewise-linearly (“kinky onset”), opposite to the canonical models, where the onset is growing exponentially (“smooth onset”). The phase plots also revealed this pattern (Fig. 6). To obtain quantitative results,  was defined as the phase plot slope at V˙ (x, ti ) > 10 mV/ms and was obtained for increasing values of cooperativity parameters (coupling strength, KJ , and cooperativity percentage, p). Table 5 Dendritic conduction velocities with respect to channel densities in AIS and NR

Fig. 4 Action potential initiation and propagation for gAI S = gN R = 3*gsoma (non-cooperative model). The dots (gray and purple) represent the numerical results. These results were fit using a linear function of latency and distance (solid lines,black and red) and the conduction velocity of each case was estimated from the slopes of the fitting functions (as given next to the solid lines). The conduction velocities for both cases were the same.gray, black : rL = 5 K.mm; purple,red : rL = 3 K.mm;. A noisy current with τI = 5 ms and an automatically adjusting μI was injected to maintain ν0 = 10 Hz. σI was set to give σV ∼ = 5 mV.

Channel Density

gN R = gsoma

gN R = 3*gsoma

gAI S = gsoma gAI S = 3*gsoma

0.58 m/s 0.67 m/s

0.63 m/s 0.56 m/s

The dendritic conduction velocities were estimated with the slope of decrease in the onset latency between x1 = -400 μm and x2 =-100 μm for the combinations of channel densities in AIS and NR. The model was non-cooperative. A noisy current with τI = 5 ms and an automatically adjusting μI was injected to maintain ν0 = 10 Hz. σI was set to give σV ∼ = 5 mV

J Comput Neurosci Fig. 5 AP waveforms in CAIS Model (gsoma = gdend ). The APs were obtained from the simulations on non-cooperative and cooperative CAIS models. The cooperativity parameters are as given in the figure. gAI S = gN R = 3*gsoma = 3*gdend . A noisy current with τI = 5 ms and an automatically adjusted μI was injected to maintain ν0 = 10 Hz. σI was set to give σV ∼ =5 mVThe arrows indicate the “kink”s.

The results from CAIS model confirmed that, APs with high action potential onset rapidness as in cortical neurons (Naundorf et al. 2006) can be generated when cooperative Na+ channel gating is implemented in a classical HodgkinHuxley type model (Fig. 5). This behavior is evident without the need of unrealistically high Na+ channel density at AIS (Mainen et al. 1995). For instance, the “kink” in Fig. 5 was obtained at gAI S = 3*gsoma with a strong (KJ = 1000 mV) coupling of a small fraction (p = 10%) of channels. The “kinks” were also evident in the comparison of phase plots for cooperative and non-cooperative models (Fig. 6).

The simulations were performed for both gsoma = gdend and gsoma = 3*gdend . In both cases, the AP waveforms and phase plots were similar (Figs. 7 and 8) and cooperative models reproduced the “kinks” that were observed for the onset of cortical APs (Naundorf et al. 2006). The results also indicated an increasing axonal onset rapidness with increasing coupling strength and the analysis of  as a function of x indicated that axonal peak point of onset rapidness moved closer to the AIS by increasing coupling strength, where the cooperative channel gating was implemented. When measured under physiological conditions, APs of cortical neurons are at least 20 ms−1 . In our CAIS model, this value was reached and/or exceeded at KJ = 600 mV for p = %20 and at KJ = 900 mV for p = %10. The spatial shift of the maximum  towards AIS was not affected by the cooperativity percentage p (Fig. 9). On the other hand, it was linearly dependent on KJ (shift rate = 0.14 μm/mV) (Fig. 10).

4 Discussion

Fig. 6 Phase Plot Diagrams (gsoma = gdend ). The phase plot diagrams of the APs given in Fig. 5. All of the APs produced in a 1 s time interval was included for phase plots. Figure 5 displays isolated APs from the same set. (Left: Phase plots for soma. Right:Phase plots for AIS. Black traces indicate the non-cooperative model, red traces indicate the cooperative model (p = 10%, KJ = 800 mV). gAI S = gN R = 3*gsoma = 3*gdend . A noisy current with τI = 5 ms and an automatically adjusted μI was injected to maintain ν0 = 10 Hz. σI was set to give σV ∼ = 5 mV

The parameters such as AP propagation and conduction velocities, AP waveform and onset rapidness and the effects of spatial variations in Na+ channel densities were first used to characterize the non-cooperative multi-compartmental model. After the completion of this step, the contribution of cooperative Na+ channel gating was investigated and was compared to the non-cooperative case. The findings of the CAIS model emphasize some intriguing points. First of all, cooperative model reproduces the target feature of cortical pyramidal neuron APs, which is

J Comput Neurosci Fig. 7 AP waveforms in CAIS Model (gsoma = 3*gdend ). The APs were obtained from the simulations on non-cooperative and cooperative CAIS models. The cooperativity parameters are as given in the figure. gAI S = gN R = 3*gsoma = 9*gdend . A noisy current with τI = 5 ms and an automatically adjusted μI was injected to maintain ν0 = 10 Hz. σI was set to give σV ∼ =5 mV The arrows indicate the “kink”s.

the frequently observed “kink” at the somatic AP onsets (Naundorf et al. 2005a, 2006; Huang et al. 2012; Volgushev et al. 2008; Baranauskas et al. 2010). This feature wasn’t successfully reproduced in most of the studies utilizing canonical Hodgkin-Huxley type models. This result was also previously reported by Naundorf et al. (2006) and it was suggested that the contribution of cooperative gating kinetics in AIS is a possible mechanism to reproduce the cortical AP onset and threshold variability. This suggestion is opposed by the argument that the rapid AP onset in cortical neurons can be explained with the back-propagation of axonal APs (Kole et al. 2008; McCormick et al. 2007; Yu et al. 2008) and the counterarguments were also discussed

Fig. 8 Phase Plot Diagrams (gsoma = 3*gdend ). The phase plot diagrams of the APs given in Fig. 7. All of the APs produced in a 1 s time interval was included for phase plots. Figure 7 displays isolated APs from the same set.Right:Phase plots for AIS. Black traces indicate the non-cooperative model, red traces indicate the cooperative model (p = 10%, KJ = 1000 mV). gAI S = gN R = 3*gsoma = 9*gdend . A noisy current with τI = 5 ms and an automatically adjusted μI was injected to maintain ν0 = 10 Hz. σI was set to give σV ∼ = 5 mV

in detail in several reports (Baranauskas et al. 2010; Naundorf et al. 2007; Brette 2013). It should be noted that the back-propagation hypothesis doesn’t account for the sharpness of spike initiation and also that the spike initiation site is quite close to soma compared to the electrotonic length of the axon. Secondly, CAIS model simulations indicated that the cooperative model reproduces the rapid AP onsets with Na channel density only 3- to 10-fold higher than the soma or without a distortion in the actual geometry of axon hillock and axon initial segment (Mainen et al. 1995; McCormick et al. 2007). Furthermore, the coupling strength and the cooperativity percentage were ineffective on the localization of action potential initiation, however, they were critically effective in increasing the onset rapidness of action potentials.

Fig. 9 The location of maximum  as a function KJ . The red dashed lines indicate the interval that corresponds to AIS. The spatial shift of maximum  is linearly dependent on KJ and the slope of the fit line is same for both p = 10 % and p = 20 %. gAI S = gN R = 3*gsoma = 3*gdend . A noisy current with τI = 5 ms and an automatically adjusted μI was injected to maintain ν0 = 10 Hz.σI was set to give σV ∼ = 5 mV

J Comput Neurosci Fig. 10 mcoop vs m at AIS and soma. Upperleft: Action potentials at AIS (red) and soma(black). Upperright: The change in m during an action potential at AIS (red) and soma(black). Lowerright: The change in mcoop during an action potential at AIS (red) and soma(black). Lowerleft: m vs mcoop during an action potential at AIS (red) and soma(black). Note that only AIS is implemented with cooperativity, therefore mcoop and m in soma is actually the same

The use of blebs as the site of axonal recordings (McCormick et al. 2007; Shu et al. 2007) present several difficulties for the accurate interpretation of Na+ channel densities or interactions at the axon. It is well-known that, in the case of such injuries, the axonal cytoskeleton rearranges itself (Schafer et al. 2009). This is a very crucial point, since the cytoskeletal proteins, especially of the ankyrin and spectrin families, is known to be important for the recruitment and localization of the ionic channels on the neuronal membrane (Grubb and Burrone 2010; Angelides et al. 1988; Dzhashiashvili et al. 2007). Indeed, it was shown that neuronal injury can cause rapid and irreversible proteolysis of the AIS cytoskeleton and loss of ion channel clusters (Schafer et al. 2009). Therefore, we suggest that the possibility of cooperative activity at AIS should be investigated with the recordings from healthy neurons, i.e. in vitro on brain slices or in vivo. Although each technique presents its own difficulties, recent progress in the patch-clamp methods and the combination with advanced optical techniques might provide more accurate information. A recent study (Brette 2013) presents a model based on the localized distribution of Nav1.2 and Nav1.6 channels in AIS geometric discontunity between soma and AIS to explain the kink at somatic AP onsets. Alhough this model utilizes a similar theoretical framework and a similar, yet simplified, compartmental structure compared to the CAIS model, the two models differ in the implementation of cooperativity at AIS. Since our model considered only one type of Na channel with varying channel densities in each subregion of the model, our non-cooperative model was canonical

Hodgkin-Huxley type. The notion of using the localization of Nav1.2 and Nav1.6 channels would be intriguing in such model and should be further investigated both theoretically and experimentally. Another recent study showed that the cable proprties of the dendritic tree may directly contribute to the sharpness of AP onsets at the initiation site (Eyal et al. 2014). It was suggested that partially acivated voltage-gated Na conductance during the AP onset makes the rise time of AP onset susceptible to effective dendritic impedance load. Therefore, the changes in the dendritic surface can strongly modulate the AP onset rapidness at AIS. This recent findings were also suggested or implied in some of the previous experimental and theoretical studies (Bekkers and Hauser 2007; Hay et al. 2013; Mainen et al. 1995). Although the contribution of dendritic load on the AP initiation and onset at AIS wasn’t included in CAIS model, this study presents important findings that should be further studied. The onset dynamics of the AP is shown to critically determine the capability of a neuron to effectively respond to temporally fast-changing inputs (Fourcaud-Trocme et al. 2003). Therefore, the effect of cooperative channel gating on the linear response properties should be also further investigated.

Acknowledgments We thank A. Neef, M. Monteforte, W.Wei and M. Gutnick for fruitful discussions.

Conflict of interests of interest.

The authors declare that they have no conflict

J Comput Neurosci

References Almers, W., & Stirling, C. (1984). Distribution of transport proteins over animal cell membranes. The Journal of Membrane Biology, 77(3), 169–186. Angelides, K., Elmer, L., Loftus, D., & Elson, E. (1988). Distribution and lateral mobility of voltage-dependent sodium channels in neurons. The Journal of Cell Biology, 106(6), 1911–1925. Baranauskas, G., & Martina, M. (2006). Sodium currents activate without a hodgkin-and-huxley-type delay in central mammalian neurons. The Journal of Neuroscience, 26, 671–684. Baranauskas, G., Mukovskiy, A., Wolf, F., & Volgushev, M. (2010). The determinants of the onset dynamics of action potentials in a computational model. Neuroscience, 167(4), 1070–1090. Bekkers, J., & Hauser, M. (2007). Targeted dendrotomy reveals active and passive contributions of the dendritic tree to synaptic integration and neuronal output. Proc Nal Acad Sci USA, 104, 11, 447–11, 452. Boucsein, C., Tetzlaff, T., Meier, R., Aertsen, A., & Naundorf, B. (2009). Dynamical response properties of neocortical neuron ensembles: multiplicative versus additive noise. The Journal of Neuroscience, 29(4), 1006–1010. Brette, R. (2013). Sharpness of spike initiation in neurons explained by compartmentalization. Plos Comput Biol, 9(12), e1003, 338. Colbert, C., & Johnston, D. (1996). Axonal action-potential initiation and na+ channel densities in the soma and axon initial segment of subicular pyramidal neurons. The Journal of Neuroscience, 16(21), 6676–6686. Colbert, C., & Pan, E. (2002). Ion channel properties underlying axonal action potential initiation in pyramidal neurons. Nature Neuroscience, 5, 533–538. Dekker, J., & Yellen, G. (2006). Cooperative gating between single hcn pacemaker channels. The Journal of General Physiology, 128(5), 561–567. Dixon, R.E., Yuan, C., Cheng, E.P., NAvedo, M.F., & Santana, L.F. (2012). Ca2+ signaling amplification by oligomerization of l-type cav1.2 channels. Proc Nal Acad Sci, 109(5), 1749–1754. Doyle, D., Cabral, J.M., Pfuetzner, R., Kuo, A., Gulbis, J., Cohen, S., Chait, B., & MacKinnon, R. (1998). The structure of the potassium channel: molecular basis of k+ conduction and selectivity. Science, 280(5360), 69–77. Dzhashiashvili, Y., Zhang, Y., Galinska, J., Lam, I., Grumet, M., & Salzer, J. (2007). Nodes of ranvier and axon initial segments are ankyrin-g-dependent domans that assemble by distinct mechanisms. The Journal of Cell Biology, 177, 857–870. Ermentrout, B., Pascal, M., & Gutkin, B. (2001). The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators. Neural Computation, 13, 1285– 1310. Eyal, G., Mansvelder, H., de Kock, C., & Segev, I. (2014). Dendrites impact the encoding capabilities of the axon. The Journal of Neuroscience, 34(24), 8063–8071. Fleidervish, I., Lasser-Ross, N., Gutnick, M., & Ross, W. (2010). Na+ imaging reveals little difference in action potential-evoked na+ influx between axon and soma. Nature Neuroscience, 13(7), 852– 860. Fourcaud-Trocme, N., & Brunel, N. (2005). Dynamics of the instantaneous firing rate in response to changes in input statistics. Journal of Computational Neuroscience, 18(3), 311–321. Fourcaud-Trocme, N., Hansel, D., van Vreeswijk, C., & Brunel, N. (2003). How spike generation mechanisms determine the neuronal response to fluctuating inputs. The Journal of Neuroscience, 23(37), 11, 628–11, 640. Gabelli, S.B., Boto, A., Kuhns, V.H., Bianchet, M.A., Farinelli, F., Aripala, S., Yoder, J., Jakoncic, J., Tomaselli, G.F., & Amzel,

L.M. (2014). Regulation of the nav1.5 cytoplasmic domain by calmodulin. Nature Communications, 5, 5126. Grubb, M., & Burrone, J. (2010). Building and maintaining the axon initial segment. Current Opinion in Neurobiology, 20(4), 481–488. H¨ausser, M., Stuart, G., Racca, C., & Sakmann, B. (1995). Axonal initiation and active dendritic propagation of action potentials in substantia nigra neurons. Neuron, 15(3), 637–647. Hay, E., Schrmann, F., Markram, H., & Segev, I. (2013). Preserving axosomatic spiking features despite diverse dendritic morphology. Journal of Neurophysiology, 109, 2972–2981. Hodgkin, A.L., & Huxley, A.F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology, 117, 500–544. Hu, W., Tian, C., Li, T., Yang, M., Hou, H., & Shu, Y. (2009). Distinct contributions of na(v)1.6 and na(v)1.2 in action potential initiation and backpropagation. Nature Neuroscience, 12(8), 996– 1002. Huang, H. (2006). Molecular mechanism of antimicrobial peptides: the origin of cooperativity. Biochimica et Biophysica Acta, 1758(9), 1292–1302. Huang, M., Volgushev, M., & Wolf, F. (2012). A small fraction of strongly cooperative sodium channels boosts neuronal encoding of high frequencies. PLoS One, 7, e37629. Kole, M., & Stuart, G. (2008). Action potential generation requires a high sodium channel density in the axon initial segment. Nature Neuroscience, 11(11), 1253–1255. Kole, M., Ilschner, S., Kampa, B., Williams, S., Ruben, P., & Stuart, G. (2008). Is action potential threshold lowest in the axon? Nature Neuroscience, 11, 178–186. K¨ondgen, H., Geisler, C., Fusi, S., Wang, X., L¨uscher, H., & Giugliano, M. (2008). The dynamical response properties of neocortical neurons to temproally modulated noisy inputs in vitro. Cerebral Cortex, 18(9), 2086–2097. Mainen, Z., Joerges, J., Huguenard, J., & Sejnowski, T.A. (1995). model of spike initiation in neocortical pyramidal neurons. Neuron, 15(6), 1427–1439. Marx, S., Ondrias, K., & Marks, A. (1998). Coupled gating between individual skeletal muscle ca2+ release channels (ryanodine receptors). Science, 281(5378), 818–821. Marx, S.O., Gaburjakova, J., Gaburjakova, M., Henrikson, C., Ondrias, K., & Marks, A.R. (2001). Coupled gating between cardiac calcium release channels (ryanodine receptors). Circulation Research, 88, 1151–1158. McCormick, D., Shu, Y., & Yu, Y. (2007). Neurophysiology: Hodgkin and huxley model - still standing?. Nature, 445, 1–2. Molina, M., Barrera, F., Fern`andez, A., Poveda, J., Renart, M., Encinar, J., Riquelme, G., & Gonz`alez-Ros, J. (2006). Clustering and coupled gating modulate the activity in kcsa, a potassium channel model. The Journal of Biological Chemistry, 281(27), 18, 837–18, 848. Naundorf, B., Geisel, T., & Wolf, F. (2005a). Action potential onset dynamics and the response speed of neuronal populations. Journal of Computational Neuroscience, 18, 297–309. Naundorf, B., Geisel, T., & Wolf, F. (2005b). Dynamical response properties of a canonical model for type-i membranes. Neurocomputing, 20, 421–428. Naundorf, B., Wolf, F., & Volgushev, M. (2006). Unique features of action potential initiation in cortical neurons. Nature, 440(7087), 1060–1063. Naundorf, B., Wolf, F., & Volgushev, M. (2007). reply to: hodgkinhuxley model – still standing?. Nature, 445, E2–E3. Neumcke, B., & St¨ampfli, R. (1983). Alteration of the conductance of na+ channels in the nodal membrane of frog nerve by holding potential and tetrodotoxin. Biochimica et Biophysica Acta, 727(1), 177–184.

J Comput Neurosci Palmer, L., & Stuart, G. (2006). Site of action potential initiation in layer 5 pyramidal neurons. The Journal of Neuroscience, 26, 1854–1863. Post, J., Leunissen-Bijvelt, J., Ruigrok, T., & Verkleij, A. (1985). Ultrastructural changes of sarcolemma and mitochondria in the isolated rabbit heart during ischemia and reperfusion. Biochimica et Biophysica Acta, 845(1), 119–123. Saito, A., Inui, M., Radermacher, M., Frank, J., & Fleischer, S. (1988). Ultrastructure of the calcium release channel of sarcoplasmic reticulum. The Journal of Cell Biology, 107(1), 211–219. Schafer, D., Jha, S., Liu, F., Akella, T., McCullogh, L., & Rasband, M. (2009). Disruption of axon initial segment cytoskeleton is a new mechanism for neuronal injury. The Journal of Neuroscience, 29(42), 13, 242–13, 254. Shu, Y., Duque, A., Yu, Y., Haider, B., & McCormick, D. (2007). Properties of action-potential initiation in neocortical pyramidal cells: evidence from whole cell axon recordings. Journal of Neurophysiology, 97(1), 746–760. Silberberg, S., & Magleby, K. (1997). Voltage-induced slow activation and deactivation of mechanosensitive channels in xenopus oocytes. The Journal of Physiology, 505(Pt 3), 551–569. Sonkusare, S.K., Dalsgaard, Bonev, A.D., Hill-Eubanks, D.C., Kotlikoff, M.I., Scott, J.D., Santana, L.F., & Nelson, M.T. (2014). Akap150-dependent cooperative trpv4 channel gating is central to endothelium-dependent vasodilation and is disrupted in hypertension. Science Signaling, 7(333).

Stuart, G., & Schiller, J. (1997). Sakmann B Action potential initiation and propagation in rat neocortical pyramidal neurons. The Journal of Physiology, 505, 617–632. Undrovinas, A., & Makielski, I.F.J. (1992). Inward sodium current at resting potentials in single cardiac myocytes induced by the ischemic metabolite lysophosphatidylcholine. Circulation Research, 71(5), 1231–1241. Ursell, T., Huang, K., Peterson, E., & Phillips, R. (2007). Cooperative gating and spatial organization of membrane proteins through elastic interactions. PLoS Computational Biology, 3(5), e81. Volgushev, M., Malyshev, A., Balaban, P., Chistiakova, M., & Volgushev, S. (2008). Wolf F Onset dynamics of action potentials in rat neocortical neurons and identified snail neurons : quantification of the difference. PLoS ONE, e1962, 3. Wang, X., & Buzs`aki, G. (1996). Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. The Journal of Neuroscience, 16, 6402–6413. Wei, W., & Wolf, F. (2011). Spike onset dynamics and response speed in neuronal populations. Physical Review Letters, 106(8), 088, 102. Wiesenfeld, K., & Moss, F. (1995). Stochastic resonance and the benefits of noise: from ice ages to crayfish and squids. Nature, 73(6509), 33–36. Yu, Y., Shu, Y., & McCormick D. (2008). Cortical action potential backpropogation explains spike threshold variability and rapidonset kinetics. The Journal of Neuroscience, 28, 7260–7272.