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

Science, Projects, and Showcases
,
10

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