2022-02-21 • FENS 2022 Abstract
Contents
2022-02-21 • FENS 2022 Abstract¶
As usual, the reusable code of the previous notebook
has been added in pkg/VoltageToMap/src/
,
and is imported below.
Setup¶
#
using Revise # Reloads src code without having to restart kernel
using MyToolbox
using VoltageToMap
[ Info: Precompiling VoltageToMap [b3b8fdc5-3c26-4000-a0c8-f17415fdf48e]
Params & sim¶
Short warm-up run. Get compilation out of the way.
p0 = ExperimentParams(sim=SimParams(inputs=previous_N_30_inputs, duration=1 * minutes));
@time sim(p0.sim);
2.311696 seconds (1.21 M allocations: 169.872 MiB, 91.71% gc time)
p = ExperimentParams(
sim = SimParams(
inputs = realistic_N_6600_inputs,
duration = 10 * minutes,
synapses = SynapseParams(Δg_multiplier = 0.066),
)
)
dumpc(p)
ExperimentParams
seed: 22022022
sim: SimParams
duration: 600.0
Δt: 0.0001
num_timesteps: 6000000
seed: 0
inputs: PoissonInputsParams
N_unconn: 100
N_exc: 5200
N_inh: 1300
N_conn: 6500
N: 6600
spike_rates: LogNormal
μ: 1.08629
σ: 0.774597
synapses: SynapseParams
Δg_exc: 4.0e-10
Δg_inh: 1.6e-9
Δg_multiplier: 0.066
E_exc: 0.0
E_inh: -0.065
g_t0: 0.0
τ: 0.007
izh_neuron: IzhikevichParams
C: 1.0e-10
k: 7.0e-7
v_rest: -0.06
v_thr: -0.04
a: 30.0
b: -2.0e-9
v_peak: 0.035
v_reset: -0.05
Δu: 1.0e-10
v_t0: -0.06
u_t0: 0.0
imaging: VoltageImagingParams
spike_SNR: 10.0
spike_SNR_dB: 20.0
spike_height: 0.095
σ_noise: 0.0095
conntest: ConnTestParams
STA_window_length: 0.1
num_shuffles: 100
seed: 22022022
evaluation: EvaluationParams
num_tested_neurons_per_group: 40
seed: 22022022
t, v, input_spikes = @time sim(p.sim);
Progress: 100%|█████████████████████████████████████████| Time: 0:02:39
159.438690 seconds (18.23 M allocations: 2.165 GiB, 10.06% gc time, 0.01% compilation time)
num_spikes = length.(input_spikes)
ComponentVector{Int64}(conn = (exc = [1228, 940, 1389, 948, 509, 997, 831, 527, 1041, 876 … 1059, 997, 1200, 1103, 350, 585, 1106, 386, 1203, 798], inh = [470, 979, 876, 183, 882, 746, 1152, 1032, 201, 661 … 1287, 677, 865, 576, 1076, 1185, 1101, 1065, 927, 1334]), unconn = [908, 983, 717, 1052, 1400, 946, 680, 1444, 759, 612 … 1435, 868, 854, 329, 1235, 1300, 1394, 725, 505, 390])
Plot¶
using Sciplotlib
""" tzoom = [200ms, 600ms] e.g. """
function plotsig(t, sig, tzoom = nothing; ax = nothing, clip_on=false, kw...)
isnothing(tzoom) && (tzoom = t[[1, end]])
izoom = first(tzoom) .≤ t .≤ last(tzoom)
if isnothing(ax)
plot(t[izoom], sig[izoom]; clip_on, kw...)
else
plot(t[izoom], sig[izoom], ax; clip_on, kw...)
end
end;
plotsig(t, v/mV, [0s, 4seconds], xlabel="Time (s)", hylabel="mV");
![../_images/2022-02-21b__FENS_abstract_17_0.png](../_images/2022-02-21b__FENS_abstract_17_0.png)
Imaging noise¶
resetrng!(p.sim.seed)
noise = randn(length(v)) * p.sim.imaging.σ_noise
vimsig = v + noise;
ax = plotsig(t, vimsig / mV, [200ms,1200ms], xlabel="Time (s)", hylabel="mV", alpha=0.7);
plotsig(t, v / mV, [200ms,1200ms], xlabel="Time (s)", hylabel="mV"; ax);
![../_images/2022-02-21b__FENS_abstract_20_0.png](../_images/2022-02-21b__FENS_abstract_20_0.png)
Window¶
const window_length = p.conntest.STA_window_length
const Δt = p.sim.Δt
const win_size = round(Int, window_length / Δt)
const t_win = linspace(zero(window_length), window_length, win_size) # for plottin
function calc_STA(presynaptic_spikes)
STA = zeros(eltype(vimsig), win_size)
win_starts = round.(Int, presynaptic_spikes / Δt)
num_wins = 0
for a in win_starts
b = a + win_size - 1
if b ≤ lastindex(vimsig)
STA .+= @view vimsig[a:b]
num_wins += 1
end
end
STA ./= num_wins
return STA
end;
function plotSTA(presynspikes)
STA = calc_STA(presynspikes)
plot(t_win/ms, STA/mV)
end;
example_presynspikes = rand(input_spikes.conn.exc)
plotSTA(example_presynspikes);
![../_images/2022-02-21b__FENS_abstract_24_0.png](../_images/2022-02-21b__FENS_abstract_24_0.png)
Test connection¶
to_ISIs(spiketimes) = [first(spiketimes); diff(spiketimes)] # copying
to_spiketimes!(ISIs) = cumsum!(ISIs, ISIs) # in place
(example_presynspikes |> to_ISIs |> to_spiketimes!) ≈ example_presynspikes # test
true
shuffle_ISIs(spiketimes) = to_spiketimes!(shuffle!(to_ISIs(spiketimes)));
test_statistic(spiketimes) = spiketimes |> calc_STA |> mean;
Note difference with 2021: there it was peak-to-peak (max - min). Here it is mean.
const num_shuffles = p.conntest.num_shuffles
function test_connection(presynspikes)
real_t = test_statistic(presynspikes)
shuffled_t = Vector{typeof(real_t)}(undef, num_shuffles)
resetrng!(p.conntest.seed)
for i in eachindex(shuffled_t)
shuffled_t[i] = test_statistic(shuffle_ISIs(presynspikes))
end
N_shuffled_larger = count(shuffled_t .> real_t)
return if N_shuffled_larger == 0
p_value = 1 / num_shuffles
else
p_value = N_shuffled_larger / num_shuffles
end
end;
Results¶
p_exc
num_trains = p.evaluation.num_tested_neurons_per_group
num_trains = 300
resetrng!(p.evaluation.seed)
tested_spike_trains_exc = rand(input_spikes.conn.exc, num_trains)
tested_spike_trains_inh = rand(input_spikes.conn.inh, num_trains)
tested_spike_trains_unconn = rand(input_spikes.unconn, min(num_trains, 100))
p_exc = Float64[]
p_inh = Float64[]
p_unconn = Float64[]
for (groupname, spiketrains, pvals) in (
("excitatory", tested_spike_trains_exc, p_exc),
("inhibitory", tested_spike_trains_inh, p_inh),
("unconnected", tested_spike_trains_unconn, p_unconn),
)
@showprogress 200ms groupname for spiketrain in spiketrains
push!(pvals, test_connection(spiketrain))
end
end
excitatory100%|█████████████████████████████████████████| Time: 0:01:38
inhibitory100%|█████████████████████████████████████████| Time: 0:01:30
unconnected100%|████████████████████████████████████████| Time: 0:00:33
best_exc = tested_spike_trains_exc[p_exc .== minimum(p_exc)]
best_inh = tested_spike_trains_inh[p_inh .== maximum(p_inh)]
length(best_exc), length(best_inh)
(8, 19)
st_exc_i, st_exc = rand(pairs(best_exc))
st_inh_i, st_inh = rand(pairs(best_inh))
plotSTA(st_exc)
plotSTA(st_inh)
st_exc_i, st_inh_i
![../_images/2022-02-21b__FENS_abstract_35_0.png](../_images/2022-02-21b__FENS_abstract_35_0.png)
(2, 9)
fig, ax = plt.subplots(figsize=(3.4,3))
function plotdot(y, x, c, jitter=0.28)
N = length(y)
x -= 0.35
plot(x*ones(N) + (rand(N).-0.5)*jitter, y, "o", color=c, ms=4.2, markerfacecolor="none", clip_on=false)
plot(x+0.35, mean(y), "k.", ms=10)
end
plotdot(p_exc, 1, "C2"); ax.text(1-0.16, -0.1, "excitatory"; color="C2", ha="center")
plotdot(p_unconn, 2, "C0"); ax.text(2-0.16, -0.1, "unconnected"; color="C0", ha="center")
plotdot(p_inh, 3, "C1"); ax.text(3-0.16, -0.1, "inhibitory"; color="C1", ha="center")
ax.boxplot([p_exc, p_unconn, p_inh], widths=0.2, medianprops=Dict("color"=>"black"))
Sciplotlib.set(ax, xlim=(0.33, 3.3), ylim=(0, 1), xaxis=:off)
hylabel(ax, L"p(\, \mathrm{shuffled\ \overline{STA}} \ > \ \mathrm{actual\ \overline{STA}}\, )"; dy=10);
![../_images/2022-02-21b__FENS_abstract_36_0.png](../_images/2022-02-21b__FENS_abstract_36_0.png)