{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Population rate model of generalized integrate-and-fire neurons\n\nThis script simulates a finite network of generalized integrate-and-fire\n(GIF) neurons directly on the mesoscopic population level using the effective\nstochastic population rate dynamics derived in the paper [1]_. The stochastic\npopulation dynamics is implemented in the NEST model gif_pop_psc_exp. We\ndemonstrate this model using the example of a Brunel network of two coupled\npopulations, one excitatory and one inhibitory population.\n\nNote that the population model represents the mesoscopic level\ndescription of the corresponding microscopic network based on the\nNEST model ``gif_psc_exp``.\n\n## References\n\n.. [1] Schwalger T, Degert M, Gerstner W (2017). Towards a theory of cortical columns: From spiking\n       neurons to interacting neural populations of finite size. PLoS Comput Biol.\n       https://doi.org/10.1371/journal.pcbi.1005507\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Loading the necessary modules:\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport nest"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We first set the parameters of the microscopic model:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# All times given in milliseconds\ndt = 0.5\ndt_rec = 1.\n\n# Simulation time\nt_end = 2000.\n\n# Parameters\nsize = 200\nN = np.array([4, 1]) * size\nM = len(N)  # number of populations\n\n# neuronal parameters\nt_ref = 4. * np.ones(M)  # absolute refractory period\ntau_m = 20 * np.ones(M)  # membrane time constant\nmu = 24. * np.ones(M)    # constant base current mu=R*(I0+Vrest)\nc = 10. * np.ones(M)     # base rate of exponential link function\nDelta_u = 2.5 * np.ones(M)   # softness of exponential link function\nV_reset = 0. * np.ones(M)    # Reset potential\nV_th = 15. * np.ones(M)      # baseline threshold (non-accumulating part)\ntau_sfa_exc = [100., 1000.]  # adaptation time constants of excitatory neurons\ntau_sfa_inh = [100., 1000.]  # adaptation time constants of inhibitory neurons\nJ_sfa_exc = [1000., 1000.]   # size of feedback kernel theta\n#                              (= area under exponential) in mV*ms\nJ_sfa_inh = [1000., 1000.]   # in mV*ms\ntau_theta = np.array([tau_sfa_exc, tau_sfa_inh])\nJ_theta = np.array([J_sfa_exc, J_sfa_inh])\n\n# connectivity\nJ = 0.3  # excitatory synaptic weight in mV if number of input connections\n#          is C0 (see below)\ng = 5.   # inhibition-to-excitation ratio\npconn = 0.2 * np.ones((M, M))\ndelay = 1. * np.ones((M, M))\n\nC0 = np.array([[800, 200], [800, 200]]) * 0.2  # constant reference matrix\nC = np.vstack((N, N)) * pconn  # numbers of input connections\n\n# final synaptic weights scaling as 1/C\nJ_syn = np.array([[J, -g * J], [J, -g * J]]) * C0 / C\n\ntaus1_ = [3., 6.]  # time constants of exc./inh. postsynaptic currents (PSC's)\ntaus1 = np.array([taus1_ for k in range(M)])\n\n\n# step current input\nstep = [[20.], [20.]]  # jump size of mu in mV\ntstep = np.array([[1500.], [1500.]])  # times of jumps\n\n# synaptic time constants of excitatory and inhibitory connections\ntau_ex = 3.  # in ms\ntau_in = 6.  # in ms"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Simulation on the mesoscopic level\n\nTo directly simulate the mesoscopic population activities (i.e. generating\nthe activity of a finite-size population without simulating single\nneurons), we can build the populations using the NEST model\n``gif_pop_psc_exp``:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.set_verbosity(\"M_WARNING\")\nnest.ResetKernel()\nnest.resolution = dt\nnest.print_time = True\nnest.local_num_threads = 1\n\nt0 = nest.biological_time\n\nnest_pops = nest.Create('gif_pop_psc_exp', M)\n\nC_m = 250.  # irrelevant value for membrane capacity, cancels out in simulation\ng_L = C_m / tau_m\n\nparams = [{\n    'C_m': C_m,\n    'I_e': mu[i] * g_L[i],\n    'lambda_0': c[i],  # in Hz!\n    'Delta_V': Delta_u[i],\n    'tau_m': tau_m[i],\n    'tau_sfa': tau_theta[i],\n    'q_sfa': J_theta[i] / tau_theta[i],  # [J_theta]= mV*ms -> [q_sfa]=mV\n    'V_T_star': V_th[i],\n    'V_reset': V_reset[i],\n    'len_kernel': -1,  # -1 triggers automatic history size\n    'N': N[i],\n    't_ref': t_ref[i],\n    'tau_syn_ex': max([tau_ex, dt]),\n    'tau_syn_in': max([tau_in, dt]),\n    'E_L': 0.\n} for i in range(M)]\nnest_pops.set(params)\n\n# connect the populations\ng_syn = np.ones_like(J_syn)  # synaptic conductance\ng_syn[:, 0] = C_m / tau_ex\ng_syn[:, 1] = C_m / tau_in\nfor i in range(M):\n    for j in range(M):\n        nest.Connect(nest_pops[j], nest_pops[i],\n                     syn_spec={'weight': J_syn[i, j] * g_syn[i, j] * pconn[i, j],\n                               'delay': delay[i, j]})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To record the instantaneous population rate `Abar(t)` we use a multimeter,\nand to get the population activity `A_N(t)` we use spike recorder:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# monitor the output using a multimeter, this only records with dt_rec!\nnest_mm = nest.Create('multimeter')\nnest_mm.set(record_from=['n_events', 'mean'], interval=dt_rec)\nnest.Connect(nest_mm, nest_pops)\n\n# monitor the output using a spike recorder\nnest_sr = []\nfor i in range(M):\n    nest_sr.append(nest.Create('spike_recorder'))\n    nest_sr[i].time_in_steps = True\n    nest.Connect(nest_pops[i], nest_sr[i], syn_spec={'weight': 1., 'delay': dt})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "All neurons in a given population will be stimulated with a step input\ncurrent:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# set initial value (at t0+dt) of step current generator to zero\ntstep = np.hstack((dt * np.ones((M, 1)), tstep))\nstep = np.hstack((np.zeros((M, 1)), step))\n\n# create the step current devices\nnest_stepcurrent = nest.Create('step_current_generator', M)\n# set the parameters for the step currents\nfor i in range(M):\n    nest_stepcurrent[i].set(amplitude_times=tstep[i] + t0,\n                            amplitude_values=step[i] * g_L[i],\n                            origin=t0,\n                            stop=t_end)\n    pop_ = nest_pops[i]\n    nest.Connect(nest_stepcurrent[i], pop_, syn_spec={'weight': 1., 'delay': dt})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now start the simulation:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.rng_seed = 1\n\nt = np.arange(0., t_end, dt_rec)\nA_N = np.ones((t.size, M)) * np.nan\nAbar = np.ones_like(A_N) * np.nan\n\n# simulate 1 step longer to make sure all t are simulated\nnest.Simulate(t_end + dt)\ndata_mm = nest_mm.events\nfor i, nest_i in enumerate(nest_pops):\n    a_i = data_mm['mean'][data_mm['senders'] == nest_i.global_id]\n    a = a_i / N[i] / dt\n    min_len = np.min([len(a), len(Abar)])\n    Abar[:min_len, i] = a[:min_len]\n\n    data_sr = nest_sr[i].get('events', 'times')\n    data_sr = data_sr * dt - t0\n    bins = np.concatenate((t, np.array([t[-1] + dt_rec])))\n    A = np.histogram(data_sr, bins=bins)[0] / float(N[i]) / dt_rec\n    A_N[:, i] = A"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "and plot the activity:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(1)\nplt.clf()\nplt.subplot(2, 1, 1)\nplt.plot(t, A_N * 1000)  # plot population activities (in Hz)\nplt.ylabel(r'$A_N$ [Hz]')\nplt.title('Population activities (mesoscopic sim.)')\nplt.subplot(2, 1, 2)\nplt.plot(t, Abar * 1000)  # plot instantaneous population rates (in Hz)\nplt.ylabel(r'$\\bar A$ [Hz]')\nplt.xlabel('time [ms]')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Microscopic (\"direct\") simulation\n\nAs mentioned above, the population model ``gif_pop_psc_exp`` directly\nsimulates the mesoscopic population activities, i.e. without the need to\nsimulate single neurons. On the other hand, if we want to know single\nneuron activities, we must simulate on the microscopic level. This is\npossible by building a corresponding network of ``gif_psc_exp`` neuron models:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.ResetKernel()\nnest.resolution = dt\nnest.print_time = True\nnest.local_num_threads = 1\n\nt0 = nest.biological_time\n\nnest_pops = []\nfor k in range(M):\n    nest_pops.append(nest.Create('gif_psc_exp', N[k]))\n\n# set single neuron properties\nfor i in range(M):\n    nest_pops[i].set(C_m=C_m,\n                     I_e=mu[i] * g_L[i],\n                     lambda_0=c[i],\n                     Delta_V=Delta_u[i],\n                     g_L=g_L[i],\n                     tau_sfa=tau_theta[i],\n                     q_sfa=J_theta[i] / tau_theta[i],\n                     V_T_star=V_th[i],\n                     V_reset=V_reset[i],\n                     t_ref=t_ref[i],\n                     tau_syn_ex=max([tau_ex, dt]),\n                     tau_syn_in=max([tau_in, dt]),\n                     E_L=0.,\n                     V_m=0.)\n\n# connect the populations\nfor i, nest_i in enumerate(nest_pops):\n    for j, nest_j in enumerate(nest_pops):\n        if np.allclose(pconn[i, j], 1.):\n            conn_spec = {'rule': 'all_to_all'}\n        else:\n            conn_spec = {\n                'rule': 'fixed_indegree', 'indegree': int(pconn[i, j] * N[j])}\n\n        nest.Connect(nest_j, nest_i,\n                     conn_spec,\n                     syn_spec={'weight': J_syn[i, j] * g_syn[i, j],\n                               'delay': delay[i, j]})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We want to record all spikes of each population in order to compute the\nmesoscopic population activities `A_N(t)` from the microscopic simulation.\nWe also record the membrane potentials of five example neurons:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# monitor the output using a multimeter and a spike recorder\nnest_sr = []\nfor i, nest_i in enumerate(nest_pops):\n    nest_sr.append(nest.Create('spike_recorder'))\n    nest_sr[i].time_in_steps = True\n\n    # record all spikes from population to compute population activity\n    nest.Connect(nest_i, nest_sr[i], syn_spec={'weight': 1., 'delay': dt})\n\nNrecord = [5, 0]    # for each population \"i\" the first Nrecord[i] neurons are recorded\nnest_mm_Vm = []\nfor i, nest_i in enumerate(nest_pops):\n    nest_mm_Vm.append(nest.Create('multimeter'))\n    nest_mm_Vm[i].set(record_from=['V_m'], interval=dt_rec)\n    if Nrecord[i] != 0:\n        nest.Connect(nest_mm_Vm[i], nest_i[:Nrecord[i]], syn_spec={'weight': 1., 'delay': dt})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As before, all neurons in a given population will be stimulated with a\nstep input current. The following code block is identical to the one for\nthe mesoscopic simulation above:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# create the step current devices if they do not exist already\nnest_stepcurrent = nest.Create('step_current_generator', M)\n# set the parameters for the step currents\nfor i in range(M):\n    nest_stepcurrent[i].set(amplitude_times=tstep[i] + t0,\n                            amplitude_values=step[i] * g_L[i],\n                            origin=t0,\n                            stop=t_end)\n    nest_stepcurrent[i].set(amplitude_times=tstep[i] + t0,\n                            amplitude_values=step[i] * g_L[i],\n                            origin=t0,\n                            stop=t_end)\n    # optionally a stopping time may be added by: 'stop': sim_T + t0\n    pop_ = nest_pops[i]\n    nest.Connect(nest_stepcurrent[i], pop_, syn_spec={'weight': 1., 'delay': dt})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now start the microscopic simulation:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.rng_seed = 1\n\nt = np.arange(0., t_end, dt_rec)\nA_N = np.ones((t.size, M)) * np.nan\n\n# simulate 1 step longer to make sure all t are simulated\nnest.Simulate(t_end + dt)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's retrieve the data of the spike recorder and plot the activity of the\nexcitatory population (in Hz):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for i in range(len(nest_pops)):\n    data_sr = nest_sr[i].get('events', 'times') * dt - t0\n    bins = np.concatenate((t, np.array([t[-1] + dt_rec])))\n    A = np.histogram(data_sr, bins=bins)[0] / float(N[i]) / dt_rec\n    A_N[:, i] = A * 1000  # in Hz\n\nt = np.arange(dt, t_end + dt, dt_rec)\nplt.figure(2)\nplt.plot(t, A_N[:, 0])\nplt.xlabel('time [ms]')\nplt.ylabel('population activity [Hz]')\nplt.title('Population activities (microscopic sim.)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This should look similar to the population activity obtained from the\nmesoscopic simulation based on the NEST model ``gif_pop_psc_exp`` (cf. figure\n1). Now we retrieve the data of the multimeter, which allows us to look at\nthe membrane potentials of single neurons. Here we plot the voltage traces\n(in mV) of five example neurons:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "voltage = []\nfor i in range(M):\n    if Nrecord[i] > 0:\n        senders = nest_mm_Vm[i].get('events', 'senders')\n        v = nest_mm_Vm[i].get('events', 'V_m')\n        voltage.append(\n            np.array([v[np.where(senders == j)] for j in set(senders)]))\n    else:\n        voltage.append(np.array([]))\n\nf, axarr = plt.subplots(Nrecord[0], sharex=True)\nfor i in range(Nrecord[0]):\n    axarr[i].plot(voltage[0][i])\n    axarr[i].set_yticks((0, 15, 30))\naxarr[i].set_xlabel('time [ms]')\naxarr[2].set_ylabel('membrane potential [mV]')\naxarr[0].set_title('5 example GIF neurons (microscopic sim.)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note that this plots only the subthreshold membrane potentials but not the\nspikes (as with every leaky integrate-and-fire model).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}