Project2_Handout_2022
Problem 1:¶
Recall the following circuit from HW4:
$h^1$ and $h^2$ refer to feedforward filters; $h^{11}$ and $h^{22}$ refer to feedback filters; $H$ refer to summation-cross-feedback filters.
Take $BSG_1, BSG_2$ to be IAF neurons with $b=1,\delta=0.01$ and $b=-1,\delta=-0.01$ respectively. Both neurons have $\kappa=1$.
Questions:¶
Encode and decode a randomly generated input stimulus using the given circuit. Plot decoding results (time plots and corresponding SNR plots).
(15 points) Encode a randomly generated signal $u(t)$ of bandwidth $50$ Hertz using the given neuron model. Plot the time-course of the membrane voltages of the two IAF neurons, and label the corresponding spikes on the same plot.
(25 points) Use the detected spikes to decode the input signal $u(t)$. Plot the time-courses of the decoded and original input signal on the same plot to illustrate the quality of recovery. Also plot the time-course of the corresponding SNR between the signals.
# import libraries
import numpy as np
import matplotlib.pyplot as plt
# utility functions
def stimulus(t, omega, s):
out = np.zeros_like(t)
omega_pi = omega / np.pi
for k in range(len(s)):
out += s[k] * omega_pi * np.sinc(omega_pi * t – k)
return out / np.max(np.abs(out))
def signal_generator(t, samples, sample_times, omega):
signal = np.zeros_like(t)
for s, st in zip(samples, sample_times):
signal += omega / np.pi * s * np.sinc(omega / np.pi * (t – st))
return signal
dt = 1e-5 # time step
ds = 10 # downsampling factor
t = np.arange(-0.1, 0.25, dt) # simulation time
omega = 2 * np.pi * 50 # cutoff frequency
N = np.floor(t[-1] / np.pi * omega).astype(int) # number of sample points
u = stimulus(t, omega, np.random.rand(N) – 0.5) # randomly generated input stimulus
# filter impulse responses
t_filt = np.arange(T_1, T_2, dt)
h1 = ( # h_1
* np.exp(-a * t_filt)
(a * t_filt) ** 3 / np.math.factorial(3)
– (a * t_filt) ** 5 / np.math.factorial(5)
t_filt = np.arange(T_1, T_2, dt)
h2 = ( # h_1
* np.exp(-a * t_filt)
(a * t_filt) ** 3 / np.math.factorial(3)
– (a * t_filt) ** 5 / np.math.factorial(5)
def own_filter(t, tk): # h^{11} = h^{22}
td = t – tk
td = td[td>0]
return np.sum(0.1 * np.exp(-td/0.01))
def cross_filter(t, tk): # H
td = t – tk
td = td[td>0]
return np.sum(0.075 * np.exp(-td/0.015))
# define model
# instantiate model
# plot results
Problem 2:¶
Consider the following circuit diagram for a neuron with two inputs $u_1, u_2$:
The linear filters $h_1, h_2$ are models for two dendrites (dendritic processing filters) incident on the neuron, while the axon-hillock (which governs spiking behavior) of the neuron is modeled by an Integrate-and-Fire (IAF) unit.
The aggregate dendritic input to the IAF is thus given by $v(t) = \sum_{m=1}^2 (u_m\ast h_m)(t)$. The IAF unit maps the input signal $v$ into the output spike sequence $(t_k)_{k=1}^n$, where $n$ denotes the total number of spikes produced on an interval $t\in[0,\,T]$.
You are each given a black-box instance of the neuron described above (with name neuron_uni.py) that takes a time vector $t$ and an input signal $u$ as input arguments.
In “Dendrite” mode, the model first feeds $u$ to the dendritic tree, then passes the output of the dendritic tree to the axon hillock, and finally returns the spike train $\sum_k \delta(t-t_k)$ generated by the axon-hillock. In this case $u$ should be a 2D array with $u[0]=u_1, u[1]=u_2$.
In “Axon” mode, $u$ is passed directly to the axon-hillock. In this case $u$ should be a 1D array of the same shape as $t$.
The model given to you in model_uni.py can be used by doing the following:
>>> from model_uni import model
>>> tk1 = model(t, np.array([u1, u2]), “Dendrite”) # inject into dendritic tree
>>> tk2 = model(t, v, “Axon”) # inject into axon-hillock
Questions:¶
Your task is to identify the circuit.
(20 points) We will start with identifying the IAF unit. Assume that $\kappa=1$ is already known. Inject the axon-hillock with appropriately chosen currents to identify the firing threshold $\delta$ and bias current $b$ of the IAF unit. Briefly explain your methodology.
(40 points) Identify the two filters $(h^1(t),h^2(t))$ on $t\in [0,0.5]$[s] by generating random non-zero bandlimited input signals ${\bf u} = (u_1, u_2)$ in the trignometric polynomial space. Choose any other parameters for your experiments as appropriate. For a range of bandwidths $\Omega^i, i=1,2,\ldots$, identify the filters $\hat{h}_{1,i}, \hat{h}_{2,i}$ by using the input with corresponding bandwidth. Calculate the mean square error $\epsilon_i = MSE(\hat{h}_{1,i-1}, \hat{h}_{1,i}) + MSE(\hat{h}_{2,i-1}, \hat{h}_{2,i})$, where $\hat{h}_{1,0}=\hat{h}_{2,0}=0$.
Plot $\epsilon_i$ v.s. $\Omega_i$. What do you find? Can you identify the effective bandwidth of the unknown dendritic processing filters?
import sys
sys.path.append(‘student’)
def getStimulusTrig(t, Omega, L, is_zero, M):
if is_zero: # returns a zero signal which can be useful for debugging
ul = np.zeros(2 * L + 1)
ul = np.random.rand(2 * L + 1)
ul[:L] = np.conj(ul[::-1][:L])
u = np.zeros_like(t) # initialize the signal
for l in range(-L, L + 1):
u = u + ul[l + L] * np.exp(1j * l * Omega * t / L)
u = np.real(u)
# normalize
max_u = np.max(np.abs(u))
if max_u > 0:
u = M * u / max_u
ul = M * ul / max_u
return u, ul
t = np.arange(0, 0.5, dt) # sec
Omega = 2 * np.pi * L / T
# generate random signals in trigonometric polynomial space
u, u_coeffs = getStimulusTrig(t, Omega, L, False, 1)
# compute q
# compute Phi
# plot results