Added neuron model aeif_cond_beta_multisynapse using GSL ODE solver (see issue 438)

[golosio]

Multisynaptic neuron model with synaptic conductance modeled by the sum of two exponential functions, with independent synaptic rise times and decay times, as described in A. Roth and M.C.W. van Rossum in Computational Modeling Methods for Neuroscientists, MIT Press 2013, Chapter 6. The model is an extension of the conductance based exponential integrate-and-fire neuron model according to Brette and Gerstner (2005), as implemented in aeif_cond_alpha_multisynapse
This pull request deal with the issue #438

[heplesser]

@golosio Thank you for your contribution! We are currently trying to properly understand and handle some numerical issues with aeif models in general (see #229). I would prefer to solve them in the simple case of the aeif_cond_alpha model first, before adding new aeif models. In that way, we can ensure that all new models will use the best numerics. I hope that this is okay for you.

Another pre-requisite for merging your model into the NEST master repository is a signed Contributor License Agreement, which you can download from our website and send me by mail.

[golosio]

@heplesser Sure, that's fine. Meanwhile I added the (equivalent) model aeif_cond_beta_multisynapse. I compared the two and they give the same results, which is not surprising as the underlying differential equations have the same solutions. aeif_cond_beta_multisynapse is closer to the original model aeif_cond_alpha_multisynapse, and it is based on the same differential equations as those proposed by dnao in "Create iaf_cond_beta" #405 . so it should be a quicker review, when you feel you can do it. By the way, if you think it can be useful I can help with the revision of iaf_cond_beta. Let me know...

[heplesser]

@golosio I am a bit confused: isn't the beta function precisely the difference-of-two-exponentials? So what precisely is the difference between the _beta_ and the _2exp_ model?

Btw, our preferred integrator for AEIF neurons is at present the GSL ODE solver used in aeif_cond_alpha. Would you consider to convert your code to that solver? In #229, we are currently looking at even better solvers that handle threshold detection as part of the integration.

Concerning contributions to #405, I'd suggest that you contact @dnao directly, since it is his/her code.

[golosio]

Right, the beta function is the difference of two exponential. The two models differ because they are based on a slightly different (but equivalent) system of differential equations, which have the same solution (i.e. the beta function) and it is probably superfluous to have both, but it was useful for me as a check of consistency. Sure, I'll have a look at the GSL ODE solver and see how difficult it would be to adapt my code for it.

[heplesser]

@golosio Good. While cross-validation of models is very useful, we want to keep the number of models in the NEST distribution within sensible limits. So I would suggest that you remove the _2exp_ variant.

[golosio]

Sure, that's easy!

[golosio]

@heplesser I removed the 2exp variant

[golosio]

@heplesser ok, I'm working on the gsl ode version, it should not be too much work if I use aeif_cond_alpha as a base.

[golosio]

@heplesser ok, I adapted my code for using the GSL ODE solver and tested it. You can find the new version in the updated pull request.

[heplesser]

@jougs Do you have any objections to merging this PR now that all tests have passed?

[jougs]

@heplesser: I was just commenting here because of the failing builds, not as a reviewer.

Maybe @ingablundell or @DimitriPlotnikov could have a look?

[DimitriPlotnikov]

@heplesser @jougs my first impression is that the model is implemented well. Before I proceed with more detailed review, I would have an acknowledgment that the model is important enough to be included into NEST repository.

[golosio]

Thank you @DimitriPlotnikov.

[heplesser]

@DimitriPlotnikov Yes, the model is relevant (otherwise I wouldn't have given it a thumb up).

[DimitriPlotnikov]

@heplesser Thank you for your replay. I still a bit confused with this 'kitschy' emojies which can appear at arbitary place in the discussion. I'll tackle it next week.

[DimitriPlotnikov]

Here I need qualified help of @heplesser:

We already discussed the issue with our interpretation of inhibitory and excitatory buffers in the current NEST state. In a model with one buffer (iaf_psc_exp) the logic to use the sigh of the weight to route spike events is at least plausible. In the case of the mutlisynapse model it leads to: a) you always have the same number of inh und exc buffers b) you have to buffers (inh/exc) at the same portnumber c) finally, this implementation is inconsistent with the implementation of the iaf_psc_alpha_multisynapse and iaf_psc_exp_multisynapse.

[golosio]

@DimitriPlotnikov I took as reference aeif_cond_alpha_multisynapse, but I agree with your observations. If you already discussed this issue for the iaf_psc_*_multisynapse models, probably there is no reason to keep two separate buffers for inhibitory and excitatory connections. We have to be careful that in conductance based models the driving force is different for inhibitory vs excitatory signals, but this should not be a big deal. I will try to modify the code and let you know.

[golosio]

@DimitriPlotnikov I modified the code according to your comment. Now there is only one buffer for both types of events, inhibitory or excitatory.
I compared the results of the test, which uses both exc and inh connections, with the old and with the updated code, and they are exactly the same.
Probably this modification should be done to all conductance-based multisynapse and multicompartment models.

[heplesser]

@golosio @DimitriPlotnikov The issue of the exc/inh buffer also occurs in #261 for the gif_cond_exp_multisynapse models. I put it on the agenda for the NEST developer VC on 19 Sept.

[golosio]

@heplesser @DimitriPlotnikov maybe we should seek a more general solution (see the issue that i just opened #477 ). My opinion is that I would not ask Hesam to review his code once more for this issue before merging it, after all it only affects performances, it can be fixed at a later stage.

[golosio]

@DimitriPlotnikov Did you have a chance to check the code with the changes that I made to meet your observation? Do you have other amendments to propose?

[DimitriPlotnikov]

I'll do it tomorrow.

[DimitriPlotnikov]

From the implementation perspective the model is well implemented. I also wrote an "brunel" simulation based on the provided model and it produced some results (Excitatory rate : 4.40 Hz, Inhibitory rate : 10.20 Hz)
e.g. non zero rates. After fixing mentioned 2 minor issues, it can be finally approved.

[DimitriPlotnikov]

Could you please add an informative comment why are you doing this computation and checks? Can it be already detected in the Parameters_::set?

[DimitriPlotnikov]

I tried to run the example that you provide in the commet. I would add

import nest 
import numpy as np

to make it an executable python program.

Another issue is that the resulting script fails at my machine with the following stacktrace:

Traceback (most recent call last):
  File "aeif_cond_alpha_multisynapse.py", line 24, in <module>
    nest.Connect(spike, neuron, syn_spec="synapse2")
  File "/home/nestml/repositories/nest_repo/bld_master_nompi/lib/python2.7/site-packages/nest/lib/hl_api_helper.py", line 230, in stack_checker_func
    return f(*args, **kwargs)
  File "/home/nestml/repositories/nest_repo/bld_master_nompi/lib/python2.7/site-packages/nest/lib/hl_api_connections.py", line 314, in Connect
    sr('Connect')
  File "/home/nestml/repositories/nest_repo/bld_master_nompi/lib/python2.7/site-packages/nest/__init__.py", line 91, in catching_sli_run
    errorname, commandname, message))
pynestkernel.NESTError: IllegalConnection in Connect_g_g_D_D: Creation of connection is not possible because:
All outgoing connections from a device must use the same synapse type.

Am I doing something wrong? Otherwise, could you fix the comment?

[golosio]

@DimitriPlotnikov Thank you for your review. I made the changes that you requested. I improved the comments and completed the example. I already noticed that the example provided in aeif_cond_alpha_multisynapse don't work, see my comments at the end of #439 (comment) and the following reply from @heplesser . On the other hand, I fixed my example and it should work.
It's interesting that you used the model in a Brunel simulation. I'm also working on a kind of balanced network, I'm trying to implement with a more realistic neuron model the network described in Lim S, Goldman M, Nature Neuroscience 16, 1306–1314 (2013) http://www.nature.com/neuro/journal/v16/n9/full/nn.3492.html

[DimitriPlotnikov]

This PR can be merged.

[golosio]

Thank you @DimitriPlotnikov

[heplesser]

@golosio Thank you for your efforts! The code is ready for merge 👍!