3A

2.Introduction to NEST Simulator

Unit Neural Coding and Brain Computing Unit (Tomoki Fukai)
Teacher Gastón Sivori
Lecture Material https://github.com/oist-ncbc/skill-pill-plus/tree/master/Tutorial-2_NEST-Introduction
Library https://github.com/iwasakishuto/University/tree/gh-pages/3A/theme/program/OIST/tutorial2.py
In [1]:
import sys
sys.path.append("path/to/library") # Like "University/3A/theme/program"

Agenda

Installation

Memo(What I did)

pip
$ pip install python-nest

※ This `nest` has lack of some module (ex. `nest.Create`)

cmake
# Install dependencies via Homebrew
$ brew install gcc cmake gsl open-mpi libtool
$ mkdir NEST
$ cd NEST
# Clone source files from Github.
$ wget --content-disposition https://github.com/nest/nest-simulat
or/archive/v2.16.0.tar.gz
$ tar zxf nest-simulator-2.16.0.tar.gz
$ mkir bld
$ cd bld
# Configure and build NEST inside the build directory 
# (replacing gcc-9 and g++-9 with the GCC compiler versions you installed with brew).
# You can check it with `g++ -v` or `gcc -v` command.
$ cmake -DCMAKE_INSTALL_PREFIX:PATH=. \
        -DCMAKE_C_COMPILER=gcc-9 \
        -DCMAKE_CXX_COMPILER=g++-9 \
        ../nest-simulator-2.16.0/
$ make -j4
$ make install
$ make installcheck

※ This `nest` also has lack of some module (ex. `nest.Create`)

Anacond
# If you want to use jupyter notebook, please include "jupyter"
$ conda create --name skill-pill-plus -c conda-forge nest-simulator python jupyter
# path/to/activate is like "/Users/iwasakishuto/.pyenv/versions/anaconda3-5.0.1/bin/activate"
$ source [path/to/activate] skill-pill-plus
# Activated!!
(skill-pill-plus) $

If you want to use this environment on the jupyter notebook, please follow this.

$ pip install environment_kernels
$ jupyter notebook --generate-config
$ vi ~/.jupyter/jupyter_notebook_config.py
# Add these line
>>> c.NotebookApp.kernel_spec_manager_class = 'environment_kernels.EnvironmentKernelSpecManager'
>>> c.EnvironmentKernelSpecManager.conda_env_dirs=['/Users/iwasakishuto/.pyenv/versions/anaconda3-5.0.1/env']
$ jupyter notebook
# If successful, you can chose your ideal environment at [Kernel → Change kernel]
# NOTE: You should install `jupyter` in the environment.
#       If you not, please run `conda install --name "Your Environment Name" jupyter`

※ This `nest` has all modules :)

In [2]:
import nest
import numpy as np
import matplotlib.pyplot as plt
# If you need to import other packages like scipy or scikit-learn
# You should import these "before" nest!!
For NEST to simulate different experiments, its kernel (core functionality) needs to be reseted (nest.ResetKernel()) for every new simulation. Please remember to run this code before running each simulation.

About NEST

  • Conceptually, neural networks consists of neurons and connections.
  • Nodes can be neurons, sub-networks or devices inside the NEST framework.
  • We will use devices to stimulate neurons and to measure their membrane potential.
  • Sub-networks, on the other hand, are arrangements of neurons whose parameters may be modified as a group. We will look at sub-networks in the following lectures.

leaky integrate-and-fire (LIF)

For now, we will focus on constructing a leaky integrate-and-fire (LIF) neuron with delta-shaped synaptic currents:

  • In this figure, a delta-shaped pulse $\delta(t-t_j^f)$ from neuron $j$ is being transmitted along its axon until it reaches neuron $i$'s dendrites through a synapse.
  • Synapses are modelled as low-pass RC circuit filters which output a current $\alpha(t-t_j^f)$.
  • These presynaptic currents reach the soma as an input current $I(t)$ which charges the capacitor $C$ (integration) while some of it leaks out through the resistance $R$. ($I(t) = I_R + I_C$)
  • The electrical components here represent a circuit modelization of the biological mechanisms happening at the synapse and the soma.
  • In this model, the membrane's potential is represented by capacitor $C$ such that when it has enough charge and reaches the threshold potential $\vartheta$, a spike $\delta(t-t_i^f)$ is generated and transmitted.
In [3]:
from OIST.tutorial2 import LIFmodel
In [4]:
model = LIFmodel(model_name="iaf_psc_alpha")

Create LIF neuron with alpha-function shaped synaptic currents.
Created a spike_generator.
Created a voltmeter.
Created a spike detector.
In [5]:
# Check the neuron params.
model.neuron_params()

C_m    : 250.0 pF
E_L    : -70.0 mV
I_e    : 0.0 pA
V_m    : -70.0 mV
V_reset: -70.0 mV
V_th   : -55.0 mV
t_ref  : 2.0 ms
tau_m  : 10.0 ms

Q1. You might have noticed that there is no parameter for the membrane resistance `R_m` as shown in the previous figure. Do you know why?

The membrane time constant tau_m corresponds to the membrane capacitance C_m and membrane resistance R_m as observed in the electrical circuit above. Therefore it's possible to calculate what is the R_m value given the following formula:

$$ \tau_m = C_m R_m $$

</span>

In [6]:
# Retrieve a particular set of parameters
model.neuron_params("V_reset", "V_th")
Out[6]:
((-70.0, -55.0),)
In [7]:
# Modify spike generation to values 10 and 50 ms
print(f"Before: {model.spike_params('spike_times')}")
model.set_params(node="spike_generator", spike_times=[10.,50.])
print(f"After : {model.spike_params('spike_times')}")

Before: ((array([], dtype=float64),),)
After : ((array([10., 50.]),),)
In [8]:
# Connect the spike generator and the voltmeter to the neuron.
model.build(weight=1e3)

Connect spike generator with a given synaptic specification ({'weight': 1000.0})
Connected voltmeter to the neuron for measurements.
Connected the neuron to spike detector for measurements.

Simulate for 100 ms and observe the results

In [9]:
ax = model.simulate(ms=100)
plt.show()

Q2. What can we observe here? Are these neuron spikes?

  • We can observe the membrane trace of the LIF neuron model excited with alpha-function shaped spikes at times 10 ms and 50 ms.
  • The observed peaks in the membrane trace are not spikes since the membrane threshold V_th is -55mV.
  • If we modify the external current or lower the membrane's threshold potential we should be able to observe spikes.
    • </span>

In [10]:
# Modifying the threshold
model = LIFmodel(model_name="iaf_psc_alpha", display=False)
model.set_params("neuron", V_th=-58.)
model.set_params("spike_generator", spike_times=[10.,50.])
model.set_params("voltmeter", interval=0.1)
model.build(weight=1e3, display=False)
ax = model.simulate(ms=100)
plt.show()
In [11]:
# modifying external current to LIF neuron
model = LIFmodel(model_name="iaf_psc_alpha", display=False)
model.set_params("neuron", I_e=300.)
model.set_params("spike_generator", spike_times=[10.,50.])
model.build(weight=1e3, display=False)
model.simulate(ms=100)
plt.show()

In this neuron model, the membrane potential does not follow the biologically observed dynamics that neurons have. Their use is limited to understand and experiment with network dynamics and spiking activity :(

HH-neuron model

For comparison, let's observe the membrane potential of a HH-neuron model (Hodgkin and Huxley, 1952).

This model was the "conclusion of a series of papers" concerned with the flow of electric current through the surface membrane of the giant nerve fibre in the squid.

Using voltage-clamp experiments and by varying extracellular sodium and potassium concentrations, Hodgkin and Huxley developed a model in which the properties of an excitable cell are described by a set of ordinary differential equations.

For simplicity, we will not focus on these equations. NEST's framework integrates these equations to obtain the instantaneous membrane potential. The following figure represents an electric circuit model of the neuron's membrane.

  • The HH-model (Hodgkin and Huxley, 1952) considers the ionic channels present on the cellular membrane.
  • These channels allows the conduction of specific types of ions.
  • Sodium channels (Na+) are scarce inside the cell, so that when sodium channel opens, positive charges rush into the cell to cause excitation.
  • Potassium ions (K+) are rich inside the cell, so that when potassium channel opens, positive charges rush out of the cell to cause inhibition.
  • The model assumes a leak current composed of all other remaining ionic currents (Cl, Mg, etc.)
  • This model is based on careful data analysis and follows the following equation: $$C_m \frac{dV}{dt} = \bar{g_{Na}}m^3h(E_{Na}-V) + \bar{g_K}n^4(E_K-V) + g_{Leak}(E_{Leak}-V) + I$$

Where the conductances $g_{Na}$ and $g_{K}$ are:

$$ g_{Na} = \bar{g_{Na}}m^3h $$$$ g_{K} = \bar{g_K}n^4 $$

Let's see these parameters in NEST by creating a hh_psc_alpha neuron. A HH-model neuron with delta-shaped synaptic currents.

In [12]:
from OIST.tutorial2 import HHmodel
In [13]:
model = HHmodel()

Create Hodgkin-Huxley neuron with alpha-function shaped synaptic currents.
Created a spike_generator.
Created a voltmeter.
Created a spike detector.
In [14]:
model.neuron_params()

C_m  : 100.0 pF
E_K  : -77.0 mV
E_L  : -54.402 mV
E_Na : 50.0 mV
g_K  : 3600.0 nS
g_L  : 30.0 nS
g_Na : 12000.0 nS
I_e  : 0.0 pA
V_m  : -65.0 mV
t_ref: 2.0 ms
In [15]:
model.set_params("voltmeter", interval=0.1)
model.set_params("spike_generator", spike_times=[10.,50.])
model.build(weight=1e3)

Connect spike generator with a given synaptic specification ({'weight': 1000.0})
Connected voltmeter to the neuron for measurements.
Connected the neuron to spike detector for measurements.
In [16]:
ax = model.simulate(ms=100)
plt.show()

Q3. Modify the previous example to observe spikes.

In [17]:
# modifying the synaptic strength (w) to observe spikes
model = HHmodel(display=False)
model.set_params("voltmeter", interval=0.1)
model.set_params("spike_generator", spike_times=[10.,50.])
model.build(weight=1e4, display=False)
ax = model.simulate(ms=100)
plt.show()
In [18]:
print(f"Number of spikes: {model.spikes}")

Number of spikes: 2

Could you observe the spikes? Just after a spike, the membrane potential resets to -77 mV and a short refractory period of 2 ms (check these parameter value t_ref from the parameter list shown above!).

During this period, the neuron is unable to react to excitation.

Q4. Can you try exciting the neuron during this period to observe this phenomena?

In [19]:
# In the following example the spike generator generates spikes at times 10ms, 11ms, 50ms and 52ms.
model = HHmodel(display=False)
model.set_params("voltmeter", interval=0.1)
model.set_params("spike_generator", spike_times=[10.,11.0,50.,52.0]) # Add!!
model.build(weight=1e4, display=False)
ax = model.simulate(ms=100)
plt.show()

The figure is the same as before due to the refractory period of the neuron.

Connecting nodes with different inputs

  • To compare neurons with different inputs, let's see how the parameters of the input affect the neurons' output.
  • In the following example, we will use
    • constant input current and belows with the same mean strength as the constant input current.
    • a Poisson spike train
  • The probability density function and firing rates of neurons will be shaped given by these inputs.

First, let us define the parameters of the simulation and input currents.

In [20]:
T = 1.2e3            # simulation time
nu_ext = 15e3        # firing rate of external Poisson source
J = 0.08             # synaptic weight
d = 0.1              # delay
C = 250.0            # membrane potential capacitance
mu = J*1e-3*nu_ext*C # mean input in pAa
I_ext = mu           # external current
In [21]:
neuron_params = {
    'C_m': C,      # (pF)
    'E_L': 0.,     # (mV)
    'I_e': 0.0,    # (pA)
    'V_m': 0.,     # (mV)
    'V_reset': 0., # (mV)
    'V_th': 15.,   # (mV)
    't_ref': 2.0,  # (ms)
    'tau_m': 10.0, # (ms)
}
In [22]:
plot_infos = {
    "color": ["blue", "red"],
    "label": ["constant input", "Poisson input"]
}

Then, make two types of models and compare them.

In [23]:
from OIST.tutorial2 import mkMultiDetectors, plotMultResults
In [24]:
const_model = LIFmodel(model_name="iaf_psc_delta", spike_generator=None, reset=True)
const_model.set_params("neuron", **neuron_params)
const_model.set_params("neuron", I_e=I_ext)

Create LIF neuron with delta-shaped synaptic currents.
Created a voltmeter.
Created a spike detector.
In [25]:
poisson_model = LIFmodel(model_name="iaf_psc_delta", spike_generator="poisson_generator", reset=False)
poisson_model.set_params("neuron", **neuron_params)
poisson_model.set_params("spike_generator", rate=nu_ext)
poisson_model.build(weight=J, delay=d)

Create LIF neuron with delta-shaped synaptic currents.
Created a poisson_generator.
Created a voltmeter.
Created a spike detector.
Connect spike generator with a given synaptic specification ({'weight': 0.08, 'delay': 0.1})
Connected voltmeter to the neuron for measurements.
Connected the neuron to spike detector for measurements.
In [26]:
neurons = const_model.neuron + poisson_model.neuron
multimeters, spikedetectors = mkMultiDetectors(2)
nest.Connect(multimeters, neurons,    'one_to_one')
nest.Connect(neurons, spikedetectors, 'one_to_one')

Created 2 multimeter to record membrane potential of the neurons.
Created 2 spike detectors to record spikes times and neuron identifiers, but not record from file.
In [27]:
nest.Simulate(T)
axL, axR = plotMultResults(multimeters, spikedetectors, rate=True, **plot_infos)
axR.legend(), axR.set_xlim(0,T)
plt.show()

What!? There are no spikes?! Well, look at the figures and consider the threshold potential value V_th. Let's try increasing the synaptic strengh J. The neuron receiving Poisson input will be in different firing rate regimes!

In [28]:
print(f"poisson_model.neuron_params('V_th') = {poisson_model.neuron_params('V_th')[0][0]}")
print(f"const_model.neuron_params('V_th') = {const_model.neuron_params('V_th')[0][0]}")

poisson_model.neuron_params('V_th') = 15.0
const_model.neuron_params('V_th') = 15.0

Q5. Try to make the neuron fire at an irregular regime.

In [29]:
T = 1.2e3            # 
nu_ext = 15e3        #
J = 0.09             # From 0.08 to 0.09.
d = 0.1              #
C = 250.0            #
mu = J*1e-3*nu_ext*C #
I_ext = mu           #
In [30]:
const_model = LIFmodel(model_name="iaf_psc_delta", spike_generator=None, reset=True, display=False)
const_model.set_params("neuron", **neuron_params)
const_model.set_params("neuron", I_e=I_ext)

poisson_model = LIFmodel(model_name="iaf_psc_delta", spike_generator="poisson_generator", reset=False, display=False)
poisson_model.set_params("neuron", **neuron_params)
poisson_model.set_params("spike_generator", rate=nu_ext)
poisson_model.build(weight=J, delay=d, display=False)

neurons = const_model.neuron + poisson_model.neuron
multimeters, spikedetectors = mkMultiDetectors(2, display=False)

nest.Connect(multimeters, neurons,    'one_to_one')
nest.Connect(neurons, spikedetectors, 'one_to_one')

nest.Simulate(T)
axL, axR = plotMultResults(multimeters, spikedetectors, rate=True, **plot_infos)
axR.legend(), axR.set_xlim(0,T)
plt.show()

Q6. Try to make the neuron spike at a regular regime.

In [31]:
T = 1.2e3            # 
nu_ext = 15e3        #
J = 0.11             # From 0.08 to 0.11.
d = 0.1              #
C = 250.0            #
mu = J*1e-3*nu_ext*C #
I_ext = mu           #
In [32]:
const_model = LIFmodel(model_name="iaf_psc_delta", spike_generator=None, reset=True, display=False)
const_model.set_params("neuron", **neuron_params)
const_model.set_params("neuron", I_e=I_ext)

poisson_model = LIFmodel(model_name="iaf_psc_delta", spike_generator="poisson_generator", reset=False, display=False)
poisson_model.set_params("neuron", **neuron_params)
poisson_model.set_params("spike_generator", rate=nu_ext)
poisson_model.build(weight=J, delay=d, display=False)

neurons = const_model.neuron + poisson_model.neuron
multimeters, spikedetectors = mkMultiDetectors(2, display=False)

nest.Connect(multimeters, neurons,    'one_to_one')
nest.Connect(neurons, spikedetectors, 'one_to_one')

nest.Simulate(T)
axL, axR = plotMultResults(multimeters, spikedetectors, rate=True, **plot_infos)
axL.legend(), axR.legend(), axR.set_xlim(0,T)
plt.show()

Since we are using the LIF neuron model with a membrane threshold potential V_th at 15mV, you will not be able to see proper spikes. Instead, you will only see the after-spike dynamics (membrane potential reset and refractory period).

Q7. Try to make HH neuron models spike at a regular regime for both input conditions.

In [33]:
T = 500.0            # simulation time
nu_ext = 15e3        # firing rate of external Poisson source
J = 0.08             # synaptic weight
d = 0.1              # delay
C = 100.0            # membrane potential capacitance
mu = J*1e-3*nu_ext*C # mean input in pAa
I_ext = mu           # external current
In [34]:
neuron_params = {
    'g_K': 3000.0, 
    'g_Na': 16000.0, 
    'g_L' : 60.0
}
In [35]:
const_model = HHmodel(model_name="hh_psc_alpha", spike_generator=None, reset=True, display=False)
const_model.set_params("neuron", **neuron_params)
const_model.set_params("neuron", I_e=I_ext)

poisson_model = HHmodel(model_name="hh_psc_alpha", spike_generator="poisson_generator", reset=False, display=False)
poisson_model.set_params("neuron", **neuron_params)
poisson_model.set_params("spike_generator", rate=nu_ext)
poisson_model.build(weight=J, delay=d, display=False)

neurons = const_model.neuron + poisson_model.neuron
multimeters, spikedetectors = mkMultiDetectors(2, display=False)

nest.Connect(multimeters, neurons,    'one_to_one')
nest.Connect(neurons, spikedetectors, 'one_to_one')

nest.Simulate(T)
axL, axR = plotMultResults(multimeters, spikedetectors, **plot_infos)
axL.legend(), axR.legend(), axR.set_xlim(0,T), axR.set_ylim(-80, 40)
plt.show()

Rate of neuron stimulated with constant input input: 56.112
Rate of neuron stimulated with Poisson input input: 48.096

F-I curve

  • As you have seen in the past figures, the firing rate of the neuron can change given different external current I_e.
  • It's possible to plot this relationship and its usually referred to as the F-I curve.
In [36]:
from OIST.tutorial2 import FiringRate
In [37]:
neuron_params = {
    'C_m': 250.0,
    'E_L': 0.,
    'V_m': 0.,
    'V_reset': 0.,
    'V_th': 15.,
    't_ref': 2.0,
    'tau_m': 10.0,
}
In [38]:
I_es = (200.0,600.0,5.0)
rates = [FiringRate("iaf_psc_delta", I_e, **neuron_params) for I_e in np.arange(*I_es)]

plt.plot(np.arange(*I_es), rates)
plt.xlabel(r'$I_e$ (pA)'), plt.ylabel('Firing rate (Hz)'), plt.grid(linestyle='-', linewidth=.25, alpha=.7)
plt.show()
  • As you can observe from this figure, given the parameter setup, the neuron in a leaky-integrate-and-fire model starts to fire spikes at around 370 pA of external input current I_e.
  • For this particular case, V_th was chosen to be low to observe the trend of the firing rate as external current is increasing.
  • Modifying the value of the capacitance C_m moves this curve in the x-axis
  • Modifying the refractory period t_ref slows down the firing rate i.e. more current is required for a given firing rate.

Q8. Observe these changes in the F-I curve by modifying the neuron model parameters.

In [39]:
neuron_params = {
    'C_m': 150.0,
    'E_L': 0.,
    'V_m': 0.,
    'V_reset': 0.,
    'V_th': 15.,
    't_ref': 2.0,
    'tau_m': 10.0,
}
In [40]:
I_es = (200.0,600.0,5.0)
rates = [FiringRate("iaf_psc_delta", I_e, **neuron_params) for I_e in np.arange(*I_es)]

plt.plot(np.arange(*I_es),rates)
plt.xlabel(r'$I_e$ (pA)'), plt.ylabel('Firing rate (Hz)'), plt.grid(linestyle='-', linewidth=.25, alpha=.7)
plt.show()
In [41]:
neuron_params = {
    'C_m': 350.0,
    'E_L': 0.,
    'V_m': 0.,
    'V_reset': 0.,
    'V_th': 15.,
    't_ref': 2.0,
    'tau_m': 10.0,
}
In [42]:
I_es = (200.0,600.0,5.0)
rates = [FiringRate("iaf_psc_delta", I_e, **neuron_params) for I_e in np.arange(*I_es)]

plt.plot(np.arange(*I_es),rates)
plt.xlabel(r'$I_e$ (pA)'), plt.ylabel('Firing rate (Hz)'), plt.grid(linestyle='-', linewidth=.25, alpha=.7)
plt.show()
In [43]:
neuron_params = {
    'C_m': 250.0,
    'E_L': 0.,
    'V_m': 0.,
    'V_reset': 0.,
    'V_th': 15.,
    't_ref': 1.0,
    'tau_m': 10.0,
}
In [44]:
I_es = (200.0,600.0,5.0)
rates = [FiringRate("iaf_psc_delta", I_e, **neuron_params) for I_e in np.arange(*I_es)]

plt.plot(np.arange(*I_es),rates)
plt.xlabel(r'$I_e$ (pA)'), plt.ylabel('Firing rate (Hz)'), plt.grid(linestyle='-', linewidth=.25, alpha=.7)
plt.show()

Q10. Repeat these simulations with the HH neuron model hh_psc_alpha by modifying channel conductances.

In [45]:
neuron_params = {
    'g_K': 3600.0, 
    'g_Na': 12000.0, 
    'g_L' : 30.0
}
In [46]:
I_es = (300.0,1000.0,5.0)
rates = [FiringRate("hh_psc_alpha", I_e, **neuron_params) for I_e in np.arange(*I_es)]

plt.plot(np.arange(*I_es),rates)
plt.grid(linestyle='-', linewidth=.25, alpha=.7), plt.xlabel(r'$I_e$ (pA)'), plt.ylabel('Firing rate (Hz)')
plt.show()
In [47]:
neuron_params = {
    'g_K': 3600.0, 
    'g_Na': 13000.0, 
    'g_L' : 30.0
}
In [48]:
I_es = (300.0,1000.0,5.0)
rates = [FiringRate("hh_psc_alpha", I_e, **neuron_params) for I_e in np.arange(*I_es)]

plt.plot(np.arange(*I_es),rates)
plt.grid(linestyle='-', linewidth=.25, alpha=.7), plt.xlabel(r'$I_e$ (pA)'), plt.ylabel('Firing rate (Hz)')
plt.show()
In [49]:
neuron_params = {
    'g_K': 3800.0, 
    'g_Na': 12000.0, 
    'g_L' : 30.0
}
In [50]:
I_es = (300.0,1000.0,5.0)
rates = [FiringRate("hh_psc_alpha", I_e, **neuron_params) for I_e in np.arange(*I_es)]

plt.plot(np.arange(*I_es),rates)
plt.grid(linestyle='-', linewidth=.25, alpha=.7), plt.xlabel(r'$I_e$ (pA)'), plt.ylabel('Firing rate (Hz)')
plt.show()
In [ ]:
hidden
Table of Contents
Published
Last Updated
Category
OIST
Tags
Contact
Other contents
  • Home
  • Blog
  • Front-End
  • Kerasy
  • Python-Charmers
  • Translation-Gummy