Modelling Pinsky and Rinzel cell: Minimum and heaviside step inside equations object

Support
, ,
1

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,
     'Cad' : 0.,
    '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
dVd/dt=(-gLd*(Vd-VL)-ICad-gKahp*qd*(Vd-VK)-gKC*cd*chid*(Vd-VK)+(gc*(Vs-Vd)+Id)/(1.0-pp))/Cm: volt
dCad/dt=  (-0.13*ICad/amp-0.075*Cad)/ms: 1 

# 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
alphaqd=  (min([0.00002*Cad,0.01]))/ms: Hz
chid=     min([Cad/250.0,1.0]): 1

# calcium current
ICad=     gCa*sd*sd*(Vd-VCa): amp
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

Expected output (if relevant)

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
2

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

3

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")
```
4

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

1 Like