BCPNN learning rule in GeNN

[bjoernapfel]

So I´m recently trying to implement BCPNN in GeNN as an alternative to STDP.



The basic idea is that the values for Pi/Pj are bigger if both pre- and postsyn neuron fire which leads to an increased synaptic weight. If only for example the presyn neuron fires and the postsyn neuron doesn´t, the P-value of the postsyn neuron gets lower which leads to a decreased weight. However, I have difficulties reproducing the latter aspect, my synaptic weights also increase (weaker, but still) if only one of the neurons fires. As an example, I plotted the synaptic weights between two pyramidal neurons over time. The synaptic weights additionally become a bit lower than expected.

Si represents presyn spiking events and is = 1 if the presyn neuron fires and = 0 if the presyn neuron doesn´t fire, same goes for Sj and the postsyn neuron. The other parameters are taken over from another paper where it worked totally fine so there has to be some issue with my implementation. I also checked it against STDP, STDP works just fine in my model. I have been sitting over my code for a bit and can´t figure out where it is off. Do you have any suggestions?

    BCPNN = genn_model.create_custom_weight_update_class(
    "BCPNN",
    param_names = ["tauZi", "tauZj", "tauP", "kappa", "fmax", "tspike", "epsilon", "wgain", "wMax", "wMin", "tau_ref_PN"],
    var_name_types = [("g", "scalar"), ("Zi", "scalar"), ("Zj", "scalar"), ("Pi", "scalar"), ("Pj", "scalar"), ("Pij", "scalar")],
    derived_params=[["ExpTC_Zi", genn_model.create_dpf_class(lambda pars, dt : np.exp((-dt * pars[3]) / pars[0]))()],
                    ["ExpTC_Zj", genn_model.create_dpf_class(lambda pars, dt : np.exp(-dt / pars[1]))()],
                    ["ExpTC_P", genn_model.create_dpf_class(lambda pars, dt : np.exp(-dt / pars[2]))()]],
    sim_code =
    """
    $(addToInSyn, $(g));
    scalar Si = 1;
    scalar alpha_Si = (Si / ($(fmax) * $(tspike))) + $(epsilon);
    $(Zi) = alpha_Si - $(ExpTC_Zi) * (alpha_Si - $(Zi));
    const scalar dt = $(t) - $(sT_post);
    if ( dt < $(tau_ref_PN)) {
        scalar Sj = 1;
        scalar alpha_Sj = (Sj / ($(fmax) * $(tspike))) + $(epsilon);
        $(Zj) = alpha_Sj - $(ExpTC_Zj) * (alpha_Sj - $(Zj));
    } else {
        scalar Sj = 0;
        scalar alpha_Sj = (Sj / ($(fmax) * $(tspike))) + $(epsilon);
        $(Zj) = alpha_Sj - $(ExpTC_Zj) * (alpha_Sj - $(Zj));
    }
    $(Pi) = $(Zi) - $(ExpTC_P) * ($(Zi) - $(Pi));
    $(Pj) = $(Zj) - $(ExpTC_P) * ($(Zj) - $(Pj));
    $(Pij) = ($(Zi)*$(Zj)) - $(ExpTC_P) * (($(Zi)*$(Zj)) - $(Pij));
    const scalar newWeight = $(g) + ($(wgain) * log($(Pij) / ($(Pi) * $(Pj))));
    $(g) = fmin($(wMax), fmax($(wMin), newWeight));
    """,
    learn_post_code=
    """
    const scalar dt = $(t) - $(sT_pre);
    if (dt > $(tau_ref_PN)) {
        scalar Si = 0;
        scalar Sj = 1;
        scalar alpha_Si = (Si / ($(fmax) * $(tspike))) + $(epsilon);
        scalar alpha_Sj = (Sj / ($(fmax) * $(tspike))) + $(epsilon);
        $(Zi) = alpha_Si - $(ExpTC_Zi) * (alpha_Si - $(Zi));
        $(Zj) = alpha_Sj - $(ExpTC_Zj) * (alpha_Sj - $(Zj));
        $(Pi) = $(Zi) - $(ExpTC_P) * ($(Zi) - $(Pi));
        $(Pj) = $(Zj) - $(ExpTC_P) * ($(Zj) - $(Pj));
        $(Pij) = ($(Zi)*$(Zj)) - $(ExpTC_P) * (($(Zi)*$(Zj)) - $(Pij));
        const scalar newWeight = $(g) + ($(wgain) * log($(Pij) / ($(Pi) * $(Pj))));
        $(g) = fmin($(wMax), fmax($(wMin), newWeight));    
    }
    """,
    is_pre_spike_time_required=True,
    is_post_spike_time_required=True)

[tnowotny]

Not sure yet what the problem is but I notice the following:

1. why is kappa involved in the derived parameter ExpTC_Zi and not ExpTC_P?
2. decorating $ are missing on some occasions - I am surprised this compiles at all, e.g. (t) should be $(t) (might recover by accident) but (sT_post) must be $(sT_post); same for (sT_pre) that should be $(sT_pre)
3. Why do you have a provision for Sj=1 in the sim_code? that code only executes upon pre-synaptic spikes
4. I think your maths for turning the decay term in the P_i / P_j equations into event-based decay calculation isn't correct. You are approximating Z_i to be constant between spikes while it is exponentially decaying
5. Similarly, the maths for P_ij does not work
   I suspect, 4 and 5 affect quantitative results more than qualitative results and your main problem is setting S_j= 1 in the sim_code.

[tnowotny]

If both pre and post-synaptic spikes occur in the same timestep, then both codes are executed (I seem to remember sim_code first then learn_post_code). At any time, the values of variables will be updated within the execution of the relevant snippet and then be available with updated values in later code snippets.

[tnowotny]

But do have a look at the paper referenced by @neworderofjamie below which has all the mathematical detail of how to do correct event-based updates for BCPNN.

[bjoernapfel]

Well today I tried the event-based implementation like this (used the analytical solution I. the "$()" are added wherever needed, I don´t know why they aren´t shown sometimes)

    BCPNN = genn_model.create_custom_weight_update_class(
    "BCPNN",
    param_names = ["tauZi", "tauZj", "tauZij", "tauP", "wGain", "wMax", "wMin"],
    var_name_types = [("g", "scalar"), ("Zi", "scalar"), ("Zj", "scalar"), ("Pi", "scalar"), ("Pj", "scalar"), ("Pij", "scalar")],
    derived_params=[["a_i", genn_model.create_dpf_class(lambda pars, dt: math.exp(pars[0] / (pars[0] - pars[3])))()],
                    ["a_j", genn_model.create_dpf_class(lambda pars, dt: math.exp(pars[1] / (pars[1] - pars[3])))()],
                    ["a_ij", genn_model.create_dpf_class(lambda pars, dt: math.exp(pars[2] / (pars[2] - pars[3])))()]],
    sim_code =
    """
    $(addToInSyn, $(g));
    const scalar last_update = fmax($(prev_sT_pre), $(prev_sT_post));
    const scalar deltaT = $(t) - last_update;
    const scalar ExpTC_Zi = exp(-deltaT / $(tauZi));
    const scalar ExpTC_Zj = exp(-deltaT / $(tauZj));
    const scalar ExpTC_Zij = exp(-deltaT / $(tauZij));
    const scalar ExpTC_P = exp(-deltaT / $(tauP));
    const scalar Si = 1;
    $(Zi) = $(Zi) * ExpTC_Zi + Si;
    $(Pi) = $(Pi) * ExpTC_P + $(Zi) * $(a_i) * (ExpTC_Zi - ExpTC_P);
    const scalar deltaT_post = $(t) - $(sT_post);
    if (deltaT_post > 0) {
        const scalar Sj = 0;
        $(Zj) = $(Zj) * ExpTC_Zj + Sj;
        $(Pj) = $(Pj) * ExpTC_P + $(Zj) * $(a_j) * (ExpTC_Zj - ExpTC_P);
        $(Pij) = $(Pij) * ExpTC_P + $(Zi) * $(Zj) * $(a_ij) * (ExpTC_Zij - ExpTC_P);
        const scalar newWeight = $(g) + ($(wGain) * log($(Pij) / ($(Pi) * $(Pj))));
        $(g) = fmin($(wMax), fmax($(wMin), newWeight));
    }
    """,
    learn_post_code=
    """
    const scalar last_update = fmax($(prev_sT_pre), $(prev_sT_post));
    const scalar deltaT = $(t) - last_update;
    const scalar ExpTC_Zi = exp(-deltaT / $(tauZi));
    const scalar ExpTC_Zj = exp(-deltaT / $(tauZj));
    const scalar ExpTC_Zij = exp(-deltaT / $(tauZij));
    const scalar ExpTC_P = exp(-deltaT / $(tauP));
    const scalar Sj = 1;
    $(Zj) = $(Zj) * ExpTC_Zj + Sj;
    $(Pj) = $(Pj) * ExpTC_P + $(Zj) * $(a_j) * (ExpTC_Zj - ExpTC_P);
    const scalar deltaT_pre = $(t) - $(sT_pre);
        if (deltaT_pre > 0) {
            const scalar Si = 0;
            $(Zi) = $(Zi) * ExpTC_Zi + Si;
            $(Pi) = $(Pi) * ExpTC_P + $(Zi) * $(a_i) * (ExpTC_Zi - ExpTC_P);
        }
    $(Pij) = $(Pij) * ExpTC_P + $(Zi) * $(Zj) * $(a_ij) * (ExpTC_Zij - ExpTC_P);
    const scalar newWeight = $(g) + ($(wGain) * log($(Pij) / ($(Pi) * $(Pj))));
    $(g) = fmin($(wMax), fmax($(wMin), newWeight));    
    """,
    is_pre_spike_time_required=True,
    is_post_spike_time_required=True,
    is_prev_pre_spike_time_required=True,
    is_prev_post_spike_time_required=True).

The results look a bit better but it feels off that the weights rise or almost immediately drop to 0.

[neworderofjamie]

const scalar deltaT_post = $(t) - $(sT_post);

and

const scalar deltaT_pre = $(t) - $(sT_pre);

will both always be zero - if you switch to $(prev_sT_pre) and $(prev_sT_post) (basically the time of the spike before the one currently being processed) this bit should work better.

Worth getting Analytical 1 working but, as I remember, Analytical II ends up being much simpler.

[tnowotny]

A couple more observations:

1. Why are there math.exp() expressions in your derived parameters - I thought those parameters were simple fractions
2. It would be much easier (and also valid) to update the pre-synaptic traces in sim_code only and the post-synaptic ones in post_learn_code. Then you don't need to find out which (pre or post) was the last spike but use last_update=deltaT= $(t) - $(prev_sT_pre) in sim_code and last_update= deltaT= $(t) - $(prev_sT_post); in post_learn_code.
3. if you do it the way you have, const scalar last_update = fmax($(prev_sT_pre), $(prev_sT_post)); isn't quite right ... should be const scalar last_update = fmax($(prev_sT_pre), $(sT_post)); in sim_code and const scalar last_update = fmax($(sT_pre), $(prev_sT_post)); in post_learn_code.

[neworderofjamie]

I believe http://www.frontiersin.org/Neuromorphic_Engineering/10.3389/fnins.2015.00002/abstract explains the event-based implementation

[neworderofjamie]

A long time ago, I did actually make a GeNN implementation of a two trace (so Z and P) BCPNN learning rule using analytical II at https://github.com/neworderofjamie/genn_examples/blob/master/common/bcpnn.h. Sadly it predates PyGeNN so is in C++ but porting it should be pretty easy and the learning rule it appears to work well enough to reproduce figure 4 of Tully 2014 using https://github.com/neworderofjamie/genn_examples/tree/master/bcpnn_simple). It also predates the $(prev_sT_pre) and $(prev_sT_post) variables which could be used to replace lastUpdate and make it more efficient.

[bjoernapfel]

thanks for the input, seems to work now!