Reviewing the code of the SPNet++ I realized that the implemented synaptic plasticity mechanism in the simulation is much simpler than I expected. It turned our that the code implements the nearest-neighbor spike model described in RELATING-STDP-TO-BCM article by Eugene M. Izhikevich and Niraj S. Desai. According to the authors of the article it may be sufficient to only consider two postsynaptic spikes -- the one that occurs before and the one that occurs after the presynaptic firing -- while determining whether a synapse should be potentiated or depressed based on recent spike activity. In addition to attempting to unify different forms of plasticity -- spike-timing-dependent plasticity (STDP) and a standard long-term potentiation and depression (LTP/LTD) -- into a single framework, this approach appears to be very computationally efficient.

update 2006/01/14: STDP-BASED-ON-LOCAL-INFORMATION paper presents an alternative rule, which is also described as computationally efficient, yet simple and reliable. This paper also provides five types of gating functions: no gating, presynaptic, postsynaptic, dual OR and dual AND gating.

There is also more information on different types of synaptic plasticity in the post on synaptic connectivity.

After reaching a dead end with neural modeling using Perl and PDL I decided to take another look at Eugene Izhikevich's implementation of a neural simulator described in one of his articles. It turned our that a group in Israel already extended his implementation and added a MATLAB interface as well as support for different kinds of input channels and recorders. Still, the current code a) implements the older model, b) has all neuron parameters hardcoded, and c) doesn't provide real-time Neural-Viewer-like visualization. Let's see how much time fixing a) and b) will take...

During the last several weeks I've been working on implementing a simulation of a simple network based on the simple model of spiking neurons. As I am very familiar with Perl I decided to try to use Perl and the Perl Data Language (PDL) module, which provides capabilities similar to those of MATLAB and IDL. The PDL module supports efficient and compact storage of multidimensional arrays and provides the fundamental operation of numerical linear algebra.

Much to my disappointment, the code that I wrote runs too slow. I'm sure that there are many ways to improve my PDL code (it's my first PDL script), which may significantly increase its performance, but still, it took more than 4 minutes (265 seconds) to simulate 400ms of model time for the network of randomly connected 1000 neurons with 28 spikes per neuron per second. It took 18 seconds to simulate 100ms of model time for 1000 neurons and 189 seconds to simulate 1s of model time for 100 neurons, which means that the implementation doesn't scale well when the number of neurons or the simulation time increases.

To compare these numbers with the performance of the simulation written in C or C++ I downloaded the source code of the program written by Eugene Izhikevich for his model and available on his website. His program compiled under CygWin calculated 1s of model time in less than a second of real time for a network with 1000 neurons with ~10 spikes per neuron per second. Even though two programs implement slightly different models (the C++ code implements the model with 4 parameters covered here and my code implements the model with 7 parameters covered here) and have different number of connections per neuron (100 for Eugene's implementation and 1000 in my implementation), the C code also implements axonal conduction delays and spike-timing-dependent plasticity (STDP), which are not included in my implementation.

This difference of two orders of magnitude means that I'll stay with C++ for now. And, BTW, here is the perl script I was running.

In addition to all good things about the Simple Model of Spiking Neurons, there are few weaknesses that I observed while working on implementation of the model:

• The model appears to be sensitive to large currents (for some types of neurons); there is probably something wrong with my implementation of the model
• The model requires very small time step for some types of neurons (0.05ms for low threshold spiking (LTS) and fast spiking (FS) neurons); there is probably not much that can be done about this, but it definitely creates challenges with simulating those types of neurons as it requires to implement various time steps for different types of neurons to speed up the simulation. What's interesting, the previous version of the model had the same time step (0.5ms) for all types of neurons.

update 2005/12/21: According to the author of the model (private email):

There is indeed a problem with stability when the injected current, I, into a thalamic cell is too strong. You may interpret this as an artifact of the model, or as the sign that such a strong current can kill the cell.

update 2005/12/22: Again, from the author of the model on the same problem (private email):

What happens is that the peak of the spike and the voltage reset variables meet, periodic spiking disappears, and the voltage variable drops. This behavior resembles the "excitation block" behavior of many neurons (no spiking when the input current is too strong), though the mechanism may be different.

update 2005/12/22: And here is the recomendation on how to fix this:

...use the following formula for the after-spike reseting of u: u <- min(u+50,530). This way, u never goes above 530, and the the spike cutoff and afterspike reset values never intersect.

Here is the result and the modified script.

After finishing playing with excitatory neurons I turned to inhibitory neurons and after several unsuccessful attempts realized that I'll probably spend more time getting them right than I was planning on. Even though I tried many different sets of parameters (based on parameters documented in Izhikevich's book Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting) I just couldn't get it to work properly. I sent an email to the author, who kindly provided his MATLAB scripts that he used to model fast spiking (FS), late spiking (LS), and thalamic interneurons (TI) neurons. The reason for most of my problems was quite simple: I was using the wrong time step. For example, this is how the diagram looks for similation of one fast spiking (FS) neuron with 0.5ms time step:

and here is how it looks with the correct time step of 0.05ms:

Even though this fixed most of the problems there were still few I didn't know what to do about. For example, the model for thalamic interneurons was unstable with large currents (even with the current just slightly larger than showed in the book: 258pA vs. 255pA):

I also ended up modifying parameters for the low threshold spiking (LTS) neurons; I took c and d parameters from the model for fast spiking neurons (based on the older version of the model). Here is how the result looks like:

And here is the script that I used to generate all these diagrams.

As I wrote in my last post on neural modeling there are several things that I like about the model of spiking neurons proposed by Eugene Izhikevich (see also WHICH-MODEL-TO-USE for comprehensive comparison with other available models):

1. It's fast: only 13 FLOPS versus 5 FLOPS for Integrate-and-Fire and 1200 FLOPS for the Hodgkin-Huxley model
2. It's biologically plausible: it exhibits the same neuro-computational properties as the most complex, Hodgkin-Huxley model
3. It allows to model various types of neurons by changing four (in the latest revision seven) parameters
4. It has no fixed threshold or absolute refractory period: these are properties rather than parameters of the model and depend on the type of a neuron. Based on the history of the membrane potential prior to the spike, the threshold potential can be as low as -55 mV or as high as -40 mV.

Izhikevich proposes classification of all neuron types as resonator/integrator (based on presence of subthreshold oscillations) and as being bistable/monostable (based on co-existance of resting and spiking states), which, in combination, make four groups. By using this classification he is able to derive several neuro-computational properties of those types of neurons:

• Inhibition impedes spiking in integrators, but can promote it in resonators (this came as news to me; I was under the impression that inhibition always inhibit spiking)
• Integrators have all-or-none spikes while resonators may not
• Integrators have well-defined threshold while resonators may not
• Integrators integrate, resonators resonate. This means integrators prefer high-frequency inputs; the higher the frequency, the sooner they fire. By contrast, the response of the resonator neuron depends on the frequency content of the input

Most cortical pyramidal neurons (including regular spiking (RS), and chattering (CH) neurons shown in the previous post) are integrators. Most cortical inhibitory neurons are resonators. According to the author of the model, a good neuronal model must reproduce not only electrophysiology but also bifurcaiton dynamics of neurons.

It seems logical to start building a neural network framework or a simulator with a neuronal model. Neuronal models and their implementations range from simple (like the model that is implemented in Neural Viewer) to very complex (like the Hodgkin-Huxley model). While there are many models to use (Which model to use for cortical spiking neurons? article provides a good coverage of the existing models) I like the model proposed by Eugene Izhikevich in Simple model of Spiking Neurons as it is computationally efficient and still biologically plausible.

Eugene's website provides MATLAB scripts to play with the models, but since I didn't have access to a MATLAB instance I decided to reproduce the same set of diagrams using Perl. After a couple of weeks of experiments and several emails to the author I was able to model various types of neurons according to the parameters described in Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Here is how the results look for the regular spiking (RS) neuron:

And here are the results for the chattering (CH) neuron: