{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Perlin spontaneous\n",
"\n",
"Author: Sebastian Spreizer\n",
"\n",
"The MIT License"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib notebook\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import nest\n",
"import noise\n",
"import numpy as np\n",
"\n",
"import lib.connectivity_map as cm\n",
"import lib.lcrn_network as lcrn\n",
"import lib.animation as animation\n",
"import lib.plot3d as pl3d\n",
"import lib.colormap as cmap\n",
"import lib.activity_sequence as seq"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Perlin landscape for connectivity"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Network size\n",
"nrowE = ncolE = 120\n",
"nrowI = ncolI = 60\n",
"npopE = nrowE * ncolE\n",
"npopI = nrowI * ncolI\n",
"\n",
"nrow = nrowE\n",
"landscape = np.round(cm.Perlin_uniform(nrow, size=3, base=1) * 7).astype(int)\n",
"move = cm.move(nrow)\n",
"\n",
"fig,ax = plt.subplots(1)\n",
"im = ax.imshow(landscape.reshape(nrow,-1), cmap=cmap.virno())\n",
"plt.colorbar(im,ax=ax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Set Kernel Status"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(0)\n",
"nest.ResetKernel()\n",
"nest.SetKernelStatus({\n",
" 'local_num_threads': 4,\n",
" 'resolution': 0.1,\n",
"})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Create nodes"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"popE = nest.Create('iaf_psc_alpha', npopE, params={\n",
" \"C_m\": 250.0,\n",
" \"E_L\": -70.0,\n",
" \"V_reset\": -70.0,\n",
" \"V_th\": -55.0,\n",
" \"t_ref\": 2.0,\n",
" \"tau_m\": 10.0,\n",
" \"tau_minus\": 20.0,\n",
" \"tau_syn_ex\": 5.0,\n",
" \"tau_syn_in\": 5.0,\n",
"})\n",
"\n",
"popI = nest.Create('iaf_psc_alpha', npopI, params={\n",
" \"C_m\": 250.0,\n",
" \"E_L\": -70.0,\n",
" \"V_reset\": -70.0,\n",
" \"V_th\": -55.0,\n",
" \"t_ref\": 2.0,\n",
" \"tau_m\": 10.0,\n",
" \"tau_minus\": 20.0,\n",
" \"tau_syn_ex\": 5.0,\n",
" \"tau_syn_in\": 5.0,\n",
"})\n",
"pop = popE + popI\n",
"\n",
"# Create devices\n",
"ngE = nest.Create('noise_generator')\n",
"ngI = nest.Create('noise_generator')\n",
"ng = ngE + ngI\n",
"\n",
"stimE = nest.Create('noise_generator')\n",
"sd = nest.Create('spike_detector')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Connect nodes"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"offsetE = popE[0]\n",
"offsetI = popI[0]\n",
"\n",
"p = 0.05\n",
"stdE = 9\n",
"stdI = 12\n",
"g = 8\n",
"shift = 1\n",
"\n",
"\n",
"for idx in range(npopE):\n",
" # E-> E\n",
" source = idx, nrowE, ncolE, nrowE, ncolE, int(p * npopE), stdE, False\n",
" targets, delay = lcrn.lcrn_gauss_targets(*source)\n",
" targets = (targets + shift * move[landscape[idx] % len(move)]) % npopE\n",
" targets = targets[targets != idx] \n",
" nest.Connect([popE[idx]], (targets + offsetE).tolist(), syn_spec={'weight': 10.0})\n",
"\n",
" # E-> I\n",
" source = idx, nrowE, ncolE, nrowI, ncolI, int(p * npopI), stdE / 2, False\n",
" targets, delay = lcrn.lcrn_gauss_targets(*source)\n",
" nest.Connect([popE[idx]], (targets + offsetI).tolist(), syn_spec={'weight': 10.0})\n",
"\n",
"for idx in range(npopI):\n",
" # I-> E\n",
" source = idx, nrowI, ncolI, nrowE, ncolE, int(p * npopE), stdI, False\n",
" targets, delay = lcrn.lcrn_gauss_targets(*source)\n",
" nest.Connect([popI[idx]], (targets + offsetE).tolist(), syn_spec={'weight': g * -10.0})\n",
"\n",
" # I-> I\n",
" source = idx, nrowI, ncolI, nrowI, ncolI, int(p * npopI), stdI / 2, False\n",
" targets, delay = lcrn.lcrn_gauss_targets(*source)\n",
" targets = targets[targets != idx]\n",
" nest.Connect([popI[idx]], (targets + offsetI).tolist(), syn_spec={'weight': g * -10.0})\n",
"\n",
"# Connect noise input device to all neurons\n",
"nest.Connect(ngE, popE, syn_spec={'weight': 10.0})\n",
"nest.Connect(ngI, popI, syn_spec={'weight': 10.0})\n",
"\n",
"centerE = nrowE * nrowE // 2 + nrowE // 2\n",
"source = centerE, nrowE, ncolE, nrowE, ncolE, 50, 4\n",
"targets, delay = lcrn.lcrn_gauss_targets(*source)\n",
"targets = np.unique(targets)\n",
"nest.Connect(stimE, (targets + offsetE).tolist(), syn_spec={'weight': 10.0})\n",
"\n",
"nest.Connect(pop, sd)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Warming up"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nest.SetStatus(ng, params={'mean': 0., 'std': 50.})\n",
"nest.Simulate(250.)\n",
"nest.SetStatus(ng, params={'mean': 15., 'std': 30.})\n",
"nest.Simulate(250.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Start simulation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for i in range(4):\n",
" nest.SetStatus(stimE, params={'mean': 50., 'std': 0.})\n",
" nest.Simulate(50.)\n",
" nest.SetStatus(stimE, params={'mean': 0., 'std': 0.})\n",
" nest.Simulate(450.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Plot spiking activity"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sdE = nest.GetStatus(sd, 'events')[0]\n",
"ts, gids = sdE['times'], sdE['senders']\n",
"fig, ax = plt.subplots(1)\n",
"ax.plot(ts, gids, '|')\n",
"ax.set_xlabel('Time [ms]')\n",
"ax.set_ylabel('Neuron')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"idx = gids - offsetE < npopE\n",
"gids, ts = gids[idx] - offsetE, ts[idx]\n",
"time = nest.GetKernelStatus('time')\n",
"\n",
"ts_bins = np.arange(time-2000., time, 10.)\n",
"h = np.histogram2d(ts, gids, bins=[ts_bins, range(npopE + 1)])[0]\n",
"hh = h.reshape(-1, nrowE, ncolE)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig,ax = plt.subplots(1)\n",
"im = ax.imshow(hh.sum(0) / 2., cmap='binary')\n",
"plt.colorbar(im,ax=ax)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1)\n",
"im1 = ax.imshow(hh[0], vmin=0, vmax=np.max(hh), cmap='binary_r')\n",
"anim = animation.imshow(fig, ax, im1, hh, ts_bins)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#animation.HTML(anim.to_jshtml())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#fig,ax = pl3d.scatter(ts,gids%nrow,gids//nrow)\n",
"#ax.set_xlim(time-1000.,time)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fig,ax = pl.subplots(1,2)\n",
"# ax[0].imshow(landscape.reshape(nrow,-1), cmap=cmap.virno(), vmin=0, vmax=2*np.pi)\n",
"# ax[1].imshow(W_EE.reshape(nrow,-1))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"td = 3.14\n",
"eps = 3.14"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"clusters1, sequences1 = seq.identify_vectors(ts, gids, nrowE, ncolE, steps=10., width=10., td=td, eps=eps)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"center = int(nrow * nrow/2 + nrow/2)\n",
"clusters2, sequences2 = seq.identify_vectors(ts, (gids + center) % npopE, nrowE, ncolE, steps=10., width=10., td=td, eps=eps)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"c1 = np.copy(clusters1[1])\n",
"c2 = clusters2[1]\n",
"\n",
"for ii in range(10):\n",
" x = np.unique(list(zip(c2,c1)), axis=0)\n",
" x2,x1 = x.T\n",
" x2_set = np.unique(x2, return_counts=True)\n",
" b2 = x2_set[0][x2_set[1] > 1][1:]\n",
" for i in b2:\n",
" dd = x1[x2 == i]\n",
" dd = dd[dd>-1]\n",
" idx = np.in1d(c1, dd)\n",
" c1[idx] = dd[0]\n",
"\n",
"c1 = np.unique(c1 + 1, return_inverse=True)[1] - 1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cmod=20\n",
"idx = c1 > -1\n",
"fig, ax = plt.subplots(1)\n",
"ax.scatter(ts[idx][::5], gids[idx][::5], c=c1[idx][::5]%cmod, s=5, cmap=cmap.virno())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cmod = 10\n",
"idx = c1 > -1\n",
"fig, ax = plt.subplots(1)\n",
"ax.scatter(ts[idx][::5] % 500., gids[idx][::5], c=c1[idx][::5]%cmod, s=5, cmap=cmap.virno())\n",
"ax.set_xlim(0,500)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fig = plt.figure(figsize=(10,8))\n",
"# ax = fig.add_subplot(111, projection='3d')\n",
"# x,y,z = ts,gids%nrow,gids//nrow\n",
"# ax.scatter(x[idx][::10], y[idx][::10], z[idx][::10], c=c[idx][::10]%cmod, s=5, cmap=cmap.virno())\n",
"# ax.set_ylim(0,nrow)\n",
"# ax.set_zlim(0,nrow)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cmod = 10\n",
"fig, ax = plt.subplots(1)\n",
"ti = ts\n",
"xi = gids % nrow\n",
"yi = gids // nrow\n",
"ci = c1 % cmod\n",
"ci[c1 == -1] = -1\n",
"\n",
"scat = ax.scatter(xi, yi, c=ci, s=5, vmin=-1, vmax=cmod, cmap=cmap.virno())\n",
"ax.set_xlim(0,nrow)\n",
"ax.set_ylim(0,nrow)\n",
"anim = animation.scatter(fig, ax, scat, ti, xi, yi, ci, ts_bins)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"animation.HTML(anim.to_jshtml())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#anim.save('EI_networks-stimulus.mp4', fps=10., extra_args=['-vcodec', 'libx264'])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}