Generalized reconfigurable memristive dynamical system (MDS) for neuromorphic applications

ORIGINAL RESEARCH published: 03 November 2015 doi: 10.3389/fnins.2015.00409 Frontiers in Neuroscience | www.frontiersin.org 1November 2015 | Volume 9 | Article 409 Edited by: Julius Georgiou, University of Cyprus, Cyprus Reviewed by: Jongkil Park, University of California, San Diego, USA Stavros G. Stavrinides, University of Thessaly, Greece *Correspondence: Bernabé Linares-Barranco [email protected] Specialty section: This article was submitted to Neuromorphic Engineering, a section of the journal Frontiers in Neuroscience Received: 22 May 2015 Accepted: 12 October 2015 Published: 03 November 2015 Citation: Bavandpour M, Soleimani H, Linares-Barranco B, Abbott D and Chua LO (2015) Generalized reconfigurable memristive dynamical system (MDS) for neuromorphic applications. Front. Neurosci. 9:409. doi: 10.3389/fnins.2015.00409 Generalized reconfigurable memristive dynamical system (MDS) for neuromorphic applications Mohammad Bavandpour1, Hamid Soleimani2, Bernabé Linares-Barranco3*, Derek Abbott4and Leon O. Chua5 1Department of Electrical and Computer Engineering, University of California, Santa Barbara, Santa Barbara, CA, USA, 2Department of Bioengineering of Imperial College, London, UK, 3Instituto de Microelectrónica de Sevilla, IMSE-CNM, Universidad de Sevilla and CSIC, Sevilla, Spain, 4School of Electrical and Electronic Engineering, The University of Adelaide, Adelaide, SA, Australia, 5Department of Electrical and Computer Engineering, University of California, Berkeley, Berkeley, CA, USA This study firstly presents (i) a novel general cellular mapping scheme for two dimensional neuromorphic dynamical systems such as bio-inspired neuron models, and (ii) an efficient mixed analog-digital circuit, which can be conveniently implemented on a hybrid memristor-crossbar/CMOS platform, for hardware implementation of the scheme. This approach employs 4nmemristors and no switch for implementing an n-cell system in comparison with 2n2memristors and 2nswitches of a Cellular Memristive Dynamical System (CMDS). Moreover, this approach allows for dynamical variables with both analog and one-hot digital values opening a wide range of choices for interconnections and networking schemes. Dynamical response analyses show that this circuit exhibits various responses based on the underlying bifurcation scenarios which determine the main characteristics of the neuromorphic dynamical systems. Due to high programmability of the circuit, it can be applied to a variety of learning systems, real-time applications, and analytically indescribable dynamical systems. We simulate the FitzHugh-Nagumo (FHN), Adaptive Exponential (AdEx) integrate and fire, and Izhikevich neuron models on our platform, and investigate the dynamical behaviors of these circuits as case studies. Moreover, error analysis shows that our approach is suitably accurate. We also develop a simple hardware prototype for experimental demonstration of our approach. Keywords: general cellular mapping, hybrid memristor-crossbar/CMOS platform, FitzHugh-Nagumo (FHN) neuron model, Adaptive Exponential (AdEx) integrate and fire neuron model, Izhikevich neuron model, dynamical behavior analysis 1. INTRODUCTION The human nervous system is an intriguing complex system capable of performing intricate tasks using an enormous number of neurons each connected via synapses to several thousand neighboring neurons. Mathematical modeling of biological neural elements such as neurons, synapses, and glial cells has been a long standing active research area (Hodgkin and Huxley, 1952; Fitzhugh, 1961; Abbott, 1999; Izhikevich, 2003; Brette and Gerstner, 2005). Recently, the hardware implementation of neural networks has been the subject of much research with three motivations: (1) exploitation of functional mechanisms of biological neural systems; (2) development of neural prosthetics (i.e., replacement of a biological neural system by an electronic circuit), and building Bavandpour et al. Generalized reconfigurable MDS an artificial brain as a long-term goal; and (3) development of artificial neural networks for use in engineering applications. Two main factors in this research area are accuracy (targeting biological plausibility) and computational complexity of the model and system (targeting efficiency in terms of the large scale simulation and implementation cost). In general, there is a tradeoff between model accuracy and its computational complexity. Different platforms have been applied for hardware realization of neuromorphic dynamical systems. Here, we investigate main approaches for this challenge: 1. Special purpose computing architectures have been developed to simulate neurobiological networks and functions using their specially designed software tools for large scale simulations (Ahmadi and Soleimani, 2011; Furber et al., 2014; Minkovich et al., 2014). Even though these systems are flexible and biologically realistic with considerably high performance due to their massively parallel architecture, the presented hardware approaches are based on bulky and power-hungry workstations with relatively high cost and development time. Hence, these approaches are often impractical for public access, and general purpose large-scale simulations. 2. Analog CMOS platform is considered to be a significant choice for direct implementation of neural dynamic functions (Linares-Barranco et al., 1991; Arthur and Boahen, 2011; Wijekoon and Dudek, 2012; Soleimani et al., 2014). However, this approach is not a general approach. This implies that the structure of the circuit must be modified if there are any changes to the starting equations. Moreover, non-linear functions in target mathematical models are directly implemented using inherent non-linearity of the circuit elements and networks. Hence, the variability and mismatch of the circuit elements drastically affect the circuit performance, and implementation of some models with special non-linear dynamics is cumbersome and intractable. Besides, this hardware platform is comparatively inflexible, and model adjustment is typically troublesome in these circuits. 3. A digital platform is used to realize bio-inspired neural cells. Most digital approaches (Weinstein et al., 2007; Soleimani et al., 2012, 2015; Cassidy et al., 2013) use digital computational units to implement the mathematical equations describing the behavior of biological neural cells. Another type of digital approach (Matsubara and Torikai, 2013) is presented based on Cellular Automata (CA) consisting logic gates, register arrays, and reconfigurable wires, and the dynamics of the target system are determined by the wire patterns. This platform can be further subcategorized into FPGA based (Weinstein et al., 2007; Soleimani et al., 2012, 2015; Matsubara and Torikai, 2013) and digital custom IC based (Cassidy et al., 2013) hardware, where FPGAs provide more configurability and lower development time, but at a higher cost in terms of area, power, and speed in comparison with the digital custom ICs. Generally, a digital platform achieves rapid development time, high reconfigurability, and immunity to device mismatch. However, its silicon area and power consumption is directly related to the mathematical complexity of the model, and it is not beneficial for complex models due to excessive utilization of high-cost computational units such as comparators, multipliers, and adders. Moreover a CA based approach (Matsubara and Torikai, 2013) uses a complicated wiring network that occupies a large switching area, decreases the maximum frequency of the overall system clock by making a long critical path, and results in a cumbersome implementation process. 4. Full custom analog/digital (mixed mode) implementations that comprise low-power, fast analog circuits and programmable, mismatch immune digital circuits. Generally, for this approach, neural computation is performed in the analog domain while the communication of spikes between nervous cells is carried out in the digital domain (Schemmel et al., 2008; Benjamin et al., 2014; Moradi and Indiveri, 2014). Recently using this approach, we have proposed a Cellular Memristive Dynamical System (CMDS; Bavandpour et al., 2014) consisting a dynamical systems cellular mapping, and a memristor-based hardware platform for implementing a wide range of neuromorphic systems. The mixed analog-digital nature of this system results in accuracy and efficiency in terms of hardware implementation cost. However, in this platform, the number of memristors employed in a twodimensional system has a quadratic relation with the number of cells where an n-cell system requires 2n2memristors. So, the area consumption is significantly increased with the number of cells. Besides, the circuit provides just one-hot digital values of the dynamical variables, which constrains the networking schemes, so an additional circuit is needed to provide the analog value of the variables. Our approach is a cellular-based system that discretizes dynamical variables resulting in a cellular phase plane, stores the equilibrium curves in a memristive crossbar-based analog memory block, evaluates the velocity and direction of the vector field in the cells, and tracks the state point in the space using VCOs and pointer registers. This system is a fully reconfigurable general approach capable of implementing a wide range of two dimensional neuromorphic dynamical systems such as FitzHugh-Nagumo (FHN; Fitzhugh, 1961), Adaptive Exponential (AdEx) integrate and fire (Brette and Gerstner, 2005), and Izhikevich neuron models (Izhikevich, 2003). One promising technology particularly suited for analog and mixed analog-digital computing is based on hybrid circuits that integrate CMOS and memristor devices (Strukov and Likharev, 2006; Bahar et al., 2007; Kavehei et al., 2012). The memristor is a two-terminal thin-film device whose resistance can be tuned in a nonvolatile and analog way (Chua, 1971; Strukov et al., 2008; Yang et al., 2008). In the context of analog and mixed circuit applications, recent advances in memristive devices and their integration with CMOS enable efficient implementations of nanoscale analog-grade resistive elements that can be fine-tuned after fabrication (Alibart et al., 2012). The proposed approach is an efficient generalized system with high reconfigurability due to using memristor-crossbar/CMOS hardware platform. This circuit is highly programmable because reconfiguration is carried out by adjusting memristance values, while the structure of the Frontiers in Neuroscience | www.frontiersin.org 2November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS system remains fixed. It implies that a wide range of target dynamical systems such as various neuron models can be easily programmed into the system and their parameters can be easily adjusted without changing the hardware structure. Furthermore, the circuit can accurately mimic different dynamic neuronal behaviors, and yield a variety of biological-like spike trains as well as single spike shapes. The implementation constraint has strongly limited the power of dynamical systems in modeling, and neuroscientists have been unable to develop an accurate model. Moreover, a number of behaviors such as special output signal shapes have no elegant analytical description. Our approach significantly alleviates the limitation of computational effort in dynamical functions. The rest of the paper is organized as follows: In Section 2, we discuss the proposed general mapping scheme. The memristorcrossbar/CMOS based hardware structure of the proposed platform is introduced in Section 3. Section 4 presents implementation of three widely used neuron models called the FHN neuron model (Fitzhugh, 1961), AdEx neuron model (Brette and Gerstner, 2005), and Izhikevich neuron model (Izhikevich, 2003) based on the proposed approach, and investigates various dynamical responses of the resultant circuits as the case studies. Section 5 introduces design steps and practical results of a simple hardware prototype of our approach. The error analysis of the proposed approach is studied in Section 6, and the networking and learning capability of the system is explained as two remarks in Section 7. Finally, the work is concluded in Section 8. 2. GENERAL CELLULAR MAPPING FOR NEUROMORPHIC DYNAMICAL SYSTEMS The first step toward developing a memristor-crossbar/CMOS hardware for implementing a reconfigurable neuromorphic dynamical system is drawing a proper mapping to represent target dynamical systems. This step plays a significant role in the whole design process because the proposed mapping ought to satisfy a number of crucial conditions: It ought to have the capability of (1) transforming the relatively complicated equations of neuromorphic dynamical systems into a number of simplified fully implementable equations, (2) obtaining a maximum accuracy considering the available hardware resources for implementation, (3) being applied to a wide range of neuromorphic dynamical systems, and (4) bringing reconfigurability into the hardware by separating the variable features of the neuromorphic dynamical systems (requisite information in phase plane for reconstructing vector field) from the shared features (computations for reproducing timedomain signals), and implementing the variable features on the reconfigurable memristive part of the circuit. In this section, we propose a cellular mapping that properly satisfies the aforementioned conditions. In Section 2.1, basic concepts of dynamical systems in neuroscience are reviewed, and a general form for two dimensional neuromorphic dynamical systems and its simplified forms are presented. One of the simplified forms, which covers most of the two dimensional neuron models, is the basis of our mapping in the next subsections. In Section 2.2, a two dimensional cellular space is defined and the mapping equations for transforming the target dynamical system from the continuous space to the cellular space is presented. Then, according to the simplified form, a technique for calculating the motion velocity of the state point in each cardinal direction in the cellular space is proposed. Applying this technique to the simplified general form in the cellular space results in a mapping which satisfies first, third and fourth out of four above-mentioned conditions. But the second condition related to the timing and asynchrony ought to be satisfied in order to maximize the accuracy of tracking state point in the cellular space. In Section 2.3, the concept of timing and asynchrony in our cellular approach is explained, the requisite timing condition is introduced, and the final satisfactory timing equations and cell change policy based on the timing parameters are presented. 2.1. Dynamical Systems in Neuroscience In computational modeling of neural cell behaviors and phenomena, neuroscientists aim to point out critical features and factors as the system variables and parameters, and draw out deterministic laws governing the evolution of these variables over time (Gerstner and Kistler, 2002). These laws are presented by mathematical functions, and the collection of resultant equations is called a dynamical system. Dynamical systems model the behavior of a given system without a thorough comprehensible knowledge on all the governing rules for its evolution. This feature significantly matches the nature of neuroscience, and turns dynamical system into a strong mathematical tool for effective modeling of individual components in the nervous system. In neuromorphic modeling, two-dimensional dynamical systems are popular because of their capability of phase portrait representation. The overall qualitative dynamics of a system can be easily investigated through the study of the phase portrait of the system including different types of equilibrium points and curves and consequently local velocity vectors in the phase space, and also a geometric representation of special trajectories determining topological behaviors of trajectories in a neighborhood in the phase space. Hence, several neuromorphic models are presented in the general form described as: (dx dt =F(x,y)+b dy dt =G(x,y)+c(1) where Fand Gare smooth functions, band care input variables and a number of auxiliary functions on x,ymay be attached. In this case, the phase plane is a two-dimensional space, (x,y)∈ Z2represents the location of state point in the plane, and (˙x,˙y) (velocity vector) determines the velocity and direction of the motion. According to Equation (1), the terms band c directly influence the velocity of the moving state point in both dimensions. Here, we temporarily ignore band c(b=c=0) in the general form to simplify developing the mapping scheme, and then apply their influences directly on the velocity vector in the last step of the hardware design. Frontiers in Neuroscience | www.frontiersin.org 3November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS In most neuromorphic dynamical systems with the general form of Equation (1) the equilibrium lines are single-valued functions. This implies that one dynamical variable is a singlevalued function of another one. Hence, the simplified general form of neuromorphic dynamical systems, ignoring input parameters band c, can be rewritten as one of these four different forms: (dx dt =α·(F(x)−y) dy dt =β·(G(x)−y)(2) (dx dt =α·(F(y)−x) dy dt =β·(G(x)−y)(3) (dx dt =α·(F(x)−y) dy dt =β·(G(y)−x)(4) (dx dt =α·(F(y)−x) dy dt =β·(G(y)−x). (5) Note that our approach supports all four different abovementioned forms, but considering this fact that most of the popular neuron models such as the FitzHugh-Nagumo (Fitzhugh, 1961), Izhikevich (Izhikevich, 2003), and AdEx (Brette and Gerstner, 2005) neural models, fit into the first simplified general form represented in Equation (2), we continue the rest of the paper based on this general form. 2.2. Cellular Space and Generalized Mapping In the first step, we map the continuous phase plane to a cellular plane. As shown in Figure 1A, in the cellular space, discrete points on the state variable axes and accordingly discrete meshlike points on the phase plane are considered. Thus, X= iwhere i∈M≡ {0,1,...,M−1}, and Y=jwhere j∈N≡ {0,1,...,N−1}are discrete variables parallel to xand yin the continuous space. The location of state point in the M×Ncellular space can be evaluated as: (X= ⌊x−xmin 1x⌋ Y= ⌊y−ymin 1y⌋(6) where 1x=xmax−xmin Mand 1y=ymax−ymin Nare cellular steps in each dimension. In the next step, which is one of the key steps in the mapping scheme, we evaluate the velocity and direction of the vector field in the cells in the way that it can be easily implemented on our proposed hardware platform. In CMDS approach (Bavandpour et al., 2014), the value of velocity and direction in each cell is directly stored in its corresponding memory cell in a cellular memory block, which requires embedding a memory block of size 2 ·N·Mmemory cells (memristors) in the hardware implementation. Here, we significantly shrink the size of memory block using a technique to evaluate the motion velocity and direction in the simplified general form. As mentioned, the first simplified general form represents the most common condition in neuromorphic dynamical systems where yis a function of xin both equations. In this condition, equilibrium curves of the system (considering dx dt =0,dy dt =0) are yeqx =F(x) and yeqy =G(x). In the cellular representation, we have: Yeqx(X)=F(xmin +X·1x) Yeqy(X)=G(xmin +X·1x)(7) where Yeqx and Yeqy are two arrays of Mreal numbers presenting the value of the yvariable for discrete points on the equilibrium curves corresponding to different discrete values of X. We refer to these arrays as the equilibrium arrays in the cellular space. Figure 1A shows the cellular phase plane and a sample equilibrium curve in the plane. As shown, the value of yvariable for cross points of equilibrium curve with X=i(i= 0,1,...,M−1) vertical lines are the elements of Yeq array. It will be substantiated that we can evaluate the velocity vector ( ˙ X,˙ Y) for all cells in the phase plane just using the equilibrium arrays. FIGURE 1 | (A) Cellular phase portrait of a sample dynamical system, one of its equilibrium curves, and evaluating the velocity of a sample point using equilibrium curve according to Equation (9), (B) absolute value of the velocity in the phase plane for a quadratic equilibrium curve used in Izhikevich neuron model (Izhikevich, 2003), (C) absolute value of the velocity in the phase plane for a cubic parabola equilibrium curve used in FHN (Fitzhugh, 1961) neuron model. Frontiers in Neuroscience | www.frontiersin.org 4November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS Consider a sample point (X=i,Y=j); the velocity vector for this point is calculated as: ˙ X(X,Y)=α·(F(xmin +X·1x)−(ymin +Y·1y)) ˙ Y(X,Y)=β·(G(xmin +X·1x)−(ymin +Y·1y)) . (8) Substituting Equation (7) in Equation (8), the velocity vector equation is given by: ˙ X(X,Y)=α·(Yeqx(X)−(ymin +Y·1y)) ˙ Y(X,Y)=β·(Yeqy(X)−(ymin +Y·1y)) . (9) As the above equation shows, the velocity vector depends on the equilibrium arrays and a number of predetermined parameters. In other words, the absolute value of the velocity at every state point is proportional to its absolute vertical distance (distance in ydirection) from the equilibrium curve, and its direction (forward/backward) is proportional to the position of the state point (above/below) to the equilibrium curve in the phase plane. The equations of cell change policy based on the timing parameters are represented in the next subsection after defining these parameters and drawing requisite timing and asynchrony conditions in our approach. Figures 1B,C are shown the absolute value of the velocity in the phase plane for two applied equilibrium curves in the neuron models. For hardware realization of our approach, we store the equilibrium arrays on the memristor crossbar and evaluate the velocity vector using hardware structure of Equation (9). Clearly, the resultant circuit is fully reconfigurable because the equilibrium arrays can be evaluated and programmed on the memristive memory block once a new dynamical system is presented. 2.3. Timing and Asynchrony in Cell Change Policy According to the previous subsection, the velocity of motion and the motion direction (forward/backward in each dimension) are variable in the cellular phase plane, and can be evaluated using Equation (9). Note that the evaluated velocities are in units of (1/s) where xand yare unitless quantities. For uniform representation in the cellular space, we ought to evaluate the cellular velocity with the unit of (cell/s). So, the velocity values are divided by the amount of cellular steps: (˙ Xcellular(X,Y)=˙ X(X,Y) 1x ˙ Ycellular(X,Y)=˙ Y(X,Y) 1y . (10) Here, with the aim to simplify the timing analysis, we convert the concept of motion velocity to the concept of motion time. Hence, the cellular motion time (the amount of time required for moving one cell) can be obtained by reversing the cellular motion velocity as:    T(X,Y) x=NF( 1 ˙ Xcellular ) T(X,Y) y=NF( 1 ˙ Ycellular )(11) where NF(·) is a function presenting restrictions imposed by hardware implementation that limits the time delay value in a specific boundary. A simple form for this boundary function can be represented as: NF(input) =   MAX input >MAX input MIN ≤input ≤MAX MIN input <MIN (12) where MIN and MAX are lower and upper boundaries. Clearly, our cellular approach is asynchronous, so the state variables are changed asynchronously. Considering this fact, it is crucial to proportionally handle the cellular motion time in one direction when other state variable is changed and consequently causes a change in motion time in both directions. In other words, time handling of the motion in each dimension is independent, but the motion in one dimension influences the motion time in another dimension by changing the vertical distance of the state point from the corresponding equilibrium curve, and it ought to be considered. Here, we illustrate the time handling procedure for a special case, and then derive a generalized rule to handle the two dimensional timing. Assume that the initial state (t=0) of the system is (i,j) and the corresponding motion times are T(i,j) xand T(i,j) ywhere |T(i,j) x|<|T(i,j) y|. At t= |T(i,j) x|the state point moves to (i+1,j) and the motion time of variable Xis updated to |T(i+1,j) x|, but this is not reasonable to update the motion time of variable Yto |T(i+1,j) y|without taking into account the elapsed time in point (i,j) resulting in no change in Ydirection. Hence, the passed time is considered as a portion of the new motion time, and deducted from this time amount: T(i+1,j) y(new) =(1 −|T(i,j) x| |T(i,j) y| )·T(i+1,j) y. (13) Regarding above explanations, an overall recursive rule can be derived to handle the dimensional motion times for moving from any given cell P∈M×Nto its neighbor cell Q∈M×N: tx r(Q)=((1 −ty r(P) tx r(P))· |T(Q) x|X(Q)=X(P) |T(Q) x|X(Q)6= X(P)(14) ty r(Q)=((1 −tx r(P) ty r(P))· |T(Q) y|Y(Q)=Y(P) |T(Q) y|Y(Q)6= Y(P)(15) where tx r(Q) and ty r(Q) are, respectively, the remaining time to change for Xand Ystate variables at Qcellular point which are decreased with passing operation time of the system. The cell change policy is given by: X(t+)=X(t)+sign(T(X(t),Y(t)) x)tx r=0 X(t) Otherwise (16) Y(t+)=Y(t)+sign(T(X(t),Y(t)) y)ty r=0 Y(t) Otherwise . (17) Thus, when the shorter remaining motion time for the faster state variable is elapsed, the state variable is changed by one Frontiers in Neuroscience | www.frontiersin.org 5November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS cell, and the remaining motion time for the other state variable is updated according to the motion time vector in the new state and using Equation (13). In the following sections, it is qualitatively and quantitatively shown that the dynamical behaviors of a given two dimensional dynamical system described generally by Equation (2) can be accurately mimicked by evaluating velocity vector given by Equation (9) using the equilibrium arrays, and following the rules of motion given by Equations (14–17). 3. MEMRISTIVE CIRCUIT IMPLEMENTATION The memristor is a nano-scale passive variable resistor with memory whose resistance changes depending on the polarity and magnitude of a voltage applied to the device terminals and the duration of this voltage application. Threshold condition is one of the key characteristics of the memristor. Based on this feature, small voltages across the memristor, below its threshold (Vth), do not cause a considerable memristance change, while larger voltages, greater than the threshold, induce much greater memristance changes (Jo et al., 2010). In this approach, we use the memristors as pre-programmed analog memories, so their memristance values ought to remain constant during operation time of the circuit. This implies that the maximum applied voltage to the memristors ought to be lower than the threshold voltage. The emerging technology of memristive circuits is one of the most promising technologies for hardware implementation of computational systems. Accordingly, developing new analog computation approaches compatible with a hybrid memristorcrossbar/CMOS platform is a significant ongoing research area (Bichler et al., 2012; Ligang et al., 2013; Merrikh-Bayat et al., 2013). In this paper, we engage this platform in the area of neuromorphic engineering and bio-modeling by developing a novel mapping scheme for neuromorphic dynamical systems, which makes them compatible for implementation on the aforementioned platform. In this section, we present the hardware structure of the system according to the proposed mapping scheme in the previous section. The proposed memristive circuit is shown in Figure 2. In this circuit, the memristors are connected in a crossbar architecture which offers flexibility, scalability, and simplicity, and also provides maximum density. Moreover, the memristor crossbar is separated from the CMOS-based part of the circuit causing the compatibility of the circuit for implementation on the promising hybrid memristor-crossbar/CMOS technology. As shown in the figure, the circuit is broken down to different parts which are explained one by one as follows: FIGURE 2 | Hardware structure of the proposed memristive dynamical circuit which is totally compatible with the hybrid memristor crossbar/CMOS platform because the memristor crossbar part of the circuit is isolated from the CMOS-based part. Frontiers in Neuroscience | www.frontiersin.org 6November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS 3.1. XR and YR In the proposed circuit, Xand Ydiscrete variables are, respectively, stored in two M-bit and N-bit bidirectional shift registers named XR and YR using one-hot encoding described as: XR(i)=1X(t)=i 0 Otherwise ,i=0,1, . . . , M−1 (18) YR(i)=1Y(t)=i 0 Otherwise ,i=0,1, . . . , N−1. (19) Clearly in this representation, increasing or decreasing one unit in the Xand Ydiscrete variables are, respectively, equivalent to one-bit forward or backward shift in the registers. Also, registers can be reset through a control signal. 3.2. X_DAC and Y_DAC These units receive the one-hot value of the discrete variables from XR and YR and produce their proportional analog values in the outputs Xa and Ya. The circuit diagram of these units contains an opamp-based multi-input analog adder where the input resistors are replaced by the memristors, so the gain of each input can be controlled by changing the value of its input memristance. Thus, as the Figure 2 shows, we set the memristances proportional to the weight array of the input bits in the one-hot encoding, which is an array of linearly increasing real numbers. Thus, the governing equation of these circuits can be given by: Xa = −6M−1 i=0XR(i)·Gmx(i)·Rf·vd Ya = −6N−1 i=0YR(i)·Gmy(i)·Rf·vd (20) where Rfis the feedback resistor of the op-amps, vdis the voltage level from digital registers representing logic “one." In this circuit, reading conductance of the memristors ought not to affect their value. On the other hand, it is mentioned that if the applied voltage to the memristor is below its threshold (Vth), it does not induce considerable change in the memristance. Hence, as an essential condition, the applied voltage to the memristors ought to be lower than the threshold value (vd< Vth). The conductance values functions Gmx(·) and Gmy(·) for the memristors are described as: Gmx(i)=(Ax.i+1) ·Gmx0 Gmy(i)=(Ay.i+1) ·Gmy0 (21) where Ax,Ayare the normalizing coefficients and Gmx0,Gmy0 are the initial conductance parameters. For uniformity, the output voltage range for both units ought to be clamped in (vop min,vop max). For satisfying lower boundary condition of this interval, we have: vop min =Gmx0 ·Rf·vd=Gmy0 ·Rf·vd ⇒Gmx0 =Gmy0 (22) and for satisfying upper boundary condition, we have: vop max =(Ax·(M−1) +1) ·Gmx0 ·Rf·vd =(Ay·(N−1) +1) ·Gmx0 ·Rf·vd ⇒Ax Ay=N−1 M−1. (23) The parameters can get various values satisfying the conditions described in Equations (22) and (23). Note that there is an undesirable negative sign multiplied in the output voltages Xa and Ya in Equation (20), but it is not the source of any difficulties in our approach, because what is calculated in the next blocks is just the difference of the voltages and the negative sign can be easily handled. 3.3. X_eqa and Y_eqa In these units, the equilibrium arrays given by Equation (7) are stored on the memristors, and the analog output voltages proportional to the equilibrium values for current YR state is produced in the Xeqa and Yeqa outputs. Thus, the governing equation of these circuits can be given by: Xeqa = −6M−1 i=0XR(i)·Geqx(i)·Rf·vd Yeqa = −6N−1 i=0XR(i)·Geqy(i)·Rf·vd (24) where Rfis the feedback resistor of the op-amps, vdis the voltage level from digital registers representing logic “one,” and the conductance values functions Geqx(·) and Geqy(·) for the memristors are described as:    Geqx(i)=(Ay·Yeqx(i)−ymin 1y+1) ·Gmy0 Geqy(i)=(Ay·Yeqy(i)−ymin 1y+1) ·Gmy0 . (25) Similar to the equations of the X_DAC, Y_DAC units presented in Equation (20), there is an undesirable negative sign in the expression for the output voltages Xeqa and Yeqa in Equation (24), but it is not the source of any difficulties. This is because voltage differences are calculated in the next blocks and the sign can be easily handled. 3.4. X_vel and Y_vel As shown in Figure 2, each of these units contains one analog voltage subtractor and one analog voltage adder which can be merged in one op-amp based circuit. The subtractor receives the proportional analog value of equilibrium arrays (Xeqa or Yeqa) and the Ystate variable (Ya), and produces the proportional analog value of the velocity (Vdx or Vdy) using Equation (9). The governing equations of the subtractor circuits are: Vdx =Gsx ·(Xeqa −Ya) Vdy =Gsy ·(Yeqa −Ya) (26) where Gsx and Gsy are the gain of subtractors. Note that the coefficients αand βin Equation (9), and the aforementioned undesirable negative sign from the previous steps ought to be considered in the gain of the subtractors in this step. In the next step, the adders add the proportional analog voltages of the b and cparameters in Equation (1), which were ignored (b= c=0) for the sake of convenience in the mathematical process and hardware design, to the output of subtractors to apply their influence directly to the velocity of moving state point. The adder equations are given by: Vx=Vdx +bb Vy=Vdy +cc . (27) Frontiers in Neuroscience | www.frontiersin.org 7November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS For adjusting the mathematical parameters band cwith their proportional analog voltage value bb and cc, they are then multiplied by the adjusting coefficients as: bb =Gb·b cc =Gc·c(28) where Gband Gcare the adjusting parameters. We will determine the essential conditions on how all the introduced parameters in X_vel and V_vel units can be chosen, after introducing the VCO units and their parameters. 3.5. X_VCO and Y_VCO These units contain VCOs and controller circuit. The VCO receives the analog voltage corresponding to the velocity (Vxor Vy) and produces the linearly proportional motion clock pulse, which is used for clocking the shift of the XR and YR registers, with a specific range of the frequency. Functional equation of the VCOs is given by: fXclk out (Vx)=   fmax |Vx|>UT Gvco · |Vx|LT <|Vx| ≤ UT fmin |Vx| ≤ LT (29) fYclk out (Vy)=   fmax |Vy|>UT Gvco · |Vy|LT <|Vy| ≤ UT fmin |Vy| ≤ LT (30) where Gvco is the VCO coefficient, the UT and LT are, respectively, the upper-hand and lower-hand input voltage threshold of the VCOs and fmax and fmin are, respectively, the maximum and minimum output frequencies of the VCOs. The VCOs produce two asynchronous clock pulse signals based on the derived equation on their output pins X_clk and Y_clk. Note that the VCOs inherently satisfy Equations (14) and (15). The controller circuit determines the direction of the motion (shift direction of the XR and YR state registers) on the pins X_SD and Y_SD. Zero logical value on these pins means that the next change in their proportional state registers is a backward shift, and the one logical value means that the next change is a forward shift: X_SD =1Vx>0 0 otherwise (31) Y_SD =1Vy>0 0 otherwise . (32) Eventually, the state registers are shifted at the positive edge of their proportional clock pulse signals X_clk and Y_clk: XR(t+ posedge)       XR ≫1 XR(M−1) 6= 1,X_SD =1 XR ≪1 XR(0) 6= 1,X_SD =0 XR XR(M−1) =1,X_SD =1 XR XR(0) =1,X_SD =0 (33) YR(t+ posedge)       YR ≫1 YR(N−1) 6= 1,Y_SD =1 YR ≪1 YR(0) 6= 1,Y_SD =0 YR YR(N−1) =1,Y_SD =1 YR YR(0) =1,Y_SD =0 (34) where t+ posedge is the time when a positive edge is appeared on the input clock pulse of the state registers. As we mentioned before, we ought to draw the essential conditions on the coefficients for proper operation of the system. According to the Equation (2), the analog velocity of the state point in xdirection is αtimes of its absolute analog vertical distance (distance in ydirection) from the equilibrium curve, and in ydirection is βtimes of its absolute analog vertical distance from the equilibrium curve. Accordingly, in the cellular space, the cellular velocity of the state point in xdirection is α·1y 1xtimes of its absolute cellular vertical distance from the equilibrium curve, and in ydirection is βtimes of its absolute cellular vertical distance from the equilibrium curve. Note that in the cellular space, the cellular velocity is equal to the frequency of VCOs. So, for a special condition where the vertical distance of the state point from the equilibrium curve is equal to one cell, the cellular output frequencies of the VCOs are fXclk out =α·1y 1xand fYclk out =β. In this condition, if we merge the Equations (21), (24–26), (29), and (30) (assuming b=c=0) to evaluate the output frequency of the VCOs, we have: Ay·Gmy0 ·Rf·vd·Gsx ·Gvco =α·1y 1x Ay·Gmy0 ·Rf·vd·Gsy ·Gvco =β. (35) By dividing these two equations, we obtain: Gsx Gsy =α·1y β·1x. (36) Now, assume the special condition where the vertical distance of the state point from the equilibrium curve is equal to zero and b=1x,c=1y. In this condition, the frequencies of the VCOs are fXclk out =fYclk out =1. According to the Equations (28–30), we have: Gb·1x·Gsx =1 Gc·1y·Gsy =1. (37) By dividing these two equations, we have: Gb Gc =β α. (38) Therefore, the parameters must satisfy the essential conditions given by Equations (35–38). Auxiliary Function: This unit is a logical block for implementing the auxiliary functions like threshold conditions and resetting functions. The proposed circuit is significantly more efficient in term of area in comparison with the previous similar approach named CMDS (Bavandpour et al., 2014) while achieving the same performance and accuracy. In the Table 1, our approach is compared with the CMDS in terms of circuit elements/blocks for implementing a two-dimensional 100 ×100 cellular dynamical system. Frontiers in Neuroscience | www.frontiersin.org 8November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS TABLE 1 | Hardware assessment of the proposed approach in comparison with the CMDS approach for implementing a two-dimensional 100 ×100 cellular dynamical system. Circuit element/block Proposed approach CMDS Memristor 400 20,000 Switch – 400 Analog adder 1 1 VCO 2 2 Auxiliary function 1 1 State variable registers 2 2 Analog output Yes No One-hot digital output Yes Yes 4. DYNAMICAL BEHAVIOR ANALYSIS The behavior of a biological neuron depends on many intrinsic and extrinsic factors, such as the morphology of its dendritic tree, the type and characteristics of ion-gated channels and voltage expressed by the neuron, the location of the stimulating input, and so on. These factors are indeed important, because they determine not only the neuronal response but also the rules that govern dynamics of the neuron. Accordingly, different neuron models are developed to mimic different responses of the neuron and their related dynamics (bifurcations). Our generalized circuit can mimic the dynamics and the various responses of a wide range of neuron models with different computational complexity. In this section, we investigate three main neuron models as the FHN (Fitzhugh, 1961), AdEx (Brette and Gerstner, 2005), and Izhikevich neuron models (Izhikevich, 2003), and show that the circuit is able to properly mimic different responses of different neuron models. 4.1. FitzHugh-Nagumo (FHN) Neuron Model The FHN is a two-dimensional neuron model derived form the simplified Hodgkin-Huxley (HH) model for biological process of spike generation in squid large axons (Fitzhugh, 1961). This model is considered as a relatively complex mathematical model due to a third power factor in its equations. FHN is represented by a simplified set of equations as: ˙v=v−v3 3−u+I ˙u=a·(v+0.7 −0.8u)(39) where vand uare, respectively, the membrane potential variable and the recovery variable, Iis the injected stimulus current, and ais an adjustable parameter. As shown, FHN model utilizes no auxiliary resetting function to reproduce spiking behaviors. Hence, it does not employ the auxiliary function block in the circuit implementation. Mathematical analysis of this model simply shows that this system behaves like a relaxation oscillator. It implies that if the stimulus input current of the model exceeds a specific threshold, the model will follow a characteristic trajectory in phase space before the state variables relax back to their resting values. Typically, this behavior is how a biological neuron generates a spike (defined as a short and fast elevation of membrane potential) after injecting a short stimulus input current. Comparing the FHN model with the general form of target systems presented in the previous section, it can be rewritten in the general form of Equation (2) and then mapped easily on the proposed platform. The general form of FHN neuron model in our mapping is given by:            x=v,y=u α=1, β =0.8a b=I,c=0 F(x)=x−v3 3 G(x)=(x+0.7)/0.8 . (40) Our FHN-MDS approach can exhibit all significant qualitative phenomena of the original FHN model and their underlying bifurcations. In this study, four main bio-inspired phenomena of FHN model are individually investigated on our approach. 4.1.1. Absence of All-or-none Spike The original FHN model is capable of producing the absence of all-or-none spikes as it is produced in HH model in response to stimulus current I. According to this behavior, the amplitude of output relaxation trajectory is directly and continuously related to the amplitude of the injected current pulse Ito the model. Weak stimulation produces a small-size relaxation trajectory as a subthreshold response. Stronger stimulation produces an intermediate-size trajectory as a partial-size spike, and strong stimulation produces a large-size trajectory as a supra-threshold firing response. Figure 3A shows time-domain and phase-plane representations of this phenomenon produced by our approach. As shown, size of trajectories in the phase plane is directly related to the amplitude of applied input pulses in the time domain. 4.1.2. Excitation Block According to this response, the neuron stops repetitive firing and goes back to a stable resting state as the amplitude of the input current increases. This type of response is based on a special bifurcation scenario in the phase plane representation. When input current Iis relatively weak, the equilibrium point (the point that the nullclines intersect) is located on the left branch of X- nullcline resulting in a stable state, and hence the model is in the resting state. As the input current increases, the nullcline shifts upward in the phase plane, so the equilibrium point moves onto the middle branch of the nullcline resulting in an unstable state. Hence, the model produces a limit cycle in the phase plane and starts a repetitive firing in the output. Further increase of the stimulus current shifts the equilibrium point to the right branch of the N-shaped nullcline resulting in a stable state again, so the previous firing state is blocked. The FHN-MDS approach can exhibit this behavior as accurate as the original FHN model. Transition from the resting state to the repetitive firing, and then the blocking state produced by the FHN-MDS approach using a ramp input current are shown in Figure 3B in both phase and time domains. Frontiers in Neuroscience | www.frontiersin.org 9November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS TABLE 2 | The relative spike timing error and relative spike energy error of our approach in comparison with the original model for tonic spiking and regular bursting responses of the FitzHugh-Nagumo, AdEx, and Izhikevich neuron models. Neuron model Response type 20-bit (%) 40-bit (%) 60-bit (%) 80-bit (%) 100-bit (%) (A) RELATIVE SPIKE TIMING ERROR FHN Tonic spiking 1.78 1.04 0.67 0.43 0.26 AdEx Tonic spiking 2.29 1.34 1.00 0.79 0.54 Regular bursting 3.52 1.73 1.08 0.81 0.65 Izhikevich Tonic spiking 2.03 1.22 0.88 0.54 0.32 Regular bursting 3.01 1.69 1.01 0.76 0.55 (B) RELATIVE SPIKE ENERGY ERROR FHN Tonic spiking 3.24 1.78 1.22 0.88 0.62 AdEx Tonic spiking 9.41 5.09 3.99 2.98 2.07 Regular bursting 17.55 8.77 5.04 4.57 3.95 Izhikevich Tonic spiking 7.85 4.08 3.12 2.01 1.44 Regular bursting 10.14 5.00 3.85 2.97 2.45 FIGURE 8 | Signal shape of regular bursting produced by MDS-AdEx approach using 20, 40, and 60 bits. explained, MDS approach produces both analog and digital one-hot outputs, and accepts analog voltage as injected input to the system. This feature makes it feasible to replace conventional hardware implementation of single neural cells with the MDS cells in a wide range of applied digital/analog weighted networking schemes with plasticity to achieve a network of optimum accurate general neural cells. Besides, it is convenient to propose novel networking schemes to connect Memristive Dynamical System (MDS) cells with variable weights, so the resultant system is flexible network of general neural cells with variable intera-cell dynamics and the capability of applying bioinspired learning schemes. Figure 9 shows the overall conceptual structure of the system which contains the network of the proposed Memristive Dynamical System (MDS) cells connecting to each other through a weighting and integrating Synaptic Unit (SU) and sharing a connection bus for programming the reconfigurable MDS cells. In this conceptual model of the system, the user is provided with a user interface to determine the equilibrium functions and the desired intervals [xmin,xmax) and [ymin,ymax) on the dynamical variables xand y. A host processor calculates the equilibrium arrays Yeqx and Yeqy and their equivalent memristances (where b=c=0), and applies a compatible programmer system to program them on the stand-alone neuro-chip using a developed communication protocol. During this process, the equilibrium memristor arrays are selected and programmed using the special pulse shapes proposed for writing into the memristors. Clearly, the reconfigurable memristors ought to be updated every time the dynamical system parameters are changed. Then, the MDS cells starts to track their deterministic trajectory in the phase plane according to the moving policy given by Equations (14–17) implemented on the neuro-chip. 7.2. Learning Applying the flexible memristors in a well-mapped general hardware brings a significant advantage of learning capability to the MDS approach. In other words, the MDS is capable of learning various intra-cell dynamics with various complexities to reproduce various signal shapes specifically spike patters in its output. One of the major advantages of our novel MDS approach over previous CMDS (Bavandpour et al., 2014) is our novel mapping which uses the vertical distance of the state point from the nullclines to calculate the velocity. In this mapping, the nullclines while in the CMDS the memristors store the value of velocity in all over the cellular phase plane. Therefore, learning process in MDS is significantly simpler and more effective than CMDS. Considering this feature, different Frontiers in Neuroscience | www.frontiersin.org 16 November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS FIGURE 9 | Conceptual network-level structure for networking the MDS cells. The user inputs required mathematical information of the model into the host computer, and the required equilibrium curves data are calculated and programmed into the neuro-chip memristors using a programmer unit. intra-cell learning schemes can be applied to the MDS to teach it different dynamical behaviors with no straight knowledge on the dynamical equations. Figure 10 shows the conceptual block diagram of a possible learning system for MDS. As shown, target system is a black box system with unknown mathematical model governing its dynamical behavior. Assume that the system can be modeled with a two-dimensional dynamical system, and the state variables of the black box are observable. In this condition, the dynamical behavior of the system can be extracted by applying different stimulus inputs in different initial states. Output signals of the black box can be applied to MDS system as teacher signals to follow its behavior. Learning controller detects the black box state, and modifies MDS behavior in the same state using teacher signal. Modification of MDS behavior during the learning phase is done by changing the memristance of the memristors using learning signals. 8. DISCUSSION This study presented an unconventional computing approach based on a novel general mapping for dynamical systems in two-dimensional cellular phase space, and then its hardware implementation on efficient hybrid memristors-crossbar/CMOS memristive circuit. The proposed approach calculates the Frontiers in Neuroscience | www.frontiersin.org 17 November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS FIGURE 10 | A possible learning system block diagram which extracts the behavior of an unknown black box system and teaches an MDS cell to imitate this behavior. velocity vector using the vertical distance of the state point from the nullclines, and applies an unconventional VCO-based asynchronous technique to track the state trajectory in the phase space. This approach employed 4nmemristors and no switch for implementing an n-cell system in comparison with the 2 ·n2memristors and 2nswitches of a Cellular Memristive Dynamical System (CMDS). In addition, our proposed circuit provides both analog and one-hot digital values of dynamical variables offering a wide range of choices for interconnections and networking schemes. The fidelity of MDS approach was investigated on widely used two-dimensional neuromorphic dynamical systems such as the FitzHugh-Nagumo model, the Adaptive Exponential (AdEx) integrate and fire model, and the Izhikevich neuron model covering a wide range of reallife neuronal dynamic behaviors and responses. The MDS-based neurons were simulated with various stimulus currents, and it was showed that they can mimic every key response produced by the original models. Moreover, we clarified different definitions of error in neuromorphic dynamical systems and investigated spike timing error and spike shape error for different responses of different neuron models. The error results showed that our approach can properly mimic the exact behavior of the dynamical system. In conclusion, the key advantages of the proposed approach can be listed as: •It is general and highly programmable. •It achieves a relatively high accuracy in a low resolution cellular space. •Its implementation cost is almost independent from the computational effort of the target mathematical model. •It is capable of implementing mathematically indescribable dynamical systems. •It is implementable on a nanoscale efficient memristor/CMOS hardware platform. •It is conveniently feasible to apply conventional analog or digital networking schemes to the MDS cells, and also propose novel networking schemes. •It is capable of learning unknown intra-cell dynamics. Toward the future roadmap, the problems which ought to be solved can be listed as: •A detailed circuit design and analysis process to achieve an optimum circuit for different blocks of the MDS system. •Developing unconventional learning methods for MDS-based dynamical systems. •A novel compatible mathematical method for stability analysis of the proposed cellular mapping. FUNDING This work was partially funded by EU HBP project under grant number FP7-ICT-2013-FET-F-604102, by EU H2020 ECOMODE project under grant agreement 604102, and by Spanish research grant (with support from the European Regional Development Fund) TEC2012-37868-C04-01 (BIOSENSE). REFERENCES Abbott, L. F. (1999). Lapique’s introduction of the integrate-and-fire model neuron. Brain Res. Bull. 50, 303–304. doi: 10.1016/S0361-9230(99)00161-6 Ahmadi, A., and Soleimani, H. (2011). “A GPU based simulation of multilayer spiking neural networks,” in 19th Iranian IEEE Conference on Electrical Engineering (ICEE) (Tehran), 1–5. Alibart, F., Gao, L., Hoskins, B. D., and Strukov, D. B. (2012). High precision tuning of state for memristive devices by adaptable variation-tolerant algorithm. Nanotechnology 23:075201. doi: 10.1088/0957-4484/23/7/075201 Arthur, J. V., and Boahen, K. (2011). Silicon neuron design: the dynamical systems approach. IEEE Trans. Circuits Syst. I 58, 1034–1043. doi: 10.1109/TCSI.2010.2089556 Bahar, R. I., Hammerstrom, D., Harlow, J., Joyner, J. W. H., Lau, D., Marculescu, C., et al. (2007). Architecturesfor silicon nanoelectronics and beyond. J. Comput. 40, 25–33. doi: 10.1109/MC.2007.7 Bavandpour, M., Soleimani, H., Bagheri-Shouraki, S., Ahmadi, A., Abbott, D., and Chua, L. O. (2014). Cellular memristive dynamical systems (CMDS). Int. J. Bifurcat. Chaos 24:1430016. doi: 10.1142/s021812741430016x Benjamin, B. V., Peiran, G., McQuinn, E., Choudhary, S., Chandrasekaran, A. R., Bussat, J. M., et al. (2014). Neurogrid: a mixed analog-digital multichip system for large-scale neural simulations. Proc. IEEE 102, 699–716. doi: 10.1109/JPROC.2014.2313565 Bichler, O., Suri, M., Querlioz, D., Vuillaume, D., DeSalvo, B., and Gamrat, C. (2012). Visual pattern extraction using energy-efficient 2-PCM synapse neuromorphic architecture. IEEE Trans. Electron Devices 59, 2206–2214. doi: 10.1109/TED.2012.2197951 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 Cassidy, A., Merolla, P., Arthur, J., Esser, S., Jackson, B., Alvarez-icaza, R., et al. (2013). “Cognitive computing building block: a versatile and efficient digital neuron model for neurosynaptic cores,” in International Joint Conference on Neural Networks (IJCNN) (Dallas, TX: IEEE). Chua, L. O. (1971). Memristor-the missing circuit element. IEEE Trans. Circuit Theory 18, 507–519. doi: 10.1109/TCT.1971.1083337 Clopath, C., Jolivet, R., Rauch, A., Luscher, H., and Gerstner, W. (2007). Predicting neuronal activity with simple models of the threshold types: adaptive Frontiers in Neuroscience | www.frontiersin.org 18 November 2015 | Volume 9 | Article 409 Bavandpour et al. Generalized reconfigurable MDS exponential integrate-and-fire model with two compartments. Neurocomputing 70, 1668–1673. doi: 10.1016/j.neucom.2006.10.047 Fitzhugh, R. (1961). Impulses and physiological states in models of nerve membrane. J. Biophys. 1, 445–466. doi: 10.1016/S0006-3495(61)86902-6 Furber, S. B., Galluppi, F., Temple, S., and Plana, L. A. (2014). The SpiNNaker project. Proc. IEEE 102, 652–665. doi: 10.1109/JPROC.2014.2304638 Gerstner, W., and Kistler, W. M. (2002). Spiking Neuron Models: Single Neurons. Cambridge, UK: Cambridge University Press. Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764 Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440 Jo, S. H., Chang, T., Ebong, I., Bhadviya, B. B., Mazumder, P., and Lu, W. (2010). Nanoscale memristor device as synapse in neuromorphic systems. Nano Lett. 10, 1297–1301. doi: 10.1021/nl904092h Kavehei, O., Al-Sarawi, S., Cho, K. R., Eshraghian, K., and Abbott, D. (2012). An analytical approach for memristive nanoarchitectures. IEEE Trans. Nanotechnol. 11, 374–385. doi: 10.1109/TNANO.2011.2174802 Ligang, G., Merrikh-Bayat, F., Alibart, F., Xinjie, G., Hoskins, B. D., Kwang-Ting, C., et al. (2013). “Digital-to-analog and analog-to-digital conversion with metal oxide memristors for ultra-low power computing,” in IEEE/ACM International Symposium on Nanoscale Architectures (NANOARCH 2013) (New York, NY), 19–22. Linares-Barranco, B., Sánchez-Sinencio, E., Rodríguez-Vázquez, A., and Huertas, J. L. (1991). A CMOS implementation of Fitzhugh-Nagumo neuron model. IEEE J. Solid State Circuits 26, 956–965. doi: 10.1109/4.92015 Matsubara, T., and Torikai, H. (2013). Asynchronous cellular automaton-based neuron: theoretical analysis and on-FPGA learning. IEEE Trans. Neural Netw. Learn. Syst. 24, 736–748. doi: 10.1109/TNNLS.2012.2230643 Merrikh-Bayat, F., Merrikh-Bayat, F., and Bagheri-Shouraki, S. (2013). Neurofuzzy computing system with the capacity of implementation on memristorcrossbar and optimization-free hardware training. IEEE Trans. Fuzzy Syst. 22, 1272–1287. doi: 10.1109/TFUZZ.2013.2290140 Minkovich, K., Thibeault, C. M., O’Brien, M. J., Nogin, A., Cho, Y., and Srinivasa, N. (2014). HRLSim: a high performance spiking neural network simulator for GPGPU clusters. IEEE Trans. Neural Netw. Learn. Syst. 25, 316–331. doi: 10.1109/TNNLS.2013.2276056 Moradi, S., and Indiveri, G. (2014). An event-based neural network architecture with an asynchronous programmable synaptic memory. IEEE Trans. Biomed. Circuits Syst. 8, 98–107. doi: 10.1109/TBCAS.2013.2255873 Schemmel, J., Fieres, J., and Meier, K. (2008). “Wafer-scale integration of analog neural networks,” in IEEE World Congress on Computational Intelligence In Neural Networks, IJCNN 2008 (Hong Kong), 431–438. Soleimani, H., Ahmadi, A., and Bavandpour, M. (2012). Biologically inspired spiking neurons: piecewise linear models and digital implementation. IEEE Trans. Circuit Syst. I 59, 2991–3004. doi: 10.1109/TCSI.2012.2206463 Soleimani, H., Ahmadi, A., Bavandpour, M., and Sharifipoor, O. (2014). A generalized analog implementation of piecewise linear neuron models using CCII building blocks. Neural Netw. 51, 26–38. doi: 10.1016/j.neunet.2013.12.004 Soleimani, H., Bavandpour, M., Ahmadi, A., and Abbott, D. (2015). Digital implementation of a biological astrocyte model and its application. IEEE Trans. Neural Netw. Learn. Syst. 26, 127–139. doi: 10.1109/TNNLS.2014. 2311839 Strukov, D., and Likharev, K. (2006). “A reconfigurable architecture for hybrid CMOS/nanodevice circuits,” in Proceedings of the 2006 ACMISIGDA 14th International Symposium on Field Programmable Gate Arrays (Monterey, CA), 131–140. Strukov, D. B., Snider, G. S., Stewart, D. R., and Williams, R. S. (2008). The missing memristor found. Nature 453, 80–83. doi: 10.1038/nature06932 Weinstein, R. K., Reid, M. S., and Lee, R. H. (2007). Methodology and design flow for assisted neural-model implementations in FPGAs. IEEE Trans. Neural Syst. Rehabil. Eng. 15, 83–93. doi: 10.1109/TNSRE.2007.891379 Wijekoon, J. H., and Dudek, P. (2012). VLSI circuits implementing computational models of neocortical circuits. J. Neurosci. Methods 210, 93–109. doi: 10.1016/j.jneumeth.2012.01.019 Yang, J. J, Pickett, M. D., Li, X., Ohlberg, D. A., Stewart, D. R., and Williams, R., S. (2008). Memristive switching mechanism for metal/oxide/metal nanodevices. Nat. Nanotechnol. 3, 429–433. doi: 10.1038/nnano.2008.160 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. Copyright © 2015 Bavandpour, Soleimani, Linares-Barranco, Abbott and Chua. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms. Frontiers in Neuroscience | www.frontiersin.org 19 November 2015 | Volume 9 | Article 409