Dynamic current leads to runaway values

Support
1

Hello,
I’m trying to inject a current which is dependent on the membrane voltage. I’ve already created a TimedArray of conductance values and multiply this by Vm to get the current. However, during the run both current and voltage quickly become too large or NaN (depending on integration method). I suspect this is due to the current and voltage amplifying each other. Is there a way to update the current based on the voltage one time-step in to the past, or an even better way?

Minimal code to reproduce the problem

from brian2 import *
import matplotlib.pyplot as plt
import numpy as np

start_scope()

#Conductances
g_exc = abs(np.random.normal(size=10000))*mS
g_exc = TimedArray(g_exc, dt=0.2*ms)

#Parameters
area = 20000*umetre**2
E_Na = 55*mV
E_K = -77*mV
E_leak = -65*mV
g_Na = (120*msiemens) 
g_K = (35*msiemens) 
g_leak = (0.3*msiemens)
C = (1*ufarad)
duration = 2000*ms 

#Make the neuron
eqs = '''
alpha_m = (0.182/mV) * (35.*mV + v) / (1 - exp(-(35.*mV + v) / (9.*mV)))/ms : Hz
alpha_h = (0.25) * exp(- (90.*mV + v) / (12.*mV))/ms : Hz
alpha_n = (0.02/mV) * (25.*mV - v) / (1 - exp(-(25.*mV - v) / (9.*mV)))/ms : Hz 

beta_m = (-0.124/mV) * (35.*mV + v) / (1 - exp((35.*mV - v) / (9.*mV)))/ms : Hz
beta_h = (0.25) * exp((62.*mV + v) / (6.*mV)) / exp((90.*mV + v) / (12.*mV))/ms : Hz
beta_n = (-0.002/mV) * (25.*mV - v) / (1 - exp((25.*mV - v) / (9.*mV)))/ms : Hz

dm/dt = alpha_m * (1-m) - beta_m * m : 1
dh/dt = alpha_h * (1-h) - beta_h * h : 1
dn/dt = alpha_n * (1-n) - beta_n * n : 1

dv/dt = (g_Na * m**3 * h * (v - E_Na) +
        g_K * n**4 * (v - E_K) +
        g_leak * (v - E_leak) + I) / C : volt
        
I = (g_exc(t) * v) : amp
'''

neuron = NeuronGroup(1, model=eqs, method='euler')
neuron.v = -65*mV
neuron.m = 0.05
neuron.h = 0.6
neuron.n = 0.32
M = StateMonitor(neuron, ['v', 'I'], record = True)

run(duration, report='text')

plt.plot(M.t/ms, M.v[0])
#plt.plot(M.t/ms, M.I[0])
plt.show()

I’ve also tried changing I (in the equations) to: dI/dt = (g_exc(t) * v)/ms : amp
Is there a way to update the current based on the voltage one time-step into the past, or an even better way?

Kind regards,
Mik

2

Hi. The problem is not due to the current, if you remove the current you see the same behaviour. The current as such is correct if you want to model an excitatory conductance. That it depends on the voltage isn’t a problem, that’s the case of all currents (in your model, the Na, K, and leak currents depend on the voltage as well). The confusing thing is maybe that you multiply the conductance directly by v, but this should be read as multiplying it by (v - E_exc) where E_exc = 0*mV. Having said all that, the real problem is that your currents have the wrong sign. You’ll have to add a minus in front of g_Na, g_K, g_leak, and I, or use e.g. (E_Na - v) instead of (v - E_Na).

Also note that the euler method is not adapted to these kind of equations, at least unless you use very short time steps. We generally recommend using exponential_euler and probably also a smaller time step, e.g. with defaultclock.dt = 0.01*ms. Also, have a look at this discussion why it would be better to reformulate some of your alpha_.../beta_... equations with the exprel function.

Finally, I think your model uses the parameters of Hodgkin & Huxley’s :squid: axon, make sure this is appropriate for what you are trying to model :grinning: .

2 Likes