Issue in multicompartmental models: 'Membrane current must be linear with respect to v'

[jacop...@hotmail.com]

Dear Brian team

I wanted to point out an issue I encountered several times (on several computers/laptops and by multiple people):

When implementing non-linear ion channels or synapses into the multicompartmental neuron, often the error 'Membrane current must be linear with respect to v' appears. However, opening a new python console usually enables the code to run normally. It seems to occur in some consoles while it is absent in others... I have included a piece of code below that shows the behaviour, and included two screenshots after running this code from two python consoles in Spyder, one throwing an error while the other runs normally.

Is there anything I am doing wrong?

Kind Regards,
Jacopo



Here is a minimal code with NMDA-type synapses that reproduces the problem:

#///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
# EXAMPLE CODE
#///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

from __future__ import division
from brian2 import *

#####################################################
# Set Parameters
#####################################################

start_scope()
defaultclock.dt = 25.*us # timeStep

#------------
# PASSIVE parameters
#------------
R_axial = 90.*ohm*cm                # axial resistance
g_leak =  0.04*msiemens/cm**2       # leak conductance
Capacit = 1.*uF/cm**2               # membrane capacitance
EL=-70.*mV                          # resting membrane voltage

#------------
# NMDA synapses
#------------
ENMDA = 0.*mV           # NMDA reversal potential
tau_NMDA = 50.*ms       # NMDA time constant
syn_weight = 0.5        # initial synaptic weight
wn_max = 1.             # maximal NMDA weight
w_min = .1            # maximal NMDA weight
nmda_cond = 2000.*1e-12*siemens  # NMDA maximal conductance


#####################################################
# Equations for the biophysical neuron
#####################################################

eqs_leak='''     
Im = g_leak*(EL-v) : amp/meter**2
'''

eqs_synapses = """
I_NMDA = gNMDA*(ENMDA-v)*Mgblock: amp (point current)
gNMDA : siemens
Mgblock = 1./(1. +  exp(-0.062*v/mV)/3.57) :1
"""

eqs = eqs_leak + eqs_synapses


######################################################
## SpatialNeuron
######################################################
morpho = Soma(diameter=30*um)
morpho.dendrite = Cylinder(length=50*um, diameter=2*um, n=5)
neuron = SpatialNeuron(morphology=morpho, model=eqs, Cm=Capacit, Ri=R_axial)


#####################################################
# Run
#####################################################
print('Start simulation ...')
run(5*ms, report=None)
print('Simulation finished!')

[Marcel Stimberg]

Hi,

this is very odd. There's no other difference between the consoles, like
different environments with different dependency versions or anything
like that? The error message means that sympy could not figure out how
to re-write your Im current in the form A*v + B, but 1) this should be
trivial in your example and 2) it should either work every time or
never. I remember a colleague had some problem with repeatedly running
sympy-related stuff in spyder, could you select "Execute in a new
dedicated Python console" in Run->Configure and see whether this fixes
the issue?

If I am right about it being a sympy issue, then it should be possible
to trigger the problem with a simpler example:

import sympy as sp
from brian2.parsing.sympytools import str_to_sympy

var = sp.Symbol('v', real=True)
wildcard = sp.Wild('_A', exclude=[var])
constant_wildcard = sp.Wild('_B', exclude=[var])
pattern = wildcard * var + constant_wildcard
s_expr = sp.collect(str_to_sympy("g_leak*(EL-v)").expand(), var)
matches = s_expr.match(pattern)
print(matches)

This should print

{_B_: EL*g_leak, _A_: -g_leak}

If it instead prints "None", then this corresponds to the error you were
seeing.

Best,

Marcel

[jacop...@hotmail.com]

Dear Marcel

Thank you for the reply. The different consoles in Brian2 are exactly the same, which makes the error very strange indeed. When encountering a console that throws the error after running the example code from my first message, Mgblock = 1./(1. +  exp(-0.062*v/mV)/3.57) :1 , is the line that causes the error to occur in some consoles. This is indeed non-linear in the voltage, and if I change it to e.g. Mgblock = 1. :1, the error disappears in that console. Is it possible that this non-linearity from synapses is reflected in the membrane current? In any case, the very odd thing is that in other consoles it runs smoothly...

- I have tried "Execute in a new dedicated Python console", but the same strange behaviour occurs here as well: the first 3 runs were without problem, the fourth one resulted in the error message. This was all with exactly the same code and within the same Spyder.

- I went back to try the simpler example as you mentioned. In a console where the error pops up when running my script, I executed your simpler example and that seems to run without issues (see screenshot below).

- Just to check whether this is a Spyder issue, I have tried to run the example code via the linux terminal and yet again, the first terminal raised the error, while a second ran the script without issue...

I hope this is helpful in finding a cause!

Best,
Jacopo

[Marcel Stimberg]

Hi Jacopo,

this is still very strange. Everything on Brian's side should be
deterministic, so I don't see why running it again in a different
console should make a difference. However, I realized that *not* raising
the error is indeed a bug, as you say your equation is non-linear in v.
You add "I_NMDA" as a point current which simply means that I_NMDA/area
will be added to the definition of "Im". Something is going wrong with
the substitution of I_NMDA and Mgblock here (it should check whether the
equation for Im is linear in v *after* substituting the expressions for
I_NMDA and Mgblock) -- possibly this bug depends on the order of the
equations in a dictionary or something similar which is not guaranteed
to be the same over runs which would explain the strange behaviour.

Now unfortunately we currently do not support non-linear Im currents (we
have plans to relax this requirement, though). The only workaround to
get the model simulated is to approximate the current by assuming that
it (or only the Mgblock variable, i.e. the non-linear part) is constant
over a time step. This introduces an error for time steps where the
voltage is changing rapidly, but with small enough time steps it should
not make a significant difference, I think. We have a new syntax for
this in the development version of Brian 2 (and it will be included in
the 2.0 version that we will release in about 2 weeks from now), but I
just realized that it does only work for normal NeuronGroups, not for
SpatialNeuronGroup. One more thing to fix...

You can do the same thing manually, though: instead of defining Mgblock
(or I_NMDA) as a subexpression, define it as a parameter:
Mgblock : 1

and update it using a run_regularly operation:
neuron.run_regularly('Mgblock = 1./(1. + exp(-0.062*v/mV)/3.57)')

Best,
Marcel

[jacop...@hotmail.com]

Dear Marcel

Thank you very much, that makes a bit more sense now. I tried the 'run_regularly' solution, and it runs without notable difference (I use small time-steps anyway since I also have several Hodgkin-Huxley type channels in the model)!

On that point, I was wondering if there is any difference in using point current as opposed to a current in amp/meter**2 regarding non-linear currents? For example, when I use Hodgkin-Huxley type conductances across the neuron (as in http://brian2.readthedocs.io/en/2.0rc/examples/compartmental.hh_with_spikes.html), I don't seem to encounter the error. The non-linear dependence on voltage is quite similar to the NMDA synapses, but the only difference seems to be that the former is not a 'point current'.

Best,
Jacopo

[Marcel Stimberg]

Hi Jacopo,

> On that point, I was wondering if there is any difference in using
> point current as opposed to a current in amp/meter**2 regarding
> non-linear currents? For example, when I use Hodgkin-Huxley type
> conductances across the neuron (as in
> http://brian2.readthedocs.io/en/2.0rc/examples/compartmental.hh_with_spikes.html),
> I don't seem to encounter the error. The non-linear dependence on
> voltage is quite similar to the NMDA synapses, but the only difference
> seems to be that the former is not a 'point current'.

the numerical integration of multi-compartmental models is a two-step
process. First, all differential equations are integrated for each
compartment independently. This can use any integration method that
Brian offers and therefore does not have any restriction on the
equations (the "exponential euler" method that is used by default for
SpatialNeuron does require that each differential equation is linear
with respect to the variable it is defining, though). This updates the
Hodgkin-Huxley gating variables n, m, h, etc. The second step updates
the membrane potential based on Im and the membrane potential of
neighbouring compartments -- this step currently requires that Im
depends linearly on "v". The equations for other variables, such as the
gating variables, do not matter, since they have already been
integrated. Your I_NMDA current on the other hand introduces a direct
non-linear dependence of Im on the membrane potential and that is not
supported. The fact that it is a point current is not the relevant
difference.

Hope that makes things clearer, best
Marcel

[jacop...@hotmail.com]

Hi Marcel

Oh yes of course, this makes a lot of sense!

Thanks for all the help!

Kind regards,
Jacopo