# Example: spike_based_homeostasis

> 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):

Following O. Breitwieser: “Towards a Neuromorphic Implementation of Spike-Based Expectation Maximization”

Two poisson stimuli are connected to a neuron. One with a varying rate and the other with a fixed rate. The synaptic weight from the varying rate stimulus to the neuron is fixed. The synaptic weight from the fixed rate stimulus to the neuron is plastic and tries to keep the neuron at a firing rate that is determined by the parameters of the plasticity rule.

Sebastian Schmitt, 2021

    import itertools
    import numpy as np
    import matplotlib.pyplot as plt

    from brian2 import TimedArray, PoissonGroup, NeuronGroup, Synapses, StateMonitor, PopulationRateMonitor
    from brian2 import defaultclock, run
    from brian2 import Hz, ms, second

    # The synaptic weight from the steady stimulus is plastic
    steady_stimulus = TimedArray([50]*Hz, dt=40*second)
    steady_poisson = PoissonGroup(1, rates='steady_stimulus(t)')

    # The synaptic weight from the varying stimulus is static
    varying_stimulus = TimedArray([25*Hz, 50*Hz, 0*Hz, 35*Hz, 0*Hz], dt=10*second)
    varying_poisson = PoissonGroup(1, rates='varying_stimulus(t)')

    # dw_plus/dw_minus determines scales the steady stimulus rate to the target firing rate, must not be larger 1
    # the magntude of dw_plus and dw_minus determines the "speed" of the homeostasis
    parameters = {
        'tau': 10*ms,  # membrane time constant
        'dw_plus': 0.05,  # weight increment on pre spike
        'dw_minus': 0.05,  # weight increment on post spike
        'w_max': 2,  # maximum plastic weight
        'w_initial': 0  # initial plastic weight
    }

    eqs = 'dv/dt = (0 - v)/tau : 1 (unless refractory)'

    neuron_with_homeostasis = NeuronGroup(1, eqs,
                                          threshold='v > 1', reset='v = -1',
                                          method='euler', refractory=1*ms,
                                          namespace=parameters)
    neuron_without_homeostasis = NeuronGroup(1, eqs,
                                             threshold='v > 1', reset='v = -1',
                                             method='euler', refractory=1*ms,
                                             namespace=parameters)

    plastic_synapse = Synapses(steady_poisson, neuron_with_homeostasis,
                               'w : 1',
                               on_pre='''
                               v_post += w
                               w = clip(w + dw_plus, 0, w_max)
                               ''',
                               on_post='''
                               w = clip(w - dw_minus, 0, w_max)
                               ''', namespace=parameters)
    plastic_synapse.connect()
    plastic_synapse.w = parameters['w_initial']

    non_plastic_synapse_neuron_without_homeostasis = Synapses(varying_poisson,
                                                              neuron_without_homeostasis,
                                                              'w : 1', on_pre='v_post += w')
    non_plastic_synapse_neuron_without_homeostasis.connect()
    non_plastic_synapse_neuron_without_homeostasis.w = 2

    non_plastic_synapse_neuron = Synapses(varying_poisson, neuron_with_homeostasis,
                                          'w : 1', on_pre='v_post += w')
    non_plastic_synapse_neuron.connect()
    non_plastic_synapse_neuron.w = 2

    M = StateMonitor(neuron_with_homeostasis, 'v', record=True)
    M2 = StateMonitor(plastic_synapse, 'w', record=True)
    M_rate_neuron_with_homeostasis = PopulationRateMonitor(neuron_with_homeostasis)
    M_rate_neuron_without_homeostasis = PopulationRateMonitor(neuron_without_homeostasis)

    duration = 40*second
    defaultclock.dt = 0.1*ms
    run(duration, report='text')

    fig, axes = plt.subplots(3, sharex=True)

    axes[0].plot(M2.t/second, M2.w[0], label="homeostatic weight")
    axes[0].set_ylabel("weight")
    axes[0].legend()

    # dt is in second
    dts = np.arange(0., len(varying_stimulus.values)*varying_stimulus.dt, varying_stimulus.dt)
    x = list(itertools.chain(*zip(dts, dts)))
    y = list(itertools.chain(*zip(varying_stimulus.values/Hz, varying_stimulus.values/Hz)))
    axes[1].plot(x, [0] + y[:-1], label="varying stimulus")
    axes[1].set_ylabel("rate [Hz]")
    axes[1].legend()

    # in ms
    smooth_width = 100*ms
    axes[2].plot(M_rate_neuron_with_homeostasis.t/second,
                 M_rate_neuron_with_homeostasis.smooth_rate(width=smooth_width)/Hz,
                 label="with homeostasis")
    axes[2].plot(M_rate_neuron_without_homeostasis.t/second,
                 M_rate_neuron_without_homeostasis.smooth_rate(width=smooth_width)/Hz,
                 label="without homeostasis")
    axes[2].set_ylabel("firing rate [Hz]")
    axes[2].legend()

    plt.xlabel('Time (s)')
    plt.show()