# NaN warnings appeared in the subexpression of SpatialNeuron

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()
``````

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
``````

Shall I wirte another individual eqs object for Cm? I would really appreciate it if someone could help .

Hi @czx. Apologies for the late reply, I completely overlooked this question (I was away for all of August).

Unfortunately, Brian does not actually support defining `Cm` by a subexpression! When you provide an expression for `Cm`, this expression will only be used to initialize `Cm`, i.e. as if you’d have written

``````neuron = SpatialNeuron(...)
neuron.Cm = "1*uF/cm**2*sqrt(1+ema)"
``````

Evaluating `ema` before `taufu` has been set to a non-zero value will lead to NaN, so `Cm` will be set to NaN… You could set it after defining the other parameters to avoid the NaN, but it would not change over time.

I don’t see any easy solution to this for now, Brian’s framework is built for a constant `Cm` (potentially varying over compartments) – some terms necessary for updating `v` at each time step are pre-calculated in the very beginning of the simulation, and some of these terms refer to `Cm`

1 Like

Hi @mstimberg
That’s all right. Thanks for your kind reply. I will check whether it is necessary to provide an expression for Cm. Since the simulation time is relatively short, and it may not lead to an obvious change.