Simple Cortical and Thalamic Neuron Models for Digital Arithmetic Circuit Implementation

Trade-off between reproducibility of neuronal activities and computational efficiency is one of crucial subjects in computational neuroscience and neuromorphic engineering. A wide variety of neuronal models have been studied from different viewpoints. The digital spiking silicon neuron (DSSN) model is a qualitative model that focuses on efficient implementation by digital arithmetic circuits. We expanded the DSSN model and found appropriate parameter sets with which it reproduces the dynamical behaviors of the ionic-conductance models of four classes of cortical and thalamic neurons. We first developed a four-variable model by reducing the number of variables in the ionic-conductance models and elucidated its mathematical structures using bifurcation analysis. Then, expanded DSSN models were constructed that reproduce these mathematical structures and capture the characteristic behavior of each neuron class. We confirmed that statistics of the neuronal spike sequences are similar in the DSSN and the ionic-conductance models. Computational cost of the DSSN model is larger than that of the recent sophisticated Integrate-and-Fire-based models, but smaller than the ionic-conductance models. This model is intended to provide another meeting point for above trade-off that satisfies the demand for large-scale neuronal network simulation with closer-to-biology models.

1. Introduction

Silicon neuronal networks have gained remarkable attention in recent years. The silicon neuronal network is composed of dedicated circuit that solves the differential equations of a neuron and synapse model. On account of their parallel and distributed structures, silicon neuronal networks can simulate neuronal activities with low power consumption and in high speed, potentially realizing an extremely large-scale network comparable to that of the human brain in future.

Their analog circuit implementation consumes ultra-low power down to several nano watts per silicon neuron (Brink et al., 2013; Kohno and Aihara, 2014; Mandloi et al., 2014), however, it includes technical hurdles of fabrication mismatch and temperature dependence to construct a large-scale network. On the other hand, digital circuit implementation solves this limitation because it is far less sensitive to these factors, though power consumption tends to be higher than the analog circuit implementations. A digital silicon neuronal network comprising 1 million spiking neurons and 256 million synapses has been implemented on a 5.4-billion transistor chip (Merolla et al., 2014). This Application Specific Integrated Circuit (ASIC) chip calculates asynchronously in real-time and consumes just 63 mW. Merolla et al. adopted the leaky integrate-and-fire (LIF) model, one of the most simple neuronal models. The LIF model describes the dynamics of the neuronal membrane potential, which is perturbed by a stimulus inputs over time and converges slowly to the resting potential. It is computationally efficient and suitable for the large-scale implementation. After refining the synaptic efficacy by a learning process, the ASIC was applied to a multi-object detection and classification.

In silicon neuronal networks, diversified neuronal models are used due to a trade-off between the reproducibility of neuronal activity and computational efficiency. For instance, ionic-conductance models can reproduce a neuronal activity accurately but demands excessive computational resources in large-scale implementations. In contrast, neuronal models that approximate a spiking process by resetting the state variable [integrate-and-fire (I&F) -based models], such as the LIF, exponential I&F (Fourcaud-Trocmé et al., 2003), adaptive exponential I&F (Brette and Gerstner, 2005), and Izhikevich (IZH) models, can be implemented at low computational cost. However, it has reduced reproducibility of complex neuronal activities. For example, these models assume fixed maximum membrane potentials during the spike process, whereas these potentials are nonuniform in the nervous system (Alle and Geiger, 2006). Reduction of computational cost benefits the power consumption, system size, and response speed of silicon neuronal networks regardless of implementation types (ASIC, Field Programmable Gate Array, and massively-parallel CPUs). These factors are particularly important for their application to neuromorphic systems (e.g., for robot control) that are required to occupy a small space and operate with restricted power supply. They are also important to implement with realistic power consumption and system size a large-scale network comparable to the human brain that is composed of about a hundred billion neurons. The DSSN model (Kohno and Aihara, 2007) is a qualitative neuronal model designed for efficient implementation in a digital arithmetic circuit. It is a non-I&F-based model that can realize several neuronal activities including the Class I and II in the Hodgkin's classification (Hodgkin, 1948). Because this model does not abbreviate the calculation of the spiking process, it can reproduce the gradient response in Class II neurons (Wang and Rinzel, 2003). Li et al. (2012) constructed an auto-associative memory with 256 fully connected digital spiking silicon neuron (DSSN) models on an FPGA. They reported that a network of Class II neurons yields higher retrieval performance than that of Class I neurons in the associative memory task. Recently, a Hebbian learning rule was applied to this network (Li et al., 2013). In Osawa and Kohno (2015), it was reported that the Class II mode of the IZH model with a standard parameter setting has discontinuous phase resetting curve (PRC) and the auto-associative memory constructed in the same way does not deliver higher performance than the Class I-mode IZH model. Increasing the value of a parameter can solve this problem, however, it considerably distorts the spiking waveform.

Pospischil et al. (2008) have found appropriate parameter sets for an ionic-conductance model that replicate experimental data in the well-known four classes of cortical and thalamic neurons; regular spiking (RS), fast spiking (FS), intrinsically bursting (IB), and low-threshold spike (LTS). A most typical cortical neuron class, RS is characterized by spike-frequency adaptation; that is, the spike frequency decreases over time in response to a constant stimulus input. Conversely, FS neurons maintain firing at a constant frequency. IB and LTS are the neuron classes with bursting ability. Neurons in the IB class generate a burst firing immediately at the onset of a stimulus, then continue spiking until its termination. The LTS class also exhibits the bursting. Neurons in the LTS class generate a burst firing just after the termination of a sufficient hyperpolarizing stimulus. Their model equations were constructed by integrating several ionic-conductance models reported in the previous researches and contains up to seven variables. For the RS and FS classes, they applied an auto fitting procedure that minimizes the error function of spike intervals by the simulated annealing method. The IZH model supports these four classes. The Pospischils' ionic-conductance model for the IB class generates a long silent phase following the burst firing (Figure 11-left), but the IZH model produced uniform silent phases (Figure 1). Moreover, the response to larger inputs evidently differed from that of Pospischils' model. Because it is not elucidated completely what properties of the neuronal activities are playing the key roles in the information processing in the brain, developing a simple non-I&F model that supports a wide variety of neuronal activities can be significant for both scientific and engineering purposes, even if it consumes computational resources than I&F-based models.

FIGURE 1

Figure 1. Behaviors of the Izhikevich model in the IB class. Stimulus step input I_stim rise at t = 0.02, where (A) I_stim = 10, (B) 15, and (C) 30, respectively. A numerical integration step dt is 0.05 ms.

In Nanami and Kohno (2015), we reduced the dimensions of Pospischils models by Keplers method and analyzed their mathematical structures by bifurcation analysis. We also reported parameter sets where the three-variable DSSN model (Kobayashi et al., 2011) produces similar activities to our target classes. While we did not evaluate their similarity, it was clear that the parameter set for the IB class could not reproduce the transient firing patterns from the first bursting phase to the later spiking phase. In this work, we cleared this limitation as follows. Firstly, we repeated the parameter search and found three-variable DSSN model parameter sets for RS, FS, and LTS modes that reproduce corresponding mathematical structures of the reduced Pospischils models and statistical characteristics of spiking patterns. Secondly, for the IB class, we incorporated an additional slow variable and found a parameter set that reproduces corresponding mathematical structures and the statistical characteristics of spiking patterns. The statistical evaluation was performed using Cv and Lv (Shinomoto et al., 2003) that are utilized to characterize the spiking activities of neuronal cells.

The remainder of this paper is organized as follows. Section 2 introduces our neuron model, its basic concepts, and details of its construction. The model is tested and evaluated by simulation in Section 3. Section 4 summarizes the work and suggests ideas for future.

2. Methods

In the nervous system, the neuronal cells are classified by their intrinsic spike properties (Llinas, 1988; Connors and Gutnick, 1990; Gupta et al., 2000). The quantitative characteristics of neuronal activity, including the degree of adaptation and duration of the burst firing, differ even among neurons of the same class. To accommodate this fact, Pospischil et al. assembled various experimental data of different nerve cells in each neuron class, and adjusted the parameters of their ionic-conductance models' equations to each dataset. We expanded the DSSN model so that it has a compatible dynamics and behaviors to a typical parameter of the Pospischil's model in each class.

2.1. DSSN Model

The DSSN model (Kohno and Aihara, 2007) is a qualitative neuron model that can simulate several classes of neuronal activities by Euler's method with fixed point operation. Given appropriate parameter sets, it can reproduce the Class I and II in the Hodgkin's classification as well as Class I* (Fujii and Tsuda 2004) which is defined by a unique mathematical structure named narrow channel. The differential equations of the model is written as follows:

$\begin{array}{cc}\frac{dv}{dt}=\frac{\varphi }{\tau }(f(v)-n+I_0+I_{\text{stim}}),Ko& \left(1\right)\end{array}$

$\begin{array}{cc}\frac{dn}{dt}=\frac{1}{\tau }(g(v)-n),& \left(2\right)\end{array}$

$\begin{array}{cc}f(v)=\{\begin{array}{l}{a}_{fn}{(v-b_{fn})}^2+c_{fn}(v<0)\\ a_{fp}{(v-b_{fp})}^2+c_{fp}(v\ge 0),\end{array}& \left(3\right)\end{array}$

$\begin{array}{cc}g(v)=\{\begin{array}{l}{a}_{gn}{(v-b_{gn})}^2+c_{gn}(v<r_g)\\ a_{gp}{(v-b_{gp})}^2+c_{gp}(v\ge r_g),\end{array}& \left(4\right)\end{array}$

And Kobayashi et al. derived three-variable DSSN model (Kobayashi et al., 2011) that realize autonomous bursting including the square-wave bursting and the elliptic bursting. As explained in the introduction, we reported parameter sets with which this model can produce firing patterns similar to the four classes of cortical and thalamic neuronal activities in Pospischil et al. (2008) but with limited affinity. The equations of this model are given by

$\begin{array}{cc}\frac{dv}{dt}=\frac{\varphi }{\tau }(f(v)-n-q+I_0+I_{\text{stim}}),& \left(5\right)\end{array}$

$\begin{array}{cc}\frac{dn}{dt}=\frac{1}{\tau }(g(v)-n),& \left(6\right)\end{array}$

$\begin{array}{cc}\frac{dq}{dt}=\frac{ϵ}{\tau }(h(v)-q),& \left(7\right)\end{array}$

$\begin{array}{cc}h(v)\text{=}\{\begin{array}{l}{a}_{hn}{(v-b_{hn})}^2+c_{hn}(v<r_h)\\ a_{hp}{(v-b_{hp})}^2+c_{hp}(v\ge r_h),\end{array}& \left(8\right)\end{array}$

where v corresponds to the membrane potential, and n and q are the fast and slow variables, respectively, that abstractly describe the activity of the ion channels. The Parameter I₀ is a bias constant and I_stim represents the input stimulus. Parameters ϕ, ϵ, and τ control the time constants of the variables. Parameters r_x, a_x, b_x, and c_x, where x = fn, fp, gn, gp, hn, or hp, are constants that adjust the nullclines of the variables. All of the variables and constants in this qualitative model are purely abstract with no physical units. Most existing qualitative neuronal models replicate the spiking dynamics by a cubed variable term (FitzHugh, 1961; Nagumo et al., 1962; Hindmarsh and Rose, 1984). Because multiplication consumes significant circuit resources in a digital arithmetic circuit, the DSSN model adopts a piecewise quadratic function so that its numerical integration step includes only one multiplication between variables.

2.2. Reduced Model

The Pospischils' models contain up to seven variables, which makes an analysis of their mathematical structure difficult. Therefore, by applying Keplers' method (Kepler et al., 1992), we derived three-variable reduced models for the RS, FS, and LTS classes, and a four-variable reduced model for the IB class. Their equations and parameter sets are listed in the Appendix of Supplementary Material. This process reduces the number of variables without critically altering the models' intrinsic structure. The reduced model of the IB class has four variables, because the time-scale differences among the variables prevent the division of the variables into three groups without altering the neuronal activity. Comparing the waveforms of the reduced and Pospischils' models, we confirmed that the dimensional reduction imposed no significant changes on their dynamical properties (Figure 2). In both models, the RS classes have the spike-frequency adaptation (Figures 2A,B). The FS classes have little adaptation. The IB classes have burst firing followed by repetitive single spikes (Figure 2C). The LTS classes generate rebound bursting after removing the long hyperpolarizing stimuli (Figure 2D).

FIGURE 2

Figure 2. Waveforms of the Pospischils' and reduced models in the cortical and thalamic neuron classes. (A) The excitatory RS class, (B) the inhibitory RS class, (C) the FS class, (D) the IB class, and (E) the LTS class, respectively. An excitatory step stimulus is applied in (A–D). An inhibitory pulse stimulus that starts at 180 ms and ends at 980 ms is applied in (E).

2.3. Analysis of Mathematical Structure

This subsection analyzes mathematical structures of the reduced models derived in the preceding section. The bifurcation diagrams of the reduced models of the excitatory RS, inhibitory RS, FS, and LTS classes are displayed on a V – V_y plane with a V_y–nullcline (Figure 3). The bifurcation analysis was conducted on the fast subsystem including V and V_x while varying the remaining slow variable as a bifurcation parameter. Their mathematical structures are projected onto the V − V_y planes. The fast subsystem of these neuron classes undergo a saddle- node on invariant circle bifurcation. Limit cycles which correspond to the periodic spiking are represented by two curves that correspond to their maximum and minimum values of V. The V_y–nullcline represents the continuous points where $\frac{d{V}_y}{dt}$ is 0. The state point on the stable limit cycle slows down as it approaches to the saddle-node point. Therefore, the slow variable V_y can control the firing frequency, and the variations in the dynamics of V_y induce the different properties of the neuron classes.

FIGURE 3

Figure 3. Mathematical structures of the reduced models. Projected mathematical structures onto the V_y − V plane. A stimulus sufficiently strong for periodic spiking is applied. The saddle-node bifurcation is observed in (A) the excitatory RS, (B) the inhibitory RS, (C) the FS, and (D) the LTS classes.

In our analysis of the IB class, we fixed the slowest variable V_z while varying the second slowest variable V_y as the parameter because V_z is sufficiently slower than V_y. The transitions in the mathematical structure of the fast subsystem induced by the change in V_z are presented in Figure 4. The fast subsystem exhibits a homoclinic loop bifurcation and the bistable range of V_y depends on the slowest variable V_z. Here, the bistable region is where the stable limit cycle overlaps the stable node with respect to the V_y-axis. In response to an appropriate step input, the system produces three qualitatively different phases (Figure 4A). They are burst firing, subsequent long silent, and monotonic single spiking phase. This phase transition proceeds as the slowest variable V_z increases (Figure 4B). When V_z is relatively small, a bursting with many spikes is followed by a long silent phase, because the system has a wide bistability region where the state point takes long time to go back to stable limit cycle (Figure 4C). When V_z gets larger, this region is small and the state point can generate only one spike while it is on the limit cycle (Figure 4D).

FIGURE 4

Figure 4. Mathematical structure of the reduced model in the IB class. (A) Periodic spike waveform of the membrane potential in response to an excitatory step input. (B) Transition of the slowest variable V_z that controls the mathematical structure of the fast system. The projected mathematical structure for the (C) V_z = –70 mV and (D) V_z = −58 mV.

2.4. Parameter Tuning and the Four-Variable DSSN Model

To support the RS, FS, and LTS classes by the DSSN model, we found appropriate parameter sets for it that produce the same structures as those in the above reduced models. Firstly, we determined appropriate parameters of the fast subsystem. They were selected so that the system undergoes a saddle-node bifurcation as observed in the reduced models for the RS, FS, and LTS classes (Figure 5). The parameters of the fast subsystem for the LTS mode was tuned to mimic the hyperpolarizing response in the reduced model. Secondly, we tuned the remaining parameters, which control the dynamics of the slow system, to precisely reproduce the distinctive neuronal activities of the corresponding neuron classes. In the RS mode, we placed the q-nullcline at relatively low area of the V_y plane so that the system state converges to the limit cycle at near the bifurcation point. It induces spike frequency adaptation because the firing frequency decreases as the slow variable q approaches the saddle-node bifurcation point. In the FS mode, the q-nullcline was adjusted in the inverse direction to suppress the adaptation. For the LTS mode, the parameters for the fast subsystem were selected to mimic the resting membrane potential in the Pospischils' model which is higher than that in the other modes. Parameters of the slow variable q were selected so that it decreases under a hyperpolarizing input. Just after the injection of the hyperpolarizing input has been terminated, because the state point is at far from the saddle-node point, it rapidly moves along the stable limit cycle. This dynamics produces the rebound bursting.

FIGURE 5

Figure 5. Mathematical structures of the DSSN models. Projected mathematical structures onto the q-v plane. A stimulus sufficiently strong for periodic spiking is applied. A saddle-node on invariant circle bifurcation is seen in (A) the excitatory RS, (B) inhibitory RS, (C) FS, and (D) LTS.

To support the IB class by the DSSN model, we added a slowest variable in order to vary the bistable area, which was explained in the last part of the previous subsection. The four-variable DSSN model is written as follows:

$\begin{array}{cc}\frac{dv}{dt}=\frac{\varphi (u)}{\tau }(f(v)-n-q+I_0+I_{\text{stim}})& \left(9\right)\end{array}$

$\begin{array}{cc}\frac{dn}{dt}=\frac{1}{\tau }(g(v)-n)& \left(10\right)\end{array}$

$\begin{array}{cc}\frac{dq}{dt}=\frac{ϵ}{\tau }(h(v)-q)& \left(11\right)\end{array}$

$\begin{array}{cc}\frac{du}{dt}=\frac{{ϵ}_u}{\tau }(v-v_0-\alpha u)& \left(12\right)\end{array}$

$\begin{array}{cc}\varphi (u)=\{\begin{array}{l}{\varphi }_0(v<r_{u0})\\ {\varphi }_1(r_{u0}\le v<r_{u1})\\ {\varphi }_2(r_{u1}<v)\end{array}& \left(13\right)\end{array}$

where, the new variable u dynamically tunes the value of ϕ. Because the equation of u does not include multiplication between variables, and ϕ is selected in the list of possible values (ϕ₀, ϕ₁, and ϕ₂), it does not increase the number of multiplications in a numerical integration step. Its mathematical structure is illustrated in Figure 6. With no stimulus, the state point resides at a point which corresponds to the resting state (Figure 6, I). This point is at the intersection of the q-nullcline and the stable node. In response to a positive sustained input, the state point jumps into the stable limit cycle and keep firing repetitively (Figure 6, II). The state point drifts rightward because most part of the limit cycle is at above the q-nullcline. Therefore, the state point arrives at the right endpoint of the stable limit cycle and falls to the stable node. Subsequently, the state point moves leftward along the stable node because it is at below the q-nullcline, which corresponds to the long silent phase (Figure 6, III). Meanwhile, the slowest variable u exceeds r_u0, and the stable limit cycle shifts leftward. Therefore, the state point goes into the limit cycle again and makes a few spikes along the limit cycle before returning to the stable node (Figure 6, IV). When u exceeds r_u1, a single spike is generated on the limit cycle before the state point goes back to the stable node (Figure 6, V). The DSSN model qualitatively reproduces the behaviors of the four neuron classes by numerical integration with the Euler's method (dt = 0.0001 s).

FIGURE 6

Figure 6. Mathematical structure of the DSSN model in the IB class. (A) Waveform of the membrane potential in response to a step input that rises at t = 0.2 s. (B) Transition of the complementary variable u. (C) Projected mathematical structure without any input stimulus. Projected mathematical structure under a sustained stimulus where (D) u < r_u0, (E) r_u0 < u < r_u1, and (F) r_u1 < u. Each roman number corresponds to a spiking phase.

The DSSN model qualitatively reproduces the behaviors of the four neuron classes by numerical integration with the Euler's method (dt = 0.0001 s). For example, the calculation of v for IB classes is conducted as follows;

$\begin{array}{cc}{v}_{\text{next}}=v\_vv+v\_v+v\_n+v\_q+v\_I+v\_c,& \left(14\right)\end{array}$

$\begin{array}{cc}v\_v{v}_{ij}=v^2·(\Delta t·a_{fj}·{\varphi }_i/\tau ),& \left(15\right)\end{array}$

$\begin{array}{cc}v\_v_{ij}=v·(\Delta t·(-2)·b_{fj}·a_{fj}·{\varphi }_i/\tau +1),& \left(16\right)\end{array}$

$\begin{array}{cc}v\_n_{ij}=n·(-\Delta t·{\varphi }_i/\tau ),& \left(17\right)\end{array}$

$\begin{array}{cc}v\_q_{ij}=q·(-\Delta t·{\varphi }_i/\tau ),& \left(18\right)\end{array}$

$\begin{array}{cc}v\_I_{ij}=I_{\text{stim}}·(\Delta t·{\varphi }_i/\tau ),& \left(19\right)\end{array}$

$\begin{array}{cc}v\_c_{ij}=(b_{fj}^2·a_{fj}+c_{fj}+I_0)·\Delta t·{\varphi }_i/\tau ,& \left(20\right)\end{array}$

$\begin{array}{cc}v\_x=\{\begin{array}{l}v\_x_{0n}\text{when}v<r_{u0}\text{and}v<0\\ v\_x_{1n}\text{when}{r}_{u0}\le v<r_{u1}\text{and}v<0\\ v\_x_{2n}\text{when}{r}_{u1}\le v\text{and}v<0\\ v\_x_{0p}\text{when}v<r_{u0}\text{and}v\ge 0\\ v\_x_{1p}\text{when}{r}_{u0}\le v<r_{u1}\text{and}v\ge 0\\ v\_x_{2p}\text{when}{r}_{u1}\le v\text{and}v\ge 0,\end{array}& \left(21\right)\end{array}$

$\text{for}x=vv,v,n,q,I,\text{and}c.$

where, i denotes 0, 1, or 2 and j denotes n or p. In these equations, we can calculate a product of parameters in advance and store it as a constant value. Parameters are a sum of at most 7 (in average 3.78) positive or negative power-of-two numbers so that the multiplication of a variable and a parameter is realized by combination of bit shift and add operations. Multiplication is only used for the calculation of v².

3. Results

3.1. Comparison of the Waveforms

We compared the waveforms generated by the DSSN and Pospischils' ionic-conductance models for several magnitudes of step inputs. The DSSN model accurately reproduce Pospischils' model in each neuron class (Figures 7–11). In the excitatory RS mode, both models exhibit the spike-frequency adaptation, and their firing frequencies increase as input stimuli become large (Figure 7). Figure 8 shows the inhibitory RS mode that has more stronger spike-frequency adaptation than the excitatory RS mode. Both models in the FS mode fire with almost constant frequency, which increases as the stimulus is strengthened. In the LTS mode (Figure 10), both models generate periodic spiking with spike-frequency adaptation in response to positive step stimulus and a rebound bursting at the end of a negative stimulus. Both models in the IB mode exhibit repetitive spiking following a burst firing in response to a weak positive input (Figure 11). The frequency of the repetitive spiking gets higher in response to a stronger positive stimulus.

FIGURE 7

Figure 7. Waveforms of Pospischil and DSSN models, corresponding to an excitatory neuron in the RS class. The Pospischils' model (left) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 1.1, (B) 1.5, (C) 2.1 nA. The DSSN model (right) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 2.9221, (B) 3.2521, (C) 3.7921.

FIGURE 8

Figure 8. Waveforms of Pospischil and DSSN models, corresponding to an inhibitory neuron in the RS class. The Pospischils' model (left) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 0.6, (B) 1.1, (C) 1.7 nA. The DSSN model (right) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 2.97, (B) 3.57, (C) 4.29.

FIGURE 9

Figure 9. Waveforms of Pospischil and DSSN models, corresponding to a neuron in the FS class. The Pospischils' model (left) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 0.78, (B) 0.98, (C) 1.22 nA. The DSSN model (right) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 2.7, (B) 2.95, (C) 3.25.

FIGURE 10

Figure 10. Waveforms of Pospischil and DSSN models, corresponding to a neuron in the LTS class. The Pospischils' model (left) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 0.34, (B) 0.79, (C) –0.6 nA. The DSSN model (right) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 0.278, (B) 0.308, (C) –1.7. (A,B) Pospischils' model (left) and DSSN model (right) generate periodic spike in response to a depolarizing step input rising at t = 0.2. (C) A hyperpolarizing step input is applied from t = 0.2 to t = 1.0.

FIGURE 11

Figure 11. Waveforms of Pospischils' model and DSSN model which correspond to a neuron in the IB class. The Pospischils' model (left) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 2.38, (B) 2.6645, (C) 3.0495 nA. The DSSN' model (right) generates periodic spikes in response to a step input rising at t = 0.2 where I_stim is (A) 0.89, (B) 1.42, (C) 1.81.

3.2. C_V − L_V

To evaluate quantitatively our DSSN models, we adopted the coefficient of variation (C_V) and local variation (L_V) (Shinomoto et al., 2003), because the spike frequency cannot capture the spiking characteristics that changes over time. These statistics were developed by Shinomoto et al. to classify the spike sequences in the cortical areas of awake monkeys (Shinomoto et al., 2003), which have been used at various situations (Koyama and Shinomoto, 2005; Miura et al., 2006; Shinomoto et al., 2009). They are defined as follows:

$\begin{array}{cc}{C}_V=\sqrt{\frac{1}{n-1}{\displaystyle \sum _{i=1}^n{(T_i-\overline{T})}^2}}/\overline{T},& \left(22\right)\end{array}$

$\begin{array}{cc}{L}_V=\frac{1}{n-1}{\displaystyle \sum _{i=1}^{n-1}\frac{3{(T_i-T_{i+1})}^2}{{(T_i+T_{i+1})}^2}},& \left(23\right)\end{array}$

where T_i denotes the ith interspike interval, $\overline{T}$ is the average of T_i, and n is the number of spikes in the sequence. The coefficient 3 in Equation (21) is chosen so that the expectation value of L_V in the Poisson spike sequence equals one. The C_V is the standard deviation divided by the mean spike interval, and L_V becomes large when the duration of the spike intervals locally fluctuates. Both of C_V and L_V are zero for a regular spike sequence with constant T_i. We measured the C_V and L_V of the spike sequences generated by Pospischils' and our DSSN models while varying the strength of the step input. As shown in the Figure 12, the results are compared by their plot on the C_V − L_V plane. The C_V − L_V characteristics of the Pospischils' and the DSSN models in each neuron class are plotted in Figures 13, 14. We confirmed that both neuron models have similar statistical properties in each neuron class.

FIGURE 12

Figure 12. Spike sequences and their representation in the C_V − L_V plane. The L_V vs. C_V for spike sequences generated by the Pospischil (red) and the DSSN (blue) models in the excitatory RS mode in response to various magnitudes of step input.

FIGURE 13

Figure 13. C_V − L_V characteristics in the RS, FS, and LTS classes. The C_V − L_V characteristics of the Pospischil (red) and the DSSN (blue) models in the excitatory and inhibitory RS, FS, and LTS classes.

FIGURE 14

Figure 14. C_V − L_V characteristics in the IB classes. The C_V − L_V characteristics of the Pospischil (red) and DSSN (blue) model.

4. Conclusion

In this paper, we expanded the DSSN model so that it can support the RS, FS, LTS, and IB classes, according to the following steps. Firstly, we reduced the number of variables of the Pospischils' ionic-conductance model utilizing the Keplers' method. The reduced models contained three variables for the RS, FS, and LTS classes and four variables for the IB class. We then elucidated their mathematical structures, including their structures on the phase portraits. The FS class is generically considered to exhibit a Hopf bifurcation; however, it yielded a saddle-node bifurcation in the Pospischils' model. Secondly, we determined the appropriate parameter sets that produce the same mathematical structures in the DSSN model as those in the reduced models. To support the IB class, the four-variable DSSN was developed by adding a slowest variable. It dynamically tunes the structure of the fast system. If ϕ₀, ϕ₁, and ϕ₂ are equal to ϕ, the four-variable DSSN model is compatible with the three-variable one. Finally, it was confirmed that the DSSN models behave very similarly to the Pospischils' models in response to several magnitudes of step input. The similarity was quantitatively confirmed by measuring the C_V and L_V statistics of the spike sequences, and confirming that both neuron models have the same statistical properties in each neuron class. The model is solved by Euler's method (Δt = 0.0001[s]), and is realized by only one multiplier that consumes significant digital circuit resources.

The DSSN model will consume more circuit resources than the IZH model because it has more variables, while the number of multiplication per a numerical integration step is the same. But it captures more aspects of the neuronal activity because it is not an I& F-based model. As is apparent from the equations, the implementation cost of the DSSN model is far lower than the ionic-conductance models. Circuit implementation of this model will be reported in our future publication. In software simulation, the computational cost of the DSSN mode is much lower than that of the ionic-conductance models, and comparable with that of the IZH per a numerical integration step. We conducted simulation of two million spikes, the DSSN and IZH models in the RS mode consumed 18.51 and 2.76 s, respectively. The DSSN model required over five times longer calculation, because the step of the DSSN model (0.0001[s]) is five times smaller than that of the IZH model (0.0005[s]). If the fixed point operation is used, it will be simulated faster. The DSSN model may be also useful for computer simulation in which more realistic model than I&F-based models are required.

By the expansion in this work, the DSSN model supports the RS, FS, LTS, and IB classes as well as the Class I and II in the Hodgkin's classification, Class I*, square-wave bursting, elliptic bursting. We expect it can be a basis for silicon- and software- based spiking neuronal networks that capture the dynamics in the nervous system approximately, which can contribute to the neuroscience from the viewpoint of analysis by construction and the neuromorphic engineering. In this work, we tuned the parameters manually. We will work on an auto-fitting method for the DSSN model with an error function based on the mathematical structure.

Author Contributions

TN performed the study and TK supervised it.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by JSPS Grant-in-Aid for scientific Exploratory Research Number 25240045.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnins.2016.00181

References

Alle, H., and Geiger, J. R. P. (2006). Combined analog and action potential coding in hippocampal mossy fibers. Science 311, 1290–1293. doi: 10.1126/science.1119055

PubMed Abstract | CrossRef Full Text | Google Scholar

Brette, R., and Gerstner, W. (2005). Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J. Neurophysiol. 94, 3637–3642. doi: 10.1152/jn.00686.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Brink, S., Nease, S., and Hasler, P. E. (2013). Computing with networks of spiking neurons on a biophysically motivated floating-gate based neuromorphic integrated circuit. Neural Netw. 45, 39–49. doi: 10.1016/j.neunet.2013.02.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Connors, B. W., and Gutnick, M. J. (1990). Intrinsic firing patterns of diverse neocortical neurons. Trends Neurosci. 13, 99–104.

PubMed Abstract | Google Scholar

FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1, 445–466.

PubMed Abstract | Google Scholar

Fourcaud-Trocmé, N., Hansel, D., van Vreeswijk, C., and Brunel, N. (2003). How spike generation mechanisms determine the neuronal response to fluctuating inputs. J. Neurosci. 23, 11628–11640.

PubMed Abstract | Google Scholar

Fujii, H., and Tsuda, I. (2004). “Itinerant dynamics of class I* neurons coupled by gap junctions,” in Computational Neuroscience: Cortical Dynamics, Vol. 3146, eds P. Érdi, A. Esposito, M. Marinaro, and S. Scarpetta (Berlin; Heidelberg: Springer), 140–160.

Google Scholar

Gupta, A., Wang, Y., and Markram, H. (2000). Organizing principles for a diversity of gabaergic interneurons and synapses in the neocortex. Science 287, 273. doi: 10.1126/science.287.5451.273

PubMed Abstract | CrossRef Full Text | Google Scholar

Hindmarsh, J. L., and Rose, R. M. (1984). A model of neuronal bursting using tree coupled first order differential equations. Philos. Trans. R. Soc. Lond. B221, 87–102.

Google Scholar

Hodgkin, A. L. (1948). The local electric changes associated with repetitive action in a non-medullated axon. J. Physiol. 107, 165–181.

PubMed Abstract | Google Scholar

Kepler, T. B., Abbott, L. F., and Marder, E. (1992). Reduction of conductance-based neuron models. Biol. Cybern. 66, 381–387.

PubMed Abstract | Google Scholar

Kobayashi, W., Kohno, T., and Aihara, K. (2011). “3-variable digital spiking silicon neuron,” in Proceedings of the 24th Workshop on Circuit and Systems (in Japanese) (Hyogo), 1–5.

Kohno, T., and Aihara, K. (2007). “Digital spiking silicon neuron: concept and behaviors in GJ-coupled network,” in Proceedings of International Symposium on Artificial Life and Robotics (Beppu).

Kohno, T., and Aihara, K. (2014). “A qualitative-modeling-based low-power silicon nerve membrane,” in Electronics, Circuits and Systems (ICECS), 2014 21st IEEE International Conference on (Kitakyushu), 199–202.

Google Scholar

Koyama, S., and Shinomoto, S. (2005). Empirical bayes interpretations of random point events. J. Phys. A Math. Gen. 38, L531–L537. doi: 10.1088/0305-4470/38/29/l04

CrossRef Full Text | Google Scholar

Li, J., Katori, Y., and Kohno, T. (2012). An FPGA-based silicon neuronal network with selectable excitability silicon neurons. Front. Neurosci. 6:183. doi: 10.3389/fnins.2012.00183

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Katori, Y., and Kohno, T. (2013). “Hebbian learning in FPGA silicon neuronal network,” in The 1st IEEE/IIAE International Conference on Intelligent Systems and Image Processing (Kitakyushu).

Google Scholar

Llinas, R. R. (1988). The intrinsic electrophysiological properties of mammalian neurons: insights into central nervous system function. Science 23, 1654–1664.

Google Scholar

Mandloi, N., Bartolozzi, C., and Indiveri, G. (2014). “Smart motion sensing for autonomous robots,” in Biomedical Circuits and Systems Conference (BioCAS) (Atlanta, GA: IEEE), 520–523.

Merolla, P., Arthur, J., Alvarez-Icaza, R., Cassidy, A., Sawada, J., Akopyan, F., et al. (2014). A million spiking-neuron integrated circuit with a scalable communication network and interface. Science 345, 668–673. doi: 10.1126/science.1254642

PubMed Abstract | CrossRef Full Text | Google Scholar

Miura, K., Okada, M., and Amari, S.-I. (2006). Estimating spiking irregularities under changing environments. Neural Comput. 18, 2359–2386. doi: 10.1162/neco.2006.18.10.2359

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagumo, J. S., Arimoto, S., and Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proc. IRE 50, 2061–2071.

Google Scholar

Nanami, T., and Kohno, T. (2015). Cortex and thalamic neuron models for digital arithmetic circuit. Seisan Kenkyu 67, 297–300. doi: 10.11188/seisankenkyu.67.297

CrossRef Full Text

Osawa, Y., and Kohno, T. (2015). “Associative memory with class I and II Izhikevich model,” in The 2015 International Conference on Artificial Life and Robotics (ICAROB 2015) (Oita).

Google Scholar

Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Fregnac, Y., et al. (2008). Minimal hodgkin-huxley type models for different classes of cortical and thalamic neurons. Biol. Cybern. 99, 427–441. doi: 10.1007/s00422-008-0263-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Shinomoto, S., Kim, H., Shimokawa, T., Matsuno, N., Funahashi, S., Shima, K., et al. (2009). Relating neuronal firing patterns to functional differentiation of cerebral cortex. PLoS Comput. Biol. 5:e1000433. doi: 10.1371/journal.pcbi.1000433

PubMed Abstract | CrossRef Full Text | Google Scholar

Shinomoto, S., Shima, K., and Tanji, J. (2003). Differences in spiking patterns among cortical neurons. Neural Comput. 15, 2823–2842. doi: 10.1162/089976603322518759

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X.-J., and Rinzel, J. (2003). “Oscillatory and bursting properties of neurons,” in The Handbook of Brain Theory and Neural Networks, 2nd Edn. ed M. A. Arbib (Cambridge, MA: MIT Press), 835–840.