2020-09-10 • STA for ≠ PSP shapes
Contents
2020-09-10 • STA for ≠ PSP shapes¶
This notebook builds on /notebooks/2020-07-29__STA.
There, we cut out spike-triggered windows of a noisy VI signal, and averaged them, to obtain one STA (spike triggered average) window per incoming neuron. The STA windows from connected neurons showed
a nice bump (at ~10 ms post spike) which made them well differentiable from those of unconnected neurons (which had flat STA windows, as expected).
The question is if this bump is due to averaging post synaptic potentials (PSPs), or postsynaptic spikes.
In this notebook we investigate that question by changing the PSP shape (by changing the synaptic integration time constant \(\tau_{syn}\)), and then checking whether the STA bump tracks this change in shape, or not.
Imports & time grid¶
from voltage_to_wiring_sim.notebook_init import *
Preloading:
- numpy … (0.33 s)
- matplotlib.pyplot … (0.47 s)
- numba … (0.80 s)
- sympy … (0.97 s)
- unyt … (2.38 s)
Imported `np`, `mpl`, `plt
Imported package `voltage_to_wiring_sim` as `v`
Imported `*` from `v.util` and from `v.units`
Setup autoreload
tg = v.TimeGrid(T=10 * minute, dt=0.1 * ms)
TimeGrid(T=Quantity(10 min, stored in s, float64), dt=Quantity(0.1 ms, stored in s, float64), N=unyt_quantity(6000000, '(dimensionless)'), t=Array([0 1.667E-06 3.333E-06 ... 10 10 10] min, stored in s, float64))
Spike trains¶
‘Network’ definition.
N_in = 30
p_connected = 0.5
N_connected = round(N_in * p_connected)
N_unconnected = N_in - N_connected
15
Have all incoming neurons spike with the same mean frequency, for now.
f_spike = 20 * Hz
Quantity(20 Hz, stored in Hz, float64)
gen_st = v.generate_Poisson_spike_train
fix_rng_seed()
%%time
spike_trains_connected = [gen_st(tg, f_spike) for _ in range(N_connected)]
spike_trains_unconnected = [gen_st(tg, f_spike) for _ in range(N_unconnected)]
Wall time: 3 s
all_spike_trains = spike_trains_connected + spike_trains_unconnected
Inspect a time excerpt..
time_slice = 1 * minute + [0, 1] * second
slice_indices = np.round(time_slice / tg.dt).astype(int)
i_slice = slice(*slice_indices)
t_slice = tg.t[i_slice].in_units(second)
..of one presynaptic neuron:
v.spike_train.plot(t_slice, all_spike_trains[0][i_slice])
![../_images/2020-09-10__STA_for_different_PSP_shapes_16_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_16_0.png)
All connected presynaptic neurons:
all_incoming_spikes = sum(spike_trains_connected)
v.spike_train.plot(t_slice, all_incoming_spikes[i_slice])
![../_images/2020-09-10__STA_for_different_PSP_shapes_18_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_18_0.png)
Synaptic conductance¶
Δg_syn = 0.16 * nS
τ_syn = 35 * ms
Quantity(35 ms, stored in s, float64)
g_syn = v.calc_synaptic_conductance(tg, all_incoming_spikes, Δg_syn, τ_syn)
plt.plot(t_slice, g_syn[i_slice])
![../_images/2020-09-10__STA_for_different_PSP_shapes_21_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_21_0.png)
Membrane voltage¶
params = v.params.cortical_RS
print(params)
IzhikevichParams
----------------
C = 100 pF
k = 0.7 nS/mV
v_r = -60 mV
v_t = -40 mV
v_peak = 35 mV
a = 0.03 1/ms
b = -2 nS
c = -50 mV
d = 100 pA
v_syn = 0 mV
%%time
sim = v.simulate_izh_neuron(tg, params, g_syn)
Wall time: 588 ms
plt.plot(t_slice, sim.V_m[i_slice])
![../_images/2020-09-10__STA_for_different_PSP_shapes_25_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_25_0.png)
Noise¶
SNR = 10
spike_height = params.v_peak - params.v_r
σ_noise = spike_height / SNR
noise = np.random.randn(tg.N) * σ_noise
Vm_noisy = (sim.V_m + noise).in_units(mV)
Vm_noisy.name = "VI signal"
plt.plot(t_slice, Vm_noisy[i_slice])
![../_images/2020-09-10__STA_for_different_PSP_shapes_27_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_27_0.png)
Spike-triggered windows¶
def get_spike_indices(spike_train):
# `nonzero` returns a tuple (one element for each array dimension).
(spike_indices,) = np.nonzero(spike_train)
return spike_indices
Test with one spike train
spike_train__example = all_spike_trains[0]
spike_indices__example = get_spike_indices(spike_train__example)
num_spikes__example = len(spike_indices__example)
num_spikes__example
12020
Spiking frequency for this presynaptic neuron:
print(Quantity(num_spikes__example / tg.T, Hz))
20.03 Hz
Extract windows from the VI signal.
window_length = 100 * ms
window_tg = v.TimeGrid(window_length, tg.dt)
window_tg.t.name = "Time after spike"
# Performance
Vm_noisy_arr = np.asarray(Vm_noisy)
window_tg_N = window_tg.N
tg_N = tg.N
import numba
@numba.njit
def _make_windows(spike_indices):
windows = []
for ix in spike_indices:
ix_end = ix + window_tg_N
if ix_end < tg_N:
windows.append(Vm_noisy_arr[ix:ix_end])
return windows
def make_windows(spike_indices):
windows = np.stack(_make_windows(spike_indices))
return Array(windows, V, name=Vm_noisy.name).in_units(mV)
windows__example = make_windows(spike_indices__example)
windows__example.shape
(12019, 1000)
An example spike-triggered window:
plt.plot(window_tg.t, windows__example[0, :])
[<matplotlib.lines.Line2D at 0x1e187eb0ca0>]
![../_images/2020-09-10__STA_for_different_PSP_shapes_38_1.png](../_images/2020-09-10__STA_for_different_PSP_shapes_38_1.png)
And some more:
fig, axes = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True)
for row, row_axes in enumerate(axes):
for col, ax in enumerate(row_axes):
i = 3 * row + col
ax.plot(window_tg.t, windows__example[i, :])
if not (row == 2 and col == 0):
ax.set(xlabel=None, ylabel=None)
![../_images/2020-09-10__STA_for_different_PSP_shapes_40_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_40_0.png)
Average windows¶
def STA(spike_train):
spike_indices = get_spike_indices(spike_train)
windows = make_windows(spike_indices)
STA = windows.mean(axis=0)
return Array(STA, V, name="Spike-triggered <VI signal>").in_units(mV)
def plot_STA(spike_train, ax=None, **kwargs):
if ax is None:
fig, ax = plt.subplots()
ax.plot(window_tg.t, STA(spike_train), **kwargs)
plot_STA(spike_trains_connected[0])
![../_images/2020-09-10__STA_for_different_PSP_shapes_42_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_42_0.png)
plot_STA(spike_trains_unconnected[0])
![../_images/2020-09-10__STA_for_different_PSP_shapes_43_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_43_0.png)
Plot average windows of all spike trains¶
%%time
fig, ax = plt.subplots()
for spike_train in spike_trains_connected:
plot_STA(spike_train, ax, alpha=0.2, color='C0')
for spike_train in spike_trains_unconnected:
plot_STA(spike_train, ax, alpha=0.2, color='C1')
plt.close()
Wall time: 20.9 s
Extreme slowness is due to unyt.
Hence I oughta use my own new yunit package, or use no unit package (plain NumPy arrays).
# Clear existing texts, for iterative positioning.
for _ in range(len(ax.texts)):
ax.texts.pop()
ax.annotate(
"Spikes from\nconnected neurons",
xy=(26.55 * ms, 0.8),
xycoords=("data", "axes fraction"),
color="C0",
ha="left",
)
ax.annotate(
"Spikes from\nunconnected neurons",
xy=(70 * ms, 0.5),
xycoords=("data", "axes fraction"),
color="C1",
ha="center",
)
fig
![../_images/2020-09-10__STA_for_different_PSP_shapes_47_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_47_0.png)
Reminder of previous STA (with τ_syn = 7 ms):
...
![../_images/2020-09-10__STA_for_different_PSP_shapes_49_0.png](../_images/2020-09-10__STA_for_different_PSP_shapes_49_0.png)
Next steps:
New:
bump could still be due to spikes. Maybe artifically clip spikes to lessen their influence on STA?
Existing:
calculate p(connected) for every presynaptic neuron
describe influence of
N_in
,p_connected
,SNR
, ..VI realism