In addition to all visual features I want, I also want the framework/simulator to be easily configurable and scriptable, and what language provides better parsing and scripting capabilities than Perl? (or Python, or Ruby, but my experience with those is very limited...)

It's (non)surprisingly easy to call perl from C as the following code shows:

#ifdef __cplusplus
  extern "C" {
#   include <EXTERN.h>
#   include <perl.h>
    static PerlInterpreter *my_perl;
  }
#endif

int main(int argc, char **argv, char **env) {

  my_perl = perl_alloc();
  perl_construct(my_perl);

  perl_parse(my_perl, NULL, argc, argv, env);
  perl_run(my_perl);

  perl_destruct(my_perl);
  perl_free(my_perl);

  return 0;
}

This code can be compiled with something like this:

g++ -c interp.c -o interp.o `perl -MExtUtils::Embed -e ccopts`
g++ interp.o -o interp `perl -MExtUtils::Embed -e ldopts`

and then interp -v should give you the same output as perl -v.

While it may be easy to compile and run this simple program, it's still a challenge to pass complex data structures between C and Perl as you need to know all those *_sv, *_pv, *_av and stack manipulation calls. Fortunately for me I own a copy of Extending and Embedding Perl. There are even libraries like PerlStream that will hide all that complexity from you.

Another problem is that it means having additional libraries and module files (if you happen to use Perl modules) in your distribution. On Windows this can be solved by using Thinstall, PEBundle, or one of many other exe packers that compress executable and dll files into a single file.

After completing all the prep work I was ready to make the jump and to modify the code to use the new model instead of the earlier proposed one that SPNet++ is based on.

The core of the change was to replace this code:

v' = 0.04 * v * v + 5 * v + 140 - u + I
u' = a * (b * v - u)

with after spike resetting:

if v >= 30mV then v = c; u = u + d

where variable v represents the membrane potential of the neuron, u represents a membrane recovery variable, I is the input current, and a, b, c, and d are model parameters.

The code above was to be replaced with this (very similar) code:

v' = (k * (v - vr) * (v - vt) - u + I) / C
u' = a * (b * (v - vr) - u)

with after spike resetting:

if v >= vpeak then v = c; u = u + d

where v is the membrane potential, u is the recovery current, C is the membrane capacitance, vr is the resting membrane potential, vt is the instantaneous threshold potential, vpeak is the spike cut-off value, a is the recovery time constant, c is the voltage reset value, and d is the total amount of outward minus inward currents activated during the spike and affecting the after-spike behavior. Though it looks like the model has ten parameters, it is equivalent to the previous version of this model shown above and hence it has only four independent parameters (DYNAMICAL-SYSTEMS-IN-NEUROSCIENCE, pp.161-3, 296)

While this change was easy, there were three other parameters, for which I didn't know what values to use: the input current and the default synaptic strength for excitatory and inhibitory neurons. After playing with various values for these parameters I finished with this result, which I'm not satisfied with as it only remotely resembles the diagram generated by the old code (also generated for the network of 1000 randomly connected neurons). Not sure which one is correct, though.

The first modification to make to implement the new model using existing SPNet++ code was to enable variable timestep, as the existing code only supported one time step (1ms) and the new model required 0.1ms time step for fast spiking (FS) neurons.

After I finished making changes and run the code for the first time the firing rate (number of spikes per neuron per second) dropped from ~7 to almost zero. After spending few minutes looking at the data I realized that by changing the time step from 1ms to 0.1ms I effectively reduced the time that input current influences the neuron and it was not enough to excite it. After a brief search on the Internet I found a paper that describes duration of spikes for excitatory and inhibitory neurons. Based on the data provided in the paper I made the input current to decay; and after applying this logic to synaptic and random thalamic input I almost got the numbers I expected.

Almost, but not quite. I ended up changing the logic that handles synaptic delays as it expected to work with delays expressed in ms and I needed delays to vary depending on the time step.

I also modified the generator of random integers from (rand()%(int)(max1)) to (int(((double)rand())/((double)(RAND_MAX)+1)*(max1))). According to the author of this page this can make a difference, depending on the type of the random number generator used.

Here is the current version of the source code and here is the diagram of the spikes it produces. It shows all spikes during the 50th second of the simulation for 1000 neurons (the top part of the picture shows inhibitory neurons):

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.

From Eugene Izhikevich's website:

On October 27, 2005 I finished simulation of a model that has the size of the human brain. The model has 100,000,000,000 neurons (hundred billion or 10¹¹) and almost 1,000,000,000,000,000 (one quadrillion or 10¹⁵) synapses. It represents 300×300 mm² of mammalian thalamo-cortical surface, specific, non-specific, and reticular thalamic nuclei, and spiking neurons with firing properties corresponding to those recorded in the mammalian brain. The model exhibited alpha and gamma rhythms, moving clusters of neurons in up- and down-states, and other interesting phenomena.
One second of simulation took 50 days on a beowulf cluster of 27 processors (3GHz each). Why did I do that?

This sounds like a great achievement, but I've always been interested in modeling a smallest, simplest, and slowest brain that exhibits intelligence. How many neurons and connections would it have? What modules and components are essential and what are not? How much of intelligence would be lost because of its size and simplicity? How fast can it be trained to do something useful? What properties would need to be wired in and what properties would emerge?