Add NeuroML implementation of PyNN's NoisyCurrentSource

[pedroernesto]

It is missing the code for 'class NoisyCurrentSource' in file 'standardmodels.py' of the neuroML module for PyNN. This is: "A Gaussian 'white' noise current source. The current amplitude changes at fixed intervals (dt), with the new value drawn from a Gaussian distribution"

Other 'source standard models' of PyNN like 'ACSource', 'StepCurrentSource' and 'DCSource' already existing in that file, have been created from class 'PulseGeneratorDL' in file 'nml.py' of LibNeuroML. Maybe the same can be done for a train of pulses with random amplitudes as 'NoiseGeneratorDL'?

(we understand that 'nml.py' is automatically created from the NeuroML Schema (NeuroML_v2beta4.xsd) and that the relevant LEMS definitions are in Inputs.xml)

[pgleeson]

Hi @pedroernesto. I assume you mean here: https://github.com/pgleeson/PyNN/blob/master/pyNN/neuroml/standardmodels.py#L272? Yes, there has been no corresponding <noisyCurrentSource> up until now. This has mainly been due to the fact that there are many possible ways to implement such a noisy current, as discussed during this recent exchange on the NeuralEnsemble mailing list.

The PyNN impl of NoisyCurrentSource is certainly a good place to start, maybe @apdavison can give some insight why that particular model was chosen.

As mentioned in the email exchange above, it should be reasonably straightforward to implement a new <noisyCurrentSource> in LEMS. I've made a start for this here. There is a PyNN script there using NoisyCurrentSource and ideally the new ComponentType should be updated to behave the same. Would you be happy to help with this? After that it would be straightforward to add this to the NML2 schema and update libNeuroML etc. so it can be exported from PyNN scripts to the nml2 equivalent.

[pgleeson]

Related to this: NeuralEnsemble/PyNN#437

[pgleeson]

@pedroernesto @apdavison, any update on whether you've had a chance to look at this?

[pedroernesto]

I have modified a little bit the LEMS file cited above by @pgleeson: OpenSourceBrain/StochasticityShowcase@aa59894#diff-4639db1e9157899cf5af89aace49c88fR14, to obtain a code for the corresponding ComponentType for such PyNN current source model. Please, have a look at: https://github.com/pedroernesto/StochasticityShowcase/blob/master/NeuroML2/LEMS_NoisyCurrentInput.xml

I have followed the hint of combining the ideas used for coding the ComponentTypes 'pulseGenerator' and (specially) 'spikeGeneratorRandom', given here: https://groups.google.com/forum/#!topic/neuralensemble/1-vCUWM7su4

[pedroernesto]

Figures for step = 0.05 ms (integ. time-step) and dt = 1 ms (interval between updates of the noisy current):

Figures for step = 0.05 ms and dt = 20 ms:

[pgleeson]

@pedroernesto, good start. Note though the mean current should be 0.55 (blue line) and a standard deviation of 1, so the current should fluctuate around the blue line mostly staying within ~ -0.5 to ~ 1.5 and going outside these bounds occasionally. You'll need a normal distribution on the value of the current, as opposed to just using the random() method in LEMS. Unfortunately there's no randn() method for this in LEMS as is used in the PyNN, but look here to see how you could generate this value from random(): https://en.wikipedia.org/wiki/Normal_distribution#Generating_values_from_normal_distribution. The Box Muller method looks appropriate.

[pedroernesto]

Thanks for the hints @pgleeson. Here are new figures obtained with the Box-Muller method implemented in LEMS. Integration time-step = 0.05 ms. Code will be uploaded soon.

dt = 1 ms:

[pedroernesto]

dt = 5 ms:

[pedroernesto]

dt=20ms:

[pedroernesto]

Code is now available at https://github.com/pedroernesto/StochasticityShowcase/blob/master/NeuroML2/LEMS_NoisyCurrentInput.xml. Hints are welcome to improve it in some directions:

(a) 'randn' might be implemented as an instance of an independent ComponentType, and then be called as a Child element inside ComponentType "noisyCurrentSource".

(b) Box-Muller's method actually can provide in each run with two normally distributed numbers instead of just one: sqrt(-2*ln(U))cos(2 * 3.14 * V) and sqrt(-2ln(U))*sin(2 * 3.14 * V). A kind of static variable that can be used as a flag for the next call is needed.

(c) Pi value used here is 3.14, but which is the built-in LEMS constant for that?

(d) Box-Muller method can be improved by using its Marsaglia's variant, which avoids the use of trigonometric functions: https://en.wikipedia.org/wiki/Normal_distribution#Generating_values_from_normal_distribution
However, the generation of uniformly distributed random numbers U and V should continue till both lie on a unit circle. Is there a way to obtain in LEMS a kind of do-while loop?

[pgleeson]

Thanks for this @pedroernesto, I've checked through it and it looks good. I've put it in a branch here, and added a notebook plotting the profile of long current injections: https://github.com/OpenSourceBrain/StochasticityShowcase/blob/test_noisy/NeuroML2/TestNoisySources.ipynb. @rgerkin's Ornstein-Uhlenbeck example is here too. I'm just waiting to resolve some issues with the PyNN example running on Travis-CI to merge this PR to master (some issue with neo+PyNN @apdavison?).

Regarding your questions:
a) The simplest way is how you've done it, just use it directly in the expression. This can be directly mapped to Neuron, etc. then (though these examples fail in mapping to mod files as Neuron doesn't like a variable named dt...)

b) True, but it would be more complicated to store/retrieve the value between events, making the componenttype less clear. It is overhead as it is, but preferable to have a simpler definition of the behaviour

c) Pi could be one of the global variables, but there's no consistent way to define those yet (the gas constant, e are others). For now it's fine to hard code it (with many more decimal places...) in the LEMS and it can be moved at a later stage.

d) Again this would be a more complex definition, any Lems compliant simulator should have access to trigonometric functions, so this is the clearest way to get the values.

[apdavison]

The PyNN issue on Travis is probably related to the release of PyNN 0.9.0, which requires Neo >= 0.5.0. The latter has backwards incompatible changes wrt Neo < 0.4