{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Weight adaptation according to the Urbanczik-Senn plasticity\n\nThis script demonstrates the learning in a compartmental neuron where the\ndendritic synapses adapt their weight according to the plasticity rule by\nUrbanczik and Senn [1]. In this simple setup, a spike pattern of 200 poisson\nspike trains is repeatedly presented to a neuron that is composed of one\nsomatic and one dendritic compartment. At the same time, the somatic\nconductances are activated to produce a time-varying matching potential.\nAfter the learning, this signal is then reproreproduced by the membrane\npotential of the neuron. This script produces Fig. 1B in [1] but uses standard\nunits instead of the unitless quantities used in the paper.\n\n[1] R. Urbanczik, W. Senn (2014): Learning by the Dendritic Prediction of\n    Somatic Spiking. Neuron, 81, 521-528.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom matplotlib import pyplot as plt\nimport nest\n\n\ndef g_inh(amplitude, t_start, t_end):\n    '''\n    returns weights for the spike generator that drives the inhibitory\n    somatic conductance.\n    '''\n    return lambda t: np.piecewise(t, [(t >= t_start) & (t < t_end)],\n                                  [amplitude, 0.0])\n\n\ndef g_exc(amplitude, freq, offset, t_start, t_end):\n    '''\n    returns weights for the spike generator that drives the excitatory\n    somatic conductance.\n    '''\n    return lambda t: np.piecewise(t, [(t >= t_start) & (t < t_end)],\n                                  [lambda t: amplitude*np.sin(freq*t) + offset, 0.0])\n\n\ndef matching_potential(g_E, g_I, nrn_params):\n    '''\n    returns the matching potential as a function of the somatic conductances.\n    '''\n    E_E = nrn_params['soma']['E_ex']\n    E_I = nrn_params['soma']['E_in']\n    return (g_E*E_E + g_I*E_I) / (g_E + g_I)\n\n\ndef V_w_star(V_w, nrn_params):\n    '''\n    returns the dendritic prediction of the somatic membrane potential.\n    '''\n    g_D = nrn_params['g_sp']\n    g_L = nrn_params['soma']['g_L']\n    E_L = nrn_params['soma']['E_L']\n    return (g_L*E_L + g_D*V_w) / (g_L + g_D)\n\n\ndef phi(U, nrn_params):\n    '''\n    rate function of the soma\n    '''\n    phi_max = nrn_params['phi_max']\n    k = nrn_params['rate_slope']\n    beta = nrn_params['beta']\n    theta = nrn_params['theta']\n    return phi_max / (1.0 + k*np.exp(beta*(theta - U)))\n\n\ndef h(U, nrn_params):\n    '''\n    derivative of the rate function phi\n    '''\n    phi_max = nrn_params['phi_max']\n    k = nrn_params['rate_slope']\n    beta = nrn_params['beta']\n    theta = nrn_params['theta']\n    return 15.0*beta / (1.0 + np.exp(-beta*(theta - U)) / k)\n\n\n'''\nsimulation params\n'''\nn_pattern_rep = 100         # number of repetitions of the spike pattern\npattern_duration = 200.0\nt_start = 2.0*pattern_duration\nt_end = n_pattern_rep*pattern_duration + t_start\nsimulation_time = t_end + 2.0*pattern_duration\nn_rep_total = int(np.around(simulation_time / pattern_duration))\nresolution = 0.1\nnest.SetKernelStatus({'resolution': resolution})\n\n'''\nneuron parameters\n'''\nnrn_model = 'pp_cond_exp_mc_urbanczik'\nnrn_params = {\n    't_ref': 3.0,        # refractory period\n    'g_sp': 600.0,       # soma-to-dendritic coupling conductance\n    'soma': {\n        'V_m': -70.0,    # initial value of V_m\n        'C_m': 300.0,    # capacitance of membrane\n        'E_L': -70.0,    # resting potential\n        'g_L': 30.0,     # somatic leak conductance\n        'E_ex': 0.0,     # resting potential for exc input\n        'E_in': -75.0,   # resting potential for inh input\n        'tau_syn_ex': 3.0,  # time constant of exc conductance\n        'tau_syn_in': 3.0,  # time constant of inh conductance\n    },\n    'dendritic': {\n        'V_m': -70.0,    # initial value of V_m\n        'C_m': 300.0,    # capacitance of membrane\n        'E_L': -70.0,    # resting potential\n        'g_L': 30.0,     # dendritic leak conductance\n        'tau_syn_ex': 3.0,  # time constant of exc input current\n        'tau_syn_in': 3.0,  # time constant of inh input current\n    },\n    # parameters of rate function\n    'phi_max': 0.15,     # max rate\n    'rate_slope': 0.5,   # called 'k' in the paper\n    'beta': 1.0 / 3.0,\n    'theta': -55.0,\n}\n\n'''\nsynapse params\n'''\nsyns = nest.GetDefaults(nrn_model)['receptor_types']\ninit_w = 0.3*nrn_params['dendritic']['C_m']\nsyn_params = {\n    'synapse_model': 'urbanczik_synapse_wr',\n    'receptor_type': syns['dendritic_exc'],\n    'tau_Delta': 100.0,  # time constant of low pass filtering of the weight change\n    'eta': 0.17,         # learning rate\n    'weight': init_w,\n    'Wmax': 4.5*nrn_params['dendritic']['C_m'],\n    'delay': resolution,\n}\n\n'''\n# in case you want to use the unitless quantities as in [1]:\n\n# neuron params:\n\nnrn_model = 'pp_cond_exp_mc_urbanczik'\nnrn_params = {\n    't_ref': 3.0,\n    'g_sp': 2.0,\n    'soma': {\n        'V_m': 0.0,\n        'C_m': 1.0,\n        'E_L': 0.0,\n        'g_L': 0.1,\n        'E_ex': 14.0 / 3.0,\n        'E_in': -1.0 / 3.0,\n        'tau_syn_ex': 3.0,\n        'tau_syn_in': 3.0,\n    },\n    'dendritic': {\n        'V_m': 0.0,\n        'C_m': 1.0,\n        'E_L': 0.0,\n        'g_L': 0.1,\n        'tau_syn_ex': 3.0,\n        'tau_syn_in': 3.0,\n    },\n    # parameters of rate function\n    'phi_max': 0.15,\n    'rate_slope': 0.5,\n    'beta': 5.0,\n    'theta': 1.0,\n}\n\n# synapse params:\n\nsyns = nest.GetDefaults(nrn_model)['receptor_types']\ninit_w = 0.2*nrn_params['dendritic']['g_L']\nsyn_params = {\n    'synapse_model': 'urbanczik_synapse_wr',\n    'receptor_type': syns['dendritic_exc'],\n    'tau_Delta': 100.0,\n    'eta': 0.0003 / (15.0*15.0*nrn_params['dendritic']['C_m']),\n    'weight': init_w,\n    'Wmax': 3.0*nrn_params['dendritic']['g_L'],\n    'delay': resolution,\n}\n'''\n\n'''\nsomatic input\n'''\nampl_exc = 0.016*nrn_params['dendritic']['C_m']\noffset = 0.018*nrn_params['dendritic']['C_m']\nampl_inh = 0.06*nrn_params['dendritic']['C_m']\nfreq = 2.0 / pattern_duration\nsoma_exc_inp = g_exc(ampl_exc, 2.0*np.pi*freq, offset, t_start, t_end)\nsoma_inh_inp = g_inh(ampl_inh, t_start, t_end)\n\n'''\ndendritic input\ncreate spike pattern by recording the spikes of a simulation of n_pg\npoisson generators. The recorded spike times are then given to spike\ngenerators.\n'''\nn_pg = 200                 # number of poisson generators\np_rate = 10.0              # rate in Hz\n\npgs = nest.Create('poisson_generator', n=n_pg, params={'rate': p_rate})\n\nprrt_nrns_pg = nest.Create('parrot_neuron', n_pg)\n\nnest.Connect(pgs, prrt_nrns_pg, {'rule': 'one_to_one'})\n\nsr = nest.Create('spike_recorder', n_pg)\nnest.Connect(prrt_nrns_pg, sr, {'rule': 'one_to_one'})\n\nnest.Simulate(pattern_duration)\nt_srs = []\nfor i, ssr in enumerate(nest.GetStatus(sr)):\n    t_sr = ssr['events']['times']\n    t_srs.append(t_sr)\n\nnest.ResetKernel()\nnest.SetKernelStatus({'resolution': resolution})\n\n'''\nneuron and devices\n'''\nnest.SetDefaults(nrn_model, nrn_params)\nnrn = nest.Create(nrn_model)\n\n# poisson generators are connected to parrot neurons which are\n# connected to the mc neuron\nprrt_nrns = nest.Create('parrot_neuron', n_pg)\n\n# excitatory input to the soma\nspike_times_soma_inp = np.arange(resolution, simulation_time, resolution)\nsg_soma_exc = nest.Create('spike_generator',\n                          params={'spike_times': spike_times_soma_inp,\n                                  'spike_weights': soma_exc_inp(spike_times_soma_inp)})\n# inhibitory input to the soma\nsg_soma_inh = nest.Create('spike_generator',\n                          params={'spike_times': spike_times_soma_inp,\n                                  'spike_weights': soma_inh_inp(spike_times_soma_inp)})\n\n# excitatory input to the dendrite\nsg_prox = nest.Create('spike_generator', n=n_pg)\n\n# for recording all parameters of the Urbanczik neuron\nrqs = nest.GetDefaults(nrn_model)['recordables']\nmm = nest.Create('multimeter', params={'record_from': rqs, 'interval': 0.1})\n\n# for recoding the synaptic weights of the Urbanczik synapses\nwr = nest.Create('weight_recorder')\n\n# for recording the spiking of the soma\nsr_soma = nest.Create('spike_recorder')\n\n'''\ncreate connections\n'''\nnest.Connect(sg_prox, prrt_nrns, {'rule': 'one_to_one'})\nnest.CopyModel('urbanczik_synapse', 'urbanczik_synapse_wr',\n               {'weight_recorder': wr[0]})\nnest.Connect(prrt_nrns, nrn, syn_spec=syn_params)\nnest.Connect(mm, nrn, syn_spec={'delay': 0.1})\nnest.Connect(sg_soma_exc, nrn,\n             syn_spec={'receptor_type': syns['soma_exc'], 'weight': 10.0*resolution, 'delay': resolution})\nnest.Connect(sg_soma_inh, nrn,\n             syn_spec={'receptor_type': syns['soma_inh'], 'weight': 10.0*resolution, 'delay': resolution})\nnest.Connect(nrn, sr_soma)\n\n'''\nsimulation divided into intervals of the pattern duration\n'''\nfor i in np.arange(n_rep_total):\n    # Set the spike times of the pattern for each spike generator\n    for (sg, t_sp) in zip(sg_prox, t_srs):\n        nest.SetStatus(\n            sg, {'spike_times': np.array(t_sp) + i*pattern_duration})\n\n    nest.Simulate(pattern_duration)\n\n'''\nread out devices\n'''\n# multimeter\nrec = nest.GetStatus(mm)[0]['events']\nt = rec['times']\nV_s = rec['V_m.s']\nV_d = rec['V_m.p']\nV_d_star = V_w_star(V_d, nrn_params)\ng_in = rec['g_in.s']\ng_ex = rec['g_ex.s']\nI_ex = rec['I_ex.p']\nI_in = rec['I_in.p']\nU_M = matching_potential(g_ex, g_in, nrn_params)\n\n# weight recorder\ndata = nest.GetStatus(wr)\nsenders = data[0]['events']['senders']\ntargets = data[0]['events']['targets']\nweights = data[0]['events']['weights']\ntimes = data[0]['events']['times']\n\n# spike recorder\ndata = nest.GetStatus(sr_soma)[0]['events']\nspike_times_soma = data['times']\n\n'''\nplot results\n'''\nfs = 22\nlw = 2.5\nfig1, (axA, axB, axC, axD) = plt.subplots(4, 1, sharex=True)\n\n# membrane potentials and matching potential\naxA.plot(t, V_s, lw=lw, label=r'$U$ (soma)', color='darkblue')\naxA.plot(t, V_d, lw=lw, label=r'$V_W$ (dendrit)', color='deepskyblue')\naxA.plot(t, V_d_star, lw=lw, label=r'$V_W^\\ast$ (dendrit)', color='b', ls='--')\naxA.plot(t, U_M, lw=lw, label=r'$U_M$ (soma)', color='r', ls='-')\naxA.set_ylabel('membrane pot [mV]', fontsize=fs)\naxA.legend(fontsize=fs)\n\n# somatic conductances\naxB.plot(t, g_in, lw=lw, label=r'$g_I$', color='r')\naxB.plot(t, g_ex, lw=lw, label=r'$g_E$', color='coral')\naxB.set_ylabel('somatic cond', fontsize=fs)\naxB.legend(fontsize=fs)\n\n# dendritic currents\naxC.plot(t, I_ex, lw=lw, label=r'$I_ex$', color='r')\naxC.plot(t, I_in, lw=lw, label=r'$I_in$', color='coral')\naxC.set_ylabel('dend current', fontsize=fs)\naxC.legend(fontsize=fs)\n\n# rates\naxD.plot(t, phi(V_s, nrn_params), lw=lw, label=r'$\\phi(U)$', color='darkblue')\naxD.plot(t, phi(V_d, nrn_params), lw=lw,\n         label=r'$\\phi(V_W)$', color='deepskyblue')\naxD.plot(t, phi(V_d_star, nrn_params), lw=lw,\n         label=r'$\\phi(V_W^\\ast)$', color='b', ls='--')\naxD.plot(t, h(V_d_star, nrn_params), lw=lw,\n         label=r'$h(V_W^\\ast)$', color='g', ls='--')\naxD.plot(t, phi(V_s, nrn_params) - phi(V_d_star, nrn_params), lw=lw,\n         label=r'$\\phi(U) - \\phi(V_W^\\ast)$', color='r', ls='-')\naxD.plot(spike_times_soma, 0.0*np.ones(len(spike_times_soma)),\n         's', color='k', markersize=2)\naxD.legend(fontsize=fs)\n\n# synaptic weights\nfig2, axA = plt.subplots(1, 1)\nfor i in np.arange(2, 200, 10):\n    index = np.intersect1d(np.where(senders == i), np.where(targets == 1))\n    if not len(index) == 0:\n        axA.step(times[index], weights[index], label='pg_{}'.format(i - 2),\n                 lw=lw)\n\naxA.set_title('Synaptic weights of Urbanczik synapses')\naxA.set_xlabel('time [ms]', fontsize=fs)\naxA.set_ylabel('weight', fontsize=fs)\naxA.legend(fontsize=fs - 4)\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}