DDA on FPGA
A modern Analog Computer

The Digital Differential Analyzer (DDA) is a device to directly compute the solution of differential equations. The DDA is a simulation of an analog computer, and is inherently parallel in operation. It can be mapped easily into an FPGA by defining mathematical integrators, adders, multipliers and other operations, then wiring them together to produce the desired dynamics. The system described here has an maximum integrator update rate of higher than 50 MHz, so it is possible to accurately simulate audio-rate signals.

For more background on analog computers see:

• Building an analog computer (Rice)
• Analog Computing (Vaxman)
  • Simulating a bouncing ball in a box
  • Workshop on Analog Computing (Hanson)
  • A practical approach to analog computers ( EAI)
  • Handbook of Analog Computation (EAI)
  • Fundamentals of the Analog Computer (Howe)
• Educational Analog Computer (Systron-Donner)

On this page you can follow these links to examples of DDAs applied to linear and nonlinear systems:

• Second order system impulse response
• Second order system driven
• Izhikevich quadratic neuron model
• FitzHugh-Nagumo (Van der Pol) cubic neuron model
• Second order system with 27-bit resolution
• Second order system controlled by a NiosII cpu
• Second order system with modularized integrator and NiosII cpu

A summary of the second order system controlled by NiosII appeared in Circuit Cellar Magazine (manuscript).

Second order system (damped spring-mass oscillator):
As an example, consider the linear, second-order differential equation resulting from a damped spring-mass system:

d2x/dt2 = -k/m*x-d/m*(dx/dt)

where k is the spring constant, d the damping coefficient, m the mass, and x the displacement. We will simulate this by converting the second-order system into a coupled first-order system. If we let v1=x and v2=dx/dt then the second order equation is equivalent to

dv1/dt = v2
dv2/dt = -k/m*v1-d/m*v2

These equations can be solved by wiring together two integrators, two multipliers and an adder as shown below. In the past this would have been done by using operational amplifiers to compute each mathematical operation. Each integrator must be supplied with an initial condition.



Converting this diagram to Verilog, the top-level module verilog code defines the 18-bit, signed, state variables and a clock divider variable (count). The clocked section resets and updates the state variables. The combinatorial statements compute the Euler approximation to the F(t,V(n)). The separate multiply module ensures that the multiplies will be instantiated as hardware multipliers. The Audio_DAC_ADC module was modifed to allow either ADC-to-DAC passthru or to connect the computation output to the DAC, depending on the position of SW17. SW17 up connects the computation.

/state variables
reg signed [17:0] v1, v2 ;
wire signed [17:0] v1new, v2new ;
//signed mult output
wire signed [17:0] v1xK_M, v2xD_M ;
// the clock divider
reg [4:0] count;


//Update state variables of simulation of spring- mass
	always @ (posedge CLOCK_50) 
	begin
		count <= count + 1; 
		if (KEY[3]==0) //reset
		begin	
			v1 <= 32'h10000 ; // 
			v2 <= 32'h00000 ;
			//count <= 0;
		end
		else if (count==0)
		begin 
			v1 <= v1new ;
			v2 <= v2new ;
		end
	end
	
	// Compute new F(t,v) with dt = 2>>9
	// v1(n+1) = v1(n) + dt*v2(n)
	assign v1new = v1 + (v2>>>9);
	// v2(n+1) = v2(n) + dt*(-k/m*v1(n) - d/m*v2(n))
	signed_mult K_M(v1xK_M, v1, 18'h10000);
	signed_mult D_M(v2xD_M, v2, 18'h00800);
	assign v2new = v2 - ((v1xK_M + v2xD_M)>>>9);
	
module signed_mult (out, a, b);
	output 		[17:0]	out;
	input 	signed	[17:0] 	a;
	input 	signed	[17:0] 	b;
	wire	signed	[17:0]	out;
	wire 	signed	[35:0]	mult_out;
	assign mult_out = a * b;
	assign out = {mult_out[35], mult_out[32:16]};
endmodule

// dt = 2>>9
// v1(n+1) = v1(n) + dt*v2(n)
assign v1new = v1 + (v2>>>9)  ;
// v2(n+1) = v2(n) + dt*(-k/m*v1(n) - d/m*v2(n))
signed_mult K_M(v1xK_M, v1, 18'h1_0000);
signed_mult D_M(v2xD_M, v2, 18'h0_0800); //light damping
//scale the input so that it does not saturate at resonance
signed_mult Sine_gain(Sinput,{sine_out[15],sine_out[15],sine_out}, 18'h0_4000);
assign v2new = v2 - ((v1xK_M + v2xD_M )>>>9) + (Sinput>>>9) ; 

Two mpegs of the scope screen show the low-damping (band pass) case for fast and slow frequency steps. The top trace is the input, and the bottom is the filter output. A sudden increase in output is seen near the resonant frequency of 485 Hz. If you stop the video at a frequency below the resonant frequency, you will see that the input and output are in phase. Above the resonant frequency, the phase is 180 degrees. The whole project is zipped here.

Izhikevich Neural Model:
One Neuron
Eugene Izhikevich developed a simple, semiempirical, model of cortical neurons. The properties of each neuron are controlled by 4 parameters. His web site includes matlab programs and detailed descriptions of the cell model. There are two state variables: v is the membrane potential and u is the membrane recovery variable. The differential equations, parameter settings and examples of model output are shown below in a figure from Izhikevich's web site. An electronic version of the figure and reproduction permissions are freely available at www.izhikevich.com.

The top-level verilog module includes the following code to build one cell. The cell parameters and state variables are reset when button KEY3 is pushed. The count variable is a clock prescaler to slow the computation down by a factor of 4096 so that it can be output through the audio codec. Note that v and u are scaled down by a factor of 100 so that the dynamic range is -1 to +1, as required by the arithmetic system and the D/A converter. Compared to the equations above, the scaled equations have all constant terms divided by 100. The a and b constants are not scaled (because they are multiplied by the scaled variables), but c and d are divided by 100. The constant (.04) on the v2 term above is scaled up by 100 because squaring (the scaled) v drops the value by 10,000. The final scaled equations are:

v1(n+1) = v1(n) + dt*(4*(v1(n)^2) + 5*v1(n) +1.40 - u1(n) + I)
u1(n+1) = u1 + dt*a*(b*v1(n) - u1(n))

where

• dt=1/16 (shift right 4 times)
• a in the range (0 to 0.1)
• b in the range (0 to 0.3)
• c in the range (-0.70 to -0.50)
• d in the range (0.0005 to 0.08)

The first pass of the code is shown below.

//analog simulation of neuron
	always @ (posedge CLOCK_50) 
	begin
		count <= count + 1; 
		if (KEY[3]==0) //reset
		begin	
			v1 <= 18'sh3_4CCD ; // -0.7
			u1 <= 18'sh3_CCCD ; // -0.2
			a <=  18'sh0_051E ; // 0.02
			b <=  18'sh0_3333 ; // 0.2
			c <=  18'sh3_8000 ; // -0.5
			d <=  18'sh0_051E ; // 0.02
			p <=  18'sh0_4CCC ; // 0.30
			c14 <= 18'sh1_6666; // 1.4
			I <= 18'sh0_2666 ;  // 0.15
		end
		else if (count==0) 
		begin
			if ((v1 > p)) 
			begin 
				v1 <= c ; 		
				u1 <= ureset ;
			end
			else
			begin
				v1 <= v1new ;
				u1 <= u1new ; 	
			end
		end 
	end
	
	// dt = 1/16 or 2>>4
	// v1(n+1) = v1(n) + dt*(4*(v1(n)^2) + 5*v1(n) +1.40 - u1(n) + I)
	// but note that what is actually coded is
	// v1(n+1) = v1(n) + (v1(n)^2) + 5/4*v1(n) +1.40/4 - u1(n)/4 + I/4)/4
	signed_mult v1sq(v1xv1, v1, v1);
	assign v1new = v1 + ((v1xv1 + v1+(v1>>>2) + (c14>>>2) - (u1>>>2) + (I>>>2))>>>2);
	
	// u1(n+1) = u1 + dt*a*(b*v1(n) - u1(n))
	signed_mult bb(v1xb, v1, b);
	signed_mult aa(du, (v1xb-u1), a);
	assign u1new = u1 + (du>>>4)  ; 
	assign ureset = u1 + d ;

Two scope images show the voltage on the top trace and the membrane recovery variable on the bottom. The traces shown use parameters suitable for a bursting neuron.

The whole project is zipped here.

Two Neurons with modularized construction
The entire neural model, including integrators can be abstracted to a module so that the top-level module has only the neuron instantiations and the interconnections. Interconnections are made by connecting the spike output of one cell to the current input of another. The bias current in this version can be set from the switches and the clock and reset signals are a little cleaner looking. Each of the neurons can have its (binary) spike train brought out to the parallel port.

  // clock divider
	always @ (posedge CLOCK_50) 
	begin
		count <= count + 1;
	end
	assign neuronClock = (count==0);	
	
	assign I1 = {2'h0,SW[15:0]} ; //the input bias currents
	assign I2 = {2'h0,SW[15:0]} ; 
	
	assign neuronReset = ~KEY[3]; //the reset pushbutton
	
	//burster parameters
	assign a1 =  18'sh0_051E ; // 0.02
	assign b =  18'sh0_3333 ; // 0.2
	assign c =  18'sh3_8000 ; // -0.5
	assign d =  18'sh0_051E ; // 0.02
	assign a2 = 18'sh0_041E ; // slightly lower than 0.02
	
	// define and connect the neurons
	// with inhib cross-connections
	Iz_neuron burst1(v1,s1, a1,b,c,d, (I1-((s2)?18'h0_8000:0)), neuronClock, neuronReset);
	Iz_neuron burst2(v2,s2, a2,b,c,d, (I2-((s1)?18'h0_8000:0)), neuronClock, neuronReset);
	
	//output spikes to parallel port
	assign GPIO_0[0] = s1;

The neuron module encapsulates all the numerical stuff. Note that the spike output is only one sample (integration step) wide. Each neuron uses 6 multipilers, so the maximum number of neurons that can run in parallel is about 11 For more that that, the system will have to be serialized or simplified. The module produces the voltage variable and a pulse output. The next 4 parameters control the dynamics of the cell. The next parameter is the current input to the cell.

//////////////////////////////////////////////////
//// Izhikevich neuron ///////////////////////////
//////////////////////////////////////////////////
module Iz_neuron(out,spike,a,b,c,d,I,clk,reset);
	output [17:0] out;
	output spike ;
	input signed [17:0] a, b, c, d, I;
	input clk, reset;
	
	wire signed	[17:0] out;
	reg spike;
	reg  signed	[17:0] v1,u1;
	wire signed	[17:0] u1reset, v1new, u1new, du1;
	wire signed	[17:0] v1xv1, v1xb;
	wire signed	[17:0] p, c14;
	
	assign p =  18'sh0_4CCC ; // 0.30
	assign c14 = 18'sh1_6666; // 1.4
	
	assign out = v1;  //copy state var to output
	
	always @ (posedge clk) 
	begin
		if (reset) 
		begin	
			v1 <= 18'sh3_4CCD ; // -0.7
			u1 <= 18'sh3_CCCD ; // -0.2
			spike <= 0;
		end	
		else 
		begin
			if ((v1 > p)) 
			begin 
				v1 <= c ; 		
				u1 <= u1reset ;
				spike <= 1;
			end
			else
			begin
				v1 <= v1new ;
				u1 <= u1new ; 
				spike <= 0;	
			end
		end 
	end
	
	// dt = 1/16 or 2>>4
	// v1(n+1) = v1(n) + dt*(4*(v1(n)^2) + 5*v1(n) +1.40 - u1(n) + I)
	// but note that what is actually coded is
	// v1(n+1) = v1(n) + (v1(n)^2) + 5/4*v1(n) +1.40/4 - u1(n)/4 + I/4)/4
	signed_mult v1sq(v1xv1, v1, v1);
	assign v1new = v1 + ((v1xv1 + v1+(v1>>>2) + (c14>>>2) - (u1>>>2) + (I>>>2))>>>2);
	
	// u1(n+1) = u1 + dt*a*(b*v1(n) - u1(n))
	signed_mult bb(v1xb, v1, b);
	signed_mult aa(du1, (v1xb-u1), a);
	assign u1new = u1 + (du1>>>4)  ; 
	assign u1reset = u1 + d ;
endmodule

The scope image below showing the voltages of the two neurons was taken with the input current to each of them set to 18'h0_2000.

Optimizing the neuron module and adding a synapse
The neural model above uses 6 out of 70 hardware multipliers for each neuron. Looking at the parameter range given in the image from Izhikevich's web site suggests that the a and b parameters don't have to be very accurate for the model to work. The two multiplies by constants were replaced by shifts to approximate the constants within a factor of two. With this change, the number of multipliers drops to two/neuron so that it is possible to put about 35 model neurons on the FPGA. The top-level module was further optimized by introducing a more realistic model of the neuron interconnect (see Ruben Guerrero-Rivera reference below) called a synapse. The new synapse module takes a pulse output from another neuron and uses the impulse (or spike) as a gate to load the initial condition of an exponential decay unit. The scope trace below shows the voltage variable from one neuron, along with the output of the synapse unit driven from the same neuron spike output. You can see that the negative weignt (18'sh3_8000) causes a negative current burst at the time of each action potential.

In the top-level module code shown below, the clock rate has been lowered to the point where individual spikes can be seen on blinking LEDs. Changing the current by flipping switches causes a variety of dynamic behaviors. A switch setting of about 16'h0180 yields alternating bursts. Higher setting cause high-rate synchronous firing.

// clock divider with: reg [15:0] count ;
	always @ (posedge CLOCK_50) 
	begin
		count <= count + 1;
	end
	assign neuronClock = (count==0);	
	
	assign I1 = {2'h0,SW[15:0]} ; //the input bias currents
	assign I2 = {2'h0,SW[15:0]}+8 ; 
	
	assign neuronReset = ~KEY[3]; //the reset pushbutton
	
	//burster parameters
	assign a1 = 6 ;  // 0.016
	assign b =  2 ;  // 0.25
	assign c =  18'sh3_8000 ; // -0.5
	assign d =  18'sh0_051E ; // 0.02
	assign a2 = 6 ; 
	
	// define and connect the neurons
	Iz_neuron neu1(v1,s1, a1,b,c,d, neu1In, neuronClock, neuronReset);
	Iz_neuron neu2(v2,s2, a2,b,c,d, neu2In, neuronClock, neuronReset);
	// with inhibitory cross-connections
	synapse   syn1(neu1In, s2,18'sh3_f000, 1,I1, 0,0, neuronClock, neuronReset);
	synapse   syn2(neu2In, s1,18'sh3_f000, 1,I2, 0,0, neuronClock, neuronReset);
	
	//output spikes to parallel port
	assign GPIO_0[1] = s1;
	assign GPIO_0[2] = s2;
	// output spikes to the leds
	assign LEDR[1] = s1;
	assign LEDR[2] = s2;

////////////////////////////////////////
//// Synapse ///////////////////////////
////////////////////////////////////////
// Acts to create an exponentially falling current at the output 
// from a spike input and a weight which can be +/-
// up to three spike inputs are defined, each with its own weight
module synapse(out,s1,w1,s2,w2,s3,w3,clk,reset);
	output [17:0] out; 				//the simulated current
	input s1,s2,s3;     			// the action potential inputs
	input signed [17:0] w1,w2,w3;   //weights
	input clk, reset;
	
	wire signed	[17:0] out;
	reg  signed	[17:0] v1 ;
	wire signed	[17:0] v1new ;
	
	assign out = v1;  //copy state var to output
	
	always @ (posedge clk) 
	begin
		if (reset) v1 <= 18'sh0; // -0.7	
		else v1 <= v1new ;		 
	end
	
	// dt = 1/16 or 2>>4 and tau=1
	// v1(n+1) = v1(n) + dt/tau*(-v1(n)+s1*w1+s2*w2+s3*w3)
	assign v1new = v1 + ((-v1)>>>4) + (s1?w1:0)+(s2?w2:0)+(s3?w3:0) ;
	
endmodule

//////////////////////////////////////////////////
//// Izhikevich neuron ///////////////////////////
//////////////////////////////////////////////////
// Modified to eliminate two of the multiplies
// Since a nd b need not be very precise 
// and because 0.02<a<0.2 and 0.1<b<0.3
// Treat the a and b inputs as an index do determine how much to 
// shift the intermediate results
// treat the a input as value>>>a with a=1 to 7
// treat the b input as value>>>b with b=1 to 3

module Iz_neuron(out,spike,a,b,c,d,I,clk,reset);
	output [17:0] out; 				//the simulated membrane voltage
	output spike ;     				// the action potential output
	input signed [17:0] c, d, I;   //a,b,c,d set cell dynamics, I is the input current
	input [3:0] a, b ;
	input clk, reset;
	
	wire signed	[17:0] out;
	reg spike;
	reg  signed	[17:0] v1,u1;
	wire signed	[17:0] u1reset, v1new, u1new, du1;
	wire signed	[17:0] v1xv1, v1xb;
	wire signed	[17:0] p, c14;
	
	assign p =  18'sh0_4CCC ; // 0.30
	assign c14 = 18'sh1_6666; // 1.4
	
	assign out = v1;  //copy state var to output
	
	always @ (posedge clk) 
	begin
		if (reset) 
		begin	
			v1 <= 18'sh3_4CCD ; // -0.7
			u1 <= 18'sh3_CCCD ; // -0.2
			spike <= 0;
		end	
		else 
		begin
			if ((v1 > p)) 
			begin 
				v1 <= c ; 		
				u1 <= u1reset ;
				spike <= 1;
			end
			else
			begin
				v1 <= v1new ;
				u1 <= u1new ; 
				spike <= 0;	
			end
		end 
	end
	
	// dt = 1/16 or 2>>4
	// v1(n+1) = v1(n) + dt*(4*(v1(n)^2) + 5*v1(n) +1.40 - u1(n) + I)
	// but note that what is actually coded is
	// v1(n+1) = v1(n) + (v1(n)^2) + 5/4*v1(n) +1.40/4 - u1(n)/4 + I/4)/4
	signed_mult v1sq(v1xv1, v1, v1);
	assign v1new = v1 + ((v1xv1 + v1+(v1>>>2) + (c14>>>2) - (u1>>>2) + (I>>>2))>>>2);
	
	// u1(n+1) = u1 + dt*a*(b*v1(n) - u1(n))
	assign v1xb = v1>>>b;         //mult (v1xb, v1, b);
	assign du1 = (v1xb-u1)>>>a ;  //mult (du1, (v1xb-u1), a);
	assign u1new = u1 + (du1>>>4) ; 
	assign u1reset = u1 + d ;
endmodule

%FitzHugh-Nagumo Euler version
clear all
clf
dt=.01;
t=0:dt:200;
v = zeros(size(t));
w = zeros(size(t));
v(1) = -0.870;
w(1) = -0.212;
current = .57;
for i=1:length(t)-1 
    v(i+1) = v(i) + dt*(v(i) - v(i)^3 - w(i) + current)  ;
    w(i+1) = w(i) + dt*(.08*(v(i) + 0.7 - 0.8*w(i)));
end
plot(t,v)

//analog simulation of FHN system
	always @ (posedge CLOCK_50) 
	begin
		count <= count + 1; 
		
		if (KEY[3]==0) //reset
		begin
			v1 <= 18'h3_2148 ; //v1(0) = -0.870;
			w1 <= 18'h3_C9BB ; //w1(0) = -0.212;
			v2 <= 18'h3_2148 ; //v2(0) = -0.870;
			w2 <= 18'h3_C9BB ; //w2(0) = -0.212;
			c <= {SW[15:0],2'h0} ;    // current = 0.35 = 0_5999
			onehalf <= 18'h0_8000;
		end
		else if (count==0)
		begin 
			v1 <= v1new ;
			w1 <= w1new ;
			v2 <= v2new ;
			w2 <= w2new ;
		end
	end
	
	//FitzHugh-Nagumo Equations
	// dt = 2>>6
	//v(i+1) = v(i) + dt*(v(i) - v(i)^3 - w(i) + current)  ;
    //w(i+1) = w(i) + dt*(1/16*(v(i) + 0.5 - 1*w(i)));
	//note that the 0.5 was omitted below
	
	//cell 1
	signed_mult v_12(v12, v1, (v1>>>1));  	// v-squared/2
	signed_mult v_13(v13, v12, (v1>>>1)); 	// v-cubed/4
	assign v1new = v1 + (((v1>>>2)-v13-(w1>>>2)+(c>>>2)-(v2>>>4))>>>4); //scale other vars to match v3
	assign w1new = w1 + (((v1>>>1)-(w1>>>1)) >>>9); // mult by a factor of 1/16 * dt
	//cell 2
	signed_mult v_22(v22, v2, (v2>>>1));  	// v-squared/2
	signed_mult v_23(v23, v22, (v2>>>1)); 	// v-cubed/4
	assign v2new = v2 + (((v2>>>2)-v23-(w2>>>2)-(v1>>>4))>>>4); //scale other vars to match v3
	assign w2new = w2 + (((v2>>>1)-(w2>>>1)) >>>9); // mult by a factor of 1/16 * dt

//////////////////////////////////////////////////
//// signed mult of 3.24 format 2'comp////////////
//////////////////////////////////////////////////
module signed_mult (out, a, b);
	output 		[26:0]	out;
	input 	signed	[26:0] 	a;
	input 	signed	[26:0] 	b;
	wire	signed	[26:0]	out;
	wire 	signed	[53:0]	mult_out;
	assign mult_out = a * b;
	assign out = {mult_out[53], mult_out[50:24]};
endmodule

Second order system controlled by NiosII.
It is often useful to be able to control an analog computer using a digital computer. Such a system is called a hybrid computer. For example you might want to sequence through a set of input frequencies applied to the simulated system to build a Bode plot. A NiosII cpu was built on the FPGA with ports to (1) control the DDS frequency, (2) control the gain of the sinewave applied to a second order system, (3) control the analog reset and start the simulation, (4) record the amplitude and phase shift of the system under test. The top-level verilog module contains the DDA simulation of the second order system and the NiosII cpu. The GCC program running on the NiosII:

1. Holds the DDA in reset, and initializes the test frequency
2. Releases the DDA reset and waits 10 cycles of the test frequency
3. Waits until the output phase and magnitude reach steady-state
4. Prints the phase and magnitude
5. Iincrements the frequency by a factor
6. Repeats the above steps for a number of frequencies

An mpeg of the frequency scan shown on an oscilloscope shows the input on the top trace and output on the bottom trace. You can easily see the dramatic increase in amplitude near resonance. You can also detect the phase shift. The results printed by the NiosII are available as text for a frequency sweep. The amplitude is in output/input ratio, phase is degrees. Note that the amplitude is not very accurate at low amplitudes because the input gain was not modulated, which causes quantization errors. The results are combined into a magnitude and phase plot with this matlab routine. The vertical red line indicates the computed resonant frequency for the linear system. the red curves are the analytical solution of a second order system with the same natural frequency and damping. The phase should be -90 degrees at the red line.

The whole project is zipped here.

Second order system with modularized integrator
A clearer version of the above project can be made by factoring the integration code into a separate module. The top-level module seems more like wiring an analog computer. A bit of the top-level module shows the wiring.

	// wire the integrators
	// time step: dt = 2>>9
	// v1(n+1) = v1(n) + dt*v2(n)
	integrator int1(v1, v2, 0,9,AnalogClock,AnalogReset);
	
	// v2(n+1) = v2(n) + dt*(-k/m*v1(n) - d/m*v2(n))
	signed_mult K_M(v1xK_M, v1, 18'h1_0000); //Mult by k/m
	signed_mult D_M(v2xD_M, v2, 18'h0_0800); //Mult by d/m
	//scale the input so that it does not saturate at resonance
	signed_mult Sine_gain(Sinput,{sine_out[15],sine_out[15],sine_out}, {2'h0,Sgain});
	integrator int2(v2, (-v1xK_M-v2xD_M+Sinput), 0,9,AnalogClock,AnalogReset);

/////////////////////////////////////////////////
//// integrator /////////////////////////////////
/////////////////////////////////////////////////
module integrator(out,funct,InitialOut,dt,clk,reset);
	output [17:0] out; 				//the state variable V
	input signed [17:0] funct;      //the dV/dt function
	input [3:0] dt ;				// in units of SHIFT-right
	input clk, reset;
	input signed [17:0] InitialOut; //the initial state variable V
	
	wire signed	[17:0] out, v1new ;
	reg signed	[17:0] v1 ;
	
	always @ (posedge clk) 
	begin
		if (reset==0) //reset	
			v1 <= InitialOut ; // 
		else 
			v1 <= v1new ;	
	end
	assign v1new = v1 + (funct>>>dt) ;
	assign out = v1 ;
endmodule
//////////////////////////////////////////////////

References

• Ruben Guerrero-Rivera, et al, Programmable Logic Construction Kits for Hyper-Real-Time Neuronal Modeling, Neural Computation, Volume 18 , Issue 11 (November 2006) Pages: 2651 - 2679
• Eugene M. Izhikevich , Simple Model of Spiking Neurons, IEEE Transactions on Neural Networks (2003) 14:1569- 1572