# Description of problem

I am trying to implement the ball-and-stick Pinsky and Rinzel cell. However, some of the state equations use heaviside step functions (example: `betacd= (1.0-heaviside(Vd+10.0*mV))*((2.0/mV)*exp((-53.5*mV-Vd)/(27.0*mV))-alphacd)/ms: Hz`), and other state equations use a minimum (example: `alphaqd= min([0.00002*Cad,0.01]): 1`). Even though these are well-defined numpy functions, I don’t think the equations object recognizes them. What is the correct way to implement either of those inside an equations object?

# Minimal code to reproduce problem

``````from brian2 import *
import numpy as np
import pandas as pd
import scipy.io
from IPython.core.debugger import set_trace

defaultclock.dt = 0.05*ms
Parameters = {
"Cm" : 3.*pF,
"gc" : 2.1*nS,
"pp" : 0.5*1,
"Is" : 0.*namp,
"Id" : -0.3*namp,
"gLs": 0.1*nS,
"gLd": 0.1*nS,
"gNa": 30.*nS,
"gKdr":15.*nS,
"gCa": 10.*nS,
"gKahp": 0.8*nS,
"gKC": 15.*nS,
"VNa": 60.*mV,
"VCa": 80.*mV,
"VK" : -75.*mV,
"VL" : -60.*mV,
"Vsyn" : 0.*mV,
"alphac": 2./ms,
"betac" : 0.1/ms,
"betaqd": 0.001/ms,
"dt" : defaultclock.dt,
"simtime": 100.*ms,
"Ibaseline": 0.*namp,
"Iscale": 0.*namp
}

InitialValues = {
'Vs' : -70.*mV,
'Vd' : -70.*mV,
'hs' : 1.,
'ns' : 0.5,
'sd' : 0.,
'cd' : 0.5,
'qd' : 0.
}

eqs = '''
dVs/dt=(-gLs*(Vs-VL)-gNa*Minfs*Minfs*hs*(Vs-VNa)-gKdr*ns*(Vs-VK)+(gc/pp)*(Vd-Vs)+Is/pp)/Cm: volt

# soma dynamics
dhs/dt= alphahs-(alphahs+betahs)*hs: 1
dns/dt= alphans-(alphans+betans)*ns: 1

# dendrite dynamics
dcd/dt= alphacd-(alphacd+betacd)*cd: 1
dqd/dt= alphaqd-(alphaqd+betaqd)*qd: 1

# calcium dynamics
dsd/dt= alphasd-(alphasd+betasd)*sd: 1

# soma state
# sodium current
alphams=  ((0.32/mV)*(-46.9*mV-Vs)/exp((-46.9*mV-Vs)/(4.0*mV)-1.0))/ms: Hz
betams=   ((0.28/mV)*(Vs+19.9*mV)/exp((Vs+19.9*mV)/(5.0*mV)-1.0))/ms: Hz
Minfs=    alphams/(alphams+betams): 1
alphahs=  (0.128*exp((-43.0*mV-Vs)/(18.0*mV)))/ms: Hz
betahs=   (4.0/(1.0+exp((-20.0*mV-Vs)/(5.0*mV))))/ms: Hz
# Kdr potassium current
alphans=  ((0.016/mV)*(-24.9*mV-Vs)/exp((-24.9*mV-Vs)/(5.0*mV)-1.0))/ms: Hz
betans=   (0.25*exp(-1.0-0.025*Vs/mV))/ms: Hz

# dendrite state
# KC potassium current
alphacd=((1.0-heaviside(Vd+10.0*mV)/mV)*exp((Vd+50.0*mV)/(11*mV)-(Vd+53.5*mV)/(27*mV))/18.975+heaviside(Vd+10.0*mV)*(2.0/mV)*exp((-53.5*mV-Vd)/(27.0*mV)))/ms: Hz
betacd=   (1.0-heaviside(Vd+10.0*mV))*((2.0/mV)*exp((-53.5*mV-Vd)/(27.0*mV))-alphacd)/ms: Hz
# AHP: calcium induced potassium current

# calcium current
alphasd=  (1.6/(1.0+exp(-0.072*(Vd-5.0*mV)/mV)))/ms: Hz
betasd=   ((0.02*(Vd+8.9*mV)/mV)/(exp((Vd+8.9*mV)/(5.0*mV))-1.0))/ms: Hz
'''

def run_simulation(Parameters, InitialValues, eqs):

start_scope()

Cm  = Parameters['Cm']
gc  = Parameters['gc']
pp  = Parameters['pp']
Is  = Parameters['Is']
Id  = Parameters['Id']
gLs = Parameters['gLs']
gLd  = Parameters['gLd']
gNa  = Parameters['gNa']
gKdr  = Parameters['gKdr']
gCa  = Parameters['gCa']
gKahp  = Parameters['gKahp']
gKC  = Parameters['gKC']
VNa  = Parameters['VNa']
VCa  = Parameters['VCa']
VK  = Parameters['VK']
VL  = Parameters['VL']
Vsyn  = Parameters['Vsyn']
alphac  = Parameters['alphac']
betac  = Parameters['betac']
betaqd  = Parameters['betaqd']

# time settings
defaultclock.dt = Parameters["dt"]

# creates a neuron and specifies its properties

PR = NeuronGroup(1, model=eqs, threshold='Minfs > 0.5', refractory=2 * ms, reset=None, dt=Parameters["dt"])
PR.set_states(InitialValues)

# these monitor the state of the variables
MVs = StateMonitor(PR, 'Vs', record = True)
S = SpikeMonitor(PR)

run(Parameters["simtime"])

return MVs, S

[MVs, S] = run_simulation(Parameters, InitialValues, eqs)
``````

# What you have aready tried

The problem is really with the dendritic and calcium equations (I have tried to run a soma only and that works well). I also tried to run it with the dendrite only having leak, ICad and compartmental currents. However, having no hyperpolarizing currents makes everything explode quickly. I tried with only the AHP current (so left the KC current out) but received the same error, making me believe it has something to do with the mininum and heaviside definitions

# Actual output (if relevant)

``````Traceback (most recent call last):

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/parsing/rendering.py:78 in render_node
return getattr(self, methname)(node)

AttributeError: 'SympyNodeRenderer' object has no attribute 'render_List'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3398 in run_code
exec(code_obj, self.user_global_ns, self.user_ns)

Input In [41] in <cell line: 1>
[MVs, S] = run_simulation(Parameters, InitialValues, eqs)

Input In [40] in run_simulation
PR = NeuronGroup(1, model=eqs, threshold='Minfs > 0.5', refractory=2 * ms, reset=None, dt=Parameters["dt"])

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/groups/neurongroup.py:499 in __init__
model = Equations(model)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/equations/equations.py:554 in __init__
self._equations = parse_string_equations(eqns)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/utils/caching.py:101 in cached_func
func._cache[cache_key] = func(*args, **kwds)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/equations/equations.py:366 in parse_string_equations
expression = Expression(p.sub(' ', expression))

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/equations/codestrings.py:107 in __init__
str_to_sympy(code)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/parsing/sympytools.py:76 in str_to_sympy
return _str_to_sympy(expr)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/utils/caching.py:101 in cached_func
func._cache[cache_key] = func(*args, **kwds)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/parsing/sympytools.py:82 in _str_to_sympy
s_expr = SympyNodeRenderer().render_expr(expr)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/parsing/rendering.py:66 in render_expr
return self.render_node(node.body)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/parsing/rendering.py:78 in render_node
return getattr(self, methname)(node)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/parsing/rendering.py:238 in render_Call
return self.render_func(node.func)(*(self.render_node(arg)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/parsing/rendering.py:238 in <genexpr>
return self.render_func(node.func)(*(self.render_node(arg)

File ~/opt/anaconda3/envs/Brian/lib/python3.10/site-packages/brian2/parsing/rendering.py:80 in render_node
raise SyntaxError(f"Unknown syntax: {nodename}")

File <string>
SyntaxError: Unknown syntax: List
``````

Sorry, I see that not all the code has been picked up as code properly

Hi @fleurzeldenrust. The right-hand-side of Brian equations uses a subset of Python syntax, in particular you cannot use indexing, lists, and other data structures. The error message you are getting is actually about the use of the list syntax `[0.00002*Cad,0.01]`. As you assumed, you cannot use all numpy functions – most importantly because the code might be translated into C++ (or CUDA C++), and we don’t implement all those functions. The list of supported functions is here: Functions — Brian 2 2.5.4 documentation
Questions about the `heaviside` function come up regularly, we should probably support it under this name for convenience. Until then, you can use `int` (that can convert a boolean operation into 0 or 1) together with a boolean operation. In your example, instead of `heaviside(Vd + 10.0*mV)`, you can use `int(Vd > -10*mV)`.
Similarly, you can replace the use of `min` by `clip`. Since the idea of `min([..., 0.01])` is to clip the value at 0.01 you can replace it by `clip(..., -inf, 0.01)`.

Hope that was clear and that the suggested solutions work for you!

PS:

The posts use markdown formatting, so to quote code within a line you can use single backticks like ``this``, and for longer code you can use triple backticks:

`````````
print("some code")
```
``````

Thank you so much for the fast response, that should do the trick.

1 Like