{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Campbell & Siegert approximation example\n\nExample script that applies Campbell's theorem and Siegert's rate\napproximation to and integrate-and-fire neuron.\n\nThis script calculates the firing rate of an integrate-and-fire neuron\nin response to a series of Poisson generators, each specified with a\nrate and a synaptic weight. The calculated rate is compared with a\nsimulation using the ``iaf_psc_alpha`` model\n\n\n\n## References:\n\n .. [1] Papoulis A (1991). Probability, Random Variables, and\n Stochastic Processes, McGraw-Hill\n .. [2] Siegert AJ (1951). On the first passage time probability problem,\n Phys Rev 81: 617-623\n\n## Authors\n\nS. Schrader, Siegert implentation by T. Tetzlaff\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we import all necessary modules for simulation and analysis. Scipy\nshould be imported before nest.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from scipy.special import erf\nfrom scipy.optimize import fmin\n\nimport numpy as np\n\nimport nest" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first set the parameters of neurons, noise and the simulation. First\nsettings are with a single Poisson source, second is with two Poisson\nsources with half the rate of the single source. Both should lead to the\nsame results.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "weights = [0.1] # (mV) psp amplitudes\nrates = [10000.] # (1/s) rate of Poisson sources\n# weights = [0.1, 0.1] # (mV) psp amplitudes\n# rates = [5000., 5000.] # (1/s) rate of Poisson sources\n\nC_m = 250.0 # (pF) capacitance\nE_L = -70.0 # (mV) resting potential\nI_e = 0.0 # (nA) external current\nV_reset = -70.0 # (mV) reset potential\nV_th = -55.0 # (mV) firing threshold\nt_ref = 2.0 # (ms) refractory period\ntau_m = 10.0 # (ms) membrane time constant\ntau_syn_ex = .5 # (ms) excitatory synaptic time constant\ntau_syn_in = 2.0 # (ms) inhibitory synaptic time constant\n\nsimtime = 20000 # (ms) duration of simulation\nn_neurons = 10 # number of simulated neurons" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For convenience we define some units.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pF = 1e-12\nms = 1e-3\npA = 1e-12\nmV = 1e-3\n\nmu = 0.0\nsigma2 = 0.0\nJ = []\n\nassert(len(weights) == len(rates))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following we analytically compute the firing rate of the neuron\nbased on Campbell's theorem [1]_ and Siegerts approximation [2]_.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for rate, weight in zip(rates, weights):\n\n if weight > 0:\n tau_syn = tau_syn_ex\n else:\n tau_syn = tau_syn_in\n\n t_psp = np.arange(0., 10. * (tau_m * ms + tau_syn * ms), 0.0001)\n\n # We define the form of a single PSP, which allows us to match the\n # maximal value to or chosen weight.\n\n def psp(x):\n return - ((C_m * pF) / (tau_syn * ms) * (1 / (C_m * pF)) *\n (np.exp(1) / (tau_syn * ms)) *\n (((-x * np.exp(-x / (tau_syn * ms))) /\n (1 / (tau_syn * ms) - 1 / (tau_m * ms))) +\n (np.exp(-x / (tau_m * ms)) - np.exp(-x / (tau_syn * ms))) /\n ((1 / (tau_syn * ms) - 1 / (tau_m * ms)) ** 2)))\n\n min_result = fmin(psp, [0], full_output=1, disp=0)\n\n # We need to calculate the PSC amplitude (i.e., the weight we set in NEST)\n # from the PSP amplitude, that we have specified above.\n\n fudge = -1. / min_result[1]\n J.append(C_m * weight / (tau_syn) * fudge)\n\n # We now use Campbell's theorem to calculate mean and variance of\n # the input due to the Poisson sources. The mean and variance add up\n # for each Poisson source.\n\n mu += (rate * (J[-1] * pA) * (tau_syn * ms) *\n np.exp(1) * (tau_m * ms) / (C_m * pF))\n\n sigma2 += rate * (2 * tau_m * ms + tau_syn * ms) * \\\n (J[-1] * pA * tau_syn * ms * np.exp(1) * tau_m * ms /\n (2 * (C_m * pF) * (tau_m * ms + tau_syn * ms))) ** 2\n\nmu += (E_L * mV)\nsigma = np.sqrt(sigma2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Having calculate mean and variance of the input, we can now employ\nSiegert's rate approximation.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "num_iterations = 100\nupper = (V_th * mV - mu) / (sigma * np.sqrt(2))\nlower = (E_L * mV - mu) / (sigma * np.sqrt(2))\ninterval = (upper - lower) / num_iterations\ntmpsum = 0.0\nfor cu in range(0, num_iterations + 1):\n u = lower + cu * interval\n f = np.exp(u ** 2) * (1 + erf(u))\n tmpsum += interval * np.sqrt(np.pi) * f\nr = 1. / (t_ref * ms + tau_m * ms * tmpsum)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now simulate neurons receiving Poisson spike trains as input,\nand compare the theoretical result to the empirical value.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nest.ResetKernel()\n\nnest.set_verbosity('M_WARNING')\nneurondict = {'V_th': V_th,\n 'tau_m': tau_m,\n 'tau_syn_ex': tau_syn_ex,\n 'tau_syn_in': tau_syn_in,\n 'C_m': C_m,\n 'E_L': E_L,\n 't_ref': t_ref,\n 'V_m': E_L,\n 'V_reset': E_L}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Neurons and devices are instantiated. We set a high threshold as we want\nfree membrane potential. In addition we choose a small resolution for\nrecording the membrane to collect good statistics.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nest.SetDefaults('iaf_psc_alpha', neurondict)\nn = nest.Create('iaf_psc_alpha', n_neurons)\nn_free = nest.Create('iaf_psc_alpha', 1, {'V_th': 1e12})\npg = nest.Create('poisson_generator', len(rates),\n [{'rate': float(rate_i)} for rate_i in rates])\nvm = nest.Create('voltmeter', 1, {'interval': .1})\nsr = nest.Create('spike_recorder')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We connect devices and neurons and start the simulation.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for indx in range(len(pg)):\n nest.Connect(pg[indx], n,\n syn_spec={'weight': float(J[indx]), 'delay': 0.1})\n nest.Connect(pg[indx], n_free, syn_spec={'weight': J[indx]})\n\nnest.Connect(vm, n_free)\nnest.Connect(n, sr)\n\nnest.Simulate(simtime)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we read out the recorded membrane potential. The first 500 steps are\nomitted so initial transients do not perturb our results. We then print the\nresults from theory and simulation.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "v_free = vm.get('events', 'V_m')[500:-1]\nprint('mean membrane potential (actual / calculated): {0} / {1}'\n .format(np.mean(v_free), mu * 1000))\nprint('variance (actual / calculated): {0} / {1}'\n .format(np.var(v_free), sigma2 * 1e6))\nprint('firing rate (actual / calculated): {0} / {1}'\n .format(nest.GetStatus(sr, 'n_events')[0] /\n (n_neurons * simtime * ms), r))" ] } ], "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 }