How to create synapse between parts of neuron

Support
1

Description of problem

HI, I am trying to create a microcircuit model of a pyramidal cell (pinsky and rinzel model) and an interneuron (wang and buszaki model). I want to create synapses between the dendrite component and the interneuron and similarly between the soma component and the interneuron. I was wondering how one could do that considering the that Synapse function only allows to connect one neuron with another.
I am forwarding my entire code which I understand is long but any input would be appreciated!

Minimal code to reproduce problem

from brian2 import *
import numpy as np
import pandas as pd
import scipy.io

defaultclock.dt = 0.05*ms
pyr_parameters = {
    "Cm" : 3.*pF,
    "gc" : 2.1*nS,
    "pp" : 0.5,
    "Is" : 0.*namp,
    "Id" : 0.*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
}

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

# commented out statements are the orignal statements
pyr_eqs = '''
#membrane potential equations
#for soma
dVs/dt=(-gLs*(Vs-VL)-gNa*Minfs*Minfs*hs*(Vs-VNa)-gKdr*ns*(Vs-VK)+(gc/pp)*(Vd-Vs)+Is/pp)/Cm: volt

#for dendrite
dVd/dt=(-gLd*(Vd-VL)-ICad-gKahp*qd*(Vd-VK)-gKC*cd*chid*(Vd-VK)+(gc*(Vs-Vd))+Id/(1.0-pp))/Cm: volt

#calcium concentration
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)*(-46.9-Vs/mV)/exp((-46.9-Vs/mV)/(4.0)-1.0))/ms: Hz
betams=   ((0.28)*(Vs/mV+19.9)/exp((Vs/mV+19.9)/(5.0)-1.0))/ms: Hz
Minfs=    alphams/(alphams+betams): 1
alphahs=  (0.128*exp((-43.0-Vs/mV)/(18.0)))/ms: Hz
betahs=   (4.0/(1.0+exp((-20.0-Vs/mV)/(5.0))))/ms: Hz
# Kdr potassium current
alphans=  ((0.016)*(-24.9-Vs/mV)/exp((-24.9-Vs/mV)/(5.0)-1.0))/ms: Hz
betans=   (0.25*exp(-1.0-0.025*Vs/mV))/ms: Hz

# dendrite state
# calcium induced potassium current
alphacd=(int(Vd/mV<=-10)*(exp((Vd/mV+50.0)/11)-(Vd/mV+53.5)/27)/18.975+int(1.0*Vd/mV>-10)*2.0*exp((-53.5-Vd/mV)/27.0))/ms  : Hz
betacd=(int(Vd/mV<=-10)*(2.0*exp((-53.5-Vd/mV)/27.0)-alphacd*ms))/ms : Hz

# AHP: calcium induced potassium current
alphaqd=  clip(0.00002*Cad, 0, 0.01)/ms: Hz
chid=     clip(Cad/250.0, 0, 1.0): 1

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

Iext : amp
'''

def run_pinsky_rinzel(pyr_parameters, pyr_initialvalues, pyr_eqs):
    
    start_scope()
    
    Cm  = pyr_parameters['Cm']
    gc  = pyr_parameters['gc']
    pp  = pyr_parameters['pp']
    Is  = pyr_parameters['Is']
    Id  = pyr_parameters['Id']
    gLs = pyr_parameters['gLs']
    gLd  = pyr_parameters['gLd']
    gNa  = pyr_parameters['gNa']
    gKdr  = pyr_parameters['gKdr']
    gCa  = pyr_parameters['gCa']
    gKahp  = pyr_parameters['gKahp']
    gKC  = pyr_parameters['gKC']
    VNa  = pyr_parameters['VNa']
    VCa  = pyr_parameters['VCa']
    VK  = pyr_parameters['VK']
    VL  = pyr_parameters['VL']
    Vsyn  = pyr_parameters['Vsyn']
    alphac  = pyr_parameters['alphac']
    betac  = pyr_parameters['betac']
    betaqd  = pyr_parameters['betaqd']

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

    #provide input current
    time = np.arange(0,100,0.1)*ms
    # current = np.sin(time) * nA
    input_current = TimedArray(sin(2*pi*100*Hz*time) * nA, dt=0.1*ms)
   
    # creates a neuron and specifies its properties
    
    # PC = NeuronGroup(1, model=eqs, threshold = 'Vs >0*mV', refractory=2 * ms, reset=None, dt=pyr_parameters["dt"])
    PC = NeuronGroup(1, model=pyr_eqs, threshold='Minfs > 0.5', refractory=2 * ms, reset=None, dt=pyr_parameters["dt"],method='euler')
    
    PC.set_states(pyr_initialvalues)
    PC.Iext = 'input_current(t)'
    
    # these monitor the state of the variables
    pyr_MVs = StateMonitor(PC, 'Vs', record = True)
    pyr_S = SpikeMonitor(PC)
    
    run(pyr_parameters["simtime"])
    
    return PC, pyr_MVs, pyr_S

int_parameters = {
    "Cmi" : 3.*pF,
    "gNai": 35.*nS,
    "gKi": 9.*nS,
    "gL"  : 0.1*nS,
    "VNa" : 55.*mV,
    "VK"  : -90.*mV,
    "VL" : -65.*mV,
    "dt" : defaultclock.dt,
    "simtime": 100.*ms,
    "Ibaseline": 0.*namp,
    "Iscale": 0.*namp,
}

int_initialvalues = {
    'v' : -70.*mV,
    'h' : 1.,
    'n' : 0.5,
}

int_eqs = '''
#Buzsaki and Wang model for interneurons
dv/dt = (-gNai*m**3*h*(v-VNa)-gKi*n**4*(v-VK)-gL*(v-VL))/Cmi : volt
m = alpha_m/(alpha_m+beta_m) : 1
alpha_m = -0.1/mV*(v+35*mV)/(exp(-0.1/mV*(v+35*mV))-1)/ms : Hz
beta_m = 4*exp(-(v+60*mV)/(18*mV))/ms : Hz
dh/dt = 5*(alpha_h*(1-h)-beta_h*h) : 1
alpha_h = 0.07*exp(-(v+58*mV)/(20*mV))/ms : Hz
beta_h = 1./(exp(-0.1/mV*(v+28*mV))+1)/ms : Hz
dn/dt = 5*(alpha_n*(1-n)-beta_n*n) : 1
alpha_n = -0.01/mV*(v+34*mV)/(exp(-0.1/mV*(v+34*mV))-1)/ms : Hz
beta_n = 0.125*exp(-(v+44*mV)/(80*mV))/ms : Hz

Iext : amp
'''

def run_buzsaki(int_parameters, int_initialvalues, int_eqs):
    start_scope()
    Cmi = int_parameters["Cmi"]
    gNai = int_parameters["gNai"]
    gKi = int_parameters["gKi"]
    gL = int_parameters["gL"]
    VNa = int_parameters["VNa"]
    VK = int_parameters["VK"]
    VL = int_parameters["VL"]

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

    #provide input current
    time = np.arange(0,100,0.1)*ms
    # current = np.sin(time) * nA
    input_current = TimedArray(sin(2*pi*100*Hz*time) * nA, dt=0.1*ms)
   
    # creates a neuron and specifies its properties
    
    # PC = NeuronGroup(1, model=eqs, threshold = 'Vs >0*mV', refractory=2 * ms, reset=None, dt=int_parameters["dt"])
    interneuron = NeuronGroup(1, model=int_eqs, threshold='m > 0.5', refractory=2 * ms, reset=None, dt=int_parameters["dt"],method='euler')
    
    interneuron.set_states(int_initialvalues)
    interneuron.Iext = 'input_current(t)'
    
    # these monitor the state of the variables
    int_MVs = StateMonitor(interneuron, 'v', record = True)
    int_S = SpikeMonitor(interneuron)

    run(int_parameters["simtime"])
    
    return interneuron, int_MVs, int_S

synapse_parameters = {
    "vsyn" : -62.*mV,
    "taurise" : 1*ms,
    "taudecay" : 5*ms,
    "thetasyn" : 0*mV,
    "sigmasyn" : 1*mV,
    "taud" : 0*ms,
    "simtime": 100.*ms,
    "Ibaseline": 0.*namp,
    "Iscale": 0.*namp,
    "dt" : defaultclock.dt,
}

synapse_eqs = '''
#synaptic equations
dIsyn/dt = gsyn*(vpost-vsyn) : amp
dSsyn/dt = alphasyn(1-Ssyn)*F(vpre(t-taud), thetasyn, sigmasyn) - betasyn*Ssyn : 1

F = 1/(1+exp(-(vpre-theta)/sigma)) : 1
alphasyn = 1/taurise : Hz
betasyn = 1/taudecay : Hz
taurise : second
taudecay : second
theta : 1
sigma : 1
thetasyn : 1
sigmasyn : 1
'''
pyramidal_cell = run_pinsky_rinzel(pyr_parameters, pyr_initialvalues, pyr_eqs)
interneuron = run_buzsaki(int_parameters, int_initialvalues, int_eqs)

#create synapse between dendrite and interneuron
synapse = Synapses(pyramidal_cell, interneuron, model=synapse_eqs, on_pre='Isyn_post+=Isyn', method='euler')
synapse.connect(i=0, j=0)
synapse.taurise = synapse_parameters["taurise"]
synapse.taudecay = synapse_parameters["taudecay"]
synapse.theta = synapse_parameters["thetasyn"]
synapse.sigma = synapse_parameters["sigmasyn"]
synapse.vsyn = synapse_parameters["vsyn"]
synapse.taud = synapse_parameters["taud"]

What you have aready tried

Expected output (if relevant)

Actual output (if relevant)

Full traceback of error (if relevant)

ERROR      Brian 2 encountered an unexpected error. If you think this is a bug in Brian 2, please report this issue either to the discourse forum at <http://brian.discourse.group/>, or to the issue tracker at <https://github.com/brian-team/brian2/issues>. Please include this file with debug information in your report: C:\Users\parih\AppData\Local\Temp\brian_debug_440v9bgp.log  Additionally, you can also include a copy of the script that was run, available at: C:\Users\parih\AppData\Local\Temp\brian_script_ds_blepc.py Thanks! [brian2]
Traceback (most recent call last):
  File "c:\Users\parih\OneDrive\Documents\GitHub\Microcircuit_Model\main.py", line 251, in <module>
    synapse = Synapses(pyramidal_cell, interneuron, model=synapse_eqs, on_pre='Isyn_post+=Isyn', method='euler')
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\synapses\synapses.py", line 833, in __init__
    self.add_dependency(source)
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\core\base.py", line 203, in add_dependency
    self._dependencies.add(obj.id)
AttributeError: 'tuple' object has no attribute 'id'
2

Hi @bindok. I am currently down with a cold, so a longer reply has to wait, but in general there is no problem in connecting such models, since from Brian’s point of view these are “normal” neurons. In the synaptic equations, you can refer to e.g. Vs_pre to refer to the pre-somatic membrane potential and Vd_post to refer to a post-synaptic dendrite. However, your synaptic equations cannot be written like you do (e.g. vpre cannot be a function of time), and you’d directly write F instead of F(…), etc.
The error message you are getting is unrelated to all that:

pyramidal_cell = run_pinsky_rinzel(pyr_parameters, pyr_initialvalues, pyr_eqs)

This variable stores a tuple of a NeuronGroup, a StateMonitor and a SpikeMonitor, so you cannot hand it over to the Synapses initializer like this.

3

Ah okay! No worries.
However, even when I only return the NeuronGroup variable in the pinsky_rinzel and interneuron functions, I still get a similar error message

raceback (most recent call last):
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\parsing\rendering.py", line 79, in render_node
    return getattr(self, methname)(node)
AttributeError: 'SympyNodeRenderer' object has no attribute 'render_Tuple'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "c:\Users\parih\OneDrive\Documents\GitHub\Microcircuit_Model\main.py", line 251, in <module>
    synapse = Synapses(pyramidal_cell, interneuron, model=synapse_eqs, on_pre='Vs', on_post='Vd', method='euler')
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\synapses\synapses.py", line 845, in __init__
    model = Equations(model)
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\equations\equations.py", line 619, in __init__
    self._equations = parse_string_equations(eqns)
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\utils\caching.py", line 107, in cached_func
    func._cache[cache_key] = func(*args, **kwds)
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\equations\equations.py", line 416, in parse_string_equations
    expression = Expression(p.sub(" ", expression))
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\equations\codestrings.py", line 108, in __init__
    str_to_sympy(code)
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\parsing\sympytools.py", line 77, in str_to_sympy
    return _str_to_sympy(expr)
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\utils\caching.py", line 107, in cached_func
    func._cache[cache_key] = func(*args, **kwds)
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\parsing\sympytools.py", line 83, in _str_to_sympy
    s_expr = SympyNodeRenderer().render_expr(expr)
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\parsing\rendering.py", line 67, in render_expr
    return self.render_node(node.body)
  File "C:\Users\parih\anaconda3\lib\site-packages\brian2\parsing\rendering.py", line 81, in render_node
    raise SyntaxError(f"Unknown syntax: {nodename}")
SyntaxError: Unknown syntax: Tuple

Also, if function F cannot be written as F(…), then would you write something like F*(…) or simply F…? Neither of the options intuitively makes sense to me as they do not denote a function of V_pre as such. Could you please explain how brian interprets functions?
Here are the synaptic equations if incase I have written them incorrectly in the first place.
image

4

Hi @bindok. Unfortunately, I do not have the time for an exhaustive response right now, maybe someone else can jump in. Brian has a system to support arbitrary functions that you can define in code, but this is a mechanism that is used/necessary only for functions that you cannot write down as simple expressions based on existing functions. E.g., while in theory you could enhance Brian’s function system by a function inverse that is defined as inverse(x) = 1/x, in practice you’d use what we call a “subexpression” (see e.g. Models and neuron groups — Brian 2 2.5.4 documentation). An equation (just for illustration, not making sense): \frac{dv}{dt} = -v \cdot \text{inverse}(\tau) : volt would be written in Brian as:

"""
dv/dt = -v * inverse : volt
inverse = 1/tau : 1/second
"""

Here, inverse isn’t really a function, but just a convenience replacement for better readability. The code is exactly equivalent to

"""
dv/dt = -v * (1/tau) : volt
"""

For your code this means to write simply F in the equation for Ssyn (instead of F(...)), and define F as

F = 1/(1+exp(-(vpre-thetasyn)/sigmasyn)) : 1

We’ve discussed having a more generic system that would be useful if you re-use the same function with different parameters several times (where with Brian’s approach you’d need to repeat the subexpression with different names), but for now we don’t have this – but then, it doesn’t seem to be necessary for your use case.

Note that the above isn’t exactly correct, since your equations state v_{\text{pre}}(t - \tau_d), i.e refer to the membrane potential from the past. This isn’t actually possible at the moment in Brian (again, there had been discussions to add this feature, but it is non-trivial to implement). Currently, these exact equations can only been implemented using network_operation or user-defined functions (see Delay in continuous connections for a lot of pointers). But in my experience, using the sigmoid function of the pre-synaptic membrane potential with a continuous delayed connection is often unnecessarily complex, replacing it with a computationally simpler spike-based connection or (if necessary) with a more abstract implementation like a “square pulse” (see "Square" signal in response to spike) will most likely not change the result in a meaningful way.

Hope that helps a bit at least!

5

Yes it does!
Thank you!