Exponential LIF with alpha shape synaptic conductance

Support
1

I am trying to connect two exponential LIF with alpha shape synapse.
following the documentation for brian2 equivalence of alpha_conductance

tau = 2.5*ms; E = 0*mV
syn_eqs = Equations('''dg_syn/dt = (s - g_syn)/tau : siemens
                       ds/dt = -s/tau : siemens''')
group = NeuronGroup(N, eqs, threshold='vm>-50*mV', reset='vm=-70*mV')
syn = Synapses(source, group, on_pre='s += 10*nS')
# ... connect synapses, etc.

I wrote this:

**Edited **

eqs_e = """
        Im = IX + 
            gL * (EL - vm) + 
            gL * DeltaT * exp((vm - VT) / DeltaT) : amp
        
        dvm/dt = (Im + I_syn_e) / C : volt
        VT : volt
        IX : amp
        I_syn_e : amp
        """
eqs_syn_e = """
        dg_syn_e/dt = (s_e - g_syn_e) / Tau_e : siemens (event-driven)
        ds_e/dt = -s_e / Tau_e : siemens                (event-driven)
        I_syn_e_post =  g_syn_e * (Erev_e - vm): amp (summed)
        """
 E_cells = b2.NeuronGroup(E_cell_params['Ncells'],
                             model=eqs_e,
                             dt=dt,
                             method=integration_method,
                             threshold="vm > 0.*mV",
                             reset="vm=-65*mV"),
                             refractory="vm > 0.*mV",
                             namespace=param_E_syn)

 
cEE = b2.Synapses(E_cells, E_cells, eqs_syn_e,
                      on_pre='g_syn_e += 0.1*nsiemens'),
                      namespace=param_E_syn)
cEE.connect(i=0, j=1)

witch give me this error:

KeyError: 'The identifier "g_syn_e" could not be resolved.'

I also want to integrate the synapse as (event-driven). I am not sure whether it is possible here.
Thanks for any guide.

2

Oh It was just a typo

dg_syn_e/dt = ...

is correct, Am I correct in the implementation of alpha shape synapse?

Selection_001
Selection_0011590Ă—412 66.9 KB

1 is connectoed to 2, I expected to see a bump on 2 as 1 spike as I was familiar from NEST or it may be just need tuning the parameters.

3

Hi. The implementation of the alpha synapse looks mostly correct, but you’ll have to increase s in the on_pre statement, instead of g_syn_e. Otherwise s will always be 0 and this will be an exponential synapse. Also, you cannot use event-driven here. The event-driven flags means that the value of the respective value is only updated whenever a spike arrives. This is not possible here, since you need the values of g_syn_e and s_e at every time step to calculate the synaptic current I_syn_e_post. Note that if you use the same time constant Tau_e for all synapses, then it is much more efficient to include all your equations in the NeuronGroup equations instead of in the Synapses object. This way, you will sum up the conductances/currents of all synapses before calculating the exponential decay only once per neuron. By including these equations in the Synapses object you are calculating the exponential decay once per synapse before summing it up. Mathematically both formulations are equivalent, but the NeuronGroup calculation is much more efficient. Of course, in your example with 2 neurons and 1 synapse this does not really matter…

Regarding your plot, I don’t really see the effect of the synapse but that might be due to the issues mentioned above. To check what is going on in detail you could record and plot g_syn_e and I_syn_e.

1 Like
4

Hi, Thank you
It seems the g_syn_e is updating properly but I_syn_e is not.

Figure_1
Figure_11000Ă—683 32.7 KB

Excuse me the code is a bit long just in case to check.

import os
import time
import pylab as plt
import numpy as np
from numpy.random import randn
import brian2 as b2
from brian2.equations.equations import Equations

b2.seed(2)
np.random.seed(2)


def simulate():

    common_params = {    # Parameters common to all neurons.
        'C': 100*b2.pfarad,
        'tau_m': 10*b2.ms,
        'EL': -60*b2.mV,
        'DeltaT': 2*b2.mV,
        'Vreset': -65,  # *b2.mV
        'VTmean': -50*b2.mV,
        'VTsd': 2*b2.mV
    }
    common_params['gL'] = common_params['C'] / common_params['tau_m']

    param_E_syn = {
        "Erev_e": 0.0*b2.mV,
        "Tau_e": 4.0*b2.ms,
        "w_e": 1.1,  # ! *b2.nsiemens 0.1
        "p_e": 1.0,
    }

    E_cell_params = {'Ncells': num_E_cells,
                     'IXmean': 120*b2.pA,
                     'IXsd': 20*b2.pA}

    eqs_e = """
        VT : volt
        IX : amp
        I_syn_e : amp
        Im = IX + 
            gL * (EL - vm) + 
            gL * DeltaT * exp((vm - VT) / DeltaT) : amp
        
        ds_e/dt = -s_e / Tau_e : siemens                
        dg_syn_e/dt = (s_e - g_syn_e) / Tau_e : siemens 
        dvm/dt = (Im + I_syn_e) / C : volt
        """
    eqs_syn_e = """
        I_syn_e_post =  g_syn_e_pre * (Erev_e - vm): amp (summed)
        """

    E_cells = b2.NeuronGroup(E_cell_params['Ncells'],
                             model=eqs_e,
                             dt=dt0,
                             method=integration_method,
                             threshold="vm > 0.*mV",
                             reset="vm={}*mV".format(common_params['Vreset']),
                             refractory="vm > 0.*mV",
                             namespace={**common_params,
                                        **param_E_syn,
                                        })

    cEE = b2.Synapses(E_cells, E_cells, eqs_syn_e,
                      on_pre='s_e += {}*nsiemens'.format(
                          param_E_syn["w_e"]),
                      dt=dt0,
                      method=integration_method,
                      namespace={**common_params,
                                 **param_E_syn,
                                 })
    cEE.connect(i=0, j=1)

    # Initialise random parameters.
    E_cells.VT = [common_params['VTmean']] * 2
    E_cells.IX = (randn(len(E_cells)) *
                  E_cell_params['IXsd'] + E_cell_params['IXmean'])

    E_cells.vm = randn(len(E_cells)) * 10 * b2.mV - 60 * b2.mV

    spike_monitor_E = b2.SpikeMonitor(E_cells)

    state_monitor_E = b2.StateMonitor(E_cells,
                                      ["vm", "g_syn_e", "I_syn_e"],
                                      record=True)

    net = b2.Network(E_cells)
    net.add(state_monitor_E)
    net.add(cEE)
    net.add(spike_monitor_E)

    b2.run(sim_duration*b2.ms)
    return spike_monitor_E, state_monitor_E


def plot(spike_monitor_E,
         state_monitor_E,
         plot_voltages=False):

    fig, ax = plt.subplots(4, figsize=(10, 5), sharex=True)

    ax[0].plot(spike_monitor_E.t/b2.ms, spike_monitor_E.i, '.k', ms=3)

    if plot_voltages:
        for i in range(num_E_cells):
            ax[1].plot(state_monitor_E.t/b2.ms,
                       state_monitor_E.vm[i]/b2.mV, label=str(i+1))
    ax[2].plot(state_monitor_E.t/b2.ms,
               state_monitor_E.g_syn_e[1]/b2.nsiemens,
               lw=1, color="b", ls="--")
    ax[3].plot(state_monitor_E.t/b2.ms,
               state_monitor_E.I_syn_e[1]/b2.pA,
               lw=1, color="b", ls="--")
    ax[2].set_ylabel(r"$g_{syn}$(nS)")
    ax[3].set_ylabel(r"$I_{syn}$(pA)")

    ax[0].set_ylabel("E cells")
    ax[1].legend(loc="upper right")
    ax[1].set_ylabel("Voltages E")
    ax[-1].set_xlabel("time (ms)")
    plt.savefig("data/E_cell.png", dpi=150)
    plt.show()


if __name__ == "__main__":

    num_E_cells = 2
    dt0 = 0.1*b2.ms

    sim_duration = 200
    stimulus = np.sin(np.arange(0, sim_duration)*2*np.pi/sim_duration)
    state = "beta"

    integration_method = "rk2"
    record_volrages = True
    plot_voltages = record_volrages

    sp_mon, st_mon = simulate()
    plot(sp_mon, st_mon, plot_voltages)
5

This calculates the synaptic current of the post-synaptic cell based on the synaptic conductance of the pre-synaptic cell and therefore gives zero in your example. But you don’t actually need a summed variable here. You can delete the equation from the synapses and define I_syn_e in the neuron equations:

I_syn_e = g_syn_e * (Erev_e - vm) : amp
1 Like
6

That solved the issue,
I still need to think about it.
thank you :slight_smile:

1 Like