# Description of problem

When I run SpatialNeuron(), Python keeps running for a long time and uses an increasing amount of memory until it finally runs out of system RAM (and crashes).

I have an equation with many ionic channels, I deleted every single one of them and started adding them one by one and found out this is due to the I_Na current but it doesn’t make much sense to me as to why.

# Minimal code to reproduce problem

``````start_scope()

eqs_IO_Ca = '''
dCa/dt = (-3*I_Ca_h*((uamp / cm**2)**-1)*mM - 0.075*Ca)/ms : mM
'''
eqs_IO_Isom = '''
I_l     = g_l*(v-V_l)               : metre**-2*amp
I_Na    = g_Na*m_inf**3*h*(v-V_Na)  : metre**-2*amp
I_Ca_l  = g_Ca_l*k*k*k*l*(v-V_Ca)   : metre**-2*amp
I_K_dr  = g_Kdr*n*n*n*n*(v-V_K)     : metre**-2*amp
I_h     = g_h*q*(v-V_h)             : metre**-2*amp
I_K_s   = g_K_s*(x_s**4)*(v-V_K)    : metre**-2*amp
'''
eqs_IO_Iden = '''
I_Ca_h  = g_Ca_h*r*r*(v-V_Ca)       : metre**-2*amp
I_K_Ca  = g_K_Ca*s*(v-V_K)          : metre**-2*amp
'''
eqs_IO_Iax = '''
I_K_a  = g_K_a *x_a**4*(v-V_K)      : metre**-2*amp
I_Na_a = g_Na_a*m_a**3*h_a*(v-V_Na) : metre**-2*amp
'''
eqs_IO_activation = '''
dh/dt = (h_inf - h)/tau_h : 1
dk/dt = (k_inf - k)/tau_k : 1
dl/dt = (l_inf - l)/tau_l : 1
dn/dt = (n_inf - n)/tau_n : 1
dq/dt = (q_inf - q)/tau_q : 1
dr/dt = (r_inf - r)/tau_r : 1
ds/dt = (s_inf - s)/tau_s : 1
m_a = m_inf_a : 1
dh_a/dt = (h_inf_a - h_a)/tau_h_a : 1
dx_a/dt = (x_inf_a - x_a)/tau_x_a : 1
dx_s/dt = (x_inf_s - x_s)/tau_x_s : 1
'''
eqs_IO_inf = '''
m_inf   = alpha_m /(alpha_m+beta_m)        : 1
h_inf   = alpha_h/(alpha_h+beta_h)         : 1
k_inf   = 1/(1+e**(-(v/mvolt+61)/4.2))    : 1
l_inf   = 1/(1+e**((v/mvolt+85.5)/8.5))   : 1
n_inf   = alpha_n/(alpha_n+beta_n)         : 1
q_inf   = 1/(1+e**((v/mvolt+75)/(5.5)))   : 1
r_inf   = alpha_r/(alpha_r + beta_r)       : 1
s_inf   = alpha_s/(alpha_s+beta_s)         : 1
m_inf_a = 1/(1+(e**((-30-v/mvolt)/ 5.5))) : 1
h_inf_a = 1/(1+(e**((-60-v/mvolt)/-5.8))) : 1
x_inf_a = alpha_x_a/(alpha_x_a+beta_x_a)   : 1
x_inf_s = alpha_x_s/(alpha_x_s + beta_x_s) : 1
'''
eqs_IO_tau = '''
tau_h   = 170*msecond/(alpha_h+beta_h)                                          : second
tau_k   = 5*msecond                                                             : second
tau_l   = 1*msecond*(35+(20*e**((v/mvolt+160)/30/(1+e**((v/mvolt+84)/7.3))))) : second
tau_n   = 5*msecond/(alpha_n+beta_n)                                            : second
tau_q   = 1*msecond/(e**((-0.086*v/mvolt-14.6))+e**((0.07*v/mvolt-1.87)))     : second
tau_r   = 5*msecond/(alpha_r + beta_r)                                          : second
tau_s   = 1*msecond/(alpha_s + beta_s)                                          : second
tau_h_a = 1.5*msecond*e**((-40-v/mvolt)/33)                                    : second
tau_x_a = 1*msecond/(alpha_x_a + beta_x_a)                                      : second
tau_x_s = 1*msecond/(alpha_x_s + beta_x_s)                                      : second
'''
eqs_IO_alpha = '''
alpha_m   = (0.1*(v/mvolt + 41))/(1-e**(-(v/mvolt+41)/10)) : 1
alpha_h   = 5.0*e**(-(v/mvolt+60)/15) : 1
alpha_n   = (v/mvolt + 41)/(1-e**(-(v/mvolt+41)/10)) : 1
alpha_r   = 1.7/(1+e**(-(v/mvolt - 5)/13.9)) : 1
alpha_s   = ((0.00002*Ca/mM)*int((0.00002*Ca/mM)<0.01) + 0.01*int((0.00002*Ca/mM)>=0.01)) : 1
alpha_x_a = 0.13*(v/mvolt + 25)/(1-e**(-(v/mvolt+25)/10)) : 1
alpha_x_s = 0.13*(v/mvolt + 25)/(1-e**(-(v/mvolt+25)/10)) : 1
'''

eqs_IO_beta = '''
beta_m = 9.0*e**(-(v/mvolt+60)/20)                        : 1
beta_h = (v/mvolt+50)/(1-e**(-(v/mvolt+50)/10))          : 1
beta_n = 12.5*e**(-(v/mvolt+51)/80)                       : 1
beta_r = 0.02*(v/mvolt + 8.5)/(e**((v/mvolt + 8.5)/5)-1) : 1
beta_s = 0.015                                             : 1
beta_x_a  = 1.69*e**(-0.0125*(v/mvolt + 35))              : 1
beta_x_s  = 1.69*e**(-0.0125*(v/mvolt+ 35))               : 1
'''

eqs_vector = '''
V_Na : volt
V_K  : volt
V_Ca : volt
V_l  : volt
V_h  : volt
g_Na   : siemens/meter**2
g_Kdr  : siemens/meter**2
g_Ca_l : siemens/meter**2
g_h    : siemens/meter**2
g_Ca_h : siemens/meter**2
g_K_Ca : siemens/meter**2
g_l : siemens/meter**2
g_Na_a   : siemens/meter**2
g_K_a   : siemens/meter**2
g_K_s   : siemens/meter**2
'''

eqs_IO = eqs_IO_beta
eqs_IO += eqs_IO_alpha
eqs_IO += eqs_IO_tau
eqs_IO += eqs_IO_inf
eqs_IO += eqs_IO_activation
eqs_IO += eqs_IO_Iax
eqs_IO += eqs_IO_Iden
eqs_IO += eqs_IO_Isom
eqs_IO += eqs_IO_Ca
eqs_IO += eqs_vector

eqs_IO_V = '''
Im = (-(I_l +I_Na +  I_Ca_l + I_K_dr + I_h + I_K_s + I_Ca_h + I_K_Ca + I_c + I_K_a + I_Na_a)) : metre**-2*amp
I_c : metre**-2*amp
Iapp : amp (point current)
'''
eqs_IO += eqs_IO_V

dt = 0.025*ms # We set the timestep for the simulation
dt_rec = 0.25*ms # We set the timestep for the recording

# We create the spatial neuron. Please note that this takes a long time (more than 5 minutes).
IO = SpatialNeuron(morphology=IO_morpho, model = eqs_IO, Cm=1 *uF*cm**-2, Ri=100.0*ohm*cm, method="exponential_euler", name = 'SchweighoferOlive', dt=dt)
``````

# What you have aready tried

Deleting all the currents and adding them one by one.

# Full traceback of error (if relevant)

``````---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
/opt/conda/lib/python3.7/site-packages/sympy/core/assumptions.py in getit(self)
261         try:
--> 262             return self._assumptions[fact]
263         except KeyError:

KeyError: 'zero'

During handling of the above exception, another exception occurred:

KeyboardInterrupt                         Traceback (most recent call last)
<timed exec> in <module>

/opt/conda/lib/python3.7/site-packages/brian2/spatialneuron/spatialneuron.py in __init__(self, morphology, model, threshold, refractory, reset, events, threshold_location, dt, clock, order, Cm, Ri, name, dtype, namespace, method, method_options)
319                 f"Cannot take the derivative of '{Im_expr.code}' with respect to v.")
320
--> 321         gtot_str = sympy_to_str(sp.simplify(-diffed))
322         I0_str = sympy_to_str(sp.simplify(Im_sympy_exp - diffed*v_sympy))
323

/opt/conda/lib/python3.7/site-packages/sympy/simplify/simplify.py in simplify(expr, ratio, measure, rational, inverse, doit, **kwargs)
598     expr = bottom_up(expr, lambda w: getattr(w, 'normal', lambda: w)())
599     expr = Mul(*powsimp(expr).as_content_primitive())
--> 600     _e = cancel(expr)
601     expr1 = shorter(_e, _mexpand(_e).cancel())  # issue 6829
602     expr2 = shorter(together(expr, deep=True), together(expr1, deep=True))

/opt/conda/lib/python3.7/site-packages/sympy/polys/polytools.py in cancel(f, *gens, **args)
6618
6619     try:
-> 6620         (F, G), opt = parallel_poly_from_expr((p, q), *gens, **args)
6621     except PolificationFailed:
6622         if not isinstance(f, (tuple, Tuple)):
``````

… follows deeper into sympy …

Hi @eliasmateo95 . Thanks for the report, I can confirm the issue on my machine. The problem seems to occur when we rearrange the equations and call sympy’s `simplify` along the way. I will try to figure out what exactly is going on next week. Until then you could try a little workaround: to simplify the equations yourself, you can mark some of the equations (not differential equations, but e.g. `m_inf`, `alpha_m`, etc.) as “(constant over dt)”. This means, that instead of plugging them into the equations and using the integration method, they are calculated once at the beginning of a time step. This makes the equations much simpler, since they can now be treated as constants (even though of course their value changes every time step). With a small simulation time step, this shouldn’t make much difference compared to the usual solution, especially for slowly changing currents. Even just marking a few of marking very few of them this way seems to change the equations enough so that sympy no longer runs into this weird endless loop.
Best,
Marcel

1 Like