In collaboration with Robert Thijs Kozma

Robert and I thought that it would be nice to finally define what the main building block of what we are talking about is! What are neurons, and how do they work? How do these relatively simple processing elements give rise to higher perceptual and cognitive functions? We are not going to answer these big questions in this post, but we have to start somewhere…. Let’s take a closer look at what is a neuron, how a simple mathematical model can capture a remarkable spectrum of neuron’s behavior, and let’s look at some simple MATLAB code that would allow neurdons to run a neuron in MATLAB at the end of this post.

Neurons are cells that make up most of the nervous system. An adult human brain weights approximately 1.3 Kg and includes 100 billion neurons. Morphologically, nerve cells have three main components: the soma, or the cell body, the axon and the dendrites. Dendrites collect impulses from other neurons, and broadcast these pulses to the soma where, if the generated potential is high enough, a spike is triggered. The neuron’s axon then communicates signals to other neurons (approximately, each neuron is connected to other 10,000 neurons). A significant portion of the brain is taken up by glia cells, which contribution to signal processing is still under debate.

Also, glia may have a great deal of functional significance, albeit it hasn’t been characterized too well.

Neurons process and transfer information in the form of electrical pulses (spikes) traveling down their axons and resulting in neurotransmitter release of excitatory, inhibitory or modulatory neurotransmitters at the synapse with another neuron, usually at the receiving cell’s dendrites or soma. There are three broad classes of neurons based on their effect on other neurons’ membrane potential: excitatory, inhibitory and modulatory neurons. Complex networks of these neurons form the computational backbones of the nervous system. Whereas excitatory inputs, summing at the cell body, bring the cell closer to its firing threshold, inhibitory cells bring the cell further from its firing threshold. Excitation and inhibition are the main ingredient of cell communication and information processing in the nervous system, where processes such as cell competition and short and long term memory depend on a delicate balance between these two components. Neuromodulation acts by modifying this balance: neurotransmitters such as dopamine, acetylcholine and serotonin modify cell’s excitability, broadening a cell’s behavior in a state-dependent way.

Research in neuroscience has progressively captured more and more details on the morphology and biophysics of neurons. At the same time, computational neuroscientists have synthesized scores of data and cell recordings in computational models that are simple enough to be implemented in personal computers. Since it is believed that neurons’ main communication “alphabet” is trains of spikes, which are caused by fluctuations in the cell membrane potential, computational neuroscientists have focused on mathematically characterizing these electrical fluctuations.

In this post, we will look at a simple and powerful model proposed by Eugene Izhikevich (http://www.izhikevich.org/publications/spikes.pdf). This mathematical model is the most recent of a long tradition of models used to study individual neurons that display spiking/bursting behavior. Examples of these models include Hodgkin-Huxley, and FitzHugh-Nagumo, just to name a few. The Izhikevic model (2003) is particularly interesting because it is a compact model that, via suitable changes of parameters, can simulate a large array of neurons’ behavior.

The Izhikevich model consists of a system of two differential equations. The two variables in the system are the following: v denotes the membrane potential of a neuron in millivolts (mV), while u denotes the generic recovery variable. There are also five free parameters in the equation:

I the external input to the neuron, such as those from synaptic inputs
a the rate of recovery of u
b the ‘sensitivity of recovery to subthreshhold fluctuations of membrane potential’
c, d are the after spike resets of v and u respectively, i.e. the values to which they get set back after a spike occurs.

The long awaited Izhikevich equation is as follows:

v’=0.04*v^2+5*v+140-u+I
u’=a*(b*v-u)

with the following reset conditions:

if v>=30 then:
v:=c
u:=u+d

(i.e. if v, the membrane potential, becomes greater than 30 mV, or equivalently, the neuron spikes, then the variables get bumped back, so the neuron has a chance to spike again). The video below shows an array of 50×50 neurons which receive as in input (the variable I) the grey level of the so-called “outstein”, the logo of the Department of Cognitive and Neural Systems at Boston Univeristy (check out the brightness illusion that this display causes here: http://cns-web.bu.edu/)

In our attached code, we will be solving these equations numerically using Euler’s method, so they appear in slightly different form. By changing the parameters of the equation, it is possible to reproduce all the spiking behaviors shown in the figures above. But be careful, if you randomly change the parameters, your model may no longer have any biological meaning!

You can download the MATLAB (.m file) code, or just copy and paste the code below in a new MATLAB document.

NOTE: download the code here for completeness (the input matrix is truncated for space issue)

——————- MATLB CODE ————–
% Simulation of a simple reduced spiking equation model.
% The formalism is based on Izhikevich (http://izhikevich.org/index.htm)
% The equation:
% v’=0.04*v^2+5*v+140-u+I
% u’=a*(b*v-u)
% with the following reset conditions:
% if v>=30 then:
% v<-c
% u<-u+d
%
% (c) 2009 Massimiliano Versace and Robert Kozma
% CNS Technology Lab, Boston University
% http://techlab.bu.edu

close all, clear all, clc

%% Initialize data and parameters
% The input image is a 50×50 pixel “outstein” image (http://cns-web.bu.edu/)
input_image = [0,0,0,0,0,0,0,0,4,10,0,0,5,9,0,0,0,4,2,0,11,0,6,5,0,13,0,0,4,10,5,0,4,1,0,0,0,1,1,1,3,3,3,2,0,0,0,1,0,0;0,0,0,0,0,0,0,0,2,0,4,15,6,0,0,7,1,0,7,5,0,5,0,0,240,244,212,14,5,0,0,0,0,0,2,4,5,4,5,5,0,0,1,0,0,0,0,1,0,0;0,0,0,0,0,0,0,0,5,0,1,5,0,0,11,8,5,3,0,0,0,11,11,195,250,252,240,0,16,16,12,8,0,0,3,3,0,0,0,0,0,1,1,1,1,1,0,0,0,0;0,0,0,0,0,0,0,0,0,1,8,12,9,6,3,0,0,4,0,9,0,0,0,197,244,245,238,8,3,0,0,4,4,2,0,0,0,0,0,1,4,3,3,3,2,1,0,0,0,0;0,0,0,0,0,0,0,0,7,8,0,0,0,0,0,11,9,0,0,12,0,4,0,8,250,244,234,6,8,1,26,0,3,0,0,0,0,4,4,4,3,1,0,0,2,2,0,0,0,0;0,0,0,0,0,0,0,0,0,9,4,0,12,21,8,0,7,0,0,2,0,1,14,0,5,8,0,9,0,19,0,4,0,0,0,0,1,2,0,0,0,0,0,0,0,2,2,0,0,0;0,0,0,0,0,0,0,0,1,0,12,13,0,0,4,0,0,9,0,6,9,0,0,20,229,253,231,3,0,0,19,0,0,2,2,0,0,0,7,13,7,3,0,0,0,2,2,0,0,0;0,0,0,0,0,0,0,0,11,14,142,235,136,11,0,3,4,4,0,2,0,6,7,0,238,255,217,9,2,4,4,9,0,2,2,0,0,9,32,51,17,11,4,1,2,3,1,0,0,0;0,0,0,0,0,0,1,1,0,98,253,249,255,118,7,0,9,0,0,16,14,0,0,4,239,229,85,1,0,11,0,2,6,0,0,7,6,118,255,247,140,22,7,13,0,4,1,4,0,0;0,0,0,0,0,0,1,1,6,106,239,255,233,128,0,9,11,7,17,0,0,26,9,0,95,251,0,0,8,0,6,4,0,20,18,0,0,151,235,255,255,89,3,0,6,11,0,3,0,0;0,0,0,0,0,0,1,1,0,39,144,247,160,146,167,13,0,0,2,11,1,5,0,30,57,254,0,16,0,3,1,0,3,0,8,0,21,119,255,241,235,85,6,0,4,0,0,0,0,0;0,0,0,0,0,0,1,1,4,22,0,3,130,254,232,108,3,19,3,0,19,13,2,0,63,223,22,2,0,15,0,7,4,17,0,109,255,130,69,70,35,1,0,0,3,8,3,17,0,0;0,0,0,0,0,0,1,1,0,3,0,0,124,247,254,252,85,0,0,0,5,0,7,15,61,240,1,0,11,0,0,7,3,0,112,247,252,255,119,0,0,8,15,4,1,3,0,0,0,0;0,0,0,0,0,0,1,1,13,15,14,17,0,153,255,245,119,0,17,16,2,6,12,0,82,243,0,14,1,0,4,0,11,2,92,255,248,118,4,0,0,0,0,2,5,8,7,5,0,0;0,0,0,0,0,0,1,1,0,0,8,0,0,5,0,159,233,96,0,0,4,0,0,8,73,215,6,13,0,12,5,0,1,81,255,132,8,15,0,2,6,1,2,7,0,0,0,0,0,0;0,0,0,0,0,0,1,1,14,0,0,0,18,2,7,6,162,255,66,0,13,2,0,5,62,255,12,0,7,0,0,6,77,251,163,9,0,0,8,0,2,0,2,6,0,0,11,0,0,0;0,18,0,0,7,13,0,0,0,3,0,12,2,0,0,9,2,179,250,52,10,0,11,0,71,240,6,0,2,4,3,66,253,157,0,0,12,0,0,6,13,0,0,1,0,6,13,0,7,6;8,0,0,28,0,0,3,24,13,17,1,0,1,11,0,0,6,0,180,255,64,3,0,8,83,230,0,19,0,0,71,247,169,7,1,7,0,2,1,0,0,23,6,7,12,0,10,18,6,5;15,2,0,0,6,1,1,0,0,0,2,0,4,0,2,18,0,8,0,185,244,22,29,0,56,236,6,0,16,62,243,194,13,0,14,1,12,3,5,22,15,0,0,8,0,0,28,0,3,2;7,0,0,27,0,16,15,6,15,6,8,3,19,0,0,0,16,3,0,10,207,252,33,2,55,255,11,8,37,240,206,0,0,14,0,0,0,9,6,0,2,9,12,9,5,0,0,10,0,0;0,8,3,0,5,0,0,0,0,0,5,0,11,4,8,4,3,8,0,5,0,222,237,28,21,0,0,28,250,190,4,7,10,4,0,4,7,20,0,10,7,0,0,8,17,12,0,0,2,1;0,5,6,0,3,1,0,7,7,0,0,0,9,0,0,12,7,0,2,4,4,8,3,3,0,12,0,58,213,29,0,0,0,5,0,0,11,0,0,6,0,7,0,94,124,106,127,15,7,5;3,5,0,12,42,244,200,0,0,14,11,3,0,4,0,1,1,20,2,5,0,0,1,0,24,11,0,0,0,0,4,12,23,7,0,12,0,0,34,0,12,36,131,147,216,255,134,101,10,7;1,0,16,26,255,255,236,176,97,238,255,133,3,3,8,0,0,2,0,3,6,10,3,0,0,4,0,9,1,28,74,56,57,62,69,68,76,77,62,181,243,255,209,224,252,255,249,115,9,5;2,13,10,32,248,255,255,179,88,249,254,240,255,244,255,254,255,254,246,255,241,116,0,4,0,0,11,4,22,87,255,253,254,241,255,244,252,253,219,255,253,250,201,242,248,253,253,68,13,5;6,0,0,65,251,249,163,8,91,241,190,44,0,1,19,0,8,0,22,1,0,9,20,3,3,15,0,15,3,11,0,8,0,10,0,10,0,0,37,0,50,188,220,237,237,242,252,144,14,6;12,19,10,6,0,0,15,0,0,11,0,1,20,0,0,13,0,12,0,15,0,0,208,14,0,0,8,0,9,18,7,2,0,21,0,0,11,0,0,10,28,55,101,141,85,62,123,108,11,5;0,0,0,11,0,3,17,1,3,10,13,7,11,12,0,0,13,10,0,0,17,201,255,32,10,1,4,40,251,170,24,8,16,0,16,14,0,4,3,5,15,18,3,1,0,0,3,14,5,1;5,10,7,11,11,0,0,0,3,0,3,3,0,0,0,24,14,0,14,1,207,255,19,0,2,50,12,0,104,255,189,0,0,0,5,11,0,2,21,0,0,13,10,0,3,20,6,0,2,0;0,5,4,0,13,2,1,12,7,4,10,2,1,16,7,0,0,18,0,180,255,27,1,9,6,241,2,16,3,84,255,175,20,0,0,0,12,14,0,9,15,0,14,24,0,0,10,5,2,0;5,3,0,0,0,6,0,0,0,0,0,2,0,0,5,1,0,5,178,255,47,11,13,0,0,239,8,0,0,13,107,255,163,20,4,13,0,7,0,7,0,13,7,0,13,9,0,1,3,1;0,0,2,12,0,9,2,0,9,5,0,4,6,0,4,0,6,174,250,47,14,1,2,4,0,255,0,14,4,1,12,105,255,147,8,1,8,0,14,0,12,0,0,7,0,0,16,0,3,1;0,0,0,0,0,0,0,0,0,8,0,8,1,9,0,18,153,241,78,12,0,0,12,0,2,244,0,3,11,0,19,12,117,252,154,0,4,4,1,4,4,0,0,7,3,0,5,2,0,0;0,0,0,0,0,0,0,0,10,4,12,0,3,15,0,136,255,97,0,0,19,9,0,10,0,245,0,1,0,0,0,5,35,109,251,176,0,0,11,7,3,7,9,12,8,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,7,23,133,240,255,95,0,17,10,0,0,0,0,0,245,9,0,5,7,0,11,0,22,126,244,165,48,0,0,2,3,0,0,5,9,1,0,0,0;0,0,0,0,0,0,0,0,10,9,0,8,124,237,255,227,73,11,3,0,6,6,22,4,0,240,8,0,11,0,3,2,10,0,117,244,254,235,127,8,6,10,9,1,7,8,0,0,0,0;0,0,0,0,0,0,0,0,0,12,36,55,68,162,255,93,1,6,2,7,5,0,0,6,0,247,2,0,7,0,3,0,0,7,19,126,245,160,39,5,0,0,10,8,6,3,0,7,0,0;0,0,0,0,0,0,0,0,6,79,255,246,152,42,0,13,2,0,1,2,0,23,10,7,0,237,0,11,0,8,2,5,3,16,0,116,175,71,156,247,99,22,0,0,3,12,0,5,0,0;0,0,0,0,0,0,0,0,3,110,242,247,255,116,6,4,16,0,6,14,3,7,0,0,13,204,14,0,0,4,0,2,0,0,11,0,2,129,250,255,246,89,11,0,5,15,0,0,0,0;0,0,0,0,0,0,0,0,12,83,255,251,160,29,12,0,0,24,6,0,3,6,0,9,57,199,68,0,21,3,0,4,4,9,1,10,3,125,252,248,251,80,9,2,0,3,0,7,0,0;0,0,0,0,0,0,0,0,0,9,33,48,37,12,0,0,1,0,1,6,3,0,0,0,251,188,85,9,0,8,4,0,0,1,5,2,2,46,153,253,92,13,0,9,3,0,3,0,0,0;0,0,0,0,0,0,0,0,5,3,6,10,6,0,0,6,0,4,5,0,0,0,4,12,227,185,98,0,13,9,0,11,2,0,0,1,1,0,6,15,30,0,0,9,0,0,4,2,0,0;0,0,0,0,0,0,0,0,3,0,0,0,1,0,0,2,0,6,0,0,4,4,0,0,6,17,0,0,3,0,12,0,8,7,3,2,5,7,0,0,0,2,7,7,2,0,1,0,0,0;0,0,0,0,0,0,0,0,0,0,0,5,8,4,0,0,0,5,0,0,11,2,0,14,54,254,46,11,4,4,1,0,0,0,0,0,0,0,6,11,7,14,6,0,6,5,0,0,0,0;0,0,0,0,0,0,0,0,4,0,0,0,0,0,0,3,0,5,1,0,0,0,9,69,245,241,247,0,0,0,9,2,13,6,4,8,7,0,1,8,4,8,0,0,8,4,0,4,0,0;0,0,0,0,0,0,0,0,1,2,2,0,0,0,0,2,0,0,0,4,2,0,7,57,241,229,236,6,7,0,0,7,6,0,0,0,0,0,0,12,0,3,1,3,10,0,0,10,0,0;0,0,0,0,0,0,0,0,0,1,5,6,5,2,0,0,3,0,0,3,6,0,0,4,244,255,229,0,9,20,0,4,0,1,11,13,7,1,1,4,0,5,7,8,8,0,0,3,0,0;0,0,0,0,0,0,0,0,1,1,0,0,0,2,2,1,0,0,1,0,0,0,2,0,9,0,17,0,0,0,6,5,1,5,1,0,0,6,5,0,0,1,0,0,3,6,1,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;];
Receptor_size = [50, 50]; % Receptors total # of neurons
int_v = 0.2500; % integration rate of variable v
int_u = 0.2500; % integration rate of variable u
inp_ratio = 0.13; % coefficient that multiplies the input RGB value (in this case, it is a single value between 0 and 255)
resting_potential = -70; % the resting potential of the cell
simulation_length = 150; % the simulated time steps
a = 0.02; % parameter a in u’=a*(b*v-u)
b = 0.02; % parameter b in u’=a*(b*v-u)
c = resting_potential; % parameter c in if v>=30 then v <= c
d = 8; % parameter d in if v>=30 then u <= -u+d

%% Receptors variables
Rec_record = resting_potential*ones(Receptor_size(1)*Receptor_size(2),simulation_length); %Initialize Receptors cell voltages
Rec_v = resting_potential*ones(Receptor_size); % Initialise Receptors membrane potential
Rec_u = b*Rec_v; % Initialise Receptors recovery variable

% Create plot
h = figure; set(h, ‘DoubleBuffer’, ‘on’, ‘color’, ‘white’, ‘Name’, ‘Simple spiking model – by www.neurdon.com’);

% Movie
mov=avifile(‘neurdon_post’,'compression’,'indeo5′); %makes test.avi w/ Indeo 5

%In order to run the avi from MATLAB command window:
!test.avi&

%% Simulate for t time steps
for t = 1:simulation_length;

%%%%%%%%%%%%% Receptors %%%%%%%%%%%%%
Rec_fired = find(Rec_v >= 30); % Indices of Receptors spikes
if ~isempty(Rec_fired) % Reseting F1 after firing
Rec_v(Rec_fired) = c;
Rec_u(Rec_fired) = Rec_u(Rec_fired)+ d;
end

Rec_v = Rec_v + int_v*((0.04*Rec_v + 5) .* Rec_v + 140 – Rec_u + input_image * inp_ratio);
Rec_u = Rec_u + int_u .*(a*(b * Rec_v – Rec_u));

% Store Receptors variables
Rec_record(:,t) = reshape(Rec_v, 1, Receptor_size(1)*Receptor_size(2));
Rec_record(Rec_record >= 30)=30; % Normalize spikes

% PLOT
axis on
subplot(2,2,1), imagesc(input_image), title(['Input image, cycle ', num2str(t), ' of ', num2str(simulation_length)]), drawnow
subplot(2,2,2), imagesc(Rec_v>30, [0 1]), drawnow, title(‘Receptor spikes’), colormap gray, drawnow
subplot(2,2,3), plot(Rec_record(1,1:t), ‘black’, ‘linewidth’, 2), drawnow, hold on, title(‘Receptor (x=5,y=23), not excited’)
subplot(2,2,4), plot(Rec_record(273,1:t), ‘linewidth’, 2), drawnow, hold on, title(‘Receptor (x=5,y=23), excited’)

% Movie
f2=getframe(gcf); % gets the gcf
mov=addframe(mov,f2); % adds the frame into mov
end;

mov=close(mov); % closes the mov

% test movie
!neurdon_post.avi&