2021-11-05 • Vary I/E synaptic conductances
Contents
2021-11-05 • Vary I/E synaptic conductances¶
Prelude¶
from voltage_to_wiring_sim.notebook_init import *
Preloading: numpy, numba, matplotlib.pyplot, seaborn.
Importing from submodules … ✔
Imported `np`, `mpl`, `plt`, `sns`, `pd`
Imported codebase (`voltage_to_wiring_sim`) as `v`
Imported `*` from `v.support.units`
Setup autoreload
v.print_reproducibility_info()
This cell was last run by lpxtf3
on DUIP74576
on Thu 11 Nov 2021, at 19:13 (UTC+0000).
Last git commit (Fri 05 Nov 2021, 17:08).
Uncommited changes to 5 files.
Calculate¶
from voltage_to_wiring_sim.experiments.N_to_1_IE import Params, simulate_and_test_connections, eval_performance
base_params = Params()
v.pprint(base_params)
Params
------
sim_duration = 600
timestep = 0.0001
τ_syn = 0.007
neuron_params = {'C': 1e-10, 'a': 30.0, 'b': -2e-09, 'c': -0.05, ...}
imaging_spike_SNR = 20
window_duration = 0.1
num_spike_trains = 30
p_inhibitory = 0.2
p_connected = 0.7
spike_rate = 20
Δg_syn_exc = 4E-10
Δg_syn_inh = 1.6E-09
v_syn_exc = 0
v_syn_inh = -0.065
rng_seed = 0
New parameters: Δg_syn _exc and _inh.
Note that other base parameters have been changed wrt previous notebooks: inh syn rev pot, and 4:1 exc:inh inputs, and 70% connected.
@v.cache_to_disk("2021-11-05__vary_Δgsyn_IE")
def sim_test_eval(params: Params):
sim_data, _, test_summaries = simulate_and_test_connections(params)
evalu_p_0_05, _, AUC, AUC_exc, AUC_inh = eval_performance(sim_data, test_summaries)
output_spike_rate = sim_data.izh_output.spike_times.size / params.sim_duration
return (output_spike_rate,
evalu_p_0_05.TPR_exc, evalu_p_0_05.TPR_inh, evalu_p_0_05.FPR,
AUC, AUC_exc, AUC_inh)
from joblib import Parallel, delayed
runner = Parallel(n_jobs=-1)
# output printing is not done here, but in jupyter terminal.
def run_parallel(f, args):
return runner(delayed(f)(arg) for arg in args);
from dataclasses import replace
from itertools import product
dgsyn_ratios = [0.5, 1, 2, 4, 8, 16, 32]
def get_dgsyn_vals(ratio, dgsyn_exc = 0.4 * nS):
return dict(Δg_syn_exc=dgsyn_exc,
Δg_syn_inh=dgsyn_exc * ratio)
seeds = [0, 1, 2];
seeds = [0, 1, 2, 3, 4, 5];
seeds = np.arange(20);
f = sim_test_eval
args = [replace(base_params,
**get_dgsyn_vals(ratio),
rng_seed=seed,
)
for (ratio, seed) in product(dgsyn_ratios, seeds)];
%%time
results = run_parallel(f, args);
Wall time: 7.59 s
M = np.reshape(results, (len(dgsyn_ratios), len(seeds), -1))
(output_spike_rate,
TPR_exc, TPR_inh, FPR,
AUC, AUC_exc, AUC_inh) = (M[:,:,i] for i in range(M.shape[-1]))
Plot¶
def plot_dots_and_mean(data, ax, x, c="0.2", label=None, ms=6, lw=5):
ax.plot(x, data, "o", alpha=0.3, c=c, ms=ms)
ax.plot(x, np.mean(data, axis=-1), "-", c=c, label=label, lw=lw)
smaller = dict(ms=4, lw=3);
def make_syncond_ax():
fig, ax = plt.subplots()
ax.set_xscale('log')
ax.set_xticks(dgsyn_ratios)
ax.set_xticklabels(f"{d:g}" for d in dgsyn_ratios)
ax.set_xlabel("Synaptic conductance ratio $\; Δg_{inh} \; / \; Δg_{exc}$")
return fig, ax
plot = partial(plot_dots_and_mean, x=dgsyn_ratios);
fig, ax = make_syncond_ax();
ax.set_ylim(0.5, 1.02)
# plot(ax, AUC);
plot(AUC_exc, ax, c=v.color_exc, label="Excitatory conn.")
plot(AUC_inh, ax, c=v.color_inh, label="Inhibitory conn.", **smaller)
ax.set_ylabel("Area under ROC curve")
ax.legend();
![../_images/2021-11-05__vary_syncond_ratio_20_0.png](../_images/2021-11-05__vary_syncond_ratio_20_0.png)
fig, ax = make_syncond_ax();
plot(TPR_exc, ax, c=v.color_exc, label="Excitatory conn.")
plot(TPR_inh, ax, c=v.color_inh, label="Inhibitory conn.", **smaller)
plot(FPR, ax, c=v.color_unconn, label="Non-connections")
ax.set_ylabel("Fraction detected as connection")
# ax.yaxis.set_major_formatter(mpl.ticker.PercentFormatter(xmax=1))
ax.legend();
![../_images/2021-11-05__vary_syncond_ratio_21_0.png](../_images/2021-11-05__vary_syncond_ratio_21_0.png)
np.mean(FPR, axis=-1)
array([0.04444, 0.06111, 0.07778, 0.03889, 0.07222, 0.05, 0.05556])
fig, ax = make_syncond_ax()
plot(output_spike_rate, ax, lw=4, ms=4)
ax.set_ylabel("Output spike rate (Hz)");
![../_images/2021-11-05__vary_syncond_ratio_23_0.png](../_images/2021-11-05__vary_syncond_ratio_23_0.png)
Investigate ROC¶
Why, at ratio = 2, is AUC perfectly 1 always, but FPR goes as high as 40% at p = 0.05?
sim_and_test = v.cache_to_disk(None)(simulate_and_test_connections)
eval_perf = v.cache_to_disk(None)(eval_performance);
d, td, ts = sim_and_test(Params(rng_seed=0, **get_dgsyn_vals(ratio=2)))
evalu_p_0_05, thr_sweep, AUC, AUC_exc, AUC_inh = eval_performance(d, ts);
evalu_p_0_05.FPR
0.2222222222222222
TPR_exc = [c.evaluation.TPR_exc for c in thr_sweep]
TPR_inh = [c.evaluation.TPR_inh for c in thr_sweep]
FPR = [c.evaluation.FPR for c in thr_sweep]
fig, ax = plt.subplots()
step_kwargs = dict(marker=".")
ax.step(FPR, TPR_inh, where="post", clip_on=False, **step_kwargs)
ax.fill_between(FPR, TPR_inh, step="post", alpha=0.1, color="grey")
ax.set_aspect("equal")
ax.set(
xlabel="#FP / #unconnected",
ylabel="#TP / #connected_inh",
xlim=(0, 1),
ylim=(0, 1),
);
![../_images/2021-11-05__vary_syncond_ratio_28_0.png](../_images/2021-11-05__vary_syncond_ratio_28_0.png)
[c.p_value_threshold for c in thr_sweep]
[0.0, 0.01, 0.03, 0.28, 0.52, 0.67, 0.72, 0.88, 0.91, 0.99]
import pandas as pd
df = pd.DataFrame(ts)
df["conn"] = d.is_connected.astype(int)
df["EI"] = ["exc" if e else "inh" for e in d.is_excitatory]
df
p_value | p_value_type | relative_STA_height | conn | EI | |
---|---|---|---|---|---|
0 | 0.01 | < | 1.946213 | 1 | inh |
1 | 0.01 | < | 2.115577 | 1 | inh |
2 | 0.01 | < | 2.096411 | 1 | inh |
3 | 0.01 | < | 2.029915 | 1 | inh |
4 | 0.99 | = | 0.843468 | 0 | inh |
5 | 0.88 | = | 0.897993 | 0 | inh |
6 | 0.01 | < | 2.721153 | 1 | exc |
7 | 0.01 | < | 3.138538 | 1 | exc |
8 | 0.01 | < | 2.724142 | 1 | exc |
9 | 0.01 | < | 3.210947 | 1 | exc |
10 | 0.01 | < | 2.833998 | 1 | exc |
11 | 0.01 | < | 3.325955 | 1 | exc |
12 | 0.01 | < | 2.849604 | 1 | exc |
13 | 0.01 | < | 3.031652 | 1 | exc |
14 | 0.01 | < | 2.975779 | 1 | exc |
15 | 0.01 | < | 3.274779 | 1 | exc |
16 | 0.01 | < | 2.737663 | 1 | exc |
17 | 0.01 | < | 3.190077 | 1 | exc |
18 | 0.01 | < | 2.983896 | 1 | exc |
19 | 0.01 | < | 3.080062 | 1 | exc |
20 | 0.01 | < | 3.150776 | 1 | exc |
21 | 0.01 | < | 2.667729 | 1 | exc |
22 | 0.01 | < | 3.062099 | 1 | exc |
23 | 0.28 | = | 1.049534 | 0 | exc |
24 | 0.72 | = | 0.935894 | 0 | exc |
25 | 0.03 | = | 1.174938 | 0 | exc |
26 | 0.91 | = | 0.900700 | 0 | exc |
27 | 0.03 | = | 1.215529 | 0 | exc |
28 | 0.52 | = | 0.986098 | 0 | exc |
29 | 0.67 | = | 0.953430 | 0 | exc |
So the culprits are rows 25 and 27: unconnected, but only 3 out of the 100 shuffled STA’s were larger than the real STA. (We’d expect it to be 50 / 100, as the reality is that they’re unconnected).
Imagine if there were only unconnected pairs. ..
So the reason for perfect AUC is..
How was it constructed? for all 30 (pre, post)-pairs, test_connection
, and a resulting p-value which gives proportion of shuffled STAs that are larger than real STA.
These p-values are aggregated and made unique, then used to actually classify each pair.
Is it to be expected that FPR = 0.05 at p < 0.05?
Yes.
How to choose p-value threshold?
“at p < x%, x out of 100 shuffled spike trains are allowed to have a higher STA than the real spike train”.
So yes, p 0.05 is weird, would make more sense to say: “NO shuffled spike trains are allwed to have higher STA”
<=> p < 1/num_shuffles i.e. 0.01.
btw is this p-value here compatible with general p-value definition?
“probability of obtaining results at least as extreme as observed if no effect”
-» “probability of obtaining STA at least as high if not connected”
So yes, very compatible.
Is it weird to do thresholding on p-values?
How is definition above expanded to threshold:
“We detect a connection if the probability of obtaining STA at least as large if they were not connected is 5%”
Inspect signals¶
from voltage_to_wiring_sim.experiments.N_to_1_IE import indices_where
fig, axes = plt.subplots(nrows=3, ncols=3, **v.figsize(width=1000))
for i, ratio in enumerate([0.5, 4, 16]):
d, td, ts = sim_and_test(replace(base_params, **get_dgsyn_vals(ratio=ratio)))
vm = d.izh_output.V_m
ax = v.plot_signal(vm.slice(t_start=10*second, duration=3*second) / mV, ax=axes[0, i])
ax.axhline(np.median(vm) / mV, c='C1')
ax.set_ylim(-68, -22);
ax.set_xlabel("Time (s)")
ax.set_title("$Δg_{inh} \; / \; Δg_{exc} = " + f"{ratio:g}$", pad=18)
ax_e = v.plot_STA(td[indices_where(d.is_excitatory & d.is_connected)[0]].original_STA, ax=axes[1, i])
ax_i = v.plot_STA(td[indices_where(d.is_inhibitory & d.is_connected)[0]].original_STA, ax=axes[2, i])
ylims = (min(ax_i.get_ylim()[0], ax_e.get_ylim()[0]),
max(ax_i.get_ylim()[1], ax_e.get_ylim()[1]))
ax_e.set_ylim(ylims)
ax_i.set_ylim(ylims)
if i == 0:
ax.set_ylabel("Membrane pot. (mV)");
ax_e.set_ylabel("STA of VI signal for an" + '\n' + r"$\bf{excitatory}$ conn. (mV)")
ax_i.set_ylabel("STA of VI signal for an" + '\n' + r"$\bf{inhibitory}$ conn. (mV)")
else:
ax_e.set_ylabel(None)
ax_i.set_ylabel(None)
fig.tight_layout()
![../_images/2021-11-05__vary_syncond_ratio_35_0.png](../_images/2021-11-05__vary_syncond_ratio_35_0.png)
d, td, ts = sim_and_test(Params(rng_seed=0, **get_dgsyn_vals(ratio=4)));
v.plot_signal(d.izh_output.V_m.slice(10*second, 1*second) / mV);
![../_images/2021-11-05__vary_syncond_ratio_37_0.png](../_images/2021-11-05__vary_syncond_ratio_37_0.png)
Reproducibility¶
v.print_reproducibility_info(verbose=True)
This cell was last run by lpxtf3
on DUIP74576
on Thu 11 Nov 2021, at 19:16 (UTC+0000).
Last git commit (Fri 05 Nov 2021, 17:08).
Uncommited changes to:
M codebase/voltage_to_wiring_sim/__init__.py
M codebase/voltage_to_wiring_sim/experiments/N_to_1_IE.py
M codebase/voltage_to_wiring_sim/support/plot_style.py
M notebooks/2021-09-16__vary_E_vs_I.ipynb
?? "notebooks/2021-11-05__vary_\316\224gsyn_IE.ipynb"
Platform:
Windows-10
CPython 3.9.6 (C:\miniforge3\python.exe)
Intel(R) Xeon(R) W-2123 CPU @ 3.60GHz
Dependencies of voltage_to_wiring_sim
and their installed versions:
numpy 1.21.1
matplotlib 3.4.2
numba 0.53.1
joblib 1.0.1
seaborn 0.11.1
scipy 1.7.0
preload 2.2
nptyping 1.4.2
Full conda list:
# packages in environment at C:\miniforge3:
#
# Name Version Build Channel
argon2-cffi 20.1.0 py39hb82d6ee_2 conda-forge
async_generator 1.10 py_0 conda-forge
attrs 21.2.0 pyhd8ed1ab_0 conda-forge
backcall 0.2.0 pyh9f0ad1d_0 conda-forge
backports 1.0 py_2 conda-forge
backports.functools_lru_cache 1.6.4 pyhd8ed1ab_0 conda-forge
black 21.9b0 pyhd8ed1ab_1 conda-forge
bleach 3.3.1 pyhd8ed1ab_0 conda-forge
brotlipy 0.7.0 py39hb82d6ee_1001 conda-forge
ca-certificates 2021.10.8 h5b45459_0 conda-forge
certifi 2021.10.8 py39hcbf5309_1 conda-forge
cffi 1.14.6 py39h0878f49_0 conda-forge
chardet 4.0.0 py39hcbf5309_1 conda-forge
charset-normalizer 2.0.0 pyhd8ed1ab_0 conda-forge
click 7.1.2 pypi_0 pypi
colorama 0.4.4 pyh9f0ad1d_0 conda-forge
colorful 0.5.4 pypi_0 pypi
conda 4.10.3 py39hcbf5309_2 conda-forge
conda-package-handling 1.7.3 py39hb3671d1_0 conda-forge
cryptography 3.4.7 py39hd8d06c1_0 conda-forge
cycler 0.10.0 pypi_0 pypi
dataclasses 0.8 pyhc8e2a94_3 conda-forge
debugpy 1.4.1 py39h415ef7b_0 conda-forge
decorator 5.0.9 pyhd8ed1ab_0 conda-forge
defusedxml 0.7.1 pyhd8ed1ab_0 conda-forge
entrypoints 0.3 pyhd8ed1ab_1003 conda-forge
icu 68.2 h0e60522_0 conda-forge
idna 3.1 pyhd3deb0d_0 conda-forge
importlib-metadata 4.6.1 py39hcbf5309_0 conda-forge
ipykernel 6.0.3 py39h832f523_0 conda-forge
ipython 7.25.0 py39h832f523_1 conda-forge
ipython_genutils 0.2.0 py_1 conda-forge
jedi 0.18.0 py39hcbf5309_2 conda-forge
jinja2 3.0.1 pyhd8ed1ab_0 conda-forge
joblib 1.0.1 pypi_0 pypi
jpeg 9d h8ffe710_0 conda-forge
jsonschema 3.2.0 pyhd8ed1ab_3 conda-forge
jupyter_client 6.1.12 pyhd8ed1ab_0 conda-forge
jupyter_contrib_core 0.3.3 py_2 conda-forge
jupyter_contrib_nbextensions 0.5.1 pyhd8ed1ab_2 conda-forge
jupyter_core 4.7.1 py39hcbf5309_0 conda-forge
jupyter_highlight_selected_word 0.2.0 py39hcbf5309_1002 conda-forge
jupyter_latex_envs 1.4.6 pyhd8ed1ab_1002 conda-forge
jupyter_nbextensions_configurator 0.4.1 py39hcbf5309_2 conda-forge
jupyterlab_pygments 0.1.2 pyh9f0ad1d_0 conda-forge
kiwisolver 1.3.1 pypi_0 pypi
libclang 11.1.0 default_h5c34c98_1 conda-forge
libiconv 1.16 he774522_0 conda-forge
libpng 1.6.37 h1d00b33_2 conda-forge
libsodium 1.0.18 h8d14728_1 conda-forge
libxml2 2.9.12 hf5bbc77_0 conda-forge
libxslt 1.1.33 h65864e5_2 conda-forge
libzlib 1.2.11 h8ffe710_1013 conda-forge
llvmlite 0.36.0 pypi_0 pypi
lxml 4.6.3 py39h4fd7cdf_0 conda-forge
markupsafe 2.0.1 py39hb82d6ee_0 conda-forge
matplotlib 3.4.2 pypi_0 pypi
matplotlib-inline 0.1.2 pyhd8ed1ab_2 conda-forge
menuinst 1.4.17 py39hcbf5309_1 conda-forge
miniforge_console_shortcut 2.0 h57928b3_0 conda-forge
mistune 0.8.4 py39hb82d6ee_1004 conda-forge
mypy_extensions 0.4.3 py39hcbf5309_4 conda-forge
nbclient 0.5.3 pyhd8ed1ab_0 conda-forge
nbconvert 6.1.0 py39hcbf5309_0 conda-forge
nbformat 5.1.3 pyhd8ed1ab_0 conda-forge
nest-asyncio 1.5.1 pyhd8ed1ab_0 conda-forge
notebook 6.4.0 pyha770c72_0 conda-forge
nptyping 1.4.2 pypi_0 pypi
numba 0.53.1 pypi_0 pypi
numpy 1.21.1 pypi_0 pypi
openssl 1.1.1l h8ffe710_0 conda-forge
packaging 21.0 pyhd8ed1ab_0 conda-forge
pandas 1.3.1 pypi_0 pypi
pandoc 2.14.1 h8ffe710_0 conda-forge
pandocfilters 1.4.2 py_1 conda-forge
parso 0.8.2 pyhd8ed1ab_0 conda-forge
pathspec 0.9.0 pyhd8ed1ab_0 conda-forge
pickleshare 0.7.5 py_1003 conda-forge
pillow 8.3.1 pypi_0 pypi
pip 21.2.1 pyhd8ed1ab_0 conda-forge
platformdirs 2.3.0 pyhd8ed1ab_0 conda-forge
preload 2.2 pypi_0 pypi
prettyprinter 0.18.0 pypi_0 pypi
prometheus_client 0.11.0 pyhd8ed1ab_0 conda-forge
prompt-toolkit 3.0.19 pyha770c72_0 conda-forge
pycosat 0.6.3 py39hb82d6ee_1006 conda-forge
pycparser 2.20 pyh9f0ad1d_2 conda-forge
pygments 2.9.0 pyhd8ed1ab_0 conda-forge
pympler 0.9 pypi_0 pypi
pyopenssl 20.0.1 pyhd8ed1ab_0 conda-forge
pyparsing 2.4.7 pyh9f0ad1d_0 conda-forge
pyqt 5.12.3 py39hcbf5309_7 conda-forge
pyqt-impl 5.12.3 py39h415ef7b_7 conda-forge
pyqt5-sip 4.19.18 py39h415ef7b_7 conda-forge
pyqtchart 5.12 py39h415ef7b_7 conda-forge
pyqtwebengine 5.12.1 py39h415ef7b_7 conda-forge
pyrsistent 0.17.3 py39hb82d6ee_2 conda-forge
pysocks 1.7.1 py39hcbf5309_3 conda-forge
python 3.9.6 h7840368_1_cpython conda-forge
python-dateutil 2.8.2 pyhd8ed1ab_0 conda-forge
python_abi 3.9 2_cp39 conda-forge
pytz 2021.1 pypi_0 pypi
pywin32 300 py39hb82d6ee_0 conda-forge
pywinpty 1.1.3 py39h99910a6_0 conda-forge
pyyaml 6.0 py39hb82d6ee_0 conda-forge
pyzmq 22.1.0 py39he46f08e_0 conda-forge
qt 5.12.9 h5909a2a_4 conda-forge
regex 2021.10.23 py39hb82d6ee_1 conda-forge
requests 2.26.0 pyhd8ed1ab_0 conda-forge
ruamel_yaml 0.15.80 py39hb82d6ee_1004 conda-forge
scipy 1.7.0 pypi_0 pypi
seaborn 0.11.1 pypi_0 pypi
send2trash 1.7.1 pyhd8ed1ab_0 conda-forge
setuptools 49.6.0 py39hcbf5309_3 conda-forge
six 1.16.0 pyh6c4a22f_0 conda-forge
sqlite 3.36.0 h8ffe710_0 conda-forge
terminado 0.10.1 py39hcbf5309_0 conda-forge
testpath 0.5.0 pyhd8ed1ab_0 conda-forge
tomli 1.2.2 pyhd8ed1ab_0 conda-forge
tornado 6.1 py39hb82d6ee_1 conda-forge
tqdm 4.61.2 pyhd8ed1ab_1 conda-forge
traitlets 5.0.5 py_0 conda-forge
typed-ast 1.4.3 py39hb82d6ee_1 conda-forge
typing_extensions 3.10.0.2 pyha770c72_0 conda-forge
typish 1.9.2 pypi_0 pypi
tzdata 2021a he74cb21_1 conda-forge
ucrt 10.0.20348.0 h57928b3_0 conda-forge
urllib3 1.26.6 pyhd8ed1ab_0 conda-forge
vc 14.2 hb210afc_5 conda-forge
voltage-to-wiring-sim 0.1 dev_0 <develop>
vs2015_runtime 14.29.30037 h902a5da_5 conda-forge
wcwidth 0.2.5 pyh9f0ad1d_2 conda-forge
webencodings 0.5.1 py_1 conda-forge
wheel 0.36.2 pyhd3deb0d_0 conda-forge
win_inet_pton 1.1.0 py39hcbf5309_2 conda-forge
wincertstore 0.2 py39hcbf5309_1006 conda-forge
winpty 0.4.3 4 conda-forge
winshell 0.6 pypi_0 pypi
yaml 0.2.5 he774522_0 conda-forge
zeromq 4.3.4 h0e60522_0 conda-forge
zipp 3.5.0 pyhd8ed1ab_0 conda-forge
zlib 1.2.11 h8ffe710_1013 conda-forge