.. currentmodule:: brian2

.. homeostatic_stdp_at_inhibitory_synapes:

Example: homeostatic_stdp_at_inhibitory_synapes
===============================================


        .. only:: html

            .. |launchbinder| image:: http://mybinder.org/badge.svg
            .. _launchbinder: https://mybinder.org/v2/gh/brian-team/brian2-binder/master?filepath=examples/synapses/homeostatic_stdp_at_inhibitory_synapes.ipynb

            .. note::
               You can launch an interactive, editable version of this
               example without installing any local files
               using the Binder service (although note that at some times this
               may be slow or fail to open): |launchbinder|_

        

This script implements a homeostatic STDP rule for inhibitory
synapses onto excitatory neurons presented in Zenke et al. (2015) and
explained in detail in the accompanying supplementary information.

In essence, the magnitude and sign of the STDP kernel at inhibitory
synapses targeting excitatory neurons is modulated by a global factor
G, which corresponds to the difference between the average firing rate
of the excitatory postsynaptic neurons and a target firing rate
:math:`\gamma`. If the average firing rate of excitatory neurons exceeds the
target firing rate, pre/post pairings (and post/pre pairings) at
inhibitory synapses connecting to excitatory neurons result in
potentiation; depression of the inhibitory synapses occurs when the
excitatory neurons are unable to match the desired target rate.

The script emulates classic STDP protocols, simulating a connected
presynaptic and postsynaptic neuron pair with precisely controlled
spike times, and a large, independently firing population of
excitatory neurons.

Paul Brodersen, 2025

Reference:
Zenke, Friedemann, Everton J. Agnes, and Wulfram Gerstner.
"Diverse Synaptic Plasticity Mechanisms Orchestrated to Form and
Retrieve Memories in Spiking Neural Networks." Nature Communications 6,
no. 1 (April 21, 2015): 6922. https://doi.org/10.1038/ncomms7922.

::

    
    import numpy as np
    import matplotlib.pyplot as plt
    import brian2 as b2
    
    net = b2.Network()
    
    # excitatory neuron population
    epoch_length = 1 * b2.second
    spike_rate = np.arange(0, 12.5, 2.5)
    total_epochs = len(spike_rate)
    population_spike_rate = b2.TimedArray(spike_rate * b2.Hz, dt=epoch_length)
    population_size = 1000
    spike_generators = b2.PoissonGroup(population_size, rates="population_spike_rate(t)")
    net.add(spike_generators)
    
    # global factor G
    gamma = 5 * b2.Hz
    tau_H = 0.1 * b2.second # NB: value in Zenke et al (2015) is 10 seconds; reduced here to allow for shorter simulations
    global_factor_model = """
    G = H - gamma    : Hz
    dH/dt = -H/tau_H : Hz
    """
    global_factor = b2.NeuronGroup(1, model=global_factor_model)
    net.add(global_factor)
    spike_aggregator = b2.Synapses(spike_generators, global_factor, on_pre="H_post += 1/(tau_H * population_size)")
    spike_aggregator.connect()
    net.add(spike_aggregator)
    
    # pre-synaptic neuron
    spike_times = epoch_length * (np.arange(total_epochs) + 0.5)
    spike_indices = np.zeros_like(spike_times / b2.second)
    presynaptic_neuron = b2.SpikeGeneratorGroup(1, spike_indices, spike_times)
    net.add(presynaptic_neuron)
    
    # postsynaptic neuron
    delay = 1 * b2.msecond
    postsynaptic_neuron = b2.SpikeGeneratorGroup(1, spike_indices, spike_times + delay)
    net.add(postsynaptic_neuron)
    
    # synapse
    tau_stdp = 20. * b2.msecond
    eta = 1. * b2.nsiemens * b2.second
    synapse_model = """
    G : Hz (linked)
    wij : siemens
    dzi/dt = -zi / tau_stdp : 1 (event-driven)
    dzj/dt = -zj / tau_stdp : 1 (event-driven)
    """
    on_pre = """
    zj += 1.
    wij += eta * G * (zi + 1)
    """
    on_post = """
    zi += 1.
    wij += eta * G * zj
    """
    synapse = b2.Synapses(
        presynaptic_neuron, postsynaptic_neuron,
        model=synapse_model, on_pre=on_pre, on_post=on_post
    )
    synapse.connect(i=0, j=0)
    synapse.G = b2.linked_var(global_factor, "G")
    synapse.wij = 0 * b2.nsiemens
    net.add(synapse)
    
    # monitor
    monitor = b2.StateMonitor(synapse, ["G", "wij"], record=True, dt=1*b2.msecond)
    net.add(monitor)
    
    # run
    net.run(total_epochs * epoch_length)
    
    # analysis
    G = np.squeeze(monitor.G / b2.Hz)
    wij = np.squeeze(monitor.wij / b2.nsiemens)
    time = np.squeeze(monitor.t / b2.ms)
    spike_times = spike_times / b2.ms
    
    # visualize experimental protocol
    fig, axes = plt.subplots(2, 1, sharex=True)
    fig.suptitle("Experimental protocol")
    axes[0].plot(time, G)
    axes[0].set_ylabel("G [Hz]")
    axes[1].plot(time, wij)
    axes[1].set_ylabel("Weight [nS]")
    for ax in axes:
        ax.vlines(spike_times, *ax.get_ylim(), color="gray", linestyle=":", label="Pre/post spike pairings")
    axes[0].legend(loc="upper left")
    axes[-1].set_xlabel("Time [ms]")
    
    # plot relationship between G and \Delta wij
    g = np.zeros_like(spike_times)
    dw = np.zeros_like(spike_times)
    for ii, spike_time in enumerate(spike_times):
        g[ii] = G[np.argmin(np.abs(time - spike_time))]
        w0 = wij[np.argmin(np.abs(time - spike_time - 10))]
        w1 = wij[np.argmin(np.abs(time - spike_time + 10))]
        dw[ii] = w1 - w0
    
    fig, ax = plt.subplots()
    ax.plot(g, dw, 'o')
    ax.set_xlabel("G [Hz]")
    ax.set_ylabel(r"$\Delta$weight [nS]")
    ax.set_title("Weight change modulation")
    
    plt.show()
    

.. image:: ../resources/examples_images/synapses.homeostatic_stdp_at_inhibitory_synapes.1.png

.. image:: ../resources/examples_images/synapses.homeostatic_stdp_at_inhibitory_synapes.2.png