BMEBW4020

Problem 1: PRC of Reduced Hodgkin-Huxley Model¶
In this problem, you are asked to do the following:

Implement Winfree’s method for approximating the Phase Response Curve (PRC) for Voltage of the Reduced Hodgkin-Huxley (RHH) Model, which is called the HodgkinHuxley3State model in the compneuro package. Use simple words to describe the procedures and your method (via #comments).
Generate spike train with RHH neuron and the I/O equivalent project-integrate-and-fire (PIF) neuron when the input is $I_{ext}(t)$. Use simple words to describe your implementation of PIF (via #comments).
Evaluate the mean $\mu$ and standard deviation $\sigma$ of the difference between the corresponding spike times of the RHH neuron and the equivalent PIF neuron for bandlimited stimulus described by equation scaled to have maximum amplitude $c > 0$ (i.e. $\max_{t}|u(t)| = c$). Plot these statistics as a function of $c$ across values of $c$. Note 1: You should generate new inputs $u(t)$ (with different random coefficients) to ensure that the error statistics are not biased by the form of the input stimulus.
Note 2: The value of $c$ should not be so large that the limit cycle of the HH model collapses to singularity. It would be instructive to look at F-I curve of the RHH model first to find a range of input current values that gives periodic spiking (this was explored in the last homework).

Note on $I_{ext} = I + u(t)$ ¶
Starting this homework, the injected current to model neurons $I_{ext}$ will have this form $I+u(t)$, where $I$ refers to injected current and $u(t)$ is the stimulus.
The ideas is as follows:

$I$: injected (bias) current is a constant value for all time that is injected to the neuron model. You can consider this value as the average/DC value of the input $I_{ext}$. Since this value is not time-varying, it contains no information (in the sense of entropy). For this reason, we don’t use refer to $I$ as stimulus. However, the bias $I$ will change the shape of the limit cycle $\Gamma$ (or the periodic solution ${\bf x}^0$) of the neuron model, where the higher the bias current $I$, the smaller the limit cycle and the faster the neuron oscillates. For this reason, we can consider the dynamic of the neuron model to be parametrized/indexed by the bias $I$.

$u(t)$: this is a time-varying signal that has information content. If we consider $I$ as the DC value of $I_{ext}$, then $u(t)$ will be a zero-mean signal that corresponds to local perturbation of the neuron dynamic around the limit cycle $\Gamma$. This is the input signal that we seeks to decode from the neuron spike times.

Input Stimulus $u(t)$ – Complex Exponential ¶
The input stimulus here is a zero-DC valued complex exponential of order $M=5$ on support $t\in[0, 200]$ ms with bandwidth $\Omega = 2\pi \cdot 20 \quad [rad\cdot s^-1]$,
u(t)= \sum^{M}_{m=-M}a_{m} \exp\left(j\frac{m\Omega t}{M}\right)

Note that for stimulus $u(t)$ to be real-valued signal, the coefficients $a_m$ need to be conjugate-symmetric. In another word, $a_m = \overline{a_{-m}} = Re(a_{-m}) – j Im(a_{-m})$,
where $Re(\cdot), Im(\cdot)$ are the real and imagninary parts of the complex number respectively. The stimulus $u(t)$ should be zero-DC valued, which means that $a_0 = 0$. This input stimulus $u(t)$ is additively coupled to a periodically spiking HH neuron.

Problem 2: Synaptic Input and reduced PIF¶
So far, the injected currents into point neuron models are chosen to be arbitrary continuous waveforms. For this problem, you will explore using synaptic current to drive postsynaptic neurons.

You are asked to do the following:

Generate 20 Poisson processes, each with constant rate $\lambda=100 Hz$. These will serve as time traces of neurotransmitter concentration for 20 GABA_A synapses we will simulate in the next step.
Simulate a single GABA_A synapse connected to a HH neuron, with the input neurotransmitter concentration for the synapse being one of the 20 Poisson processes generated in Step 1, scaled by a constant $c$. Generate the following plots Synaptic current in the postsynaptic neuron $I_{syn}^{i} = g^{i}(V^{i} – E^i_{rev})$
Membrane Voltage of the postsynaptic neuron $V^{i}$
Do you see any spiking behavior in the postsynaptic neuron for any value of $c$? Why or why not?

Simulate 20 GABAA synapses connected to a HH neuron, with neurotransmitter concentrations given by the Poisson processes generated in Step 1, scaled up by $c=500$. Simulate this convergent synapse-neuron circuit where the synaptic current is the sum $I{syn} = \sum{i}g^{i}(V-E^i{rev}), i={1,2,\ldots,20}$. Plot the synaptic current and the neuron membrane voltage and check if the neuron spikes now.
Assuming that the model neuron is biased at a current level $I = mean(I_{syn})$ ($I_{syn}$ from step 3), implement reduced PIF neuron under this bias level $I$. Drive the reduced PIF neuron with input current $u(t) = I_{syn}(t)-I$. Plot the spike times of the model neuron and the reduced PIF. Is the reduced PIF a good approximation in this case? Why or why not?

Poisson Spiking $\lambda$¶
In time bin $[t,t+\Delta t]$, the probability of spike is $\lambda \cdot \Delta t$. Therefore, the spike state is a binary number indicating if a in this time bin is $s(t) = \mathbb{1}_{x<\lambda \cdot \Delta t}$, where $x \sim Unif([0,1])$. Here $\mathbb{1}$ denotes an indicator random variable, Synaptic Current $I_{syn}(t)= g_{syn}(t)\cdot(V_{post}(t)-E_{rev})$¶ Currents are injected into neurons through opening of ion-channels that may be controlled by membrane voltage of ligand bindings. In the case of chemical synapses (which is the case here), the ion-channels on post-synaptic neurons are open when neuro-transmitters relased by the pre-synaptic neurons bind to the receptors. There are many different kinds of such ligand-controlled ion-channels and they all have different dynamics, you are encouraged to look up concepts like Ionotropic vs Metabotropic. We can model all synaptic current going through ligand-gated ion channels as follows: $$I_{syn}(t)= g_{syn}(t)\cdot(V_{post}(t)-E_{rev})$$ $g_{syn}(t)$ is a time varying function that controls the synaptic conductance $V_{post}(t)$ is the membrane voltage of the postsynaptic neuron $E_{rev}$ is the reverse potential of the synapse Note: By convention, synaptic current $I_{syn}$ has negative value because it denotes inward current to the neuron. Therefore, to use simulation kernels (including hodgkin_huxley), the input current $I_{ext}$ will be $I_{ext} = -I_{syn}$. Note: Use the GABA_A synapse as defined in the compneuro # Import libraries %matplotlib inline import numpy as np import matplotlib.pyplot as plt from tqdm.notebook import tqdm from compneuro.neurons.hodgkin_huxley_3state import HodgkinHuxley3State # For Problem 1 from compneuro.neurons.hodgkin_huxley import HodgkinHuxley # For Problem 2 from compneuro.utils.signal import spike_detect np.random.seed(0) # fix random seed plt.rcParams["figure.figsize"] = [10, 5] plt.rcParams["figure.dpi"] = 80 Problem 1 - PRC of RHH model¶ Initialization¶ # TODO: Implement signal generator that takes in [time, coefficients, order, Bandwidth] as input def bandlmt_sig(t,a,M,omega): return np.nan # TODO: Specify the time resolution for computing PRC. # The smaller the value the better the Winfree method performs # NOTE: start with a smaller value to ensure the code works before making it larger dt = np.nan # TODO: Specify the bias of HH Model (I in I_ext = I + u(t)). # This value will be used in both the neuron model itself and its PIF. bias = np.nan # unit is in pico amp # we suggest refering the injection current amplitude (aka without stimulus) # that you've experimented with in HW1 that gives robust firing Item (i): Extract the PRC¶ # TODO: Implement Winfree's Method # As a suggestion, the function should return # 1. period: the period of oscillation, maybe in number of time-steps # 2. limitCycle: the output of the state variables [V, n, m, h, a, b] on a limit cycle # 3. PRC: the phase response curve of the Voltage along the limit cycle # rhh = HodgkinHuxley() # def winfree(model, I, dt, ...): # return ... # TODO: Simulate your model period, limitCycle, PRC = winfree(np.nan) # TODO: Visualize the limit cycle and the PRC along it Item (ii): Compare Spikes of PIF and RHH¶ We first generate a bandlimited signal with random coefficients # temporal support of the signal [0, 200] ms t = np.arange(0, 200, dt) # you can use [0,0.2]s as well if it suits you # bandwidth: 20Hz Omega = 2*pi*20 # order: 5 # Amplitude of the signal, you can use other values. Sig_Amp = 1 # TODO: generate coefficients am = np.nan # generate input signal. Make sure that the signal's DC value (mean) is 0 u1 = bandlmt_sig(t, am, M, Omega) u1 = Sig_Amp*u1/max(u1) Simulate the RHH neuron with bias and bandlimited input added. # TODO: simulate model and find spikes # rhh = HodgkinHuxley3State() # ... = rhh.solve(np.nan) # TODO: find spike time indices of RHH # tk_idx_RHH = np.nan Output of PIF¶ The t-transform of the PIF neuron can be written as, \int_{t_k}^{t_{k+1}}(1+\psi(s+\tau(s))u(s))ds = \delta_{t_{k+1}}-\delta_{t_{k}} \approx T, \\ \tau(t) = \psi(t+\tau(t))u(t),~\tau(0)=0. We assume that $\tau(s)=0$ , and hence the $t$-transform is reduced to \int_{t_k}^{t_{k+1}}(1+\psi(s)u(s))ds = T. # TODO: Implement the PIF model using the input defined above and the PRC curve computed # def pif(np.nan): # return V # TODO: execute the pif and find spike times # t_pif = np.nan # TODO: Plot RHH, PIF outputs and compare the Inter-spike intervals (time_diff) Item (iii): Record PIF Error Statistics¶ # TODO: Initalize experiment # Find a range of amplitude $C$ that is within the range of permissable input current to the RHH model. # Too large and the limit cycle could collapse C = np.nan # amplitude of u(t) <-- Edit this into an array of c vals avg_diff = np.zeros(len(C)) var_diff = np.zeros(len(C)) # TODO: Simulate the RHH and PIF for inputs I_ext = I + u(t) for different amplitude of u(t): max|u(t)| = c # TODO: Compute the difference between interspike interval of the RHH and PIF, # and calculate the mean/standard deviation of the error. # TODO: visualize the error across stimulus amplitude. Problem 2: Synaptic Input and reduced PIF¶ Step 1. Generate Spike Trains¶ # TODO: generate 20 spike trains and plot them spike_trains = [] Step 2. Simulate Synapse connected to HH neuron¶ from compneuro.synapses.gaba_a import GABA_A spike_state = spike_trains[0] c = np.nan NT = c * spike_state # TODO: link the synapse to HH neuron, simulate over t using Euler's method # syn = GABA_A() # syn_res = syn.solve(...) #use "Euler" solver # hh = HodgkinHuxley() # TODO: Plot (pls also include the chosed spike train) TODO: Does it spike? Why or why not? (Answer in markdown) Step 3. Run Convergent Circuit¶ All synapses provided input to the same output neuron # TODO: base your code on the code from Step 2. # TODO: Plot Step 4. Reduced PIF¶ # TODO: Use synaptic current from the previous step to find the bias current and signal. # TODO: Find the period and PRC of the HH model at the given bias current # TODO: Simulate PIF using the PRC found above and the signal deteremined from the synaptic current in Step 3 # (use the PIF function you defined in Q1) # TODO: Plot TODO: Is PIF a good approximation? Why and why not? (Answer in markdown) HINT: think about the difference in the assumption about the input between the standard PIF and the synaptic input used here.