Spiking models¶
SRM 0 Spiking Model¶
(1)\[u_i(t) = \eta(t - \hat{t}_i) + \sum_j w_{ij} \sum_{t_j^f} \epsilon_0 (t - t_j^f)\]
TODO: description of w, and \(\hat{t}_i\).
Eta function¶
The \(\eta(s)\) function describes the after-potential of a neuron after a spike.
We approximate it like this (x-axis: time, y-axis: potential):
(2)\[\eta(s) = \delta(s) - \eta_0 \exp (-\frac{s}{\tau_m})\]
Reset eta_0:
Membrane time constant
Epsilon function¶
\(\epsilon\) describes the time course of a postsynaptic potential caused by spike arrival.
We approximated this function like this:
(3)\[\epsilon (s) = \frac{1}{1 - \frac{\tau_c}{\tau_m}} (\exp(\frac{-s}{\tau_m}) - \exp(\frac{-s}{\tau_c}))\]
Where
- \(s\) is the time after the arriving spike
- \(\tau_c\) is the current time constant.
- \(\tau_m\) is the membrane time constant.
When plotted, it looks like this (x-axis: time, y-axis: potential):
Current time constant:
Membrane time constant
Vectorization¶
- Because we wanted to have a fast implementation of our model, we rewrote the original SRM equation (1) as
- a operations on matrices. Those operations (dot-product, row-sum, element-wise product) are fast to compute.
We calculate the membrane potential for every neuron at time t \(u(t)\):
(4)\[ u(t) = Z(t) + rowsum((W^T \cdot S ) \circ \mathcal{E}_t)\]
Where
- \(n\) is the number of neurons
- \(W^T\) the transposed Weight matrix
- \(S\) is the Spiketrain previous to time \(t\) in binary notation
- \(\mathcal{E}_t\) is a helper matrix \(\begin{pmatrix} \epsilon_1(t) & \epsilon_1(t-1) & \dots & \epsilon_1(0) \\ \colon & & & \colon \\ \epsilon_n(t) & \epsilon_n(t-1) & \dots & \epsilon_n(0) \end{pmatrix}\)
- \(\epsilon_i(s)\) is the epsilon function for neuron i at time s
- \(Z(t) = \begin{pmatrix}\eta(t - \hat{t}_1 & \dots & \eta(t - \hat{t}_n \end{pmatrix}^T\)
- \(\circ\) is element-wise product
- rowsum(M) is the sum over the rows in a matrix M
To understand this equation better, let’s have a look at it’s components:
- \(W^T \cdot S\) is a \(n \times t\) matrix, that says us how many weighted spikes arrive at a neuron at each time
- \(rowsum((W^T \cdot S) \circ \mathcal{E}_t)\) is a \(n\)-dimensional vector, that gives us the membrane potential of a neuron caused by incoming spikes
- \(Z(t) + \dots\) in the last step we add the after-potential of each neuron
You can find those dot-products implemented in spiking.py:
incoming_spikes = np.dot(weights.T, spiketrain_window)
incoming_potential = np.sum(incoming_spikes * epsilon_matrix, axis=1)
total_potential = self.eta(np.ones(neurons)*t - self.last_spike) + incoming_potential
# Any new spikes?
neurons_high_current = np.where(total_potential > self.threshold)
spiketrain[neurons_high_current, t] = True
# Update last_spike
spiking_neurons = np.where(spiketrain[:, t])
self.last_spike[spiking_neurons] = t
if self.verbose:
print("SRM Time step", t)
print("Incoming current", incoming_potential)
print("Total potential", total_potential)
print("Last spike", self.last_spike)
print("")
return total_potential
class SRM_X(SRM):
def __init__(self, neurons, threshold, t_current, t_membrane, eta_reset, ax_delay, simulation_window_size=100, verbose=False):
"""
Like the SRM model, but additionally it supports axonal delays.
:param neurons: Number of neurons
:param threshold: Spiking threshold
:param t_current: Current-time-constant (:math:`t_s`)
:type t_current: Float or Numpy Float Array
:param t_membrane: Membrane-time-constant (t_m)
:param eta_reset: Reset constant
:param ax_delay: Axonal delays
:param simulation_window_size: Only look at the n last spikes
:param verbose:
:return: ``None``
"""
# Check user input
# TODO
SRM.__init__(self, neurons, threshold, t_current, t_membrane, eta_reset, simulation_window_size=simulation_window_size,
verbose=verbose)
self.ax_delay = ax_delay
def eps(self, s):
r"""
Evaluate the Epsilon function with an axonal delay :math:`\tau_d`.
.. math:: \epsilon (s) = \frac{1}{1 - \frac{\tau_c}{\tau_m}} (\exp(\frac{-(s-\tau_d)}{\tau_m}) - \exp(\frac{-(s - \tau_d)}{\tau_c}))
:label: epsilon_axdelay
Returns a single Float Value if the time constants (current, membrane) are the same for each neuron.
Returns a Float Vector with eps(s) for each neuron, if the time constants are different for each neuron.
:param s: Time s
:return: Function eps(s) at time s
:rtype: Float or Vector of Floats
"""
eps = (1/(1-self.t_current/self.t_membrane))*(np.exp(-(s - self.ax_delay)/self.t_membrane)
- np.exp(- (s - self.ax_delay)/self.t_current))
eps[np.where(eps<0)] = 0
return eps
class Izhikevich:
"""
Izhikevich model
http://www.izhikevich.org/publications/spikes.htm
"""
# TODO this model looks wrong, because we receive far to less output spikes!
def __init__(self, neurons, a, b, c, d, v0=-65, threshold=30, verbose=False):
"""
:param neurons: total number of neurons
:param a: decay rate
:param b: sensitivity
:param c: reset
:param d: reset
:param v0: Initial voltage
:param threshold:
:param verbose: Verbose output to console. Default: False.
:type verbose: Boolean
:return:
"""
print("The Izhikevich model implementation hasn't been tested yet -- use with care")
self.a = a * np.ones((neurons, 1))
self.b = b * np.ones((neurons, 1))
self.c = c * np.ones((neurons, 1))
self.d = d * np.ones((neurons, 1))
self.threshold = threshold
self.verbose = verbose
# Initial u and v
self.v = v0 * np.ones((neurons, 1))
self.u = self.b*self.v
self.v_plot = np.empty((neurons, 0))
def check_spikes(self, spikes, weights, t):
# Input
# TODO: is this correct, or should I use another variable instead threshold?
I = self.threshold * spikes[:, t]
I = I[:, np.newaxis]
fired = np.nonzero(self.v >= self.threshold)[0] # Indices of neurons that spike
print(fired)
t = np.array([t+0*fired, fired])
spikes[fired, t] = True
self.v_plot = np.hstack((self.v_plot, self.v))
self.v[fired] = self.c[fired]
self.u[fired] = self.u[fired] + self.d[fired]
I = I + np.sum(weights[fired, :], axis=0, keepdims=True)
self.v = self.v + 0.5 * (0.04 * self.v ** 2 + 5 * self.v + 140 - self.u + I) # step 0.5 ms
self.v = self.v + 0.5 * (0.04 * self.v ** 2 + 5 * self.v + 140 - self.u + I) # for numerical
self.u = self.u + self.a * (np.multiply(self.b, self.v) - self.u) # stability
if self.verbose:
print("Izhikevich Time step", t)
print("Total current", self.v)
print("Fired: ", fired)
print("")
return self.v
if __name__ == "__main__":
srm_model = SRM(neurons=3, threshold=1, t_current=0.3, t_membrane=20, eta_reset=5, verbose=True)
izhikevich_model = Izhikevich(neurons=3, a=0.02, b=0.2, c=-65, d=8, threshold=30, verbose=True)
models = [srm_model, izhikevich_model]
for model in models:
print("-"*10)
if isinstance(model, SRM):
print('Demonstration of the SRM Model')
elif isinstance(model, Izhikevich):
print('Demonstration of the Izhikevich Model')
s = np.array([[0, 0, 1, 0, 0, 0, 1, 1, 0, 0],
[1, 0, 0, 0, 0, 0, 1, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
w = np.array([[0, 0, 1], [0, 0, 1], [0, 0, 0]])
neurons, timesteps = s.shape
for t in range(timesteps):
total_current = model.simulate(s, w, t)
print("Spiketrain:")
print(s)
Note
Of course, those matrices could grow very large with time \(t\). So in practice, we use an approximation on a smaller time-slice. We don’t look at the whole spiketrain, but only on a small fraction of it (maybe the last 500ms).
Other Spiking models¶
Other spiking models can be easily integrated later.
All that you need is to implement a