Description of problem
Hi everyone. I’m working with Brian2.7 and studying the example hh_with_spikes , and trying to change the Cm of the neuron from a constant to a subexpression. The subexpression contains a square root operation. However, when running the code, the variables raised NaN warnings just like below. I think it may be the problem with the SQRT operation, since I also tried other simple subexpressions and they successfully worked. I wonder what the reason is and whether this issue can be avoided.
Minimal code to reproduce problem
from brian2 import *
from scipy import stats
defaultclock.dt = 0.01*ms
morpho = Cylinder(length=10*cm, diameter=2*238*um, n=1000, type='axon')
El = 10.613*mV
ENa = 115*mV
EK = -12*mV
gl = 0.3*msiemens/cm**2
gNa0 = 120*msiemens/cm**2
gK = 36*msiemens/cm**2
# Typical equations
eqs = '''
Im = gl * (El-v) + gNa * m**3 * h * (ENa-v) + gK * n**4 * (EK-v) : amp/meter**2
I : amp (point current) # applied current
dm/dt = alpham * (1-m) - betam * m : 1
dn/dt = alphan * (1-n) - betan * n : 1
dh/dt = alphah * (1-h) - betah * h : 1
alpham = (0.1/mV) * 10*mV/exprel((-v+25*mV)/(10*mV))/ms : Hz
betam = 4 * exp(-v/(18*mV))/ms : Hz
alphah = 0.07 * exp(-v/(20*mV))/ms : Hz
betah = 1/(exp((-v+30*mV) / (10*mV)) + 1)/ms : Hz
alphan = (0.01/mV) * 10*mV/exprel((-v+10*mV)/(10*mV))/ms : Hz
betan = 0.125*exp(-v/(80*mV))/ms : Hz
gNa : siemens/meter**2
# strain evolution, the membraine capacitance Cm is related to it
ema = (ema0-eDa0)*exp(-t/taufu) + eDa0 : 1
ema0 : 1
eDa0 : 1
taufu : second
'''
neuron = SpatialNeuron(morphology=morpho, model=eqs, method="exponential_euler",
refractory="m > 0.4", threshold="m > 0.5",
Cm="1*uF/cm**2*sqrt(1+ema)", Ri=35.4*ohm*cm) #Cm is related to ema
neuron.v = 0*mV
neuron.h = 1
neuron.m = 0
neuron.n = .5
neuron.I = 0*amp
neuron.gNa = gNa0
#initialize strain parameters
neuron.taufu = 111.5*second
neuron.ema0 = 0.05
neuron.eDa0 = 0
M = StateMonitor(neuron, ['v','ema'], record=True)
run(50*ms, report='text')
neuron.I[0] = 1*uA # current injection at one end
run(3*ms)
neuron.I = 0*amp
run(50*ms, report='text')
subplot(211)
for i in range(10):
plot(M.t/ms, M.v.T[:, i*100]/mV)
ylabel('v')
subplot(212)
plot(M.t/ms, M.ema[0])
ylabel('ema')
xlabel('Time (ms)')
show()
What you have aready tried
Tried other simple subexpressions such as Cm = 2* 1* uF/cm** 2, and it worked.
Expected output (if relevant)
The variables in the eqs should be normally calculated , since the ema in the square root is positive.
Actual output (if relevant)
NaN warnings like below. Electrical variables such as v is not available.
Full traceback of error (if relevant)
Starting simulation at t=0. s for a duration of 50. ms
WARNING spatialneuron's variable 'm' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING spatialneuron's variable 'h' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING spatialneuron's variable 'n' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
50. ms (100%) simulated in < 1s
Starting simulation at t=53. ms for a duration of 50. ms
50. ms (100%) simulated in < 1s
