Implementation of stochastic noise in a Hodgkin-Huxley model but as a current

Hi all! As the title suggests, I am curious if anyone has coded a model using stochastic noise as a current instead of a voltage.

The provided example G = NeuronGroup(10, 'dv/dt = -v/tau + sigma*sqrt(2/tau)*xi : volt') is a very nice representation for when the noise is added as a voltage. However, what if in a Hodgkin-Huxley model we have to model noise as a Gaussian with (let’s say) a mean mu=0 and a variance of sigma=10*pA? What would be the proper syntax?

The reason I am stumped is the tau in the sqrt; tau should be equal to Cm/gL correct? It’s not trivial for me to see how a current source would be modeled… Thanks in advance!

P.S. I have checked the relevant links here about modeling noise using randn(). It would be great to add to the discussion the current noise too!

Hi @nvar
My first approach would be to add PoissonInput (see PoissonInput) object but control a current variable in HH model or even better to add ‘noise synapses’ and control the dynamic variable for these synapses. Although there are limitations, for example you can’t use double exponential synapses, as noise synapses, this will be the most neuron-like noise.

Hi,
You can try this:

‘dI_noise/dt = (mu_noise-I_noise) / tau_noise + sigma_noise * (sqrt(2/tau_noise) * xi) :amp’

Indicative values:

tau_noise=20ms, sigma_noise=3pA, mu_noise=0*pA

xi is just a unitless stochastic variable that in every dt draws a number from a normal distribution.

Also see the Methods from Bono & Clopath, 2017.

Good luck

1 Like

Hi all and thanks for your replies.

  1. @rth The PoissonInput would just model noise as random spikes, however, it follows a Poisson distribution and the noise would not be white. Also, it depends on the model of the synapses instead of adding a current source inside the neuron equations.

  2. @mpagkalos Yes, that would work, and thanks for the reference to the Clopath paper. However, I was looking for an implementation in a Hodgkin-Huxley model, where I do not have tau defined like so, but it’s rather more complicated. Furthermore, xi is not unitless, from the Brian2 docs it’s unit is 1/sqrt(sec).

What I was looking for is something along the lines of the following:

'''
dV/dt = (- I_L - I_noise) / Cm : volt
    I_L = g_L * (V - E_L) : amp
    noise = ... : amp
'''

What would the noise term have to be in the above to work properly as a random variable such that X~N(0,10pA)? Also, in my implementation I avoid using xi but rather randn(), as per this post.

In short, the noise as posted by Nina is given by:

'''...
noise = sigma*(2*gl/Cm)**.5*randn()/sqrt(dt) : volt/second (constant over dt)
'''

which is just voltage.

Thanks again!

1 Like

Hi Nvar, I am relatively new to all this. I was curious along the same line but concerning adding stochastic noise to a leaky integrate and fire model. Below is a picture of the ODEs implemented in my eqn for my NeuronGroup object. I would like your input on this; I do have a tau variable. Looking forward to your response. Thanks

Hi and welcome. While your example follows along the lines of the provided example in the documentation here, you are adding the noise outside the differential equation for V. In your model you explicitly define g_l and C_m, and tau in most cases is defined as tau=RmCm. But if you have a tau already specified for your noise, the definition you give is perfectly fine from my point of view. Otherwise you could maybe use tau_noise = Cm/g_l.

Thanks for the response. However, I seem to get an error when I run the code with this definition. I have tried writing it in a different way but to no avail. It describes an inconsistent unit error when I try to run it this way. I have tried writing it as “+ sigma_noisexitau**-0.5”.

This is the kind of error I get.

Hi everyone :wave: I wasn’t following this thread in detail, but just regarding the error message posted by @amanforthepeople : you seem to be using a dimensionless value for sigma, where it should be a voltage value (e.g. sigma = 1*mV).

2 Likes

Thanks for the stopping by! If it’s not too much of trouble, my original question is kind of unanswered as of now. Would you kindly chime in? In short, I was curious of how I should impement the stochastic noise as a current instead of a voltage in a Hodgkin-Huxley formalism. The provided equations in the documentation always seem to add stochastic noise as an extra voltage term and it is quite difficult for me to figure it out.

In your opinion, how would you add noise as a current with a mean of 1nA and a sigma of 10nA? Thank you in advance!

Hi @nvar . This topic that is much more complex than one might think when one reads papers that “add a Gaussian noise current” to the model :blush: It has been on my list for a while to write a longer document on this topic, but well…

Here are a few remarks:

There are different reasons to model “noise” in a model, and therefore different mathematical descriptions. Two commonly considered sources of noise is synaptic bombardment, where one does not want to model the individual spikes: see @rth 's comment above regarding PoissonInput. If one assumes a very simple synapse model, low individual weights and the limit of “infinitely many incoming spikes” then this leads to an Ornstein-Uhlenbeck process with the time constant \tau_m added to the membrane potential, which is what we use in many examples. There are similar approximations when considering e.g. ion channel noise.

Now, you want to add “noise as a current with a mean of 1nA and a sigma of 10nA”. Taken literally, this is easy: you can include I_noise in the dv/dt equation, as you did above, and then define

I_noise = mu_noise + sigma_noise*randn() : amp (constant over dt)

Something to this effect is actually done in a number of papers. The problem with this formulation is that this describes an implementation in a step-based simulation with a certain time step, and not an “ideal” mathematical model. You could think of this as similar to describing the membrane potential not by a differential equation, but by something like v_new = v_old + v_update. Or to change the perspective a bit: in electrophysiology, it wouldn’t be enough to say “we inject a \cal{N}(\mu, \sigma) white noise current”, but you’d either say something like “a \cal{N}(\mu, \sigma) white noise current filtered in this way”, or “we inject current, with its strength drawn every 1ms from \cal{N}(\mu, \sigma)”.

Here’s a simple simulation showing that this matters: it uses the equation from above to choose a new current every time step:

Code example 1
from brian2 import *

mu_noise = .1*nA
sigma_noise = 1*nA
E_L = -70*mV
g_L = 10*nS
C_m = 100*pF

def run_sim(dt):
    defaultclock.dt = dt
    group = NeuronGroup(1,
                         '''dv/dt = (g_L*(E_L - v) + I_noise)/C_m : volt
                            I_noise = mu_noise + sigma_noise*randn() : amp (constant over dt)''',
                        method='exact')
    group.v = E_L
    mon = StateMonitor(group, ['v', 'I_noise'], record=0)

    run(200*ms)
    return mon.t[:], mon.v[0], mon.I_noise[0]

t1, v1, I_noise1 = run_sim(0.5*ms)
t2, v2, I_noise2 = run_sim(0.01*ms)

fig, axs = plt.subplots(4, 2, sharey='row',
                       gridspec_kw={'width_ratios': [3, 1]},
                       constrained_layout=True)
axs[0, 0].plot(t1/ms, I_noise1/nA)
axs[0, 0].set_ylabel('I_noise (nA)')
axs[0, 1].hist(I_noise1/nA, bins=np.arange(-5, 5, 0.5), orientation='horizontal')
axs[1, 0].plot(t1/ms, v1/mV)
axs[1, 0].set_ylim(-80, -20)
axs[1, 0].set_ylabel('v (mV)')
axs[1, 1].axis('off')

axs[2, 0].plot(t2/ms, I_noise2/nA)
axs[2, 0].set_ylabel('I_noise (nA)')
axs[2, 1].hist(I_noise2/nA, bins=np.arange(-5, 5, 0.5), orientation='horizontal')
axs[3, 0].plot(t2/ms, v2/mV)
axs[3, 0].set_ylim(-80, -20)
axs[3, 0].set_ylabel('v (mV)')
axs[3, 1].axis('off')
axs[3, 0].set_xlabel('t (ms)')

plt.show()
For demonstration purposes, I run this with `dt=0.5*ms` (top) and `dt=0.01*ms` (bottom):

noise_1

As you can see in the histograms on the right, in both cases we have currents distributed according to a Gaussian with the same mean and variance. But as you can clearly see, the membrane potential reacts quite differently in the two cases. An alternative would be to use a run_regularly operation to only update the I_noise variable at regular intervals (e.g. once very ms), independent of the time step.

To abstract from the implementation and use a mathematical description, you can use differential equations and the variable \xi to describe the input current as an OU-process that reverts back to the mean:

    group = NeuronGroup(1,
                         '''dv/dt = (g_L*(E_L - v) + I_noise)/C_m : volt
                            dI_noise/dt = (mu_noise - I_noise)/tau_noise + sigma_noise*sqrt(2/tau_noise)*xi : amp''',
                        method='euler')

Note that tau_noise is a parameter of the noise (how fast the noise changes) and not necessarily Cm/gL – using the membrane time constant for the noise is something that comes out of certain approximations, but it depends on what kind of noise you want to model specifically. This gives you noise that has a comparable effect independent of the time step, e.g. with the same time steps as above and tau_noise=1*ms it gives

noise_2

The implementation in the Bono & Clopath paper is very similar – it might even be equivalent, I did not look into this in detail.

As you noted in the other thread, using xi does not work with certain algorithms like exponential_euler, so in this case unfortunately you have to again leave the beautiful realm of mathematics and come back into the ugly reality of simulation. Replacing the above by:

'''....
dI_noise/dt = (mu_noise - I_noise)/tau_noise + noise : amp
noise = sigma_noise*sqrt(2/tau_noise)*randn()/sqrt(dt) : amp/second (constant over dt)
'''

will give you the Euler approximation to the noise term (so in the above code where it uses euler already, this will not change anything), which is then compatible with exponential_euler and all other algorithms that do not support stochastic differential equations.

I hope that made things a bit clearer?

2 Likes

Dear Marcel,

Thank you for taking the time to explain things this thoroughly. It is very clear and I see the differences of the two, plus it makes more sense now that the sampling time affects the evolution of the membrane potential.

On a final note, in your last piece of code you used a variable tau_noise. I assume that this would be more accurately defined as tau_noise = Cm/g_noise? I am trying to make sense of the math behind it, but I’ll get there eventually.

Thanks again for taking the time and providing such a brilliant response!

It’s quite difficult to say in general, it depends on what the current noise is supposed to represent. The tau_noise variable defines the time constant of the noise, which is related to the cut-off frequency of a filter. Here’s a quote from a somewhat random paper on the topic:

It is the same spectrum of a first order RC filter with a cut-off frequency f_c = \frac{1}{2πτ} and this justifies our interpretation of such a stochastic process as a model for a filtered white noise. For low frequencies, in fact, the spectrum is flat, however high frequencies are cutted-off with slope −2.

I don’t think it should depend on the neuron’s Cm, since its effect on the membrane potential will already depend on the value of Cm, including it in the calculation of the time constant seems redundant.

2 Likes