genn-team / genn

GeNN is a GPU-enhanced Neuronal Network simulation environment based on code generation for Nvidia CUDA.
http://genn-team.github.io/
GNU Lesser General Public License v2.1
231 stars 57 forks source link

GeNN only supports writing to a single postsynaptic variable ? #120

Closed Vafa-Andalibi closed 7 years ago

Vafa-Andalibi commented 7 years ago

Hello,

Working with brian2genn, we noticed a limitation on writing to multiple postsynaptic variable. Do you have any timeline for adding this feature?

Best,

tnowotny commented 7 years ago

Hi, what kind of scenario is it that you are after?

Vafa-Andalibi commented 7 years ago

Hi @tnowotny, It's about implementing homeostatic synaptic scaling (HS). The HS is a common property of neurons, useful for simulations in general, and thus probably useful for plastic models. We are implementing the HS into the dynamic equations. The HS requires that the firing rate of the cell scales the strength of the synapses; the excitatory and inhibitory synapses are scaled with mutually inverted (1/f(x) )functions. The scaling is dependent on detecting the postsynaptic firing rate (dynamic spike sensor), which determines the dynamic evolution of synaptic scaling factor. All this happens at the cell level. At the synapse level, we first determine whether the synapse is excitatory or inhibitory, and if inhibitory, we invert the function. Second, we clip the function to known limits, and third, we multiply the weight with the scaling factor to get the conductance change which is passed on. All these synaptic calculations take place on_pre , i.e. only when the presynaptic spike arrives.

tnowotny commented 7 years ago

Hi Vafa,as I understand it this means your modifications (the write variable) are the synaptic conductances. This is not a problem either in GeNN or brian2genn. Conductances are synapse variables, not variables of post-synpatic neurons. Read access to post-synaptic variables is unrestricted (unless something went wrong in our code that I am not aware of). So , I wonder where your implementation is running into problems? Does it have to do with that you want to calculate some (post-synaptic) neuron properties on pre-spikes only? In principle, these calculations can take place as part of the synaptic update ... but I am not sure I am yet getting the full picture. Do you have sample code how you would do it in plain Brian 2?

Vafa-Andalibi commented 7 years ago

Brian2genn checks whether the synaptic model is changing more than one postsynaptic variable and if so, it will raise an error stating GeNN only supports writing to a single postsynaptic variable

This is the synaptic equation that we have :

def STDP_with_scaling(self):
        '''
        The method for implementing the STDP synaptic connection.
        '''
        self.output_synapse['equation'] = Equations('''
            wght : siemens
            wght0 : siemens
            dapre/dt = -apre/taupre : siemens (event-driven)
            dapost/dt = -apost/taupost : siemens (event-driven)
            ''')
        self.output_synapse['pre_eq'] = ''
        if self.output_synapse['receptor'] in ['gi']:
            self.output_synapse['pre_eq'] += '''
            synaptic_scaling_factor = 1./synaptic_scaling_factor
            '''
        self.output_synapse['pre_eq'] += '''synaptic_scaling_factor = clip(synaptic_scaling_factor, 0.66, 1.5)'''
        if self.output_synapse['namespace']['Apre'] >= 0:
            self.output_synapse['pre_eq'] += '''
                        %s += synaptic_scaling_factor * wght
                        apre += Apre * wght0 * Cp
                        wght = clip(wght + apost, 0, wght_max)
                        ''' % (self.output_synapse['receptor'] + self.output_synapse['post_comp_name'] + '_post')
        else:
            self.output_synapse['pre_eq'] += '''
                        %s += synaptic_scaling_factor * wght
                        apre += Apre * wght * Cd
                        wght = clip(wght + apost, 0, wght_max)
                        ''' % (self.output_synapse['receptor'] + self.output_synapse['post_comp_name'] + '_post')
        if self.output_synapse['namespace']['Apost'] <= 0:
            self.output_synapse['post_eq'] = '''
                        apost += Apost * wght * Cd
                        wght = clip(wght + apre, 0, wght_max)
                        '''
        else:
            self.output_synapse['post_eq'] = '''
                        apost += Apost * wght0 * Cp
                        wght = clip(wght + apre, 0, wght_max)

This is an addition to our all neuron group models:

            self.output_neuron['equation'] = Equations('''
            dsynaptic_scaling_factor/dt = scaling_speed  * (1 - (spike_sensor / (ap_target_frequency*tau_synaptic_scaling))) : 1
            dspike_sensor/dt = -spike_sensor/tau_synaptic_scaling : 1
            ''')
            self.output_neuron['reset'] += '; spike_sensor +=1'

and this is the neuron group model generator:

    def PC(self):
        '''
        This method build up the equation for PC neurons based on the number of compartments. The final equation is then saved in output_neuron['equation'].
        Main internal variables:
        * eq_template_soma: Contains template somatic equation, the variables in side the equation could be replaced later using "Equation" function in brian2. :
            ::
                dvm/dt = (gL*(EL-vm) + gealpha * (Ee-vm) + gealphaX * (Ee-vm) + gialpha * (Ei-vm) + gL * DeltaT * exp((vm-VT) / DeltaT) +I_dendr) / C : volt (unless refractory)
                dge/dt = -ge/tau_e : siemens
                dgealpha/dt = (ge-gealpha)/tau_e : siemens
                dgeX/dt = -geX/tau_eX : siemens
                dgealphaX/dt = (geX-gealphaX)/tau_eX : siemens
                dgi/dt = -gi/tau_i : siemens
                dgialpha/dt = (gi-gialpha)/tau_i : siemens
                x : meter
                y : meter
        * eq_template_dend: Contains template dendritic equation:
            ::
                dvm/dt = (gL*(EL-vm) + gealpha * (Ee-vm) + gealphaX * (Ee-vm) + gialpha * (Ei-vm) + gL * DeltaT * exp((vm-VT) / DeltaT) +I_dendr) / C : volt (unless refractory)
                dge/dt = -ge/tau_e : siemens
                dgealpha/dt = (ge-gealpha)/tau_e : siemens
                dgeX/dt = -geX/tau_eX : siemens
                dgealphaX/dt = (geX-gealphaX)/tau_eX : siemens
                dgi/dt = -gi/tau_i : siemens
                dgialpha/dt = (gi-gialpha)/tau_i : siemens
        '''
        #: The template for the somatic equations used in multi compartmental neurons, the inside values could be replaced later using "Equation" function in brian2.
        # eq_template_soma = self.value_extractor(self.cropped_df_for_current_type,'eq_template_soma')
        # eq_template_dend = self.value_extractor(self.cropped_df_for_current_type,'eq_template_dend')
        # print '\nAP target is ', str(self.output_neuron['namespace']['ap_target_frequency'])
        # print '\ntau syn sc is ', str(self.output_neuron['namespace']['tau_synaptic_scaling'])
        # print '\nSpike sensor target value is ', str(self.output_neuron['namespace']['ap_target_frequency'] *
        #                                              self.output_neuron['namespace']['tau_synaptic_scaling'])
        eq_template_soma = '''
        dvm/dt = ((gL*(EL-vm) + gealpha * (Ee-vm) + gealphaX * (Ee-vm) + gialpha * (Ei-vm) + gL * DeltaT * exp((vm-VT) / DeltaT) +I_dendr) / C) +  noise_sigma*xi*taum_soma**-0.5 : volt (unless refractory)
        dge/dt = -ge/tau_e : siemens
        dgealpha/dt = (ge-gealpha)/tau_e : siemens
        dgeX/dt = -geX/tau_eX : siemens
        dgealphaX/dt = (geX-gealphaX)/tau_eX : siemens
        dgi/dt = -gi/tau_i : siemens
        dgialpha/dt = (gi-gialpha)/tau_i : siemens
        '''
        #: The template for the dendritic equations used in multi compartmental neurons, the inside values could be replaced later using "Equation" function in brian2.
        eq_template_dend = '''
        dvm/dt = (gL*(EL-vm) + gealpha * (Ee-vm) + gealphaX * (Ee-vm) + gialpha * (Ei-vm) +I_dendr) / C : volt
        dge/dt = -ge/tau_e : siemens
        dgealpha/dt = (ge-gealpha)/tau_e : siemens
        dgeX/dt = -geX/tau_eX : siemens
        dgealphaX/dt = (geX-gealphaX)/tau_eX : siemens
        dgi/dt = -gi/tau_i : siemens
        dgialpha/dt = (gi-gialpha)/tau_i : siemens
        '''
        self.output_neuron['equation'] += Equations(eq_template_dend, vm="vm_basal", ge="ge_basal",
                                                   gealpha="gealpha_basal",
                                                   C=self.output_neuron['namespace']['C'][0],
                                                   gL=self.output_neuron['namespace']['gL'][0],
                                                   gi="gi_basal", geX="geX_basal", gialpha="gialpha_basal",
                                                   gealphaX="gealphaX_basal", I_dendr="Idendr_basal")
        self.output_neuron['equation'] += Equations(eq_template_soma, gL=self.output_neuron['namespace']['gL'][1],
                                                    ge='ge_soma', geX='geX_soma', gi='gi_soma', gealpha='gealpha_soma',
                                                    gealphaX='gealphaX_soma',
                                                    gialpha='gialpha_soma', C=self.output_neuron['namespace']['C'][1],
                                                    I_dendr='Idendr_soma',taum_soma=self.output_neuron['namespace']['taum_soma'])
        for _ii in range(self.output_neuron['dend_comp_num'] + 1):  # extra dendritic compartment in the same level of soma
            self.output_neuron['equation'] += Equations(eq_template_dend, vm="vm_a%d" % _ii,
                                                        C=self.output_neuron['namespace']['C'][_ii],
                                                        gL=self.output_neuron['namespace']['gL'][_ii],
                                                        ge="ge_a%d" % _ii,
                                                        gi="gi_a%d" % _ii, geX="geX_a%d" % _ii,
                                                        gealpha="gealpha_a%d" % _ii, gialpha="gialpha_a%d" % _ii,
                                                        gealphaX="gealphaX_a%d" % _ii, I_dendr="Idendr_a%d" % _ii)
        # basal self connection
        self.output_neuron['equation'] += Equations('I_dendr = gapre*(vmpre-vmself)  : amp',
                                                    gapre=1 / (self.output_neuron['namespace']['Ra'][0]),
                                                    I_dendr="Idendr_basal", vmself="vm_basal", vmpre="vm")
        self.output_neuron['equation'] += Equations('I_dendr = gapre*(vmpre-vmself)  + gapost*(vmpost-vmself) : amp',
                                                    gapre=1 / (self.output_neuron['namespace']['Ra'][1]),
                                                    gapost=1 / (self.output_neuron['namespace']['Ra'][0]),
                                                    I_dendr="Idendr_soma", vmself="vm",
                                                    vmpre="vm_a0", vmpost="vm_basal")
        self.output_neuron['equation'] += Equations('I_dendr = gapre*(vmpre-vmself) + gapost*(vmpost-vmself) : amp',
                                                    gapre=1 / (self.output_neuron['namespace']['Ra'][2]),
                                                    gapost=1 / (self.output_neuron['namespace']['Ra'][1]),
                                                    I_dendr="Idendr_a0", vmself="vm_a0", vmpre="vm_a1", vmpost="vm")
        for _ii in arange(1, self.output_neuron['dend_comp_num']):
            self.output_neuron['equation'] += Equations('I_dendr = gapre*(vmpre-vmself) + gapost*(vmpost-vmself) : amp',
                                                        gapre=1 / (self.output_neuron['namespace']['Ra'][_ii]),
                                                        gapost=1 / (self.output_neuron['namespace']['Ra'][_ii - 1]),
                                                        I_dendr="Idendr_a%d" % _ii, vmself="vm_a%d" % _ii,
                                                        vmpre="vm_a%d" % (_ii + 1), vmpost="vm_a%d" % (_ii - 1))
        self.output_neuron['equation'] += Equations('I_dendr = gapost*(vmpost-vmself) : amp',
                                                    I_dendr="Idendr_a%d" % self.output_neuron['dend_comp_num'],
                                                    gapost=1 / (self.output_neuron['namespace']['Ra'][-1]),
                                                    vmself="vm_a%d" % self.output_neuron['dend_comp_num'],
                                                    vmpost="vm_a%d" % (self.output_neuron['dend_comp_num'] - 1))
        self.output_neuron['equation'] += Equations('''x : meter
                            y : meter''')

In particular the source of error is this line in the first snippet: %s += synaptic_scaling_factor * wght where %s corresponds to % (self.output_synapse['receptor'] + self.output_synapse['post_comp_name'] + '_post' which is a postsynaptic variable (it's hard to immediately spot it in the third code snippet as it's being generated dynamically) and synaptic_scaling_factor is another postsynaptic variable which is defined in the snippet as: dsynaptic_scaling_factor/dt = scaling_speed * (1 - (spike_sensor / (ap_target_frequency*tau_synaptic_scaling))) : 1 According to brian2genn , it is not allowed to use two post-synaptic variables in the synaptic model.

tnowotny commented 7 years ago

Hi again - thanks for clarifying. Your problem is that you are inverting synaptic_scaling_factor "in place", i.e. synaptic_scaling_factor = 1./synaptic_scaling_factor which means you are writing to a post-synaptic variable. Then the manipulation of the (self.output_synapse['receptor'] + self.output_synapse['post_comp_name'] + '_post' variable is a second output variable. Without knowing your exact model, are you sure you actually want to invert synaptic_scaling_factor and store the inverted value in the post-synaptic neuron? or did you only want to use a local variable sfactor that is either synaptic_scaling_factor or 1.0/synaptic_scaling_factor? ... the latter makes more sense to me and in that case there is no problem ...?

tnowotny commented 7 years ago

One more detail, I thought Brian 2 now requires a _post extension when accessing post-synaptic variables in synaptic equations? So synaptic_scaling_factor_post ...? And just to clarify, brian2genn does not allow writing to two post-synaptic variables.

Vafa-Andalibi commented 7 years ago

Apologies for late reply. Adding the _post suffix to the variable actually solved the problem. Many thanks for your support