I was re-reading my reply to Chris Chatman's post about biologically accurate modeling and realized that I was focusing on the wrong aspect of his message. Chris is right. There are so many different features that may need to be modeled to achive the result we want; yet every new feature complicates the model, makes it run slower, and makes it more difficult to analyze the behavior of the model and tweak it further. So, what can we ignore?

  1. Biologically plausible model(s) for different types of neurons (like the model proposed by Eugene Izhikevich)
    1. Number of neuron types to model; Can we simply use two types: regular spiking (RS) for all excitatory and fast spiking (FS) for all inhibitory neurons. Probably not.
    2. Model time step; Is 1ms enough? Some of the types modeled with Izhikevich model require 0.1ms time step.
  2. Axonal delays
  3. Synaptic plasticity and restructuring
  4. Synaptic connectivity
  5. Synaptic noise
  6. Synaptic currents
  7. Brain modularity
  8. Size of the model; Number of neurons? Number of connections? Ratio of excitatory and inhibitory neurons? Different in different modules?
  9. Sensory input / motor output; How may different types of sensors? How sensitive/complex? Attention mechanism?
  10. Developmental aspects
  11. Emotions
  12. Various types of neurotransmitters and neuromodulators; neuron geometry and localized effects of some neurotransmitters/modulators
  13. Real-time/on-line vs. static/off-line
  14. Sleep
  15. Dendritic geometry (including dendrite-to-dendrite and axon-to-axon connections)
  16. Glial and other supporting cells
  17. Wired-in vs. emergent behavior
  18. Reward/punishment

Mind Hacks posted a link to an interesting paper (GRAY-AND-WHITE-MATTER-IN-THE-BRAIN) that aims to answer a question "What's the reason for brain segregation into white and gray matter?". According to the authors of the paper:

... the optimal design [of the brain] depends on the number of neurons, interneuronal connectivity, and axon diameter. In particular, the requirement to connect neurons with many fast axons drives the segregation of the brain into white and gray matter.

One of the authors, Dmitri Chklovskii, co-authored several other papers on the same subject (Network Motifs: Simple Building Blocks of Complex Networks and Geometry and Structural Plasticity of Synaptic Connectivity). The second paper has a brief summary of three categories of synaptic plasticity and different timeframes associated with those categories (the paper itself focuses on the third category and its implications for the structure of the brain):

  • Changes in pre-existing synapses without alterations of interneuronal connectivity that can be realized within a minute
  • Formation of new spines that takes tens of minutes, and
  • Major remodeling of dendritic and axonal branches, usually occuring on the time scale of days.

Dmitri Chklovskii's page also has information on how optimization theory can be applied to brain design with wiring economy principle being one example of a successful optimization hypothesis.

There is also a paper that presents the map of the whole cortical circuit; the connectivity model from this paper was used in the simulation of a model with the number of neurons there is in a human brain.

Chris Chatman on his weblog portraits a gloomy picture awaiting biologically accurate neural modeling:

By most estimates, there are 100 billion neurons in the brain. Some neurons are known to have more than 1,000 dendrites, and up to about 1,000 different branchings of their axons. There are some 50 known neurotransmitters, and who knows how many other neuromodulators may exist (hormones, neural growth factors, neurosteroids). There are also many different receptor types for each of the neurotransmitters. A conservative estimate of the number of interactions you'd have to model to be biologically accurate is somewhere around 225,000,000,000,000,000 (225 million billion).

While this is definitely true, I'm almost equally split between two positions. On one side, the situation may be even more hopeless if you agree with Eugene Izhikevich's polychronization theory:

...spiking networks with delays have more groups than neurons. Thus, the system has potentially enormous memory capacity and will never run out of groups, which could explain how networks of mere 1011 neurons (the size of the human neocortex) could have such a diversity of behavior. (p.270)

At the same time, I would not try to set a goal of building/modeling/simulating something that can exhibit general purpose intelligence; after all, you don't expect to have a conversation about weather with one of the participants of the Grand Challenge. I'm more interested in building brains that demonstrate special purpose intelligence. After all, there are "only" 20000 nerve cells in Aplysia, which is a much more manageable size to model.

I'm keenly interested in modeling the smallest and simplest brain that shows signs of intelligence (defined in Cotterill's terms).

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.

Out of several options I looked at for visualizing results of a simulation OpenGL/GLUT combination looks the most promising. In additional to being able to see the simulation in real-time, I'd also like to have the following:

  • Network activity; firings per second; separate for inh and exc neurons? separate for different types of neurons?
  • Neuron connections and delays; would show incoming/outgoing connections and their delays; spike propagation
  • Neuron activity; would include the input current (I), the recovery current (u), and the membrane potential (v); would generate something similar to this digram only in real-time
  • Synaptic weights; would show current and previous value of all synaptic connections (incoming, outgoing, or both?) of a neuron
  • Neuron properties; to view/edit properties of individual neurons
  • Neural complexes; to research polychronous groups

Neural Viewer already provides most of what I need (it also has help, comments, console, and network selection screens), but, unfortunately, the source code for it is not available.

To implement this in OpenGL/GLUT in addition to some basic functions I would need to at least know the following:

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):