Contents

2020-09-18 • Clip spikes before STA

2020-09-18 • Clip spikes before STA

This notebook is similar to /notebooks/2020-07-29__STA and /notebooks/2020-09-10__STA_for_different_PSP_shapes.\

We still wonder whether the bump in the STA is due to averaging post synaptic potentials (PSPs), or postsynaptic spikes.

Here, we try to remove the influence of spikes in the STA by clipping them before averaging windows.

Reproducibility

%reload_ext watermark
%watermark -ntz -v -hm -rgb -iv

Fri Sep 18 2020 16:08:54 GMT Daylight Time 

CPython 3.8.3
IPython 7.13.0

compiler   : MSC v.1916 64 bit (AMD64)
system     : Windows
release    : 10
machine    : AMD64
processor  : Intel64 Family 6 Model 142 Stepping 9, GenuineIntel
CPU cores  : 4
interpreter: 64bit
host name  : LABTOP
Git hash   : 5c2958efe77ec99f89399e38845283c44798eadc
Git repo   : git@github.com:tfiers/voltage-to-wiring-sim.git
Git branch : main

Nothing new here

All sections nested under this are the same as in the two previous notebooks.

Imports & time grid

from voltage_to_wiring_sim.notebook_init import *

Preloading:
 - numpy … (0.20 s)
 - matplotlib.pyplot … (0.35 s)
 - numba … (0.68 s)
 - sympy … (0.89 s)
 - unyt … (3.99 s)

Imported `np`, `mpl`, `plt
Imported package `voltage_to_wiring_sim` as `v`
Imported `*` from `v.util` and from `v.units`
Setup autoreload

Reproducibility info. This notebook was last run on/with:

Fri Sep 18 2020 16:09:01 GMT Daylight Time 

CPython 3.8.3
IPython 7.13.0

compiler   : MSC v.1916 64 bit (AMD64)
system     : Windows
release    : 10
machine    : AMD64
processor  : Intel64 Family 6 Model 142 Stepping 9, GenuineIntel
CPU cores  : 4
interpreter: 64bit
host name  : LABTOP
Git hash   : 5c2958efe77ec99f89399e38845283c44798eadc
Git repo   : git@github.com:tfiers/voltage-to-wiring-sim.git
Git branch : main

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: 5.26 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-18__Clip_spikes_before_STA_19_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-18__Clip_spikes_before_STA_21_0.png

Synaptic conductance

Δg_syn = 0.8 * nS
τ_syn = 7 * ms;

g_syn = v.calc_synaptic_conductance(tg, all_incoming_spikes, Δg_syn, τ_syn)

plt.plot(t_slice, g_syn[i_slice]);
../_images/2020-09-18__Clip_spikes_before_STA_24_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: 1.05 s

plt.plot(t_slice, sim.V_m[i_slice]);
../_images/2020-09-18__Clip_spikes_before_STA_28_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-18__Clip_spikes_before_STA_30_0.png

Clip spikes

import seaborn as sns

%%time
ax = sns.distplot(np.array(Vm_noisy.in_display_units))
ax.set_xlabel("VI signal (mV)");

Wall time: 41.8 s

Text(0.5, 0, 'VI signal (mV)')
../_images/2020-09-18__Clip_spikes_before_STA_33_2.png

We’ll arbitrarily put a threshold here:

clip_threshold = -20 * mV;

Vm_noisy_clipped = np.clip(Vm_noisy, a_min=None, a_max= clip_threshold).in_units(mV);

fig, ax = plt.subplots()
ax.plot(t_slice, Vm_noisy[i_slice])
ax.plot(t_slice, Vm_noisy_clipped[i_slice]);
../_images/2020-09-18__Clip_spikes_before_STA_37_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

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_clipped_arr = np.asarray(Vm_noisy_clipped)
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_clipped_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_clipped.name).in_units(mV)

spike_train__example = all_spike_trains[0]
spike_indices__example = get_spike_indices(spike_train__example)
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, :]);
../_images/2020-09-18__Clip_spikes_before_STA_44_0.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-18__Clip_spikes_before_STA_46_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-18__Clip_spikes_before_STA_48_0.png 
plot_STA(spike_trains_unconnected[0])
../_images/2020-09-18__Clip_spikes_before_STA_49_0.png

Plot STAs 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: 31.9 s

# 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-18__Clip_spikes_before_STA_52_0.png

This was the first STAs plot, with unclipped spikes (and also with τ_syn = 7 ms):

../_images/2020-07-29__STA_58_0.png

Next steps:

  • calculate p(connected) for every presynaptic neuron

  • describe influence of N_in, p_connected, SNR, ..

  • VI realism

previous

2020-10-23 • Delete signal around spikes before STA

next

2020-09-10 • STA for ≠ PSP shapes