.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/CampbellSiegert.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_CampbellSiegert.py: Campbell & Siegert approximation example ---------------------------------------- .. only:: html ---- Run this example as a Jupyter notebook: .. card:: :width: 25% :margin: 2 :text-align: center :link: https://lab.ebrains.eu/hub/user-redirect/git-pull?repo=https%3A%2F%2Fgithub.com%2Fnest%2Fnest-simulator-examples&urlpath=lab%2Ftree%2Fnest-simulator-examples%2Fnotebooks%2Fnotebooks%2FCampbellSiegert.ipynb&branch=main :link-alt: JupyterHub service .. image:: https://nest-simulator.org/TryItOnEBRAINS.png .. grid:: 1 1 1 1 :padding: 0 0 2 0 .. grid-item:: :class: sd-text-muted :margin: 0 0 3 0 :padding: 0 0 3 0 :columns: 4 See :ref:`our guide ` for more information and troubleshooting. ---- Example script that applies Campbell's theorem and Siegert's rate approximation to and integrate-and-fire neuron. This script calculates the firing rate of an integrate-and-fire neuron in response to a series of Poisson generators, each specified with a rate and a synaptic weight. The calculated rate is compared with a simulation using the ``iaf_psc_alpha`` model References ~~~~~~~~~~ .. [1] Papoulis A (1991). Probability, Random Variables, and Stochastic Processes, McGraw-Hill .. [2] Siegert AJ (1951). On the first passage time probability problem, Phys Rev 81: 617-623 Authors ~~~~~~~ S. Schrader, Siegert implementation by T. Tetzlaff .. GENERATED FROM PYTHON SOURCE LINES 51-53 First, we import all necessary modules for simulation and analysis. Scipy should be imported before nest. .. GENERATED FROM PYTHON SOURCE LINES 53-59 .. code-block:: Python import nest import numpy as np from scipy.optimize import fmin from scipy.special import erf .. GENERATED FROM PYTHON SOURCE LINES 60-64 We first set the parameters of neurons, noise and the simulation. First settings are with a single Poisson source, second is with two Poisson sources with half the rate of the single source. Both should lead to the same results. .. GENERATED FROM PYTHON SOURCE LINES 64-83 .. code-block:: Python weights = [0.1] # (mV) psp amplitudes rates = [10000.0] # (1/s) rate of Poisson sources # weights = [0.1, 0.1] # (mV) psp amplitudes # rates = [5000., 5000.] # (1/s) rate of Poisson sources C_m = 250.0 # (pF) capacitance E_L = -70.0 # (mV) resting potential I_e = 0.0 # (nA) external current V_reset = -70.0 # (mV) reset potential V_th = -55.0 # (mV) firing threshold t_ref = 2.0 # (ms) refractory period tau_m = 10.0 # (ms) membrane time constant tau_syn_ex = 0.5 # (ms) excitatory synaptic time constant tau_syn_in = 2.0 # (ms) inhibitory synaptic time constant simtime = 20000 # (ms) duration of simulation n_neurons = 10 # number of simulated neurons .. GENERATED FROM PYTHON SOURCE LINES 84-85 For convenience we define some units. .. GENERATED FROM PYTHON SOURCE LINES 85-97 .. code-block:: Python pF = 1e-12 ms = 1e-3 pA = 1e-12 mV = 1e-3 mu = 0.0 sigma2 = 0.0 J = [] assert len(weights) == len(rates) .. GENERATED FROM PYTHON SOURCE LINES 98-100 In the following we analytically compute the firing rate of the neuron based on Campbell's theorem [1]_ and Siegerts approximation [2]_. .. GENERATED FROM PYTHON SOURCE LINES 100-146 .. code-block:: Python for rate, weight in zip(rates, weights): if weight > 0: tau_syn = tau_syn_ex else: tau_syn = tau_syn_in # We define the form of a single PSP, which allows us to match the # maximal value to or chosen weight. def psp(x): return -( (C_m * pF) / (tau_syn * ms) * (1 / (C_m * pF)) * (np.exp(1) / (tau_syn * ms)) * ( ((-x * np.exp(-x / (tau_syn * ms))) / (1 / (tau_syn * ms) - 1 / (tau_m * ms))) + (np.exp(-x / (tau_m * ms)) - np.exp(-x / (tau_syn * ms))) / ((1 / (tau_syn * ms) - 1 / (tau_m * ms)) ** 2) ) ) min_result = fmin(psp, [0], full_output=1, disp=0) # We need to calculate the PSC amplitude (i.e., the weight we set in NEST) # from the PSP amplitude, that we have specified above. fudge = -1.0 / min_result[1] J.append(C_m * weight / (tau_syn) * fudge) # We now use Campbell's theorem to calculate mean and variance of # the input due to the Poisson sources. The mean and variance add up # for each Poisson source. mu += rate * (J[-1] * pA) * (tau_syn * ms) * np.exp(1) * (tau_m * ms) / (C_m * pF) sigma2 += ( rate * (2 * tau_m * ms + tau_syn * ms) * (J[-1] * pA * tau_syn * ms * np.exp(1) * tau_m * ms / (2 * (C_m * pF) * (tau_m * ms + tau_syn * ms))) ** 2 ) mu += E_L * mV sigma = np.sqrt(sigma2) .. GENERATED FROM PYTHON SOURCE LINES 147-149 Having calculate mean and variance of the input, we can now employ Siegert's rate approximation. .. GENERATED FROM PYTHON SOURCE LINES 149-161 .. code-block:: Python num_iterations = 100 upper = (V_th * mV - mu) / (sigma * np.sqrt(2)) lower = (E_L * mV - mu) / (sigma * np.sqrt(2)) interval = (upper - lower) / num_iterations tmpsum = 0.0 for cu in range(0, num_iterations + 1): u = lower + cu * interval f = np.exp(u**2) * (1 + erf(u)) tmpsum += interval * np.sqrt(np.pi) * f r = 1.0 / (t_ref * ms + tau_m * ms * tmpsum) .. GENERATED FROM PYTHON SOURCE LINES 162-164 We now simulate neurons receiving Poisson spike trains as input, and compare the theoretical result to the empirical value. .. GENERATED FROM PYTHON SOURCE LINES 164-180 .. code-block:: Python nest.ResetKernel() nest.set_verbosity("M_WARNING") neurondict = { "V_th": V_th, "tau_m": tau_m, "tau_syn_ex": tau_syn_ex, "tau_syn_in": tau_syn_in, "C_m": C_m, "E_L": E_L, "t_ref": t_ref, "V_m": E_L, "V_reset": E_L, } .. GENERATED FROM PYTHON SOURCE LINES 181-184 Neurons and devices are instantiated. We set a high threshold as we want free membrane potential. In addition we choose a small resolution for recording the membrane to collect good statistics. .. GENERATED FROM PYTHON SOURCE LINES 184-191 .. code-block:: Python n = nest.Create("iaf_psc_alpha", n_neurons, params=neurondict) n_free = nest.Create("iaf_psc_alpha", params=dict(neurondict, V_th=1e12)) pg = nest.Create("poisson_generator", len(rates), {"rate": rates}) vm = nest.Create("voltmeter", params={"interval": 0.1}) sr = nest.Create("spike_recorder") .. GENERATED FROM PYTHON SOURCE LINES 192-193 We connect devices and neurons and start the simulation. .. GENERATED FROM PYTHON SOURCE LINES 193-202 .. code-block:: Python pg_n_synspec = {"weight": np.tile(J, ((n_neurons), 1)), "delay": 0.1} nest.Connect(pg, n, syn_spec=pg_n_synspec) nest.Connect(pg, n_free, syn_spec={"weight": [J]}) nest.Connect(vm, n_free) nest.Connect(n, sr) nest.Simulate(simtime) .. GENERATED FROM PYTHON SOURCE LINES 203-206 Here we read out the recorded membrane potential. The first 500 steps are omitted so initial transients do not perturb our results. We then print the results from theory and simulation. .. GENERATED FROM PYTHON SOURCE LINES 206-212 .. code-block:: Python v_free = vm.events["V_m"] Nskip = 500 print(f"mean membrane potential (actual / calculated): {np.mean(v_free[Nskip:])} / {mu * 1000}") print(f"variance (actual / calculated): {np.var(v_free[Nskip:])} / {sigma2 * 1e6}") print(f"firing rate (actual / calculated): {sr.n_events / (n_neurons * simtime * ms)} / {r}") .. _sphx_glr_download_auto_examples_CampbellSiegert.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: CampbellSiegert.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: CampbellSiegert.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_