In [1]:
%pylab inline --no-import-all
import numpy as np
# Gaussian function
def sum_of_gaussians(a, positions, c):
f = lambda x: sum(a * np.exp( -(x-b)**2 / (2 * c**2)) for b in positions)
return f
def poisson_firing(ratefn, time):
"""
ratefn describes the prescribed rate for each point in time t.
(Try to use a lambda function for the ratefn input)
Time can be in seconds or milliseconds
Assumption:
Assuming only one spike can happen in a second (reasonable assumption)
Effective refractive period for continuous firing is 1 ms
(i.e max firing = 1000 Hz)
Good enough assumption as long as a rate function with an insane firing
rate is used
"""
dt = time[1]-time[0]
spiketimes = []
for t in time:
P = np.exp(-1 * dt * ratefn(t))
rand = np.random.random()
if np.random.random() > P:
spiketimes.append(t)
return spiketimes
In [2]:
# 1-D line
vel = 1 # 1 unit per sec
t = np.arange(0, 30, 0.001)
x = lambda t: vel * t
# Preference distribution of neurons (tuning curve) - max firing - 50 Hz
neuron1_preference = sum_of_gaussians(50, [10], 5)
neuron2_preference = sum_of_gaussians(50, [15], 5)
neuron3_preference = sum_of_gaussians(50, [20], 5)
# Plot preference
pylab.plot(x(t), neuron1_preference(x(t)))
pylab.plot(x(t), neuron2_preference(x(t)))
pylab.plot(x(t), neuron3_preference(x(t)))
Out[2]:
In [5]:
# Obtain poisson based firing rates of each (place cell) neuron
trials = 10
neuron1_firing = []
neuron2_firing = []
neuron3_firing = []
for _ in range(trials):
neuron1_firing.append(poisson_firing(lambda t: neuron1_preference(vel * t), t))
neuron2_firing.append(poisson_firing(lambda t: neuron2_preference(vel * t), t))
neuron3_firing.append(poisson_firing(lambda t: neuron3_preference(vel * t), t))
In [6]:
# Raster plot each of the three neurons
for i in range(trials):
pylab.plot(neuron1_firing[i], (5 * i + 1) * np.ones(len(neuron1_firing[i])), '|b')
pylab.plot(neuron2_firing[i], (5 * i + 2) * np.ones(len(neuron2_firing[i])), '|g')
pylab.plot(neuron3_firing[i], (5 * i + 3) * np.ones(len(neuron3_firing[i])), '|r')
# beautify
pylab.axes().set_xticks(np.arange(0,30+1,5))
pylab.axes().set_yticks(np.arange(0,trials*2,3))
Out[6]: