In [ ]:
try:
import pypet
except ImportError:
# You can defintely ignore this line, this is solely for me, Robert,
# because I did not install pypet and want to use my own developer version ;-)
import sys
sys.path.append('/media/data/PYTHON_WORKSPACE/pypet-project')
# For plotting into the notebook
%matplotlib inline
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
We will numerically integrate the linear differential equation:
$$\frac{dV}{dt} = -\frac{1}{\tau_V} V + I$$with a non-linear reset rule $V \leftarrow 0$ if $V \geq 1$ and additional refractory period of length $\tau_{ref}$.
Our neuron simulator is based on the simple Euler integration method.
V_init - initial value of VI - applied currenttau_V - membrane time constant $\tau_V$tau_ref - the refractory period $\tau_{ref}$dt - Simulation timestepduration - Simulation durationOur function returns 3 data items:
V_array - the development of V over timespiketimes - list of times the neuron spikedtimes - numpy array of simulation times
In [ ]:
def euler_neuron(V_init, I, tau_V, tau_ref, dt, duration):
""" Simulation of a leaky integrate and fire neuron.
Simulates the equation
dV/dT = -1/tau_V * V + I
with a simple Euler scheme.
"""
# Createa a times array
times = np.arange(0, duration, dt)
steps = len(times)
# Create the V array for the membrane development over time
V_array = np.zeros(steps)
V_array[0] = V_init
spiketimes = [] # List to collect all times of action potentials
# Do the Euler integration:
for step in range(1, steps):
if V_array[step-1] >= 1:
# The membrane potential crossed the threshold and we mark this as
# an action potential
V_array[step] = 0
spiketimes.append(times[step-1])
elif spiketimes and times[step] - spiketimes[-1] <= tau_ref:
# We are in the refractory period, so we simply clamp the voltage
# to 0
V_array[step] = 0
else:
# Euler Integration step:
dV = -1.0 / tau_V * V_array[step-1] + I
V_array[step] = V_array[step-1] + dV * dt
return V_array, spiketimes, times
Running and plotting the results gives:
In [ ]:
V_init = 0.0 # unitless
I = 0.1 # unitless
tau_V = 12.0 # ms
tau_ref = 5.0 # ms
dt = 0.1 # ms
duration = 250.0 # ms
V_array, spiketimes, times = euler_neuron(V_init, I, tau_V, tau_ref, dt, duration)
plt.plot(times, V_array)
plt.xlabel('t[ms]')
plt.ylabel('V')
Let's use pypet to manage our simulation. We will use a trajectory container to store our data.
Let's assume we already have such a container, and we just want to put the parameters of our simulation into it.
In [ ]:
def add_parameters(traj):
"""Adds all parameters to `traj`"""
traj.f_add_parameter('neuron.V_init', 0.0,
comment='The initial condition for the membrane potential')
traj.parameters = Parameter('neuron.I', 0.0,
comment='The externally applied current.')
traj.par.neuron.tau_V = 12.0, 'The membrane time constant in milliseconds'
traj.par.neuron.tau_ref = 5.0, 'The refractory period in milliseconds '
traj.f_add_parameter_group('simulation', comment='Group containing dt and the duration')
traj.simulation.duration = 1000.0, 'The duration of the experiment in milliseconds.'
traj.simulation.dt = 0.1, 'The step size of an Euler integration step.'
Now we don't simply want to run this parameter set, but explore the parameter space.
We will try different values of I and tau_ref.
We will use a full grid search, using the cartesian_prodcut function of pypet.
In [ ]:
def add_exploration(traj):
"""Explores different values of `I` and `tau_ref`."""
explore_dict = {'neuron.I': np.arange(0, 1.01, 0.05).tolist(),
'neuron.tau_ref': [5.0, 7.5, 10.0]}
explore_dict = cartesian_product(explore_dict, ('neuron.tau_ref', 'neuron.I'))
# The second argument, the tuple, specifies the order of the cartesian product,
# The variable on the right most side changes fastest and defines the
# 'inner for-loop' of the cartesian product
traj.f_explore(explore_dict)
We still need a wrapper function that passes our data from the trajectory to the simulator and back.
In [ ]:
def pypet_neuron(traj):
""" Wraps the `euler_neuron` function and returns an estimate of the firing rate in Hz."""
parameter_dict = traj.parameters.f_to_dict(short_names=True, fast_access=True)
V_array, spiketimes, times = euler_neuron(**parameter_dict)
# Calculate the spikes
nspikes = len(spiketimes)
traj.f_add_result('neuron.$', V=V_array, times=times, nspikes=nspikes,
comment='Contains the development of the membrane potential over time, '
'the corresponding timesteps, '
'as well as the number of spikes.')
# Return an estimate of the firing rate.
# * 1000 because we assume units of milliseconds
return nspikes/traj.duration * 1000.0
Finally we need to start our simulation, for this we also need an environment.
In [ ]:
from pypet import Parameter, cartesian_product, Environment
filename = os.path.join('hdf5', 'FiringRate.hdf5')
env = Environment(trajectory='FiringRate',
comment='Experiment to measure the firing rate '
'of a leaky integrate and fire neuron. '
'Exploring different input currents, '
'as well as refractory periods',
add_time=False, # We don't want to add the current time to the name,
log_stdout=False, # Important since we operate in Ipython notebook
multiproc=True,
ncores=4, #My laptop has 4 cores ;-)
wrap_mode='LOCK',
filename=filename,
overwrite_file=True)
traj = env.v_trajectory
# Add parameters
add_parameters(traj)
# Let's explore
add_exploration(traj)
In [ ]:
# Run the experiment
env.f_run(pypet_neuron)
# Finally disable logging and close all log-files
env.f_disable_logging()
So we can inspect the HDF5 file and all our data is there.
Let's do a bit more, our pypet_neuron function currently returns an
estimate of the firing rate. But this is lost at the moment.
We can define custom post-processing, that will be automatically called after the execution of the simulation.
We will collect all results and put these into a table - aka pandas DataFrame - to have a nice overview of input parameters and the resulting firing rate.
In [ ]:
def neuron_postproc(traj, result_list):
"""Postprocessing, sorts computed firing rates into a table"""
# Let's create a pandas DataFrame to sort the computed firing rate according to the
# parameters. We could have also used a 2D numpy array.
# But a pandas DataFrame has the advantage that we can index directly with
# the parameter values without translating these into integer indices.
I_range = traj.par.neuron.f_get('I').f_get_range()
ref_range = traj.par.neuron.f_get('tau_ref').f_get_range()
I_index = sorted(set(I_range))
ref_index = sorted(set(ref_range))
rates_frame = pd.DataFrame(columns=ref_index, index=I_index)
# This frame is basically a two dimensional table
# that we can index with our parameters.
#
# The table looks a bit like this:
#
# | tau_ref 5.0 | 7.5 | 10.0 |
# I |---------------|-------------|-------------|
# 0.0 | f[Hz] | | |
# 0.05 | | | |
# ... | | | |
# Now iterate over the results. The `result_list` is a list of tuples, with the
# run index at first position and our result at the second.
for result_tuple in result_list:
run_idx = result_tuple[0]
firing_rates = result_tuple[1]
I_val = I_range[run_idx]
ref_val = ref_range[run_idx]
rates_frame.loc[I_val, ref_val] = firing_rates # Put the firing rate into the
# data frame
# Finally we store our new firing rate table into the trajectory
traj.f_add_result('summary.firing_rates', rates_frame=rates_frame,
comment='Contains a pandas data frame with all firing rates.')
Ok let's run everything again, but also add the post-processing to the environment.
In [ ]:
filename = os.path.join('hdf5', 'FiringRate.hdf5')
env = Environment(trajectory='FiringRate',
comment='Experiment to measure the firing rate '
'of a leaky integrate and fire neuron. '
'Exploring different input currents, '
'as well as refractory periods',
add_time=False, # We don't want to add the current time to the name,
log_stdout=False, # Important since we operate in Ipython notebook
multiproc=True,
ncores=4, #My laptop has 4 cores ;-)
wrap_mode='LOCK',
filename=filename,
overwrite_file=True)
traj = env.v_trajectory
# Add parameters
add_parameters(traj)
# Let's explore
add_exploration(traj)
# NEW LINE! ADDING POST-PROCESSING
env.f_add_postprocessing(neuron_postproc)
In [ ]:
# Run the experiment
env.f_run(pypet_neuron)
# Finally disable logging and close all log-files
env.f_disable_logging()
In [ ]:
from pypet import Trajectory
# This time we don't need an environment since we are just going to look
# at the data in the trajectory
traj = Trajectory('FiringRate', add_time=False)
# Let's load the trajectory from the file
# Only load the parameters, we will load the results on the fly as we need them
filename = os.path.join('hdf5', 'FiringRate.hdf5')
traj.f_load(load_parameters=2,
load_derived_parameters=0,
load_results=0,
load_other_data=0,
filename=filename)
# We'll simply use auto loading so all data will be loaded when needed.
traj.v_auto_load = True
# Here we load the data automatically on the fly
rates_frame = traj.res.summary.firing_rates.rates_frame
plt.figure(figsize=(8,6))
plt.subplot(2,1,1)
# Let's iterate through the columns and plot the different firing rates:
for tau_ref, I_col in rates_frame.iteritems():
plt.plot(I_col.index, I_col, label='Avg. Rate for tau_ref=%s' % str(tau_ref))
# Label the plot
plt.xlabel('I')
plt.ylabel('f[Hz]')
plt.title('Firing as a function of input current `I`')
plt.legend(loc='best')
# Also let's plot an example run, how about run 13 ?
example_run = 13
traj.v_idx = example_run # Set the trajectory to a particular run!
# Get the example data
example_I = traj.I
example_tau_ref = traj.tau_ref
example_V = traj.results.neuron.crun.V # Here crun stands for run_00000013
# Get the sampled timesteps
dt_array = traj.results.neuron.crun.times
# And plot the development of V over time,
# Since this is rather repetitive, we only
# plot the first eighth of it.
plt.subplot(2,1,2)
plt.plot(dt_array, example_V)
plt.xlim((0, dt_array[1000]))
# Label the axis
plt.xlabel('t[ms]')
plt.ylabel('V')
plt.title('Example of development of V for I=%s, tau_ref=%s in run %d' %
(str(example_I), str(example_tau_ref), traj.v_idx))
plt.tight_layout()
# Finally revoke the `traj.v_idx=13` statement and set everything back to normal.
# Since our analysis is done here, we could skip that, but it is always a good idea
# to do that.
traj.f_restore_default()
In [ ]: