diff --git a/examples/tvb_netpyne/__init__.py b/examples/tvb_netpyne/__init__.py
index e69de29bb..47564932c 100644
--- a/examples/tvb_netpyne/__init__.py
+++ b/examples/tvb_netpyne/__init__.py
@@ -0,0 +1 @@
+pip3 import tvb
\ No newline at end of file
diff --git a/examples/tvb_netpyne/models/ET_wilson_cowan.py b/examples/tvb_netpyne/models/ET_wilson_cowan.py
new file mode 100644
index 000000000..413032555
--- /dev/null
+++ b/examples/tvb_netpyne/models/ET_wilson_cowan.py
@@ -0,0 +1,167 @@
+# -*- coding: utf-8 -*-
+
+from tvb_multiscale.tvb_netpyne.interfaces.models.wilson_cowan import WilsonCowanTVBNetpyneInterfaceBuilder
+from tvb_multiscale.tvb_netpyne.netpyne_models.models.wilson_cowan import WilsonCowanBuilder
+from tvb_multiscale.tvb_netpyne.netpyne_models.models.thalamic_VIM_exc_io_inh_i import WilsonCowanThalamicVIMBuilder
+
+from examples.tvb_netpyne.example import main_example
+from examples.models.wilson_cowan import wilson_cowan_example as wilson_cowan_example_base
+# 0------
+from IPython.display import Image, display
+import os
+from collections import OrderedDict
+import time
+import numpy as np
+from tvb.basic.profile import TvbProfile
+TvbProfile.set_profile(TvbProfile.LIBRARY_PROFILE)
+# 1 ------
+from tvb.datatypes.connectivity import Connectivity
+from tvb.simulator.integrators import HeunStochastic
+from tvb.simulator.monitors import Raw # , Bold, EEG
+
+from tvb_multiscale.tvb_netpyne.config import *
+
+def wilson_cowan_example(**kwargs):
+ params = {
+ "simulation_length": 500,
+ # "model_params": {"lamda": 0.5}
+ }
+
+ #(0) Setup stuff
+ # For a minimal example, select:
+ N_REGIONS = None # total TVB brain regions
+ NETPYNE_NODES_INDS = np.array([0, 1]) # the brain region nodes to place spiking networks from [0, N_REGIONS-1] interval
+ N_NEURONS = 100 # number of neurons per spiking population
+
+ # Interface basic configurations:
+ INTERFACE_MODEL = "RATE" # The only available option for NetPyNE so far
+ INTERFACE_COUPLING_MODE = "spikeNet" # "spikeNet" # "TVB"
+ # -----------------------------------------------
+ work_path = os.getcwd()
+ outputs_path = os.path.join(work_path, "outputs/NetPyNE_wilsoncowan_%s_%s" % (INTERFACE_MODEL, INTERFACE_COUPLING_MODE))
+
+ config = Config(output_base=outputs_path)
+ config.figures.SHOW_FLAG = True
+ config.figures.SAVE_FLAG = True
+ config.figures.FIG_FORMAT = 'png'
+ config.figures.DEFAULT_SIZE= config.figures.NOTEBOOK_SIZE
+ FIGSIZE = config.figures.DEFAULT_SIZE
+
+ from tvb_multiscale.core.plot.plotter import Plotter
+ plotter = Plotter(config.figures)
+
+ #1. Load structural data (minimally a TVB connectivity) & prepare TVB simulator (region mean field model, integrator, monitors etc)
+ # This would run on TVB only before creating any multiscale cosimulation interface connections.
+ from tvb_multiscale.core.tvb.cosimulator.models.reduced_wong_wang_exc_io_inh_i import ReducedWongWangExcIOInhI
+
+ # Load full TVB connectome connectivity: Glasser parcellation
+ data_path = os.path.expanduser("~/packages/tvb-multiscale/examples/tvb_netpyne/models")
+ conn_path = os.path.join(data_path, "temp_connectivity_Glasser")
+
+ w = np.loadtxt(os.path.join(conn_path, "weights.txt")) #np.zeros((11,11))
+ c = np.loadtxt(os.path.join(conn_path, "centres.txt"), usecols=range(1,3))
+ rl = np.loadtxt(os.path.join(conn_path, "region_labels.txt"), dtype="str", usecols = 1)
+ t = np.loadtxt(os.path.join(conn_path, "tract_lengths.txt"))
+
+ ######## PETERSEN ATLAS ############
+ #doesn't do anything for now, except change variable name used from this point onwards
+ # later, if we add additional connections from e.g. another atlas, we can add/modify/merge the values/labels and specify the new filenames here
+ wTVB = w
+ cTVB = c
+ rlTVB = rl
+ tTVB = t
+ #####################################
+
+ number_of_regions = len(rlTVB)
+ speed = np.array([4.0])
+
+ connTVB = Connectivity(region_labels=rlTVB, weights=wTVB, centres=cTVB, tract_lengths=tTVB, speed=speed)
+
+ # -------------- Pick a minimal brain of only the first N_REGIONS regions: ----------------
+ if N_REGIONS is not None:
+ connTVB.number_of_regions = N_REGIONS
+ connTVB.region_labels = connTVB.region_labels[:N_REGIONS]
+ connTVB.centres = connTVB.centres[:N_REGIONS]
+ connTVB.areas = connTVB.areas[:N_REGIONS]
+ connTVB.orientations = connTVB.orientations[:N_REGIONS]
+ connTVB.hemispheres = connTVB.hemispheres[:N_REGIONS]
+ connTVB.cortical = connTVB.cortical[:N_REGIONS]
+ connTVB.weights = connTVB.weights[:N_REGIONS][:, :N_REGIONS]
+ connTVB.tract_lengths = connTVB.tract_lengths[:N_REGIONS][:, :N_REGIONS]
+ # -----------------------------------------------------------------------------------------
+
+ # Remove diagonal self-connections:
+ np.fill_diagonal(connTVB.weights, 0.0)
+
+ # Normalize connectivity weights
+ connTVB.weights = connTVB.scaled_weights(mode="region")
+ connTVB.weights /= np.percentile(connTVB.weights, 99)
+
+ # CoSimulator builder--------------------------------
+ from tvb_multiscale.core.tvb.cosimulator.cosimulator_builder import CoSimulatorSerialBuilder
+
+ simulator_builder = CoSimulatorSerialBuilder()
+ simulator_builder.config = config
+ simulator_builder.model = ThalamicVIMBuilder()
+ simulator_builder.connectivity = connTVB
+ simulator_builder.model_params = model_params
+ simulator_builder.initial_conditions = np.array([0.0])
+
+ simulator_builder.configure()
+ simulator_builder.print_summary_info_details(recursive=1)
+
+ simulator = simulator_builder.build()
+ # -----
+
+ simulator.configure()
+
+
+ simulator.print_summary_info_details(recursive=1)
+
+ # Plot TVB connectome:
+ plotter.plot_tvb_connectivity(simulator.connectivity)
+
+
+ #(2) NetPyNE Network Builder
+ netpyne_network_builder = WilsonCowanThalamicVIMBuilder()
+
+ #change connectivity to Glasser
+ kwargs['connectivity']='/home/docker/packages/tvb_data/tvb_data/connectivity/connectivity_Glasser.zip'
+ kwargs.update(params)
+
+ netpyne_network_builder.tvb_to_spiking_dt_ratio = 2 # 2 NetPyNE integration steps for 1 TVB integration step
+ netpyne_network_builder.monitor_period = 1.0
+
+ #TVB (Glasser parcellation) regions to be substituted with spiking populations (netpyne pops)
+ #M1R (8), PMC (8), brainstem (379), thalamusR (371), cerebellumL (361), muscle (NA) ###DEFAULT[50, 58,] #50:precentral_R, 58:superiofrontal_R (taken as temp proxy for premotor cortex)
+
+ netpyne_network_builder.primary_motor_cortex_R = 8
+ netpyne_network_builder.cerebellar_cortex_L = 361
+ netpyne_network_builder.thalamus_R = 371
+ netpyne_network_builder.brainstem = 379
+
+ params["spiking_proxy_inds"] = [
+ netpyne_network_builder.primary_motor_cortex_R,
+ netpyne_network_builder.cerebellar_cortex_L,
+ netpyne_network_builder.thalamus_R,
+ netpyne_network_builder.brainstem,
+ ]
+ kwargs.update(params)
+
+ tvb_netpyne_interface_builder = WilsonCowanTVBNetpyneInterfaceBuilder()
+
+ # input/output interfaces
+ import numpy as np
+ tvb_netpyne_interface_builder.default_coupling_mode == "spikeNet"
+ if tvb_netpyne_interface_builder.default_coupling_mode == "spikeNet":
+ pass
+ # coupling is computing in NetPyNE, we need as many TVB proxies
+ # as TVB regions coupling to STN and Striatum
+ #??? proxy_inds = np.unique(np.concatenate([CTXtoSTNinds, CTXtoSNinds]))
+
+
+ return main_example(wilson_cowan_example_base, netpyne_network_builder, tvb_netpyne_interface_builder, **kwargs)
+
+
+if __name__ == "__main__":
+ wilson_cowan_example()
diff --git a/examples/tvb_netpyne/models/red_wong_wang.py b/examples/tvb_netpyne/models/red_wong_wang.py
index 0b1493d5c..4906834b5 100644
--- a/examples/tvb_netpyne/models/red_wong_wang.py
+++ b/examples/tvb_netpyne/models/red_wong_wang.py
@@ -1,4 +1,5 @@
from tvb_multiscale.tvb_netpyne.netpyne_models.models.default_exc_io_inh_i import DefaultExcIOInhIBuilder
+from tvb_multiscale.tvb_netpyne.netpyne_models.models.thalamic_VIM_exc_io_inh_i import ThalamicVIMBuilder
from tvb_multiscale.tvb_netpyne.interfaces.models.default import RedWongWangExcIOInhITVBNetpyneInterfaceBuilder
from examples.tvb_netpyne.example import main_example
@@ -10,8 +11,12 @@ def excio_inhi_example(**kwargs):
# "model_params": {"lamda": 0.5}
}
- netpyne_network_builder = DefaultExcIOInhIBuilder()
- params["spiking_proxy_inds"] = [0, 1]
+ # netpyne_network_builder = DefaultExcIOInhIBuilder()
+ # params["spiking_proxy_inds"] = [0, 1]
+
+ netpyne_network_builder = ThalamicVIMBuilder()
+ params["spiking_proxy_inds"] = [46]
+
kwargs.update(params)
diff --git a/examples/tvb_netpyne/models/temp_connectivity_Glasser.zip b/examples/tvb_netpyne/models/temp_connectivity_Glasser.zip
new file mode 100644
index 000000000..a630b645c
Binary files /dev/null and b/examples/tvb_netpyne/models/temp_connectivity_Glasser.zip differ
diff --git a/examples/tvb_netpyne/models/wilson_cowan.py b/examples/tvb_netpyne/models/wilson_cowan.py
deleted file mode 100644
index a43f21388..000000000
--- a/examples/tvb_netpyne/models/wilson_cowan.py
+++ /dev/null
@@ -1,21 +0,0 @@
-# -*- coding: utf-8 -*-
-
-from tvb_multiscale.tvb_netpyne.interfaces.models.wilson_cowan import WilsonCowanTVBNetpyneInterfaceBuilder
-from tvb_multiscale.tvb_netpyne.netpyne_models.models.wilson_cowan import WilsonCowanBuilder
-
-from examples.tvb_netpyne.example import main_example
-from examples.models.wilson_cowan import wilson_cowan_example as wilson_cowan_example_base
-
-
-def wilson_cowan_example(**kwargs):
- # model_params = {"model_params": {"lamda": 0.5}}
- # kwargs.update(model_params)
-
- netpyne_model_builder = WilsonCowanBuilder()
- tvb_netpyne_interface_builder = WilsonCowanTVBNetpyneInterfaceBuilder()
-
- return main_example(wilson_cowan_example_base, netpyne_model_builder, tvb_netpyne_interface_builder, **kwargs)
-
-
-if __name__ == "__main__":
- wilson_cowan_example()
diff --git a/examples/tvb_netpyne/notebooks/Roopa-TVB-corr-cortex-opt.ipynb b/examples/tvb_netpyne/notebooks/Roopa-TVB-corr-cortex-opt.ipynb
new file mode 100644
index 000000000..ade54004a
--- /dev/null
+++ b/examples/tvb_netpyne/notebooks/Roopa-TVB-corr-cortex-opt.ipynb
@@ -0,0 +1,2369 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# TVB-NetPyNE: Bridging multiscale activity by co-simulation\n",
+ "\n",
+ "## Step-by-step learn how to perform a co-simulation embedding spiking neural networks into large-scale brain networks using TVB."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from IPython.display import Image, display"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Reduced Wong-Wang TVB mean field model\n",
+ "\n",
+ "For every region node $n\\prime$ modelled as a mean-field node in TVB:\n",
+ "\n",
+ "(Post)Synaptic gating dynamics (i.e., proportion of synapse channels open at any given time):\n",
+ "\n",
+ "$\\dot{S_{n\\prime}} = - \\frac{1}{\\tau}{S_{n\\prime}}(t) + (1-{S_{n\\prime}}(t))\\gamma {R_{n\\prime}}(t)$\n",
+ "\n",
+ "and $ {R_{n\\prime}}(t) $ is the postsynaptic firing rate given by:\n",
+ "\n",
+ "$ {R_{n\\prime}}(t) = H({I_{syn_{n\\prime}}}(t), a, b, d) $\n",
+ "\n",
+ "where\n",
+ "\n",
+ "$ H({I_{syn_{n\\prime}}}(t), a, b, d) = \\frac{aI_{syn_{n\\prime}}(t)-b}{1-e^{-d(a{I_{syn_{n\\prime}}}(t)-b)}}$ \n",
+ "\n",
+ "is a sigmoidal activation function of the input presynaptic current.\n",
+ "\n",
+ "The total input presynaptic current to excitatory populations is given by: \n",
+ "\n",
+ "$ {I_{syn_{n\\prime}}}(t) = I_o + w_+J_{N}{S_{n\\prime}}(t) + GJ_{N}\\sum_{{m\\prime}\\neq {n\\prime}}C_{{m\\prime}{n\\prime}}S_{m\\prime}(t-\\tau_{{m\\prime}{n\\prime}})$"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Reduced Wong-Wang mean field model\n",
+ "\n",
+ "### Parameters following Deco et al 2013:\n",
+ "\n",
+ "- structural TVB connectivity weights $w_{{m\\prime}{n\\prime}}$ (${m\\prime}->{n\\prime}$)\n",
+ "- structural TVB connectivity delays $\\tau_{{m\\prime}{n\\prime}}$ (${m\\prime}->{n\\prime}$)\n",
+ "- global structural brain connectivity coupling constant $G$\n",
+ "- overall effective external input current $I_o = 0.3nA$ \n",
+ "- excitatory synaptic coupling $J_{N} = 0.2609nA$ \n",
+ "- local excitatory recurrence $w_+ = 0.9$\n",
+ "- excitatory kinetic parameter $\\gamma = 0.641 s$\n",
+ "- excitatory sigmoidal functions parameters $a = 2710nC^{-1}$, $b = 108Hz$, $d = 0.154s$"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Izhikevich Spiking network model in ANNarchy\n",
+ "\n",
+ "For every neuron $i$ in region node $n$ modelled in ANNarchy as a spiking network:\n",
+ "\n",
+ "Membrane potential:\n",
+ "\n",
+ "$ \\dot{V}_m = n_2V_m^2 + n_1V_m + n_0140 - U_m/C $\n",
+ "\n",
+ "$\\;\\;\\;\\;\\;\\;\\;- g_{AMPA}(V_m-E_{AMPA}) - g_{GABA}(V_m-E_{GABA}) - g_{BASE}V_m + I_e $\n",
+ "\n",
+ "where the conductances follow the equations:\n",
+ "\n",
+ "$ \\dot{g}_{AMPA} = - g_{AMPA} / \\tau_{AMPA} + \\left[\\sum_k \\delta(t-t_k) \\right]_{Exc}$\n",
+ "\n",
+ "$ \\dot{g}_{GABA} = - g_{GABA} / \\tau_{GABA} + \\left[\\sum_k \\delta(t-t_k) \\right]_{Inh}$\n",
+ "\n",
+ "$ \\dot{g}_{BASE} = - g_{BASE} / \\tau_{BASE} + \\left[\\sum_k \\delta(t-t_k) \\right]_{BASE}$\n",
+ "\n",
+ "and recovery variable:\n",
+ "\n",
+ "$ \\dot{U}_m = a(bV_m - U_m)$\n",
+ "\n",
+ "\n",
+ "When $ V_m > V_{th} $ , $ V_m $ is set to $ c $, and $ U_m $ is incremented by $ d $."
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# WORKFLOW:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-07-12T20:35:57.561354Z",
+ "start_time": "2019-07-12T20:35:52.475653Z"
+ },
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ },
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "from collections import OrderedDict\n",
+ "import time\n",
+ "import numpy as np\n",
+ "\n",
+ "from tvb.basic.profile import TvbProfile\n",
+ "TvbProfile.set_profile(TvbProfile.LIBRARY_PROFILE)\n",
+ "\n",
+ "from tvb_multiscale.tvb_netpyne.config import *\n",
+ "\n",
+ "work_path = os.getcwd()\n",
+ "data_path = os.path.expanduser(\"~/packages/tvb-multiscale/examples/data/basal_ganglia\")\n",
+ "fit_data_path = os.path.join(data_path, \"ANNarchyFittedModels/dataFits_2020_02_05/databestfits\", )\n",
+ "data_mode = \"patient\" # \"control\", or \"patient\"\n",
+ "control_data = os.path.join(fit_data_path, \"controlleft/OutputSim_Patient08.mat\")\n",
+ "patient_data = os.path.join(fit_data_path, \"patientleft/OutputSim_Patient09.mat\")\n",
+ "\n",
+ "INTERFACE_COUPLING_MODE = \"spikeNet\" # \"TVB\" or \"spikeNet\"\n",
+ "\n",
+ "if data_mode == \"patient\":\n",
+ " subject_data = patient_data\n",
+ " if INTERFACE_COUPLING_MODE == \"TVB\":\n",
+ " STN_factor = 20.5 * 1e-2 # 0.205 \n",
+ " dSN_factor = 22.2 * 1e-2 # 0.222 \n",
+ " iSN_factor = 20.2 * 1e-2\n",
+ " else:\n",
+ " STN_factor = 55.0 * 1e-3 # 0.055 - 29.5\n",
+ " dSN_factor = 35.5 * 1e-3 # 0.0355 - 19.7\n",
+ " iSN_factor = 21.3 * 1e-3 # 0.0213 - 9.5\n",
+ "else:\n",
+ " subject_data = control_data\n",
+ " if INTERFACE_COUPLING_MODE == \"TVB\":\n",
+ " STN_factor = 20.2 * 1e-2 # 0.202 \n",
+ " dSN_factor = 23.4 * 1e-2 # 0.2340 \n",
+ " iSN_factor = 21.1 * 1e-2 # 0.211 \n",
+ " else:\n",
+ " STN_factor = 45.6 * 1e-3 \n",
+ " dSN_factor = 30.0 * 1e-3 \n",
+ " iSN_factor = 19.5 * 1e-3 \n",
+ " \n",
+ "simulation_length = 150.0\n",
+ "transient = 25.0\n",
+ "start_stimulus = 100.0\n",
+ "init_cond_jitter = 0.0\n",
+ "\n",
+ "SPIKING_NODES_DELAYS = False\n",
+ "\n",
+ "simulation_mode = \"rs\" # \"stim\" or \"rs\"\n",
+ "stim_target = \"GPi\" # \"STN\", \"GPi\"\n",
+ "stim_mode = \"simple\" # \"bi\" | \"mono\" | \"simple\"\n",
+ " # -------------------------\n",
+ "stim_freq = 0.0 # 130.0 | 120.0 | 0.0 \n",
+ "stim_ampl = -10.0 # 20.0 | -35.0 | -10.0 \n",
+ "stim_duration = 0.0 # 0.3 | 0.3 | 0.0\n",
+ "if simulation_mode == \"stim\":\n",
+ " simulation_mode = simulation_mode + \"_%s_%s\" % (stim_target, stim_mode)\n",
+ "\n",
+ "outputs_path = os.path.join(work_path, \"outputs\")\n",
+ "sim_mode_path = os.path.join(outputs_path, \"TVBcortex_%s_coupl\" % INTERFACE_COUPLING_MODE, \n",
+ " data_mode, simulation_mode)\n",
+ "config = Config(output_base=sim_mode_path)\n",
+ "config.figures.SHOW_FLAG = True \n",
+ "config.figures.SAVE_FLAG = True\n",
+ "config.figures.FIG_FORMAT = 'png'\n",
+ "config.figures.DEFAULT_SIZE= config.figures.NOTEBOOK_SIZE\n",
+ "FIGSIZE = config.figures.DEFAULT_SIZE\n",
+ "\n",
+ "from tvb_multiscale.core.plot.plotter import Plotter\n",
+ "plotter = Plotter(config.figures)\n",
+ "\n",
+ "# For interactive plotting:\n",
+ "# %matplotlib notebook \n",
+ "\n",
+ "# Otherwise:\n",
+ "%matplotlib inline \n"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## 1. Load structural data
(minimally a TVB connectivity)
& prepare TVB simulator
(region mean field model, integrator, monitors etc)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-07-12T20:35:59.690799Z",
+ "start_time": "2019-07-12T20:35:57.571529Z"
+ },
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ },
+ "scrolled": true,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from tvb_multiscale.core.tvb.cosimulator.models.reduced_wong_wang_exc_io import ReducedWongWangExcIO\n",
+ "\n",
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "# ----Uncomment below to modify the simulator by changing the default options:--------------------------------------\n",
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "\n",
+ "from tvb.datatypes.connectivity import Connectivity\n",
+ "from tvb.simulator.integrators import HeunStochastic\n",
+ "from tvb.simulator.monitors import Raw # , Bold, EEG\n",
+ " \n",
+ "\n",
+ "# 0. GPe_Left, 1. GPi_Left, 2. STN_Left, 3. Striatum_Left, 4. Thal_Left\n",
+ "BG_opt_matrix_weights = np.zeros((5, 5))\n",
+ "conn_mode = \"subject\" # subject\" or \"average\"\n",
+ "if conn_mode == \"average\":\n",
+ " weights_maith = np.array([1.93, 3.56, 1.46, 4.51, 3.52, 2.30, 2.34, 3.78, 1.98, \n",
+ " 1.30, 1.82, 3.56, 3.02, 1.78, 1.36, 2.27, 4.13, 2.74, 3.27])*1e-3 # controls\n",
+ "# weights_maith = np.array([3.27, 3.80, 2.65, 3.66, 3.06, 3.06, 3.25, 4.02, 3.32, \n",
+ "# 2.98, 3.45, 3.64, 2.50, 2.12, 2.86, 2.79, 3.96, 3.69, 3.87])*1e-3 # patients\n",
+ " # probs_maith = ????\n",
+ "else:\n",
+ " import scipy.io as sio\n",
+ " weights=sio.loadmat(subject_data) # weights start from index 19\n",
+ " weights_maith = weights[\"X\"][0, 19:] # these are indices 19 till 37\n",
+ " probs_maith = weights[\"X\"][0, :19] # these are indices 0 till 18\n",
+ "\n",
+ "wdSNGPi = BG_opt_matrix_weights[3, 1] = weights_maith[0].item()\n",
+ "wiSNGPe = BG_opt_matrix_weights[3, 0] = weights_maith[1].item()\n",
+ "wGPeSTN = BG_opt_matrix_weights[0, 2] = weights_maith[2].item()\n",
+ "wSTNGPe = BG_opt_matrix_weights[2, 0] = weights_maith[3].item()\n",
+ "wSTNGPi = BG_opt_matrix_weights[2, 1] = weights_maith[4].item()\n",
+ "wGPeGPi = BG_opt_matrix_weights[0, 1] = weights_maith[5].item() \n",
+ "wGPiTh = BG_opt_matrix_weights[1, 4] = weights_maith[8].item()\n",
+ "wThdSN = BG_opt_matrix_weights[4, 3] = weights_maith[10].item() # Th -> dSN\n",
+ " \n",
+ "sliceBGnet = slice(0,5)\n",
+ "\n",
+ "wGPeGPe = weights_maith[6].item() # \"GPe\" -> \"GPe\" \n",
+ "wGPiGPi = weights_maith[7].item() # \"GPi\" -> \"GPi\" \n",
+ "wThiSN = weights_maith[9].item() # \"Eth\" -> \"IiSN\" \n",
+ "\n",
+ "wdSNdSN = weights_maith[11].item() # \"IdSN\" -> \"IdSN\" \n",
+ "wiSNiSN = weights_maith[12].item() # \"IiSN\" -> \"IiSN\" \n",
+ "wCtxdSN = weights_maith[13].item() # \"CxE\" -> \"IdSN\" \n",
+ "wCtxiSN = weights_maith[14].item() # \"CxE\" -> \"IiSN\" \n",
+ "wCtxSTN = weights_maith[15].item() # \"CxE\" -> \"Estn\"\n",
+ "wCtxEtoI = weights_maith[16].item() # \"CxE\" -> \"CxI\"\n",
+ "wCtxItoE = weights_maith[17].item() # \"CxI\" -> \"CxE\"\n",
+ "wCtxItoI = weights_maith[18].item() # \"CxI\" -> \"CxI\"\n",
+ "\n",
+ "pdSNGPi = probs_maith[0].item()\n",
+ "piSNGPe = probs_maith[1].item()\n",
+ "pGPeSTN = probs_maith[2].item()\n",
+ "pSTNGPe = probs_maith[3].item()\n",
+ "pSTNGPi = probs_maith[4].item()\n",
+ "pGPeGPi = probs_maith[5].item()\n",
+ "pGPeGPe = probs_maith[6].item() # \"GPe\" -> \"GPe\" \n",
+ "pGPiGPi = probs_maith[7].item() # \"GPi\" -> \"GPi\" \n",
+ "pGPiTh = probs_maith[8].item()\n",
+ "pThiSN = probs_maith[9].item() # \"Eth\" -> \"IiSN\n",
+ "pThdSN = probs_maith[10].item() # Th --> dSN\n",
+ "pdSNdSN = probs_maith[11].item() # \"IdSN\" -> \"IdSN\" \n",
+ "piSNiSN = probs_maith[12].item() # \"IiSN\" -> \"IiSN\" \n",
+ "pCtxdSN = probs_maith[13].item() # \"CxE\" -> \"IdSN\" \n",
+ "pCtxiSN = probs_maith[14].item() # \"CxE\" -> \"IiSN\" \n",
+ "pCtxSTN = probs_maith[15].item() # \"CxE\" -> \"Estn\"\n",
+ "pCtxEtoI = probs_maith[16].item() # \"CxE\" -> \"CxI\"\n",
+ "pCtxItoE = probs_maith[17].item() # \"CxI\" -> \"CxE\"\n",
+ "pCtxItoI = probs_maith[18].item() # \"CxI\" -> \"CxI\"\n",
+ "pCtxCtx = probs_maith[16:19].mean() # \"Ctx\" -> \"Ctx\"\n",
+ "\n",
+ "loadedParams ={'dSNGPi_probs': probs_maith[0],\n",
+ " \t'dSNGPi_weights' : weights_maith[0],\n",
+ " \t'iSNGPe_probs' : probs_maith[1],\n",
+ " \t'iSNGPe_weights' : weights_maith[1],\n",
+ " \t'GPeSTN_probs' : probs_maith[2],\n",
+ " \t'GPeSTN_weights' : weights_maith[2],\n",
+ " \t'STNGPe_probs' : probs_maith[3],\n",
+ " \t'STNGPe_weights' : weights_maith[3],\n",
+ " \t'STNGPi_probs' : probs_maith[4],\n",
+ " \t'STNGPi_weights' : weights_maith[4],\n",
+ " \t'GPeGPi_probs' : probs_maith[5],\n",
+ " \t'GPeGPi_weights' : weights_maith[5],\n",
+ " \t'GPeGPe_probs' : probs_maith[6],\n",
+ " \t'GPeGPe_weights' : weights_maith[6],\n",
+ " \t'GPiGPi_probs' : probs_maith[7],\n",
+ " \t'GPiGPi_weights' : weights_maith[7],\n",
+ " \t'GPiThal_probs' : probs_maith[8],\n",
+ " \t'GPiThal_weights' : weights_maith[8],\n",
+ " \t'ThaliSN_probs' : probs_maith[9],\n",
+ " \t'ThaliSN_weights' : weights_maith[9],\n",
+ " \t'ThaldSN_probs' : probs_maith[10],\n",
+ " \t'ThaldSN_weights' : weights_maith[10],\n",
+ " \t'dSNdSN_probs' : probs_maith[11],\n",
+ " \t'dSNdSN_weights' : weights_maith[11],\n",
+ " \t'iSNiSN_probs' : probs_maith[12],\n",
+ " \t'iSNiSN_weights' : weights_maith[12],\n",
+ " \t'CtxdSN_probs' : probs_maith[13],\n",
+ " \t'CtxdSN_weights' : weights_maith[13],\n",
+ " \t'CtxiSN_probs' : probs_maith[14],\n",
+ " \t'CtxiSN_weights' : weights_maith[14],\n",
+ " \t'CtxSTN_probs' : probs_maith[15],\n",
+ " \t'CtxSTN_weights' : weights_maith[15],\n",
+ " \t'CtxECtxI_probs' : probs_maith[16],\n",
+ " \t'CtxECtxI_weights' : weights_maith[16],\n",
+ " \t'CtxICtxE_probs' : probs_maith[17],\n",
+ " \t'CtxICtxE_weights' : weights_maith[17],\n",
+ " \t'CtxICtxI_probs' : probs_maith[18],\n",
+ " \t'CtxICtxI_weights' : weights_maith[18],\n",
+ " 'CtxThal_weights': 0.0,\n",
+ " 'CtxThal_probs': 1.0}\n",
+ "print(loadedParams)\n",
+ "\n",
+ "assert_loadedParams = dict(zip(loadedParams.values(), loadedParams.keys()))\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "display(Image(filename='pics/Connectome.png', width=1000, unconfined=False))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true,
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "# Load full TVB connectome connectivity\n",
+ "\n",
+ "conn_path = os.path.join(data_path, \"conn\")\n",
+ "\n",
+ "#Load AAL atlas normative connectome including the Basal Ganglia regions from Petersen et al. atlas\n",
+ "wTVB = np.loadtxt(os.path.join(conn_path, \"conn_denis_weights.txt\"))\n",
+ "cTVB = np.loadtxt(os.path.join(conn_path, \"aal_plus_BG_centers.txt\"),usecols=range(1,4))\n",
+ "rlTVB = np.loadtxt(os.path.join(conn_path, \"aal_plus_BG_centers.txt\"),dtype=\"str\", usecols=(0,))\n",
+ "tlTVB = np.loadtxt(os.path.join(conn_path, \"BGplusAAL_tract_lengths.txt\"))\n",
+ "\n",
+ "# Remove the second Thalamus, Pallidum (GPe/i), Putamen and Caudate (Striatum):\n",
+ "inds_Th = (rlTVB.tolist().index(\"Thalamus_L\"), rlTVB.tolist().index(\"Thalamus_R\"))\n",
+ "inds_Pall = (rlTVB.tolist().index(\"Pallidum_L\"), rlTVB.tolist().index(\"Pallidum_R\"))\n",
+ "inds_Put = (rlTVB.tolist().index(\"Putamen_L\"), rlTVB.tolist().index(\"Putamen_R\"))\n",
+ "inds_Caud = (rlTVB.tolist().index(\"Caudate_L\"), rlTVB.tolist().index(\"Caudate_R\"))\n",
+ "inds_rm = inds_Th + inds_Pall + inds_Put + inds_Caud\n",
+ "print(\"Connections of Thalami, Pallidum (GPe/i), Putamen and Caudate (Striatum) removed!:\\n\", \n",
+ " wTVB[inds_rm, :][:, inds_rm])\n",
+ "wTVB = np.delete(wTVB, inds_rm, axis=0)\n",
+ "wTVB = np.delete(wTVB, inds_rm, axis=1)\n",
+ "tlTVB = np.delete(tlTVB, inds_rm, axis=0)\n",
+ "tlTVB = np.delete(tlTVB, inds_rm, axis=1)\n",
+ "rlTVB = np.delete(rlTVB, inds_rm, axis=0)\n",
+ "cTVB = np.delete(cTVB, inds_rm, axis=0)\n",
+ "\n",
+ "number_of_regions = len(rlTVB)\n",
+ "speed = np.array([4.0])\n",
+ "min_tt = speed.item() * 0.1\n",
+ "sliceBG = [0, 1, 2, 3, 6, 7]\n",
+ "sliceCortex = slice(10, number_of_regions)\n",
+ "\n",
+ "# Remove BG -> Cortex connections\n",
+ "print(\"Removing BG -> Cortex connections with max:\")\n",
+ "print(wTVB[sliceCortex, :][:, sliceBG].max())\n",
+ "wTVB[sliceCortex, sliceBG] = 0.0\n",
+ "tlTVB[sliceCortex, sliceBG] = min_tt\n",
+ "\n",
+ "# Remove Cortex -> Thalamus connections\n",
+ "sliceThal = [8, 9]\n",
+ "print(\"Removing Cortex -> Thalamus connections with summed weight:\")\n",
+ "print(wTVB[sliceThal, sliceCortex].sum())\n",
+ "wTVB[sliceThal, sliceCortex] = 0.0\n",
+ "tlTVB[sliceThal, sliceCortex] = min_tt\n",
+ "\n",
+ "# Remove Cortex -> GPe/i connections\n",
+ "sliceGP = [0, 1, 2, 3]\n",
+ "print(\"Removing Cortex -> GPe/i connections with max:\")\n",
+ "print(wTVB[sliceGP, sliceCortex].max())\n",
+ "wTVB[sliceGP, sliceCortex] = 0.0\n",
+ "tlTVB[sliceGP, sliceCortex] = min_tt\n",
+ "\n",
+ "# # Minimize all delays for the optimized network\n",
+ "# tlTVB[:7][:, :7] = min_tt\n",
+ "\n",
+ "connTVB = Connectivity(region_labels=rlTVB, weights=wTVB, centres=cTVB, tract_lengths=tlTVB, speed=speed)\n",
+ "\n",
+ "# Normalize connectivity weights\n",
+ "connTVB.weights = connTVB.scaled_weights(mode=\"region\")\n",
+ "connTVB.weights /= np.percentile(connTVB.weights, 99)\n",
+ "\n",
+ "# Keep only left hemisphere and remove Vermis:\n",
+ "sliceLeft = slice(0, connTVB.number_of_regions -8, 2)\n",
+ "\n",
+ "connLeft = Connectivity(region_labels=connTVB.region_labels[sliceLeft], \n",
+ " centres=connTVB.centres[sliceLeft],\n",
+ " weights=connTVB.weights[sliceLeft][:, sliceLeft],\n",
+ " tract_lengths=connTVB.tract_lengths[sliceLeft][:, sliceLeft], \n",
+ " speed=connTVB.speed)\n",
+ "connLeft.configure()\n",
+ "\n",
+ "print(\"\\nLeft cortex connectome, after removing direct BG -> Cortex and intehemispheric BG <-> BG connections:\")\n",
+ "plotter.plot_tvb_connectivity(connLeft);"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true,
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "sliceBGnet = slice(0,5)\n",
+ "connTVBleftBG = Connectivity(region_labels=connLeft.region_labels[sliceBGnet], \n",
+ " centres=connLeft.centres[sliceBGnet],\n",
+ " weights=connLeft.weights[sliceBGnet][:, sliceBGnet],\n",
+ " tract_lengths=connLeft.tract_lengths[sliceBGnet][:, sliceBGnet], \n",
+ " speed=connLeft.speed)\n",
+ "connTVBleftBG.configure()\n",
+ "\n",
+ "print(\"\\nLeft BG TVB network:\")\n",
+ "plotter.plot_tvb_connectivity(connTVBleftBG);\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "scaleBGoptTOtvb = np.percentile(BG_opt_matrix_weights, 95) /\\\n",
+ " np.percentile(connTVBleftBG.weights, 95)\n",
+ " \n",
+ "print(\"Scaling factor of TVB BG network connectome to optimal one = %g\" % scaleBGoptTOtvb)\n",
+ "# confitmed: 0.000771577\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true,
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "# Construct the final connectivity to use for simulation:\n",
+ "# Rescale the \n",
+ "ww = scaleBGoptTOtvb * np.array(connLeft.weights)\n",
+ "ww[sliceBGnet, sliceBGnet] = BG_opt_matrix_weights.T # !!!NOTE TVB indices convention Wi<-j!!!\n",
+ "\n",
+ "connectivity = Connectivity(region_labels=connLeft.region_labels, \n",
+ " centres=connLeft.centres,\n",
+ " weights=ww, tract_lengths=connLeft.tract_lengths, \n",
+ " speed=connLeft.speed)\n",
+ "connectivity.configure()\n",
+ "\n",
+ "# Construct only the optimized BG connectivity only for plotting:\n",
+ "connBGopt = Connectivity(region_labels=connectivity.region_labels[sliceBGnet], \n",
+ " centres=connectivity.centres[sliceBGnet],\n",
+ " weights=connectivity.weights[sliceBGnet][:, sliceBGnet],\n",
+ " tract_lengths=connectivity.tract_lengths[sliceBGnet][:, sliceBGnet], \n",
+ " speed=connectivity.speed)\n",
+ "connBGopt.configure()\n",
+ "\n",
+ "print(\"\\nLeft BG optimized network:\")\n",
+ "plotter.plot_tvb_connectivity(connBGopt);\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "from tvb_multiscale.core.tvb.cosimulator.cosimulator_serial import CoSimulatorSerial as CoSimulator\n",
+ "\n",
+ "#white_matter_coupling = coupling.Linear(a=0.014)\n",
+ "# Create a TVB simulator and set all desired inputs\n",
+ "# (connectivity, model, surface, stimuli etc)\n",
+ "# We choose all defaults in this example\n",
+ "simulator = CoSimulator()\n",
+ "#simulator.use_numba = False\n",
+ "model_params = {\"G\": np.array([15.0/scaleBGoptTOtvb])}\n",
+ "simulator.model = ReducedWongWangExcIO(**model_params)\n",
+ "\n",
+ "simulator.connectivity = connectivity\n",
+ "\n",
+ "simulator.integrator = HeunStochastic()\n",
+ "simulator.integrator.dt = 0.1\n",
+ "simulator.integrator.noise.nsig = np.array([1e-4]) # 1e-5\n",
+ "\n",
+ "mon_raw = Raw(period=1.0) # ms\n",
+ "simulator.monitors = (mon_raw, )\n",
+ "\n",
+ "init_cond_filepath = os.path.join(data_path, \"tvb_init_cond_left.npy\")\n",
+ "init_cond = np.load(init_cond_filepath) # \n",
+ "init_cond = np.abs(init_cond *(1 + init_cond_jitter * np.random.normal(size=init_cond.shape)))\n",
+ "simulator.connectivity.set_idelays(simulator.integrator.dt)\n",
+ "simulator.initial_conditions = init_cond * np.ones((simulator.connectivity.idelays.max(),\n",
+ " simulator.model.nvar,\n",
+ " simulator.connectivity.number_of_regions,\n",
+ " simulator.model.number_of_modes))\n",
+ "\n",
+ "\n",
+ "print(\"\\nConnectome used for simulations:\")\n",
+ "plotter.plot_tvb_connectivity(simulator.connectivity);\n",
+ "\n",
+ "\n",
+ "# # Serializing TVB cosimulator is necessary for parallel cosimulation:\n",
+ "# from tvb_multiscale.core.utils.file_utils import dump_pickled_dict\n",
+ "# from tvb_multiscale.core.tvb.cosimulator.cosimulator_serialization import serialize_tvb_cosimulator\n",
+ "# sim_serial_filepath = os.path.join(config.out.FOLDER_RES, \"tvb_serial_cosimulator.pkl\")\n",
+ "# sim_serial = serialize_tvb_cosimulator(simulator)\n",
+ "# display(sim_serial)\n",
+ "\n",
+ "# # Dumping the serialized TVB cosimulator to a file will be necessary for parallel cosimulation.\n",
+ "# dump_pickled_dict(sim_serial, sim_serial_filepath)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "display(Image(filename='pics/Network.png', width=1000, unconfined=False))"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2. Build and connect the ANNarchy network model
(networks of spiking neural populations for fine-scale
regions, stimulation devices, spike detectors etc)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-07-12T20:36:10.862262Z",
+ "start_time": "2019-07-12T20:36:10.000332Z"
+ },
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ },
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "from tvb_multiscale.tvb_netpyne.netpyne_models.builders.base import NetpyneNetworkBuilder\n",
+ "\n",
+ "from ANNarchy import HomogeneousCorrelatedSpikeTrains # PoissonPopulation\n",
+ "\n",
+ "# Select the regions for the fine scale modeling with ANNarchy spiking networks\n",
+ "spiking_nodes_ids = [0, 1, 2, 3, 4] # the indices of fine scale regions modeled with NetPyNE\n",
+ "\n",
+ "# Build a ANNarchy network model with the corresponding builder\n",
+ "from tvb_multiscale.tvb_netpyne.netpyne_models.builders.netpyne_factory import load_netpyne\n",
+ "netpyne = load_netpyne(config=config)\n",
+ "\n",
+ "spiking_network_builder = NetpyneNetworkBuilder(simulator, spiking_nodes_ids, netpyne_instance=netpyne, config=config)\n",
+ "\n",
+ "\n",
+ "# # ----------------------------------------------------------------------------------------------------------------\n",
+ "# # ----Uncomment below to modify the builder by changing the default options:--------------------------------------\n",
+ "# # ----------------------------------------------------------------------------------------------------------------\n",
+ "from copy import deepcopy\n",
+ "\n",
+ "spiking_network_builder.population_order = 200 # reduce for speed\n",
+ "\n",
+ "# When any of the properties model, params and scale below depends on regions,\n",
+ "# set a handle to a function with\n",
+ "# arguments (region_index=None) returning the corresponding property\n",
+ "\n",
+ "# was: below\n",
+ "# ann_model_builder.params_common = \\\n",
+ "# {\"E_ampa\": 0.0, \"E_gaba\": -90.0, \"v_th\": 30.0, \"Vr\": 0.0, \"c\": -65.0,\n",
+ "# \"C\": 1.0, \"I\": 0.0,\n",
+ "# \"tau_syn\": 1.0, \"tau_ampa\": 10.0, \"tau_gaba\": 10.0,\n",
+ "# \"n0\": 140.0, \"n1\": 5.0, \"n2\": 0.04, \n",
+ "# \"v\": -72.0, \"u\": -14.0, \n",
+ "# \"noise\": 0.0}\n",
+ "\n",
+ "# ann_model_builder._paramsI = deepcopy(ann_model_builder.params_common)\n",
+ "# ann_model_builder._paramsI.update({\"a\": 0.005, \"b\": 0.585, \"d\": 4.0, \n",
+ "# \"v\": -70.0, \"u\": -18.55})\n",
+ "# ann_model_builder._paramsE = deepcopy(ann_model_builder.params_common)\n",
+ "# ann_model_builder._paramsE.update({\"v\": -70.0, \"u\": -18.55})\n",
+ "# ann_model_builder.paramsStr = deepcopy(ann_model_builder.params_common)\n",
+ "# ann_model_builder.paramsStr.update({\"v_th\": 40.0, \"C\": 50.0, \"Vr\": -80.0,\n",
+ "# \"n0\": 61.65119, \"n1\": 2.594639, \"n2\": 0.022799, \n",
+ "# \"a\": 0.05, \"b\": -20.0, \"c\": -55.0, \"d\": 377.0, \n",
+ "# \"v\": -70.0, \"u\": -18.55})\n",
+ "\n",
+ "# ann_model_builder.Igpe_nodes_ids = [0]\n",
+ "# ann_model_builder.Igpi_nodes_ids = [1]\n",
+ "# ann_model_builder.Estn_nodes_ids = [2]\n",
+ "# ann_model_builder.Eth_nodes_ids = [4]\n",
+ "# ann_model_builder.Istr_nodes_ids = [3]\n",
+ "\n",
+ "from netpyne import specs\n",
+ "netParams = specs.NetParams()\n",
+ "cfg = specs.SimConfig()\n",
+ "cfg.recordTraces = {'V_soma':{'sec':'soma','pointp':'izhi','var':'V'}}\n",
+ "\n",
+ "def izhiCellParams(a, b, c, d, C, Ie, n0, n1, n2, uInit=-18.55, vInit=-70, vR=0, vThresh=30):\n",
+ " params = {\n",
+ " 'mod': 'Izhi2003c', # model name (corresponds to what described in .mod file)\n",
+ " 'vref': 'V', # custom variable for membrane potential in Izhi2003c\n",
+ " 'a': a, 'b': b, 'c': c, 'd': d,\n",
+ " 'C': C,\n",
+ " 'Ie': Ie,\n",
+ " 'n0': n0, 'n1': n1, 'n2': n2,\n",
+ " 'uInit': uInit, 'vInit': vInit,\n",
+ " 'Vr': vR,\n",
+ " 'thresh': vThresh\n",
+ " }\n",
+ " return {\n",
+ " 'secs': {\n",
+ " 'soma': {'pointps': {'izhi': params}}\n",
+ " }\n",
+ " }\n",
+ "\n",
+ "# a b c d C Ie n0 n1 n2\n",
+ "netParams.cellParams['Str'] = izhiCellParams(0.05, -20, -55, 377, 50, 0, 61.65119, 2.594639, 0.022799, vR=-80, vThresh=40)\n",
+ "netParams.cellParams['GPi'] = izhiCellParams(0.005, 0.585, -65, 4, 1, 30, 140, 5, 0.04)\n",
+ "netParams.cellParams['GPe'] = izhiCellParams(0.005, 0.585, -65, 4, 1, 12, 140, 5, 0.04)\n",
+ "netParams.cellParams['STN'] = izhiCellParams(0.005, 0.265, -65, 2, 1, 3, 140, 5, 0.04)\n",
+ "netParams.cellParams['Thl'] = izhiCellParams(0.02, 0.25, -65, 0.05, 1, 3.5, 140, 5, 0.04)\n",
+ "netParams.cellParams['CxE'] = izhiCellParams(0.02, 0.2, -72, 6, 1, 7, 140, 5, 0.04, uInit=-14.0, vInit=-72) # I = 7 in resting, 50 when fitting\n",
+ "netParams.cellParams['CxI'] = izhiCellParams(0.02, 0.2, -72, 6, 1, 0, 140, 5, 0.04, uInit=-14.0, vInit=-72)\n",
+ "\n",
+ "Igpe_nodes_ids = [0]\n",
+ "Igpi_nodes_ids = [1]\n",
+ "Estn_nodes_ids = [2]\n",
+ "Eth_nodes_ids = [4]\n",
+ "Istr_nodes_ids = [3]\n",
+ "\n",
+ "I_nodes_ids = Igpe_nodes_ids + Igpi_nodes_ids\n",
+ "E_nodes_ids = Estn_nodes_ids + Eth_nodes_ids\n",
+ "\n",
+ "\n",
+ "def paramsE_fun(node_id):\n",
+ " paramsE = {} # deepcopy(ann_model_builder._paramsE)\n",
+ " if node_id in Estn_nodes_ids:\n",
+ " paramsE.update({\"a\": 0.005, \"b\": 0.265, \"d\": 2.0, \"I\": 3.0}) # dictionary of params for Estn\n",
+ " elif node_id in Eth_nodes_ids:\n",
+ " paramsE.update({\"a\": 0.02, \"b\": 0.25, \"d\": 0.05, \"I\": 3.5}) # dictionary of params for Eth\n",
+ " return paramsE\n",
+ " \n",
+ "def paramsI_fun(node_id):\n",
+ " # For the moment they are identical, unless you differentiate the noise parameters\n",
+ " paramsI = {} # deepcopy(ann_model_builder._paramsI)\n",
+ " if node_id in Igpe_nodes_ids:\n",
+ " paramsI.update({\"I\": 12.0})\n",
+ " elif node_id in Igpi_nodes_ids:\n",
+ " paramsI.update({\"I\": 30.0})\n",
+ " return paramsI\n",
+ "\n",
+ "paramsStr = {}\n",
+ " \n",
+ "# Populations' configurations\n",
+ "# When any of the properties model, params and scale below depends on regions,\n",
+ "# set a handle to a function with\n",
+ "# arguments (region_index=None) returning the corresponding property\n",
+ "# TODO: allow for having several nodes in single entry\n",
+ "spiking_network_builder.populations = [\n",
+ " {\"label\": \"E\", \"model\": \"STN\", \n",
+ " \"params\": paramsE_fun, \n",
+ " \"nodes\": Estn_nodes_ids, # E_nodes_ids, # Estn in [2] # \n",
+ " \"scale\": 1.0},\n",
+ " {\"label\": \"E\", \"model\": \"Thl\", \n",
+ " \"params\": paramsE_fun, \n",
+ " \"nodes\": Eth_nodes_ids, # Eth in [4]\n",
+ " \"scale\": 1.0},\n",
+ " {\"label\": \"I\", \"model\": \"GPe\", \n",
+ " \"params\": paramsI_fun, \n",
+ " \"nodes\": Igpe_nodes_ids, # I_nodes_ids, # Igpe in [0] # \n",
+ " \"scale\": 1.0},\n",
+ " {\"label\": \"I\", \"model\": \"GPi\", \n",
+ " \"params\": paramsI_fun, \n",
+ " \"nodes\": Igpi_nodes_ids, # Igpi in [1]\n",
+ " \"scale\": 1.0},\n",
+ " {\"label\": \"IdSN\", \"model\": \"Str\", \n",
+ " \"params\": paramsStr, \n",
+ " \"nodes\": Istr_nodes_ids, # IdSN in [3]\n",
+ " \"scale\": 1.0},\n",
+ " {\"label\": \"IiSN\", \"model\": \"Str\", # IiSN in [3]\n",
+ " \"params\": paramsStr, \n",
+ " \"nodes\": Istr_nodes_ids, # None means \"all\"\n",
+ " \"scale\": 1.0}\n",
+ "]\n",
+ "\n",
+ "# Within region-node connections\n",
+ "# When any of the properties model, conn_spec, weight, delay, receptor_type below\n",
+ "# set a handle to a function with\n",
+ "# arguments (region_index=None) returning the corresponding property\n",
+ "\n",
+ "synapse_model = None\n",
+ "# conn_spec = {'rule': \"all_to_all\", \n",
+ "# \"allow_self_connections\": True, \"force_multiple_weights\": False}\n",
+ "# conn_spec_fixed_probability = conn_spec.copy()\n",
+ "# conn_spec_fixed_probability.update({'rule': \"fixed_probability\", \"probability\": 0.1})\n",
+ "\n",
+ "def conn_spec_fixed_prob(prob=None):\n",
+ " return {\"rule\": {\"prob\": prob}}\n",
+ "\n",
+ "within_node_delay = 1.0 # ms\n",
+ " \n",
+ "\n",
+ "# for each connection, we have a different probability\n",
+ "spiking_network_builder.populations_connections = [\n",
+ " # source -> target\n",
+ " {\"source\": \"I\", \"target\": \"I\", # I -> I This is a self-connection for population \"Igpe\"\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pGPeGPe), # conn_spec\n",
+ " \"weight\": np.abs(wGPeGPe).item(), \"delay\": within_node_delay,\n",
+ " \"receptor_type\": \"gaba\", \"nodes\": Igpe_nodes_ids}, # None means apply to all\n",
+ " {\"source\": \"I\", \"target\": \"I\", # I -> I This is a self-connection for population \"Igpi\"\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pGPiGPi), # conn_spec\n",
+ " \"weight\": np.abs(wGPiGPi).item(), \"delay\": within_node_delay,\n",
+ " \"receptor_type\": \"gaba\", \"nodes\": Igpi_nodes_ids}, # None means apply to all\n",
+ " {\"source\": \"IdSN\", \"target\": \"IdSN\", # IdSN -> IdSN This is a self-connection for population \"IdSN\"\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pdSNdSN), # conn_spec\n",
+ " \"weight\": np.abs(wdSNdSN).item(), \"delay\": within_node_delay,\n",
+ " \"receptor_type\": \"gaba\", \"nodes\": Istr_nodes_ids},\n",
+ " {\"source\": \"IiSN\", \"target\": \"IiSN\", # IiSN -> IiSN This is a self-connection for population \"IiSN\"\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(piSNiSN), # conn_spec\n",
+ " \"weight\": np.abs(wiSNiSN).item(), \"delay\": within_node_delay,\n",
+ " \"receptor_type\": \"gaba\", \"nodes\": Istr_nodes_ids},\n",
+ " ]\n",
+ "\n",
+ "\n",
+ "# Among/Between region-node connections\n",
+ "# Given that only the AMPA population of one region-node couples to\n",
+ "# all populations of another region-node,\n",
+ "# we need only one connection type\n",
+ " \n",
+ "# When any of the properties model, conn_spec, weight, delay, receptor_type below\n",
+ "# depends on regions, set a handle to a function with\n",
+ "# arguments (source_region_index=None, target_region_index=None)\n",
+ "\n",
+ "from tvb_multiscale.core.spiking_models.builders.templates import scale_tvb_weight, tvb_delay\n",
+ "\n",
+ "# We set global coupling scaling to 1.0,\n",
+ "# because we need the Maith et al optimized weights without any scaling:\n",
+ "spiking_network_builder.global_coupling_scaling = 1.0 \n",
+ " \n",
+ "# Function that will return the TVB weight with optional scaling:\n",
+ "class TVBWeightFun(object):\n",
+ " \n",
+ " def __init__(self, \n",
+ " global_coupling_scaling=spiking_network_builder.global_coupling_scaling, \n",
+ " tvb_weights = spiking_network_builder.tvb_weights):\n",
+ " self.global_coupling_scaling = float(global_coupling_scaling)\n",
+ " self.tvb_weights = tvb_weights.copy()\n",
+ " \n",
+ " def __call__(self, source_node, target_node):\n",
+ " return scale_tvb_weight(source_node, target_node, self.tvb_weights,\n",
+ " scale=self.global_coupling_scaling)\n",
+ "\n",
+ "# Function that will return the TVB delay unless SPIKING_NODES_DELAYS == False:\n",
+ "tvb_delay_fun = \\\n",
+ " lambda source_node, target_node: \\\n",
+ " np.maximum(spiking_network_builder.tvb_dt, \n",
+ " tvb_delay(source_node, target_node, spiking_network_builder.tvb_delays)) \\\n",
+ " if SPIKING_NODES_DELAYS else within_node_delay\n",
+ " \n",
+ " \n",
+ "spiking_network_builder.nodes_connections = [\n",
+ " # source -> target\n",
+ " {\"source\": \"IdSN\", \"target\": \"I\", # \"IdSN\" -> \"Igpi\"\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pdSNGPi), # conn_spec\n",
+ " \"weight\": TVBWeightFun(),\n",
+ " \"delay\": lambda source_node, target_node: tvb_delay_fun(source_node, target_node), \n",
+ " \"receptor_type\": \"gaba\", \n",
+ " \"source_nodes\": Istr_nodes_ids, \n",
+ " \"target_nodes\": Igpi_nodes_ids}, # None means apply to all\n",
+ " {\"source\": \"IiSN\", \"target\": \"I\", # \"IiSN\" -> \"Igpe\"\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(piSNGPe), # conn_spec\n",
+ " \"weight\": TVBWeightFun(),\n",
+ " \"delay\": lambda source_node, target_node: tvb_delay_fun(source_node, target_node), \n",
+ " \"receptor_type\": \"gaba\", \n",
+ " \"source_nodes\": Istr_nodes_ids, \n",
+ " \"target_nodes\": Igpe_nodes_ids}, # None means apply to all\n",
+ " {\"source\": \"I\", \"target\": \"I\", # \"Igpe\" -> \"Igpi\"\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pGPeGPi), # conn_spec\n",
+ " \"weight\": TVBWeightFun(),\n",
+ " \"delay\": lambda source_node, target_node: tvb_delay_fun(source_node, target_node), \n",
+ " \"receptor_type\": \"gaba\", \n",
+ " \"source_nodes\": Igpe_nodes_ids, \n",
+ " \"target_nodes\": Igpi_nodes_ids}, # None means apply to all\n",
+ " {\"source\": \"I\", \"target\": \"E\", # \"Igpi\" -> \"Eth\"\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pGPiTh), # conn_spec\n",
+ " \"weight\": TVBWeightFun(),\n",
+ " \"delay\": lambda source_node, target_node: tvb_delay_fun(source_node, target_node), \n",
+ " \"receptor_type\": \"gaba\", \n",
+ " \"source_nodes\": Igpi_nodes_ids, \n",
+ " \"target_nodes\": Eth_nodes_ids}, # None means apply to all\n",
+ " {\"source\": \"I\", \"target\": \"E\", # \"Igpe\" -> \"Estn\"\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pGPeSTN), # conn_spec\n",
+ " \"weight\": TVBWeightFun(),\n",
+ " \"delay\": lambda source_node, target_node: tvb_delay_fun(source_node, target_node), \n",
+ " \"receptor_type\": \"gaba\", \n",
+ " \"source_nodes\": Igpe_nodes_ids, \n",
+ " \"target_nodes\": Estn_nodes_ids}, # None means apply to all\n",
+ " {\"source\": \"E\", \"target\": \"IdSN\", # \"Eth\" -> [\"IdSN\"] \n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pThdSN), # conn_spec\n",
+ " \"weight\": wThdSN, # TVBWeightFun(),\n",
+ " \"delay\": lambda source_node, target_node: tvb_delay_fun(source_node, target_node), \n",
+ " \"receptor_type\": \"ampa\", \n",
+ " \"source_nodes\": Eth_nodes_ids, \n",
+ " \"target_nodes\": Istr_nodes_ids}, # None means apply to all\n",
+ " {\"source\": \"E\", \"target\": \"IiSN\", # \"Eth\" -> [\"IiSN\"] \n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pThiSN), # conn_spec\n",
+ " \"weight\": wThiSN,\n",
+ " \"delay\": lambda source_node, target_node: tvb_delay_fun(source_node, target_node), \n",
+ " \"receptor_type\": \"ampa\", \n",
+ " \"source_nodes\": Eth_nodes_ids, \n",
+ " \"target_nodes\": Istr_nodes_ids}, # No\n",
+ " {\"source\": \"E\", \"target\": \"I\", # \"Estn\" -> [\"Igpe\"]\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pSTNGPe), # conn_spec\n",
+ " \"weight\": TVBWeightFun(),\n",
+ " \"delay\": lambda source_node, target_node: tvb_delay_fun(source_node, target_node), \n",
+ " \"receptor_type\": \"ampa\", \n",
+ " \"source_nodes\": Estn_nodes_ids, \n",
+ " \"target_nodes\": Igpe_nodes_ids},\n",
+ " {\"source\": \"E\", \"target\": \"I\", # \"Estn\" -> [\"Igpi\"]\n",
+ " \"synapse_model\": synapse_model, \"conn_spec\": conn_spec_fixed_prob(pSTNGPi), # conn_spec\n",
+ " \"weight\": TVBWeightFun(),\n",
+ " \"delay\": lambda source_node, target_node: tvb_delay_fun(source_node, target_node), \n",
+ " \"receptor_type\": \"ampa\", \n",
+ " \"source_nodes\": Estn_nodes_ids, \n",
+ " \"target_nodes\": Igpi_nodes_ids},\n",
+ " ]\n",
+ "\n",
+ "# Creating devices to be able to observe ANNarchy activity:\n",
+ "\n",
+ "spiking_network_builder.output_devices = []\n",
+ "\n",
+ "period = 1.0\n",
+ "\n",
+ "# Creating devices to be able to observe ANNarchy activity:\n",
+ "params = {} # deepcopy(spiking_network_builder.config.ANNARCHY_OUTPUT_DEVICES_PARAMS_DEF[\"SpikeMonitor\"])\n",
+ "for pop in spiking_network_builder.populations:\n",
+ " connections = OrderedDict({})\n",
+ " # label <- target population\n",
+ " connections[pop[\"label\"]] = pop[\"label\"]\n",
+ " spiking_network_builder.output_devices.append(\n",
+ " {\"model\": \"spike_recorder\", \"params\": deepcopy(params),\n",
+ " \"connections\": connections, \"nodes\": pop[\"nodes\"]}) # None means apply to \"all\"\n",
+ "\n",
+ "# Labels have to be different for every connection to every distinct population\n",
+ "# params for baladron implementation commented out for the moment\n",
+ "# TODO: use baladron neurons\n",
+ "# was:\n",
+ "# params = {} # deepcopy(spiking_network_builder.config.ANNARCHY_OUTPUT_DEVICES_PARAMS_DEF[\"Monitor\"])\n",
+ "# params.update({\"period\": period, \n",
+ "# 'variables': [\"v\", \"u\", \"I_syn\", \"I_syn_ex\", \"I_syn_in\", \"g_ampa\", \"g_gaba\", \"g_base\"]})\n",
+ "# for pop in spiking_network_builder.populations:\n",
+ "# connections = OrderedDict({})\n",
+ "# # label <- target population\n",
+ "# connections[pop[\"label\"] + \"_ts\"] = pop[\"label\"]\n",
+ "# spiking_network_builder.output_devices.append(\n",
+ "# {\"model\": \"Monitor\", \"params\": deepcopy(params),\n",
+ "# \"connections\": connections, \"nodes\": pop[\"nodes\"]}) # None means apply to all\n",
+ " \n",
+ " \n",
+ " \n",
+ "# # Create a spike stimulus input device\n",
+ "spiking_network_builder.input_devices = [] #\n",
+ "\n",
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "\n",
+ "# This will be transferred to NetPyNE\n",
+ "config.simulation_length = simulator.simulation_length\n",
+ "\n",
+ "\n",
+ "spiking_network_builder.configure(netParams=netParams, simConfig=cfg)\n",
+ "# confirmed: now spiking_network_builder.global_coupling_scaling = 0.00390625\n",
+ "\n",
+ "\n",
+ "def synaptic_weight_scale_func(is_coupling_mode_tvb):\n",
+ " # TODO: to be tuned or removed?\n",
+ " if is_coupling_mode_tvb: # \"TVB\"\n",
+ " return 1 # 1e-2\n",
+ " else: # \"spikeNet\"\n",
+ " return 1 # 5\n",
+ "\n",
+ "\n",
+ "# spiking_network_builder.global_coupling_scaling *= simulator.model.G\n",
+ "spiking_network_builder.netpyne_synaptic_weight_scale = synaptic_weight_scale_func(is_coupling_mode_tvb=INTERFACE_COUPLING_MODE==\"TVB\")\n",
+ "\n",
+ "\n",
+ "netpyne_network = spiking_network_builder.build() # set_defaults=False\n",
+ "# confirmed: 6 spiking pops of 200 neurs\n",
+ "\n",
+ "for i, (connId, conn) in enumerate(spiking_network_builder.netpyne_instance.netParams.connParams.items()):\n",
+ " print(f\"{i}. {connId}: {conn['weight']} {conn['probability']}\")\n",
+ "# confirmed: 13 connections between pops (weights and probs confirmed)\n",
+ "\n",
+ "# confirmed: 6 spike recorders, as in ANNarchy when \"_ts\" stuff commented out\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# TODO: 6 pops vs 4 in annarchy (but should be okay, because of different way of setting spiking_network_builder.populations)\n",
+ "populations_sizes = []\n",
+ "print(\"Population sizes: \")\n",
+ "for pop in spiking_network_builder.populations:\n",
+ " populations_sizes.append(int(np.round(pop[\"scale\"] * spiking_network_builder.population_order)))\n",
+ " print(\"%s: %d\" % (pop[\"label\"], populations_sizes[-1]))\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "\n",
+ "# The stimuli models:\n",
+ "\n",
+ "if simulation_mode != \"rs\":\n",
+ " annarchy_instance = annarchy_network.annarchy_instance\n",
+ "\n",
+ "\n",
+ " # amplitude in mV or V? Was V in Michmizos paper, but had value 5.\n",
+ " # frequency in Hz and later divided by 1000 because time scale is ms\n",
+ " # duration in ms!\n",
+ " #dt inside annarchy needs to be <=0.01 to realize a 60µs pulse width, but setting it to that leads to errors. \n",
+ " # realistic \"duration\" parameter for DBSInput would be 0.06\n",
+ "\n",
+ " if stim_mode.find(\"mono\") > -1:\n",
+ " DBSInput = annarchy_instance.Neuron(\n",
+ " parameters=\"\"\"\n",
+ " amplitude = 5.0\n",
+ " kappa = 8\n",
+ " frequency = 130\n",
+ " duration = 0.1\n",
+ " \"\"\",\n",
+ " equations=\"\"\"\n",
+ " h1 = if sin(2*pi*(frequency/1000)*t) > 0 : 1 else : 0\n",
+ " h2 = if sin(2*pi*(frequency/1000)*(t+duration)) > 0 : 1 else : 0\n",
+ " r = amplitude*kappa*h1*(1-h2)\n",
+ " \"\"\"\n",
+ " )\n",
+ " elif stim_mode.find(\"bi\") > -1:\n",
+ " DBSInput = annarchy_instance.Neuron(\n",
+ " parameters=\"\"\"\n",
+ " amplitude = 5.0\n",
+ " ampfactor = 10\n",
+ " kappa = 8\n",
+ " frequency = 130\n",
+ " duration = 0.3\n",
+ " \"\"\",\n",
+ " equations=\"\"\"\n",
+ " h1 = if sin(2*pi*(frequency/1000)*t) > 0 : 1 else : 0\n",
+ " h2 = if sin(2*pi*(frequency/1000)*(t-duration)) > 0 : 1 else : 0\n",
+ " h4 = if sin(2*pi*(frequency/1000)*(t-((ampfactor+1)*duration))) > 0 : 1 else : 0\n",
+ " i1 = amplitude*(1+ 1/ampfactor)*kappa*h1*(1-h2)\n",
+ " i2 = (-amplitude/ampfactor) *kappa *h1*(1-h4)\n",
+ " r = i1 + i2\n",
+ " \"\"\"\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "display(Image(filename='pics/Stimuli.png', width=1000, unconfined=False))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Prepare the stimuli:\n",
+ "\n",
+ "if simulation_mode != \"rs\":\n",
+ " \n",
+ " if stim_target.find(\"STN\") > -1:\n",
+ " stim_target = \"STN_Left\" \n",
+ " else:\n",
+ " stim_target = \"GPi_Left\" \n",
+ " \n",
+ " if stim_mode != \"simple\":\n",
+ " # now add DBS population invisible to TVB for now\n",
+ " target_pop = annarchy_network.brain_regions[stim_target][0].population\n",
+ " dbs_pop = annarchy_instance.Population(target_pop.size, DBSInput)\n",
+ " dbs_proj = annarchy_instance.CurrentInjection(dbs_pop, target_pop, 'dbs')\n",
+ " dbs_proj.connect_current()\n",
+ " # switch the stimulation off for the first time segment\n",
+ " dbs_pop.amplitude = 0\n",
+ " dbs_pop.frequency = stim_freq\n",
+ " dbs_pop.duration = stim_duration\n",
+ "\n",
+ " # just to check, add some annarchy monitors -> delete later \n",
+ " m1 = annarchy_instance.Monitor(dbs_pop, \"r\")\n",
+ " m2 = annarchy_instance.Monitor(target_pop, [\"g_dbs\", \"spike\"])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 3. Build the TVB-ANNarchy interface"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "display(Image(filename='pics/TVB-ANNarchy.png', width=1000, unconfined=False))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "display(Image(filename='pics/ANNarchy-TVB.png', width=1000, unconfined=False))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "from tvb_multiscale.tvb_netpyne.interfaces.models.basal_ganglia_izhikevich import BasalGangliaIzhikevichTVBNetpyneInterfaceBuilder\n",
+ "\n",
+ " \n",
+ "# Build a TVB-ANNarchy interface with all the appropriate connections between the\n",
+ "# TVB and ANNarchy modelled regions\n",
+ "tvb_spikeNet_model_builder = BasalGangliaIzhikevichTVBNetpyneInterfaceBuilder()\n",
+ "\n",
+ "tvb_spikeNet_model_builder.config = config\n",
+ "tvb_spikeNet_model_builder.tvb_cosimulator = simulator \n",
+ "tvb_spikeNet_model_builder.spiking_network = netpyne_network\n",
+ "# This can be used to set default tranformer and proxy models:\n",
+ "tvb_spikeNet_model_builder.model = \"RATE\" # \"RATE\" (or \"SPIKES\", \"CURRENT\") TVB->ANNarchy interface\n",
+ "tvb_spikeNet_model_builder.input_flag = True # If True, ANNarchy->TVB update will be implemented\n",
+ "tvb_spikeNet_model_builder.output_flag = True # If True, TVB->ANNarchy coupling will be implemented\n",
+ "# If default_coupling_mode = \"TVB\", large scale coupling towards spiking regions is computed in TVB\n",
+ "# and then applied with no time delay via a single \"TVB proxy node\" / ANNarchy device for each spiking region,\n",
+ "# \"1-to-1\" TVB->ANNarchy coupling.\n",
+ "# If any other value, we need 1 \"TVB proxy node\" / ANNarchy device for each TVB sender region node, and\n",
+ "# large-scale coupling for spiking regions is computed in ANNarchy, \n",
+ "# taking into consideration the TVB connectome weights and delays, \n",
+ "# in this \"1-to-many\" TVB->ANNarchy coupling.\n",
+ "tvb_spikeNet_model_builder.default_coupling_mode = INTERFACE_COUPLING_MODE # \"spikeNet\" # \"TVB\" \n",
+ "# Number of neurons per population to be used to compute population mean instantaneous firing rates:\n",
+ "tvb_spikeNet_model_builder.proxy_inds = np.array(spiking_nodes_ids)\n",
+ "tvb_spikeNet_model_builder.N_E = spiking_network_builder.population_order\n",
+ "tvb_spikeNet_model_builder.GPe_proxy_inds = np.array(Igpe_nodes_ids)\n",
+ "tvb_spikeNet_model_builder.GPi_proxy_inds = np.array(Igpi_nodes_ids)\n",
+ "tvb_spikeNet_model_builder.STN_proxy_inds = np.array(Estn_nodes_ids)\n",
+ "tvb_spikeNet_model_builder.Striatum_proxy_inds = np.array(Istr_nodes_ids)\n",
+ "tvb_spikeNet_model_builder.Thal_proxy_inds = np.array(Eth_nodes_ids)\n",
+ "\n",
+ "# Set exclusive_nodes = True (Default) if the spiking regions substitute for the TVB ones:\n",
+ "tvb_spikeNet_model_builder.exclusive_nodes = True \n",
+ "\n",
+ "tvb_spikeNet_model_builder.output_interfaces = []\n",
+ "tvb_spikeNet_model_builder.input_interfaces = []\n",
+ " \n",
+ "# was:\n",
+ "# # options for a nonopinionated builder:\n",
+ "# from tvb_multiscale.core.interfaces.base.transformers.models.models import Transformers\n",
+ "# from tvb_multiscale.core.interfaces.base.transformers.builders import \\\n",
+ "# DefaultTVBtoSpikeNetTransformers, DefaultSpikeNetToTVBTransformers, \\\n",
+ "# DefaultTVBtoSpikeNetModels, DefaultSpikeNetToTVBModels\n",
+ "# from tvb_multiscale.tvb_annarchy.interfaces.builders import \\\n",
+ "# TVBtoANNarchyModels, ANNarchyInputProxyModels, DefaultTVBtoANNarchyModels, \\\n",
+ "# ANNarchyToTVBModels, ANNarchyOutputProxyModels, DefaultANNarchytoTVBModels\n",
+ "\n",
+ "from tvb_multiscale.tvb_netpyne.interfaces.builders import NetpyneInputProxyModels\n",
+ " # TVBtoNetpyneModels, \n",
+ "\n",
+ "# , DefaultTVBtoNetpyneModels, \\\n",
+ "\n",
+ "\n",
+ " \n",
+ " \n",
+ "# def print_enum(enum):\n",
+ "# print(\"\\n\", enum)\n",
+ "# for name, member in enum.__members__.items():\n",
+ "# print(name,\"= \", member.value)\n",
+ " \n",
+ " \n",
+ "# print(\"Available input (NEST->TVB update) / output (TVB->NEST coupling) interface models:\")\n",
+ "# print_enum(TVBtoANNarchyModels)\n",
+ "# print_enum(ANNarchyToTVBModels)\n",
+ " \n",
+ " \n",
+ "# print(\"\\n\\nAvailable input (spikeNet->TVB update) / output (TVB->spikeNet coupling) transformer models:\")\n",
+ "\n",
+ "# print_enum(DefaultTVBtoSpikeNetModels)\n",
+ "# print_enum(DefaultTVBtoSpikeNetTransformers)\n",
+ " \n",
+ "# print_enum(DefaultSpikeNetToTVBModels)\n",
+ "# print_enum(DefaultSpikeNetToTVBTransformers) \n",
+ " \n",
+ " \n",
+ "# print(\"\\n\\nAvailable input (NEST->TVB update) / output (TVB->NEST coupling) proxy models:\")\n",
+ "\n",
+ "# print_enum(DefaultTVBtoANNarchyModels)\n",
+ "# print_enum(ANNarchyInputProxyModels)\n",
+ " \n",
+ "# print_enum(ANNarchyOutputProxyModels)\n",
+ "# print_enum(DefaultANNarchytoTVBModels)\n",
+ " \n",
+ "# print(\"\\n\\nAll basic transformer models:\")\n",
+ "# print_enum(Transformers)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "\n",
+ "# TVB applies a global coupling scaling of coupling.a * model.G\n",
+ "tvb_spikeNet_model_builder.global_coupling_scaling = \\\n",
+ " tvb_spikeNet_model_builder.tvb_cosimulator.coupling.a[0].item() * tvb_spikeNet_model_builder.G\n",
+ "print(\"global_coupling_scaling = %g\" % tvb_spikeNet_model_builder.global_coupling_scaling)\n",
+ "\n",
+ "# Total TVB indegree weight to STN:\n",
+ "wTVBSTNs = simulator.connectivity.weights[Estn_nodes_ids, 5:].squeeze()\n",
+ "wTVBSTN = wTVBSTNs.sum().item()\n",
+ "print(\"wTVBSTN = %g\" % wTVBSTN)\n",
+ "CTXtoSTNinds = 5 + np.where(wTVBSTNs > 0.0)[0] # indices of TVB regions coupling to STN\n",
+ "\n",
+ "# Total TVB indegree weight to Striatum:\n",
+ "wTVBSNs = simulator.connectivity.weights[Istr_nodes_ids, 5:].squeeze()\n",
+ "wTVBSN = wTVBSNs.sum().item()\n",
+ "print(\"wTVBSN = %g\" % wTVBSN)\n",
+ "CTXtoSNinds = 5 + np.where(wTVBSNs > 0.0)[0] # indices of TVB regions coupling to Striatum\n",
+ "\n",
+ "# Approximate effective scaling of TVB coupling to STN\n",
+ "# after normalizing with STN indegree and \n",
+ "# multiplying with the Maith et al. optimized CTX -> STN weight\n",
+ "iwCtxSTN = STN_factor*wCtxSTN / wTVBSTN\n",
+ "print(\"iwCtxSTN = %g\" % iwCtxSTN)\n",
+ "\n",
+ "# Approximate effective scaling of TVB coupling to dSN\n",
+ "# after normalizing with Striatum indegree and \n",
+ "# multiplying with the Maith et al. optimized CTX -> dSN weight\n",
+ "iwCtxdSN = dSN_factor*wCtxdSN / wTVBSN \n",
+ "print(\"iwCtxdSN = %g\" % iwCtxdSN)\n",
+ "\n",
+ "# Approximate effective scaling of TVB coupling to iSN\n",
+ "# after normalizing with Striatum indegree and \n",
+ "# multiplying with the Maith et al. optimized CTX -> iSN weight\n",
+ "iwCtxiSN = iSN_factor*wCtxiSN / wTVBSN\n",
+ "print(\"iwCtxiSN = %g\" % iwCtxiSN)\n",
+ "\n",
+ "# confirmed:\n",
+ "# global_coupling_scaling = 75.9402\n",
+ "# wTVBSTN = 0.013739\n",
+ "# wTVBSN = 0.0503774\n",
+ "# iwCtxSTN = 0.0530743\n",
+ "# iwCtxdSN = 0.00666963\n",
+ "# iwCtxiSN = 0.0045056\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "# ----Uncomment below to modify the builder by changing the default options:--------------------------------------\n",
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "\n",
+ "from tvb_multiscale.core.interfaces.base.transformers.models.red_wong_wang import RedWongWangExc\n",
+ "\n",
+ "\n",
+ "# --------For spike transmission from TVB to ANNarchy devices acting as TVB proxy nodes with TVB delays:--------\n",
+ "\n",
+ "\n",
+ "# TVB -> ANNarchy\n",
+ "\n",
+ "\n",
+ "if tvb_spikeNet_model_builder.default_coupling_mode == \"spikeNet\":\n",
+ " \n",
+ " # If coupling is computing in ANNarchy, we need as many TVB proxies \n",
+ " # as TVB regions coupling to STN and Striatum\n",
+ " proxy_inds = np.unique(np.concatenate([CTXtoSTNinds, CTXtoSNinds]))\n",
+ " \n",
+ " # This is the TVB coupling weight function that will determine the connections' weights \n",
+ " # from TVB proxies to the target STN and dSN/iSN populations:\n",
+ " class TVBWeightFunInterface(object):\n",
+ " \n",
+ " def __init__(self, scaling):\n",
+ " self.scaling = float(scaling)\n",
+ "\n",
+ " def __call__(self, source_node, target_node, tvb_weights):\n",
+ " return (scale_tvb_weight(source_node, target_node, tvb_weights, scale=self.scaling))\n",
+ "\n",
+ " # A similar function for TVB coupling delays is also applied in the background \n",
+ " # without need to be explicitly defined by the user\n",
+ " \n",
+ "# Optionally adjust interface scale factors here \n",
+ "# to have the same result counter act possible changes to G and coupling.a:\n",
+ "# STN_factor /= tvb_spikeNet_model_builder.global_coupling_scaling\n",
+ "# dSN_factor /= tvb_spikeNet_model_builder.global_coupling_scaling\n",
+ "# iSN_factor /= tvb_spikeNet_model_builder.global_coupling_scaling\n",
+ "\n",
+ "dSN_factor *= 40\n",
+ "iSN_factor *= 40\n",
+ "\n",
+ "# \n",
+ "tvb_spikeNet_model_builder.synaptic_weight_scale_func = synaptic_weight_scale_func\n",
+ " \n",
+ "tvb_spikeNet_model_builder.output_interfaces = []\n",
+ "# Mean spike rates are applied in parallel to all target neurons\n",
+ "for trg_pop, target_nodes, conn_scaling, this_conn_spec, scale_factor in \\\n",
+ " zip([\"E\", \"IdSN\", \"IiSN\"], # ANNarchy target populations\n",
+ " # Target region indices in ANNarchy: \n",
+ " [tvb_spikeNet_model_builder.STN_proxy_inds, \n",
+ " tvb_spikeNet_model_builder.Striatum_proxy_inds, \n",
+ " tvb_spikeNet_model_builder.Striatum_proxy_inds], \n",
+ " # Maith et al optimized... \n",
+ " [wCtxSTN, wCtxdSN, wCtxiSN], # ...weights \n",
+ " # ...and probabilities for CTX -> STN/Striatum connections\n",
+ " [conn_spec_fixed_prob(prob=pCtxSTN), # pCtxSTN \n",
+ " conn_spec_fixed_prob(prob=pCtxdSN), # pCtxdSN\n",
+ " conn_spec_fixed_prob(prob=pCtxiSN)], # pCtxiSN\n",
+ " # Interface scaling factors scaled by TVB weights' indegree to STN/Striatum:\n",
+ " [STN_factor/wTVBSTN, dSN_factor/wTVBSN, iSN_factor/wTVBSN]): \n",
+ " tvb_spikeNet_model_builder.output_interfaces.append(\n",
+ " {\"voi\": np.array([\"R\"]), # Source TVB state variable\n",
+ " \"populations\": np.array([trg_pop]), # ANNarchy target population\n",
+ " \"model\": \"RATE\",\n",
+ " \"spiking_proxy_inds\": target_nodes, # Target region indices in ANNarchy\n",
+ " # This spike generator device generates spike trains \n",
+ " # with autocorrelation corr at a time scale tau\n",
+ " \"proxy_model\": NetpyneInputProxyModels.RATE, # ANNarchyInputProxyModels.RATE_TO_SPIKES, # \n",
+ " \n",
+ " # TODO: above needed?\n",
+ " # 'transformer_params': {'scale_factor': np.array([1.0])}, # due to the way Netpyne generates spikes, no scaling by population size is needed\n",
+ "\n",
+ " \"proxy_params\": {\"geometry\": 600, \"record\": [\"spike\"],\n",
+ " \"corr\": 0.3, \"tau\": 10.0, # comment for RATE_TO_SPIKES\n",
+ " \"number_of_neurons\": tvb_spikeNet_model_builder.N_E, # todo: de-hardcode\n",
+ " },\n",
+ " 'conn_spec': this_conn_spec, # dictionary of connection properties\n",
+ " 'coupling_mode': tvb_spikeNet_model_builder.default_coupling_mode\n",
+ " }) # None means all here\n",
+ " \n",
+ " # For both coupling modes, we scale the TVB rate already at the TVB -> ANNarchy transformer\n",
+ " # with the interface scale factor (normalized by TVB indegree to STN/Striatum)\n",
+ " # and the global coupling scaling.\n",
+ " if tvb_spikeNet_model_builder.output_interfaces[-1][\"coupling_mode\"] == \"spikeNet\":\n",
+ "# if trg_pop == \"E\":\n",
+ "# tvb_spikeNet_model_builder.output_interfaces[-1][\"proxy_inds\"] = CTXtoSTNinds\n",
+ "# else:\n",
+ "# tvb_spikeNet_model_builder.output_interfaces[-1][\"proxy_inds\"] = CTXtoSNinds\n",
+ " tvb_spikeNet_model_builder.output_interfaces[-1][\"proxy_inds\"] = proxy_inds\n",
+ " # In this case connections from each TVB proxy to STN/Striatum \n",
+ " # are scaled additionally with the Maith et al. optimized weights\n",
+ " tvb_spikeNet_model_builder.output_interfaces[-1][\"weights\"] = TVBWeightFunInterface(conn_scaling)\n",
+ " tvb_spikeNet_model_builder.output_interfaces[-1][\"transformer_params\"] = \\\n",
+ " {\"scale_factor\": scale_factor * tvb_spikeNet_model_builder.global_coupling_scaling}\n",
+ " # In total:\n",
+ " # From each TVB proxy node we get a rate scaled as (\n",
+ " # (coupling.a * G * STN_factor/wTVBSTN) * R_i, (i for all TVB regions)\n",
+ " # Then, spikes generated from each TVB proxy are transferred via connections \n",
+ " # with weights TVB_w_ji * wCtxSTN or wCtxiSN or wCtxdSN (j for STN or Striatum) \n",
+ " # and probabilities pCtxSTN or pCtxiSN or pCtxdSN, respectively\n",
+ " else:\n",
+ " # In this case connections from each TVB proxy to STN/Striatum \n",
+ " # are equal to the Maith et al. optimized weights\n",
+ " tvb_spikeNet_model_builder.output_interfaces[-1][\"weights\"] = conn_scaling\n",
+ " # In this case coupling.a is already applied during computing TVB coupling.\n",
+ " # Therefore we scale only with model.G\n",
+ " tvb_spikeNet_model_builder.output_interfaces[-1][\"transformer_params\"] = \\\n",
+ " {\"scale_factor\": scale_factor * tvb_spikeNet_model_builder.G}\n",
+ " # In total:\n",
+ " # From each TVB proxy node we get a total coupling rate scaled \n",
+ " # as (coupling.a * G STN_factor/wTVBSTN) * R_j, (j for STN or Striatum)\n",
+ " # Then, spikes generated from each TVB proxy are transferred via connections \n",
+ " # with weights wCtxSTN or wCtxiSN or wCtxdSN and \n",
+ " # probabilities pCtxSTN or pCtxiSN or pCtxdSN, respectively\n",
+ " \n",
+ " \n",
+ "from tvb_multiscale.core.interfaces.base.transformers.models.red_wong_wang import ElephantSpikesRateRedWongWangExc\n",
+ "\n",
+ "tvb_spikeNet_model_builder.input_interfaces = []\n",
+ "# TVB <-- ANNarchy:\n",
+ "for src_pop, nodes, in zip(\n",
+ " # Source populations in ANNarchy:\n",
+ " [np.array([\"I\"]), np.array([\"E\"]), np.array([\"IdSN\", \"IiSN\"])],\n",
+ " # Source regions indices in ANNarchy:\n",
+ " [tvb_spikeNet_model_builder.I_proxy_inds, # GPe and GPi\n",
+ " tvb_spikeNet_model_builder.E_proxy_inds, # STN and Thalamus\n",
+ " tvb_spikeNet_model_builder.Striatum_proxy_inds]): # Striatum\n",
+ " # TVB <- ANNarchy\n",
+ " tvb_spikeNet_model_builder.input_interfaces.append(\n",
+ " {\"voi\": np.array([\"S\", \"R\"]), # Target state variables in TVB\n",
+ " \"populations\": src_pop, # Source populations in ANNarchy\n",
+ " # This transformer not only converts spike counts to rates for state variable R,\n",
+ " # but also integrates the dS/dt to compute the respective S!:\n",
+ " \"transformer\": ElephantSpikesRateRedWongWangExc,\n",
+ " \"transformer_params\": \n",
+ " # Spike counts are converted to rates via:\n",
+ " # number_of_spikes / number_of_neurons_per_population / number_of_populations\n",
+ " # (mind that there are 2 populations in Striatum)\n",
+ " {\"scale_factor\": np.array([1.0]) / tvb_spikeNet_model_builder.N_E / len(src_pop),\n",
+ " # The integrator used to integrate dS/dt\n",
+ " \"integrator\":CONFIGURED.DEFAULT_TRANSFORMER_INTEGRATOR_MODEL(\n",
+ " dt=simulator.integrator.dt),\n",
+ " \"state\": np.zeros((2, len(nodes))), # initial condition\n",
+ " # Parameters of the dS/dt differential equation:\n",
+ " \"tau_s\": simulator.model.tau_s, # time constant of integrating S\n",
+ " \"tau_r\": np.array([10.0]), # time constant of integrating R to low pass filter it\n",
+ " \"gamma\": simulator.model.gamma}, \n",
+ " \"proxy_inds\": np.array(nodes)}) # None means all here\n",
+ " \n",
+ "\n",
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "# ----------------------------------------------------------------------------------------------------------------\n",
+ "\n",
+ "# Configure and build:\n",
+ "tvb_spikeNet_model_builder.configure()\n",
+ "# tvb_spikeNet_model_builder.print_summary_info_details(recursive=1)\n",
+ " \n",
+ "print(\"\\noutput (TVB->NEST coupling) interfaces' configurations:\\n\")\n",
+ "display(tvb_spikeNet_model_builder.output_interfaces)\n",
+ " \n",
+ "print(\"\\ninput (NEST->TVB update) interfaces' configurations:\\n\")\n",
+ "display(tvb_spikeNet_model_builder.input_interfaces)\n",
+ " \n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "\n",
+ "simulator = tvb_spikeNet_model_builder.build()\n",
+ "\n",
+ "# NetPyNE model is built in two steps. First need to create declarative-style specification for both spiking network itself and TVB-Netpyne proxy devides (interfaces).\n",
+ "# Once it's done above using builders, network can be instantiated based on the specification:\n",
+ "\n",
+ "netpyne.simConfig.duration = simulation_length # TODO: do it properly! (find out why it gets lost)\n",
+ "netpyne.simConfig.dt = spiking_network_builder.spiking_dt\n",
+ "netpyne.allowSelfConns = True\n",
+ "\n",
+ "netpyne_network.netpyne_instance.instantiateNetwork()\n",
+ "\n",
+ "\n",
+ "simulator.simulate_spiking_simulator = netpyne_network.netpyne_instance.run # set the method to run NetPyNE\n",
+ " \n",
+ "# simulator.print_summary_info(recursive=3)\n",
+ "# simulator.print_summary_info_details(recursive=3)\n",
+ " \n",
+ "print(\"\\n\\noutput (TVB->NetPyNE coupling) interfaces:\\n\")\n",
+ "simulator.output_interfaces.print_summary_info_details(recursive=2)\n",
+ " \n",
+ "print(\"\\n\\ninput (NetPyNE->TVB update) interfaces:\\n\")\n",
+ "simulator.input_interfaces.print_summary_info_details(recursive=2)\n"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 3. Configure simulator, simulate, gather results"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-07-12T20:36:18.879872Z",
+ "start_time": "2019-07-12T20:36:11.148945Z"
+ },
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ },
+ "scrolled": true,
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "# -----------------------------------4. Compile network ---------------------------------------------------------\n",
+ "# Compile the ANNarchy network...\n",
+ "# tic_compile = time.time()\n",
+ "# netpyne_network.configure()\n",
+ "# print(\"Compiled! in %f min\" % ((time.time() - tic_compile) / 60))\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "netpyne_network.print_summary_info_details(recursive=2, connectivity=False)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ " # -----------------------------------5. Simulate and gather results-------------------------------------------------\n",
+ "# simulation_length = 1500.0\n",
+ "# transient = 500.0 # simulation_length/11\n",
+ "# ...and simulate!\n",
+ "\n",
+ "simulator.configure() # \n",
+ "\n",
+ "simulation_length = \\\n",
+ " np.ceil(simulation_length / simulator.synchronization_time) * simulator.synchronization_time\n",
+ "\n",
+ "advance_simulation_for_delayed_monitors_output = True\n",
+ "if simulation_mode == \"rs\":\n",
+ " simulation_length1 = simulation_length\n",
+ "else:\n",
+ " start_stimulus = np.ceil(start_stimulus / simulator.synchronization_time) * simulator.synchronization_time\n",
+ " simulation_length1 = start_stimulus\n",
+ " advance_simulation_for_delayed_monitors_output = False\n",
+ "\n",
+ "t_start = time.time()\n",
+ "\n",
+ "results = simulator.run(simulation_length=simulation_length1, \n",
+ " advance_simulation_for_delayed_monitors_output=advance_simulation_for_delayed_monitors_output\n",
+ " ) # 35.0 with stimulus application\n",
+ "\n",
+ "if simulation_mode != \"rs\":\n",
+ " \n",
+ " if simulation_mode.find(\"simpl\") > -1:\n",
+ " # for stimulus application:\n",
+ " if stim_target.find(\"STN\") > -1:\n",
+ " pop = \"E\" \n",
+ " else:\n",
+ " pop = \"I\" \n",
+ " annarchy_network.brain_regions[stim_target][pop].Set(\n",
+ " {\"I\": stim_ampl + annarchy_network.brain_regions[stim_target][pop].Get(\"I\")[\"I\"]})\n",
+ "\n",
+ " else:\n",
+ " \n",
+ " # for stimulus application:\n",
+ " dbs_pop.amplitude = stim_ampl\n",
+ " \n",
+ " results2 = simulator.run(simulation_length=simulation_length-simulation_length1, \n",
+ " advance_simulation_for_delayed_monitors_output=True)\n",
+ " \n",
+ " results[0] = list(results[0])\n",
+ " results[0][0] = np.concatenate([results[0][0], results2[0][0]], axis=0)\n",
+ " results[0][1] = np.concatenate([results[0][1], results2[0][1]], axis=0)\n",
+ " \n",
+ "print(\"\\nSimulated in %f secs!\" % (time.time() - t_start))\n",
+ "\n",
+ "netpyne.finalize()\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from netpyne import sim\n",
+ "popIds = [id for id in sim.net.pops.keys()][:6]\n",
+ "sa = sim.analysis\n",
+ "%matplotlib inline\n",
+ "sim.analysis.plotConn(showFig=True, includePre=popIds, includePost=popIds, feature='weight');\n",
+ "sim.analysis.plotConn(showFig=True, includePre=popIds, includePost=popIds, feature='prob');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sim.analysis.plotConn(showFig=True, includePre=popIds, includePost=popIds, feature='probability');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 4. Plot results and write them to HDF5 files"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# set to False for faster plotting of only mean field variables and dates, apart from spikes\" rasters:\n",
+ "from scipy.io import savemat\n",
+ "plot_per_neuron = False \n",
+ "MAX_VARS_IN_COLS = 3\n",
+ "MAX_REGIONS_IN_ROWS = 10\n",
+ "MIN_REGIONS_FOR_RASTER_PLOT = 9\n"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### TVB plots"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# If you want to see what the function above does, take the steps, one by one\n",
+ "try:\n",
+ " # We need framework_tvb for writing and reading from HDF5 files\n",
+ " from tvb_multiscale.core.tvb.io.h5_writer import H5Writer\n",
+ " from examples.plot_write_results import write_RegionTimeSeriesXarray_to_h5\n",
+ " writer = H5Writer()\n",
+ "except:\n",
+ " writer = False\n",
+ " \n",
+ "# Put the results in a Timeseries instance\n",
+ "from tvb.contrib.scripts.datatypes.time_series_xarray import TimeSeriesRegion as TimeSeriesXarray\n",
+ "\n",
+ "source_ts = TimeSeriesXarray( # substitute with TimeSeriesRegion fot TVB like functionality\n",
+ " data=results[0][1], time=results[0][0] - results[0][0][0],\n",
+ " connectivity=simulator.connectivity,\n",
+ " labels_ordering=[\"Time\", \"State Variable\", \"Region\", \"Neurons\"],\n",
+ " labels_dimensions={\"State Variable\": list(simulator.model.variables_of_interest),\n",
+ " \"Region\": simulator.connectivity.region_labels.tolist()},\n",
+ " sample_period=simulator.integrator.dt)\n",
+ "source_ts.configure()\n",
+ "\n",
+ "t = source_ts.time\n",
+ "\n",
+ "datadict = source_ts._data.to_dict()\n",
+ "datadict[\"transient\"] = transient\n",
+ "savemat(os.path.join(config.out.FOLDER_RES, \"tvb_timeseries.mat\"), datadict)\n",
+ "del datadict\n",
+ "\n",
+ "# Write to file\n",
+ "if writer:\n",
+ " write_RegionTimeSeriesXarray_to_h5(source_ts, writer,\n",
+ " os.path.join(config.out.FOLDER_RES, source_ts.title)+\".h5\")\n",
+ " \n",
+ "source_ts\n",
+ "\n",
+ "\n",
+ "# init_cond = source_ts[-100:]._data.values.mean(axis=0)\n",
+ "# init_cond[2:] = 0.0\n",
+ "# init_cond[1:, 5:] = 0.0\n",
+ "# np.save(init_cond_filepath, init_cond)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "TVBrate = source_ts[transient:][:, \"R\", 5:]._data.values.mean(axis=0).squeeze()\n",
+ "TVBrateSTNeff = TVBrate[wTVBSTNs > 0.0].mean()\n",
+ "TVBrateSNeff = TVBrate[wTVBSNs > 0.0].mean()\n",
+ "TVBrateSTN = (TVBrate * wTVBSTNs * STN_factor).sum()\n",
+ "TVBratedSN = (TVBrate * wTVBSNs * dSN_factor).sum()\n",
+ "TVBrateiSN = (TVBrate * wTVBSNs * iSN_factor).sum() \n",
+ "mean_tvb_rate = TVBrate.mean()\n",
+ "total_tvb_rate = TVBrate.sum()\n",
+ "print(\"Mean TVB rate = %g\" % mean_tvb_rate)\n",
+ "print(\"Total TVB rate = %g\" % total_tvb_rate)\n",
+ "print(\"Total weighted TVB rate to STN = %g ~= %g (expected) ~= %g (Maith et al)\" % \n",
+ " (TVBrateSTN, TVBrateSTNeff*STN_factor*wTVBSTN, TVBrateSTNeff*STN_factor*wCtxSTN))\n",
+ "print(\"Total weighted TVB rate to dSN = %g ~= %g (expected) ~= %g (Maith et al)\" \n",
+ " % (TVBratedSN, TVBrateSNeff*dSN_factor* wTVBSN, TVBrateSNeff*dSN_factor*wCtxdSN))\n",
+ "print(\"Total weighted TVB rate to iSN = %g ~= %g (expected) ~= %g (Maith et al)\" \n",
+ " % (TVBrateiSN, TVBrateSNeff*iSN_factor*wTVBSN, TVBrateSNeff*iSN_factor*wCtxiSN))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# Plot TVB time series\n",
+ "source_ts[:, 0].plot_timeseries(plotter_config=plotter.config, \n",
+ " hue=\"Region\" if source_ts.shape[2] > MAX_REGIONS_IN_ROWS else None, \n",
+ " per_variable=source_ts.shape[1] > MAX_VARS_IN_COLS, \n",
+ " figsize=FIGSIZE);\n",
+ "\n",
+ "source_ts[:, 1].plot_timeseries(plotter_config=plotter.config, \n",
+ " hue=\"Region\" if source_ts.shape[2] > MAX_REGIONS_IN_ROWS else None, \n",
+ " per_variable=source_ts.shape[1] > MAX_VARS_IN_COLS, \n",
+ " figsize=FIGSIZE);"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# # TVB time series raster plot:\n",
+ "# if source_ts.number_of_labels > MIN_REGIONS_FOR_RASTER_PLOT:\n",
+ "# source_ts.plot_raster(plotter_config=plotter.config, \n",
+ "# per_variable=source_ts.shape[1] > MAX_VARS_IN_COLS,\n",
+ "# figsize=FIGSIZE);"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# Focus on the nodes modelled in ANNarchy: \n",
+ "n_spiking_nodes = len(spiking_nodes_ids)\n",
+ "source_ts_ann = source_ts[:, :, spiking_nodes_ids]\n",
+ "source_ts_ann[:, 0].plot_timeseries(plotter_config=plotter.config, \n",
+ " hue=\"Region\" if source_ts_ann.shape[2] > MAX_REGIONS_IN_ROWS else None, \n",
+ " per_variable=source_ts_ann.shape[1] > MAX_VARS_IN_COLS, \n",
+ " figsize=FIGSIZE, figname=\"Spiking nodes TVB Time Series\");\n",
+ "\n",
+ "source_ts_ann[:, 1].plot_timeseries(plotter_config=plotter.config, \n",
+ " hue=\"Region\" if source_ts_ann.shape[2] > MAX_REGIONS_IN_ROWS else None, \n",
+ " per_variable=source_ts_ann.shape[1] > MAX_VARS_IN_COLS, \n",
+ " figsize=FIGSIZE, figname=\"Spiking nodes TVB Time Series\");"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# # Focus on the nodes modelled in ANNarchy: raster plot\n",
+ "# if source_ts_ann.number_of_labels > MIN_REGIONS_FOR_RASTER_PLOT:\n",
+ "# source_ts_ann.plot_raster(plotter_config=plotter.config, \n",
+ "# per_variable=source_ts_ann.shape[1] > MAX_VARS_IN_COLS,\n",
+ "# figsize=FIGSIZE, figname=\"Spiking nodes TVB Time Series Raster\");"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Interactive time series plot"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# # ...interactively as well\n",
+ "# # For interactive plotting:\n",
+ "# %matplotlib notebook \n",
+ "# plotter.plot_timeseries_interactive(TimeSeriesRegion().from_xarray_DataArray(source_ts._data, \n",
+ "# connectivity=source_ts.connectivity))"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Spiking Network plots"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from tvb_multiscale.core.data_analysis.spiking_network_analyser import SpikingNetworkAnalyser\n",
+ "# Create a SpikingNetworkAnalyzer:\n",
+ "spikeNet_analyzer = \\\n",
+ " SpikingNetworkAnalyser(spikeNet=annarchy_network,\n",
+ " start_time=period, end_time=simulation_length, \n",
+ " period=period, transient=transient, # transient,\n",
+ " time_series_output_type=\"TVB\", return_data=True,\n",
+ " # elephant_mean_firing_rate=False,\n",
+ " force_homogeneous_results=True, connectivity=simulator.connectivity)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Plot spikes' raster and mean spike rates and correlations"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# Spikes rates and correlations per Population and Region\n",
+ "spikes_res = \\\n",
+ " spikeNet_analyzer.\\\n",
+ " compute_spikeNet_spikes_rates_and_correlations(\n",
+ " populations_devices=None, regions=None,\n",
+ " rates_methods=[], rates_kwargs=[{}],rate_results_names=[],\n",
+ " corrs_methods=[], corrs_kwargs=[{}], corrs_results_names=[], bin_kwargs={},\n",
+ " data_method=spikeNet_analyzer.get_spikes_from_device, data_kwargs={},\n",
+ " return_devices=False\n",
+ " );"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "if spikes_res:\n",
+ " print(spikes_res[\"mean_rate\"])\n",
+ " print(spikes_res[\"spikes_correlation_coefficient\"])\n",
+ " # Plot spikes' rasters together with mean population's spikes' rates' time series\n",
+ " if plotter:\n",
+ " plotter.config.FONTSIZE = 20 # plotter.config.LARGE_FONTSIZE # LARGE = 12, default = 10\n",
+ " plotter.plot_spike_events(spikes_res[\"spikes\"], time_series=spikes_res[\"mean_rate_time_series\"], \n",
+ " mean_results=spikes_res[\"mean_rate\"], #\n",
+ " stimulus=[start_stimulus] if simulation_mode != \"rs\" else None,\n",
+ " mean_results_units=\"Hz\", stimulus_linewidth=5.0,\n",
+ " spikes_markersize=1.0, figsize=(20, 40), \n",
+ " n_y_ticks=3, n_time_ticks=4, show_time_axis=True, \n",
+ " time_axis_min=0.0, time_axis_max=simulation_length\n",
+ " ) # \n",
+ " from tvb_multiscale.core.plot.correlations_plot import plot_correlations\n",
+ " plot_correlations(spikes_res[\"spikes_correlation_coefficient\"], plotter)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "if spikes_res:\n",
+ " print(\"Mean spike rates:\")\n",
+ " for pop in spikes_res[\"mean_rate\"].coords[\"Population\"]:\n",
+ " for reg in spikes_res[\"mean_rate\"].coords[\"Region\"]:\n",
+ " if not np.isnan(spikes_res[\"mean_rate\"].loc[pop, reg]):\n",
+ " print(\"%s - %s: %0.1f\" % (pop.values.item().split(\"_spikes\")[0], reg.values.item(), \n",
+ " spikes_res[\"mean_rate\"].loc[pop, reg].values.item()))\n",
+ "\n",
+ " savemat(os.path.join(config.out.FOLDER_RES, \"spikes_mean_rates.mat\"), spikes_res[\"mean_rate\"].to_dict())\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "spikeNet_analyzer.resample = True\n",
+ "spikes_sync = None\n",
+ "spikes_sync = \\\n",
+ " spikeNet_analyzer.compute_spikeNet_synchronization(populations_devices=None, regions=None,\n",
+ " comp_methods=[spikeNet_analyzer.compute_spikes_sync, \n",
+ " spikeNet_analyzer.compute_spikes_sync_time_series, \n",
+ " spikeNet_analyzer.compute_spikes_distance, \n",
+ " spikeNet_analyzer.compute_spikes_distance_time_series,\n",
+ " spikeNet_analyzer.compute_spikes_isi_distance, \n",
+ " spikeNet_analyzer.compute_spikes_isi_distance_time_series],\n",
+ " computations_kwargs=[{}], data_kwargs={},\n",
+ " return_spikes_trains=False, return_devices=False)\n",
+ "# print(spikes_sync)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "if spikes_sync:\n",
+ " plotter.config.FONTSIZE = 20 # plotter.config.LARGE_FONTSIZE # LARGE = 12, default = 10\n",
+ " plotter.plot_spike_events(spikes_res[\"spikes\"], \n",
+ " time_series=spikes_sync[\"spikes_sync_time_series\"], \n",
+ " mean_results=spikes_sync[\"spikes_sync\"], \n",
+ " stimulus=[start_stimulus] if simulation_mode != \"rs\" else None,\n",
+ " plot_spikes=True, spikes_alpha=0.25,\n",
+ " spikes_markersize=1.0, stimulus_linewidth=5.0, time_series_marker=\"*\", \n",
+ " figsize=(20, 40), n_y_ticks=3, n_time_ticks=4, show_time_axis=True,\n",
+ " time_axis_min=0.0, time_axis_max=simulation_length\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "if spikes_sync:\n",
+ " plotter.config.FONTSIZE = 20 # plotter.config.LARGE_FONTSIZE # LARGE = 12, default = 10\n",
+ " plotter.plot_spike_events(spikes_res[\"spikes\"], \n",
+ " time_series=spikes_sync[\"spikes_distance_time_series\"], \n",
+ " mean_results=spikes_sync[\"spikes_distance\"], \n",
+ " stimulus=[start_stimulus] if simulation_mode != \"rs\" else None,\n",
+ " plot_spikes=True, spikes_alpha=0.25,\n",
+ " spikes_markersize=1.0, stimulus_linewidth=5.0, time_series_marker=\"*\", \n",
+ " figsize=(20, 40), n_time_ticks=4, show_time_axis=True, n_y_ticks=4,\n",
+ " time_axis_min=0.0, time_axis_max=simulation_length\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "if spikes_sync:\n",
+ " plotter.config.FONTSIZE = 20 # plotter.config.LARGE_FONTSIZE # LARGE = 12, default = 10\n",
+ " plotter.plot_spike_events(spikes_res[\"spikes\"], \n",
+ " time_series=spikes_sync[\"spikes_isi_distance_time_series\"], \n",
+ " mean_results=spikes_sync[\"spikes_isi_distance\"], \n",
+ " stimulus=[start_stimulus] if simulation_mode != \"rs\" else None,\n",
+ " plot_spikes=True, spikes_alpha=0.25,\n",
+ " spikes_markersize=1.0, stimulus_linewidth=5.0, time_series_marker=\"*\", \n",
+ " figsize=(20, 40), n_y_ticks=3, n_time_ticks=4, show_time_axis=True,\n",
+ " time_axis_min=0.0, time_axis_max=simulation_length\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "if spikes_sync:\n",
+ " print(\"Spike synchronization:\")\n",
+ " for pop in spikes_sync[\"spikes_sync\"].coords[\"Population\"]:\n",
+ " for reg in spikes_sync[\"spikes_sync\"].coords[\"Region\"]:\n",
+ " if not np.isnan(spikes_sync[\"spikes_sync\"].loc[pop, reg]):\n",
+ " print(\"%s - %s: %g\" % (pop.values.item().split(\"_spikes\")[0], reg.values.item(), \n",
+ " spikes_sync[\"spikes_sync\"].loc[pop, reg].values.item()))\n",
+ "\n",
+ " savemat(os.path.join(config.out.FOLDER_RES, \"spikes_sync.mat\"), spikes_sync[\"spikes_sync\"].to_dict())\n",
+ " savemat(os.path.join(config.out.FOLDER_RES, \"spikes_sync_time_series.mat\"), spikes_sync[\"spikes_sync_time_series\"].to_dict())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "if spikes_sync:\n",
+ " print(\"Spike distance:\")\n",
+ " for pop in spikes_sync[\"spikes_distance\"].coords[\"Population\"]:\n",
+ " for reg in spikes_sync[\"spikes_distance\"].coords[\"Region\"]:\n",
+ " if not np.isnan(spikes_sync[\"spikes_distance\"].loc[pop, reg]):\n",
+ " print(\"%s - %s: %g\" % (pop.values.item().split(\"_spikes\")[0], reg.values.item(), \n",
+ " spikes_sync[\"spikes_distance\"].loc[pop, reg].values.item()))\n",
+ "\n",
+ " savemat(os.path.join(config.out.FOLDER_RES, \"spikes_distance.mat\"), spikes_sync[\"spikes_distance\"].to_dict())\n",
+ " savemat(os.path.join(config.out.FOLDER_RES, \"spikes_distance_time_series.mat\"), spikes_sync[\"spikes_distance_time_series\"].to_dict())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "if spikes_sync:\n",
+ " print(\"Spike ISI distance:\")\n",
+ " for pop in spikes_sync[\"spikes_isi_distance\"].coords[\"Population\"]:\n",
+ " for reg in spikes_sync[\"spikes_isi_distance\"].coords[\"Region\"]:\n",
+ " if not np.isnan(spikes_sync[\"spikes_isi_distance\"].loc[pop, reg]):\n",
+ " print(\"%s - %s: %g\" % (pop.values.item().split(\"_spikes\")[0], reg.values.item(), \n",
+ " spikes_sync[\"spikes_isi_distance\"].loc[pop, reg].values.item()))\n",
+ "\n",
+ " savemat(os.path.join(config.out.FOLDER_RES, \"spikes_isi_distance.mat\"), spikes_sync[\"spikes_isi_distance\"].to_dict())\n",
+ " savemat(os.path.join(config.out.FOLDER_RES, \"spikes_isi_distance_time_series.mat\"), spikes_sync[\"spikes_isi_distance_time_series\"].to_dict())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "if spikes_res and writer:\n",
+ " writer.write_object(spikes_res[\"spikes\"].to_dict(), \n",
+ " path=os.path.join(config.out.FOLDER_RES, \"Spikes\") + \".h5\");\n",
+ " writer.write_object(spikes_res[\"mean_rate\"].to_dict(),\n",
+ " path=os.path.join(config.out.FOLDER_RES,\n",
+ " spikes_res[\"mean_rate\"].name) + \".h5\");\n",
+ " write_RegionTimeSeriesXarray_to_h5(spikes_res[\"mean_rate_time_series\"], writer,\n",
+ " os.path.join(config.out.FOLDER_RES,\n",
+ " spikes_res[\"mean_rate_time_series\"].title) + \".h5\",\n",
+ " recursive=False);\n",
+ " writer.write_object(spikes_res[\"spikes_correlation_coefficient\"].to_dict(),\n",
+ " path=os.path.join(config.out.FOLDER_RES,\n",
+ " spikes_res[\"spikes_correlation_coefficient\"].name) + \".h5\");"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Get SpikingNetwork mean field variable time series and plot them"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# Continuous time variables' data of spiking neurons\n",
+ "if plot_per_neuron:\n",
+ " spikeNet_analyzer.return_data = True\n",
+ "else:\n",
+ " spikeNet_analyzer.return_data = False\n",
+ "spikeNet_ts = \\\n",
+ " spikeNet_analyzer. \\\n",
+ " compute_spikeNet_mean_field_time_series(populations_devices=None, regions=None, variables=None,\n",
+ " computations_kwargs={}, data_kwargs={}, return_devices=False)\n",
+ "if spikeNet_ts:\n",
+ " if plot_per_neuron:\n",
+ " mean_field_ts = spikeNet_ts[\"mean_field_time_series\"] # mean field\n",
+ " spikeNet_ts = spikeNet_ts[\"data_by_neuron\"] # per neuron data\n",
+ " else:\n",
+ " mean_field_ts = spikeNet_ts\n",
+ " if mean_field_ts and mean_field_ts.size > 0:\n",
+ " mean_field_ts.plot_timeseries(plotter_config=plotter.config, \n",
+ " per_variable=mean_field_ts.shape[1] > MAX_VARS_IN_COLS)\n",
+ " if mean_field_ts.number_of_labels > MIN_REGIONS_FOR_RASTER_PLOT:\n",
+ " mean_field_ts.plot_raster(plotter_config=plotter.config, \n",
+ " per_variable=mean_field_ts.shape[1] > MAX_VARS_IN_COLS,\n",
+ " linestyle=\"--\", alpha=0.5, linewidth=0.5)\n",
+ "else:\n",
+ " mean_field_ts = None\n",
+ "\n",
+ "datadict = mean_field_ts._data.to_dict() \n",
+ "datadict[\"transient\"] = transient\n",
+ "savemat(os.path.join(config.out.FOLDER_RES, \"ANNarchy_timeseries.mat\"), datadict)\n",
+ "del datadict"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Write results to file:\n",
+ "if mean_field_ts and writer:\n",
+ " write_RegionTimeSeriesXarray_to_h5(mean_field_ts, writer,\n",
+ " os.path.join(config.out.FOLDER_RES, mean_field_ts.title) + \".h5\", \n",
+ " recursive=False)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Compute per neuron spikes' rates times series and plot them"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "if spikes_res and plot_per_neuron:\n",
+ " from tvb.simulator.plot.base_plotter import pyplot\n",
+ " spikeNet_analyzer.return_data = False\n",
+ " rates_ts_per_neuron = \\\n",
+ " spikeNet_analyzer. \\\n",
+ " compute_spikeNet_rates_time_series(populations_devices=None, regions=None,\n",
+ " computations_kwargs={}, data_kwargs={},\n",
+ " return_spikes_trains=False, return_devices=False);\n",
+ " if rates_ts_per_neuron is not None and rates_ts_per_neuron.size:\n",
+ " # Regions in rows\n",
+ " row = rates_ts_per_neuron.dims[2] if rates_ts_per_neuron.shape[2] > 1 else None\n",
+ " if row is None:\n",
+ " # Populations in rows\n",
+ " row = rates_ts_per_neuron.dims[1] if rates_ts_per_neuron.shape[1] > 1 else None\n",
+ " col = None\n",
+ " else:\n",
+ " # Populations in columns\n",
+ " col = rates_ts_per_neuron.dims[1] if rates_ts_per_neuron.shape[1] > 1 else None\n",
+ " pyplot.figure()\n",
+ " rates_ts_per_neuron.plot(y=rates_ts_per_neuron.dims[3], row=row, col=col, cmap=\"jet\")\n",
+ " plotter.base._save_figure(figure_name=\"Spike rates per neuron\")\n",
+ " # del rates_ts_per_neuron # to free memory"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Plot per neuron SpikingNetwork time series"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Regions in rows\n",
+ "if plot_per_neuron and spikeNet_ts.size:\n",
+ " row = spikeNet_ts.dims[2] if spikeNet_ts.shape[2] > 1 else None\n",
+ " if row is None:\n",
+ " # Populations in rows\n",
+ " row = spikeNet_ts.dims[3] if spikeNet_ts.shape[3] > 1 else None\n",
+ " col = None\n",
+ " else:\n",
+ " # Populations in cols\n",
+ " col = spikeNet_ts.dims[3] if spikeNet_ts.shape[3] > 1 else None\n",
+ " for var in spikeNet_ts.coords[spikeNet_ts.dims[1]]:\n",
+ " this_var_ts = spikeNet_ts.loc[:, var, :, :, :]\n",
+ " this_var_ts.name = var.item()\n",
+ " pyplot.figure()\n",
+ " this_var_ts.plot(y=spikeNet_ts.dims[4], row=row, col=col, cmap=\"jet\", figsize=FIGSIZE)\n",
+ " plotter.base._save_figure(\n",
+ " figure_name=\"Spiking Network variables' time series per neuron: %s\" % this_var_ts.name)\n",
+ " del spikeNet_ts # to free memory"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {
+ "run_control": {
+ "frozen": false,
+ "read_only": false
+ }
+ },
+ "source": [
+ "# References\n",
+ "\n",
+ "1 Sanz Leon P, Knock SA, Woodman MM, Domide L,
\n",
+ " Mersmann J, McIntosh AR, Jirsa VK (2013)
\n",
+ " The Virtual Brain: a simulator of primate brain network dynamics.
\n",
+ " Frontiers in Neuroinformatics 7:10. doi: 10.3389/fninf.2013.00010
\n",
+ " https://www.thevirtualbrain.org/tvb/zwei
\n",
+ " https://github.com/the-virtual-brain
\n",
+ "\n",
+ "2 Ritter P, Schirner M, McIntosh AR, Jirsa VK (2013).
\n",
+ " The Virtual Brain integrates computational modeling
\n",
+ " and multimodal neuroimaging. Brain Connectivity 3:121–145.
\n",
+ "\n",
+ "3 Vitay J, Dinkelbach HÜ and Hamker FH (2015).
\n",
+ " ANNarchy: a code generation approach to neural simulations on parallel hardware.
\n",
+ " Frontiers in Neuroinformatics 9:19. doi:10.3389/fninf.2015.00019
\n",
+ " For more details see https://annarchy.readthedocs.io/en/latest/
\n",
+ "\n",
+ "4 Baladron, J., Nambu, A., & Hamker, F. H. (2019).
\n",
+ " The subthalamic nucleus‐external globus pallidus loop biases
\n",
+ " exploratory decisions towards known alternatives: A neuro‐computational study.
\n",
+ " European Journal of Neuroscience, 49:754–767. https://doi.org/10.1111/ejn.13666
\n",
+ " \n",
+ "5 Maith O, Villagrasa Escudero F, Ülo Dinkelbach H, Baladron J,
\n",
+ " Horn, A, Irmen F, Kühn AA, Hamker FH (2020).
\n",
+ " A computational model‐based analysis of basal ganglia pathway changes
\n",
+ " in Parkinson’s disease inferred from resting‐state fMRI
\n",
+ " European Journal of Neuroscience, 00:1–18. https://doi.org/10.1111/ejn.14868"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(\"Configured connections:\\n\")\n",
+ "\n",
+ "print(\"Within node's connections:\")\n",
+ "for iP, pop in enumerate(ann_model_builder.populations_connections):\n",
+ " if hasattr(pop[\"weight\"], \"__call__\"):\n",
+ " weight = pop[\"weight\"](pop[\"nodes\"])\n",
+ " else:\n",
+ " weight = pop[\"weight\"]\n",
+ " try:\n",
+ " p = pop[\"conn_spec\"][\"probability\"]\n",
+ " except:\n",
+ " p = 1.0\n",
+ " w_conn = assert_loadedParams[np.abs(weight)]\n",
+ " conn = w_conn.split(\"_weights\")[0]\n",
+ " p_conn = loadedParams[\"%s_probs\" % conn]\n",
+ " w_conn = loadedParams[w_conn]\n",
+ " print(\"%d. %s -> %s (%s) = %g (p=%g): %s (%g, %g)\" \n",
+ " % (iP+1, pop[\"source\"], pop[\"target\"], pop[\"receptor_type\"], weight, p,\n",
+ " conn, w_conn, p_conn))\n",
+ " \n",
+ "print(\"\\nAmong node's connections:\")\n",
+ "for iP, pop in enumerate(ann_model_builder.nodes_connections):\n",
+ " if hasattr(pop[\"weight\"], \"__call__\"):\n",
+ " weight = pop[\"weight\"](pop[\"source_nodes\"], pop[\"target_nodes\"])\n",
+ " else:\n",
+ " weight = pop[\"weight\"]\n",
+ " try:\n",
+ " p = pop[\"conn_spec\"][\"probability\"]\n",
+ " except:\n",
+ " p = 1.0\n",
+ " w_conn = assert_loadedParams[np.abs(weight)]\n",
+ " conn = w_conn.split(\"_weights\")[0]\n",
+ " p_conn = loadedParams[\"%s_probs\" % conn]\n",
+ " w_conn = loadedParams[w_conn]\n",
+ " print(\"%d. %s -> %s (%s) = %g (p=%g): %s (%g, %g)\" \n",
+ " % (iP + 5, pop[\"source\"], pop[\"target\"], pop[\"receptor_type\"], weight, p,\n",
+ " conn, w_conn, p_conn))\n",
+ " \n",
+ " \n",
+ "print(\"\\nEffective connections:\")\n",
+ "conns = [\"GPeGPe\", \"GPiGPi\", \"dSNdSN\", \"iSNiSN\", # \"CtxICtxI\", \"CtxECtxI\", \"CtxICtxE\", \n",
+ " \"dSNGPi\" , \"iSNGPe\", \"GPeGPi\", \"GPiThal\", \"GPeSTN\", \"ThaldSN\", \"ThaliSN\", \n",
+ " \"STNGPe\", \"STNGPi\"] # ,\"CtxThal\", \"CtxSTN\", \"CtxdSN\", \"CtxiSN\"\n",
+ "for iC, (name, proj) in enumerate(zip(conns, \n",
+ " annarchy_network.annarchy_instance.Global._network[0][\"projections\"])):\n",
+ " meanNconns = np.mean([len(d.pre_ranks) for d in proj.dendrites])\n",
+ " p = meanNconns / proj.pre.size\n",
+ " print(\"%d. %s: %s w = %g (%g) (%s), effective_probability = %g (%g)\" % \n",
+ " (iC+1, proj.name, name, proj.w, loadedParams[\"%s_weights\" % name], proj.target,\n",
+ " p, loadedParams[\"%s_probs\" % name])) # meanNconns = %g, meanNconns,"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "names = [\"CtxSTN\", \"CtxdSN\", \"CtxiSN\"]\n",
+ "weights = np.array([wCtxSTN, wCtxdSN, wCtxiSN])\n",
+ "# iweights = np.array([iwCtxSTN, iwCtxdSN, iwCtxiSN])\n",
+ "probs = np.array([pCtxSTN, pCtxdSN, pCtxiSN])\n",
+ "# wTVBs = tvb_spikeNet_model_builder.global_coupling_scaling[0] * np.array([iwCtxSTN * wTVBSTNs, iwCtxdSN * wTVBSNs, iwCtxiSN * wTVBSNs])\n",
+ "eff_weights = []\n",
+ "eff_probs = []\n",
+ "n_conn = -1\n",
+ "# G = tvb_spikeNet_model_builder.global_coupling_scaling[0]\n",
+ "for iC, (proj, interface) in enumerate(\n",
+ " zip(annarchy_network.annarchy_instance.Global._network[0][\"projections\"][13:], \n",
+ " simulator.output_interfaces.interfaces)):\n",
+ " w = proj.w # * interface.transformer.scale_factor[0] * G * np.sum(wTVBs[iC])\n",
+ " name = names[iC]\n",
+ " weight = weights[iC]\n",
+ " prob = probs[iC]\n",
+ " eff_weights.append(w)\n",
+ " meanNconns = np.mean([len(d.pre_ranks) for d in proj.dendrites])\n",
+ " p = meanNconns / proj.pre.size\n",
+ " eff_probs.append(p)\n",
+ " print(\"%d. %s: eff w = %g == %g (exp), eff p = %g ~= %g (exp)\" % \n",
+ " (iC+1, proj.name, w, \n",
+ " weight, # G*np.sum(wTVBs[iC]), \n",
+ " p, prob))\n",
+ "\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rates = []\n",
+ "spikes = {}\n",
+ "TVBrates = [TVBrateSTN, TVBratedSN, TVBrateiSN]\n",
+ "for iR, (interface, TVBrate) in enumerate(zip(simulator.output_interfaces.interfaces, TVBrates)):\n",
+ " region = interface.proxy.target.keys()[0]\n",
+ " monitor = interface.proxy.target[0]._record\n",
+ " if spikes.get(region, False) == False:\n",
+ " spikes[region] = monitor.get(\"spike\").copy()\n",
+ " spike_times, ranks = monitor.raster_plot(spikes[region].copy())\n",
+ " if len(ranks):\n",
+ " eff_spike_times = np.array(spike_times)[np.where(spike_times > transient)]\n",
+ " rate = len(eff_spike_times) / np.unique(ranks).size / (simulation_length - transient) * 1000.0\n",
+ " else:\n",
+ " rate = 0.0\n",
+ " rates.append(rate)\n",
+ " print(\"%s rate = %g ~= %g (exp)\" % (region, rate, TVBrate))\n",
+ "total_rate = np.sum(rates)\n",
+ "print(\"mean rate = %g (%g) ~= %g (exp)\" % (np.mean(rates), total_rate/3, np.array(TVBrates).mean()))"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$ \\frac{dx}{dt} = \\frac{\\mu-x}{\\tau} + \\sigma\\frac{\\xi}{\\sqrt{\\tau}} $$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "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.9.2"
+ },
+ "latex_envs": {
+ "bibliofile": "biblio.bib",
+ "cite_by": "apalike",
+ "current_citInitial": 1,
+ "eqLabelWithNumbers": true,
+ "eqNumInitial": 0
+ },
+ "nav_menu": {},
+ "toc": {
+ "navigate_menu": true,
+ "number_sections": true,
+ "sideBar": true,
+ "threshold": 6,
+ "toc_cell": false,
+ "toc_section_display": "block",
+ "toc_window_display": false
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/examples/tvb_netpyne/notebooks/Roopa-TVB-corr-cortex-opt.py b/examples/tvb_netpyne/notebooks/Roopa-TVB-corr-cortex-opt.py
new file mode 100644
index 000000000..d6bb40e62
--- /dev/null
+++ b/examples/tvb_netpyne/notebooks/Roopa-TVB-corr-cortex-opt.py
@@ -0,0 +1,862 @@
+from IPython.display import Image, display
+import os
+from collections import OrderedDict
+import time
+import numpy as np
+
+from tvb.basic.profile import TvbProfile
+TvbProfile.set_profile(TvbProfile.LIBRARY_PROFILE)
+
+from tvb_multiscale.tvb_netpyne.config import *
+
+work_path = "/home/docker/packages/tvb-multiscale/examples/tvb_netpyne/notebooks" #os.getcwd()
+data_path = os.path.expanduser("~/packages/tvb-multiscale/examples/data/basal_ganglia")
+fit_data_path = os.path.join(data_path, "ANNarchyFittedModels/dataFits_2020_02_05/databestfits", )
+data_mode = "patient" # "control", or "patient"
+control_data = os.path.join(fit_data_path, "controlleft/OutputSim_Patient08.mat")
+patient_data = os.path.join(fit_data_path, "patientleft/OutputSim_Patient09.mat")
+
+INTERFACE_COUPLING_MODE = "spikeNet" # "TVB" or "spikeNet"
+
+if data_mode == "patient":
+ subject_data = patient_data
+ if INTERFACE_COUPLING_MODE == "TVB":
+ TC_factor = 20.5 * 1e-2 # 0.205
+ else:
+ TC_factor = 20.5 * 1e-2 # 0.205
+else:
+ subject_data = control_data
+ if INTERFACE_COUPLING_MODE == "TVB":
+ TC_factor = 20.5 * 1e-2 # 0.205
+
+ else:
+ TC_factor = 20.5 * 1e-2 # 0.205
+
+
+simulation_length = 150.0
+transient = 25.0
+start_stimulus = 100.0
+init_cond_jitter = 0.0
+
+SPIKING_NODES_DELAYS = False
+
+simulation_mode = "rs" # "stim" or "rs"
+stim_target = "PY_pop" # "PY_pop", "TC_pop"
+stim_mode = "simple" # "bi" | "mono" | "simple"
+ # -------------------------
+stim_freq = 0.0 # 130.0 | 120.0 | 0.0
+stim_ampl = -10.0 # 20.0 | -35.0 | -10.0
+stim_duration = 0.0 # 0.3 | 0.3 | 0.0
+if simulation_mode == "stim":
+ simulation_mode = simulation_mode + "_%s_%s" % (stim_target, stim_mode)
+
+outputs_path = os.path.join(work_path, "outputs")
+sim_mode_path = os.path.join(outputs_path, "TVBcortex_%s_coupl" % INTERFACE_COUPLING_MODE,
+ data_mode, simulation_mode)
+config = Config(output_base=sim_mode_path)
+config.figures.SHOW_FLAG = True
+config.figures.SAVE_FLAG = True
+config.figures.FIG_FORMAT = 'png'
+config.figures.DEFAULT_SIZE= config.figures.NOTEBOOK_SIZE
+FIGSIZE = config.figures.DEFAULT_SIZE
+
+from tvb_multiscale.core.plot.plotter import Plotter
+plotter = Plotter(config.figures)
+
+# For interactive plotting:
+# %matplotlib notebook
+
+# Otherwise:
+#%matplotlib inline
+
+from tvb_multiscale.core.tvb.cosimulator.models.reduced_wong_wang_exc_io import ReducedWongWangExcIO
+
+# ----------------------------------------------------------------------------------------------------------------
+# ----Uncomment below to modify the simulator by changing the default options:--------------------------------------
+# ----------------------------------------------------------------------------------------------------------------
+
+from tvb.datatypes.connectivity import Connectivity
+from tvb.simulator.integrators import HeunStochastic
+from tvb.simulator.monitors import Raw # , Bold, EEG
+
+
+# 0. GPe_Left, 1. GPi_Left, 2. STN_Left, 3. Striatum_Left, 4. Thal_Left
+BG_opt_matrix_weights = np.zeros((5, 5))
+conn_mode = "subject" # subject" or "average"
+if conn_mode == "average":
+ weights_maith = np.array([1.93, 3.56, 1.46, 4.51, 3.52, 2.30, 2.34, 3.78, 1.98,
+ 1.30, 1.82, 3.56, 3.02, 1.78, 1.36, 2.27, 4.13, 2.74, 3.27])*1e-3 # controls
+# weights_maith = np.array([3.27, 3.80, 2.65, 3.66, 3.06, 3.06, 3.25, 4.02, 3.32,
+# 2.98, 3.45, 3.64, 2.50, 2.12, 2.86, 2.79, 3.96, 3.69, 3.87])*1e-3 # patients
+ # probs_maith = ????
+else:
+ import scipy.io as sio
+ weights=sio.loadmat(subject_data) # weights start from index 19
+ weights_maith = weights["X"][0, 19:] # these are indices 19 till 37
+ probs_maith = weights["X"][0, :19] # these are indices 0 till 18
+
+wdSNGPi = BG_opt_matrix_weights[3, 1] = weights_maith[0].item()
+wiSNGPe = BG_opt_matrix_weights[3, 0] = weights_maith[1].item()
+wGPeSTN = BG_opt_matrix_weights[0, 2] = weights_maith[2].item()
+wSTNGPe = BG_opt_matrix_weights[2, 0] = weights_maith[3].item()
+wSTNGPi = BG_opt_matrix_weights[2, 1] = weights_maith[4].item()
+wGPeGPi = BG_opt_matrix_weights[0, 1] = weights_maith[5].item()
+wGPiTh = BG_opt_matrix_weights[1, 4] = weights_maith[8].item()
+wThdSN = BG_opt_matrix_weights[4, 3] = weights_maith[10].item() # Th -> dSN
+
+sliceBGnet = slice(0,5)
+
+wGPeGPe = weights_maith[6].item() # "GPe" -> "GPe"
+wGPiGPi = weights_maith[7].item() # "GPi" -> "GPi"
+wThiSN = weights_maith[9].item() # "Eth" -> "IiSN"
+
+wdSNdSN = weights_maith[11].item() # "IdSN" -> "IdSN"
+wiSNiSN = weights_maith[12].item() # "IiSN" -> "IiSN"
+wCtxdSN = weights_maith[13].item() # "CxE" -> "IdSN"
+wCtxiSN = weights_maith[14].item() # "CxE" -> "IiSN"
+wCtxSTN = weights_maith[15].item() # "CxE" -> "Estn"
+wCtxEtoI = weights_maith[16].item() # "CxE" -> "CxI"
+wCtxItoE = weights_maith[17].item() # "CxI" -> "CxE"
+wCtxItoI = weights_maith[18].item() # "CxI" -> "CxI"
+
+pdSNGPi = probs_maith[0].item()
+piSNGPe = probs_maith[1].item()
+pGPeSTN = probs_maith[2].item()
+pSTNGPe = probs_maith[3].item()
+pSTNGPi = probs_maith[4].item()
+pGPeGPi = probs_maith[5].item()
+pGPeGPe = probs_maith[6].item() # "GPe" -> "GPe"
+pGPiGPi = probs_maith[7].item() # "GPi" -> "GPi"
+pGPiTh = probs_maith[8].item()
+pThiSN = probs_maith[9].item() # "Eth" -> "IiSN
+pThdSN = probs_maith[10].item() # Th --> dSN
+pdSNdSN = probs_maith[11].item() # "IdSN" -> "IdSN"
+piSNiSN = probs_maith[12].item() # "IiSN" -> "IiSN"
+pCtxdSN = probs_maith[13].item() # "CxE" -> "IdSN"
+pCtxiSN = probs_maith[14].item() # "CxE" -> "IiSN"
+pCtxSTN = probs_maith[15].item() # "CxE" -> "Estn"
+pCtxEtoI = probs_maith[16].item() # "CxE" -> "CxI"
+pCtxItoE = probs_maith[17].item() # "CxI" -> "CxE"
+pCtxItoI = probs_maith[18].item() # "CxI" -> "CxI"
+pCtxCtx = probs_maith[16:19].mean() # "Ctx" -> "Ctx"
+
+loadedParams ={'dSNGPi_probs': probs_maith[0],
+ 'dSNGPi_weights' : weights_maith[0],
+ 'iSNGPe_probs' : probs_maith[1],
+ 'iSNGPe_weights' : weights_maith[1],
+ 'GPeSTN_probs' : probs_maith[2],
+ 'GPeSTN_weights' : weights_maith[2],
+ 'STNGPe_probs' : probs_maith[3],
+ 'STNGPe_weights' : weights_maith[3],
+ 'STNGPi_probs' : probs_maith[4],
+ 'STNGPi_weights' : weights_maith[4],
+ 'GPeGPi_probs' : probs_maith[5],
+ 'GPeGPi_weights' : weights_maith[5],
+ 'GPeGPe_probs' : probs_maith[6],
+ 'GPeGPe_weights' : weights_maith[6],
+ 'GPiGPi_probs' : probs_maith[7],
+ 'GPiGPi_weights' : weights_maith[7],
+ 'GPiThal_probs' : probs_maith[8],
+ 'GPiThal_weights' : weights_maith[8],
+ 'ThaliSN_probs' : probs_maith[9],
+ 'ThaliSN_weights' : weights_maith[9],
+ 'ThaldSN_probs' : probs_maith[10],
+ 'ThaldSN_weights' : weights_maith[10],
+ 'dSNdSN_probs' : probs_maith[11],
+ 'dSNdSN_weights' : weights_maith[11],
+ 'iSNiSN_probs' : probs_maith[12],
+ 'iSNiSN_weights' : weights_maith[12],
+ 'CtxdSN_probs' : probs_maith[13],
+ 'CtxdSN_weights' : weights_maith[13],
+ 'CtxiSN_probs' : probs_maith[14],
+ 'CtxiSN_weights' : weights_maith[14],
+ 'CtxSTN_probs' : probs_maith[15],
+ 'CtxSTN_weights' : weights_maith[15],
+ 'CtxECtxI_probs' : probs_maith[16],
+ 'CtxECtxI_weights' : weights_maith[16],
+ 'CtxICtxE_probs' : probs_maith[17],
+ 'CtxICtxE_weights' : weights_maith[17],
+ 'CtxICtxI_probs' : probs_maith[18],
+ 'CtxICtxI_weights' : weights_maith[18],
+ 'CtxThal_weights': 0.0,
+ 'CtxThal_probs': 1.0}
+#print(loadedParams)
+
+assert_loadedParams = dict(zip(loadedParams.values(), loadedParams.keys()))
+
+# Load full TVB connectome connectivity
+
+conn_path = os.path.join(data_path, "conn")
+
+#Load AAL atlas normative connectome including the Basal Ganglia regions from Petersen et al. atlas
+wTVB = np.loadtxt(os.path.join(conn_path, "conn_denis_weights.txt"))
+cTVB = np.loadtxt(os.path.join(conn_path, "aal_plus_BG_centers.txt"),usecols=range(1,4))
+rlTVB = np.loadtxt(os.path.join(conn_path, "aal_plus_BG_centers.txt"),dtype="str", usecols=(0,))
+tlTVB = np.loadtxt(os.path.join(conn_path, "BGplusAAL_tract_lengths.txt"))
+
+# Remove the second Thalamus, Pallidum (GPe/i), Putamen and Caudate (Striatum):
+inds_Th = (rlTVB.tolist().index("Thalamus_L"), rlTVB.tolist().index("Thalamus_R"))
+inds_Pall = (rlTVB.tolist().index("Pallidum_L"), rlTVB.tolist().index("Pallidum_R"))
+inds_Put = (rlTVB.tolist().index("Putamen_L"), rlTVB.tolist().index("Putamen_R"))
+inds_Caud = (rlTVB.tolist().index("Caudate_L"), rlTVB.tolist().index("Caudate_R"))
+inds_rm = inds_Th + inds_Pall + inds_Put + inds_Caud
+#print("Connections of Thalami, Pallidum (GPe/i), Putamen and Caudate (Striatum) removed!:\n",
+wTVB[inds_rm, :][:, inds_rm]
+wTVB = np.delete(wTVB, inds_rm, axis=0)
+wTVB = np.delete(wTVB, inds_rm, axis=1)
+tlTVB = np.delete(tlTVB, inds_rm, axis=0)
+tlTVB = np.delete(tlTVB, inds_rm, axis=1)
+rlTVB = np.delete(rlTVB, inds_rm, axis=0)
+cTVB = np.delete(cTVB, inds_rm, axis=0)
+
+number_of_regions = len(rlTVB)
+speed = np.array([4.0])
+min_tt = speed.item() * 0.1
+sliceBG = [0, 1, 2, 3, 6, 7]
+sliceCortex = slice(10, number_of_regions)
+
+# Remove BG -> Cortex connections
+#print("Removing BG -> Cortex connections with max:")
+#print(wTVB[sliceCortex, :][:, sliceBG].max())
+wTVB[sliceCortex, sliceBG] = 0.0
+tlTVB[sliceCortex, sliceBG] = min_tt
+
+"""
+# Remove Cortex -> Thalamus connections
+# They need to be removed because I'm creating an interface between the TVB motor cortex and the spiking thalamus (TC_pop).
+# That will be the only/"replacement connection" between the motor cortex and the thalamus
+# TODO: add this code back. I'm keeping it for now because otherwise I get a "divide by zero" error downstream
+# TODO: why is sliceThal = 8 & 9 here when thalamus node number is 4?
+sliceThal = [8, 9]
+#print("Removing Cortex -> Thalamus connections with summed weight:")
+#print(wTVB[sliceThal, sliceCortex].sum())
+wTVB[sliceThal, sliceCortex] = 0.0
+tlTVB[sliceThal, sliceCortex] = min_tt
+"""
+
+# Remove Cortex -> GPe/i connections
+sliceGP = [0, 1, 2, 3]
+#print("Removing Cortex -> GPe/i connections with max:")
+#print(wTVB[sliceGP, sliceCortex].max())
+wTVB[sliceGP, sliceCortex] = 0.0
+tlTVB[sliceGP, sliceCortex] = min_tt
+
+# # Minimize all delays for the optimized network
+# tlTVB[:7][:, :7] = min_tt
+
+connTVB = Connectivity(region_labels=rlTVB, weights=wTVB, centres=cTVB, tract_lengths=tlTVB, speed=speed)
+
+# Normalize connectivity weights
+connTVB.weights = connTVB.scaled_weights(mode="region")
+connTVB.weights /= np.percentile(connTVB.weights, 99)
+
+# Keep only left hemisphere and remove Vermis:
+sliceLeft = slice(0, connTVB.number_of_regions -8, 2)
+
+connLeft = Connectivity(region_labels=connTVB.region_labels[sliceLeft],
+ centres=connTVB.centres[sliceLeft],
+ weights=connTVB.weights[sliceLeft][:, sliceLeft],
+ tract_lengths=connTVB.tract_lengths[sliceLeft][:, sliceLeft],
+ speed=connTVB.speed)
+connLeft.configure()
+
+#print("\nLeft cortex connectome, after removing direct BG -> Cortex and interhemispheric BG <-> BG connections:")
+plotter.plot_tvb_connectivity(connLeft);
+
+sliceBGnet = slice(0,5)
+connTVBleftBG = Connectivity(region_labels=connLeft.region_labels[sliceBGnet],
+ centres=connLeft.centres[sliceBGnet],
+ weights=connLeft.weights[sliceBGnet][:, sliceBGnet],
+ tract_lengths=connLeft.tract_lengths[sliceBGnet][:, sliceBGnet],
+ speed=connLeft.speed)
+connTVBleftBG.configure()
+
+#print("\nLeft BG TVB network:")
+plotter.plot_tvb_connectivity(connTVBleftBG);
+
+scaleBGoptTOtvb = np.percentile(BG_opt_matrix_weights, 95) /\
+ np.percentile(connTVBleftBG.weights, 95)
+
+#print("Scaling factor of TVB BG network connectome to optimal one = %g" % scaleBGoptTOtvb)
+# confirmed: 0.000771577
+
+# Construct the final connectivity to use for simulation:
+# Rescale the
+ww = scaleBGoptTOtvb * np.array(connLeft.weights)
+ww[sliceBGnet, sliceBGnet] = BG_opt_matrix_weights.T # !!!NOTE TVB indices convention Wi<-j!!!
+
+connectivity = Connectivity(region_labels=connLeft.region_labels,
+ centres=connLeft.centres,
+ weights=ww, tract_lengths=connLeft.tract_lengths,
+ speed=connLeft.speed)
+connectivity.configure()
+
+#print("connectivity info")
+#print(connLeft.region_labels[4])
+#print(ww[4,:])
+#print(ww[:,4])
+
+print("\nFinal connectivity:")
+plotter.plot_tvb_connectivity(connectivity);
+
+import matplotlib.pyplot as plt
+
+# plot connectome weights & tract lengths
+f = plt.matshow(wTVB)
+plt.savefig(sim_mode_path+"/figs/"+"connectivity.png")
+f = plt.matshow(connLeft.tract_lengths)
+plt.savefig(sim_mode_path+"/figs/"+"tracts_lengths.png")
+
+# Construct only the optimized BG connectivity only for plotting:
+connBGopt = Connectivity(region_labels=connectivity.region_labels[sliceBGnet],
+ centres=connectivity.centres[sliceBGnet],
+ weights=connectivity.weights[sliceBGnet][:, sliceBGnet],
+ tract_lengths=connectivity.tract_lengths[sliceBGnet][:, sliceBGnet],
+ speed=connectivity.speed)
+connBGopt.configure()
+
+from tvb_multiscale.core.tvb.cosimulator.cosimulator_serial import CoSimulatorSerial as CoSimulator
+
+#white_matter_coupling = coupling.Linear(a=0.014)
+# Create a TVB simulator and set all desired inputs
+# (connectivity, model, surface, stimuli etc)
+# We choose all defaults in this example
+simulator = CoSimulator()
+#simulator.use_numba = False
+model_params = {"G": np.array([15.0/scaleBGoptTOtvb])}
+simulator.model = ReducedWongWangExcIO(**model_params)
+
+simulator.connectivity = connectivity
+
+simulator.integrator = HeunStochastic()
+simulator.integrator.dt = 0.1
+simulator.integrator.noise.nsig = np.array([1e-4]) # 1e-5
+
+mon_raw = Raw(period=1.0) # ms
+simulator.monitors = (mon_raw, )
+
+init_cond_filepath = os.path.join(data_path, "tvb_init_cond_left.npy")
+init_cond = np.load(init_cond_filepath) #
+init_cond = np.abs(init_cond *(1 + init_cond_jitter * np.random.normal(size=init_cond.shape)))
+simulator.connectivity.set_idelays(simulator.integrator.dt)
+simulator.initial_conditions = init_cond * np.ones((simulator.connectivity.idelays.max(),
+ simulator.model.nvar,
+ simulator.connectivity.number_of_regions,
+ simulator.model.number_of_modes))
+
+
+print("\nConnectome used for simulations:")
+plotter.plot_tvb_connectivity(simulator.connectivity);
+
+
+# # Serializing TVB cosimulator is necessary for parallel cosimulation:
+# from tvb_multiscale.core.utils.file_utils import dump_pickled_dict
+# from tvb_multiscale.core.tvb.cosimulator.cosimulator_serialization import serialize_tvb_cosimulator
+# sim_serial_filepath = os.path.join(config.out.FOLDER_RES, "tvb_serial_cosimulator.pkl")
+# sim_serial = serialize_tvb_cosimulator(simulator)
+# display(sim_serial)
+
+# # Dumping the serialized TVB cosimulator to a file will be necessary for parallel cosimulation.
+# dump_pickled_dict(sim_serial, sim_serial_filepath)
+
+###################################################################################################################
+# START MODIFYING SCRIPT FROM HERE
+###################################################################################################################
+
+from tvb_multiscale.tvb_netpyne.netpyne_models.builders.base import NetpyneNetworkBuilder
+
+# Select the regions for the fine scale modeling with NetPyNE spiking networks
+# numbers come from overall connectome, so find where CCTC regions are. Add more for e.g. DCN and NO if they're the same region? or just do thalamus for now
+spiking_nodes_ids = [4] # the indices of fine scale regions modeled with NetPyNE
+
+# Build a ANNarchy network model with the corresponding builder
+from tvb_multiscale.tvb_netpyne.netpyne_models.builders.netpyne_factory import load_netpyne
+netpyne = load_netpyne(config=config)
+
+spiking_network_builder = NetpyneNetworkBuilder(simulator, spiking_nodes_ids, netpyne_instance=netpyne, config=config)
+
+
+# # ----------------------------------------------------------------------------------------------------------------
+# # ----Uncomment below to modify the builder by changing the default options:--------------------------------------
+# # ----------------------------------------------------------------------------------------------------------------
+from copy import deepcopy
+
+spiking_network_builder.population_order = 200 # reduce for speed
+
+##TODO: add rest of cfg/netParams function, or is this enough?
+from tvb_multiscale.tvb_netpyne.netpyne_models.models.thalamic_VIM_ET.src.netParams import netParams
+from tvb_multiscale.tvb_netpyne.netpyne_models.models.thalamic_VIM_ET.src.cfg import cfg
+
+# Populations' configurations
+# When any of the properties model, params and scale below depends on regions,
+# set a handle to a function with
+# arguments (region_index=None) returning the corresponding property
+
+# TODO: allow for having several nodes in single entry
+vim_node_id = [4]
+spiking_network_builder.populations = [
+ {"label": "TC_pop",
+ "params": {"global_label": "TC_pop"},
+ "nodes": vim_node_id, # Eth in [4]
+ "scale": 1.0}
+]
+
+# Within region-node connections
+# When any of the properties model, conn_spec, weight, delay, receptor_type below
+# set a handle to a function with
+# arguments (region_index=None) returning the corresponding property
+
+synapse_model = None
+# conn_spec = {'rule': "all_to_all",
+# "allow_self_connections": True, "force_multiple_weights": False}
+# conn_spec_fixed_probability = conn_spec.copy()
+# conn_spec_fixed_probability.update({'rule': "fixed_probability", "probability": 0.1})
+
+def conn_spec_fixed_prob(prob=None):
+ return {"rule": {"prob": prob}}
+
+within_node_delay = 1.0 # ms
+
+"""
+# for each connection, we have a different probability
+spiking_network_builder.populations_connections = [
+ # source -> target
+ {"source": "I", "target": "I", # I -> I This is a self-connection for population "Igpe"
+ "synapse_model": synapse_model, "conn_spec": conn_spec_fixed_prob(pGPeGPe), # conn_spec
+ "weight": np.abs(wGPeGPe).item(), "delay": within_node_delay,
+ "receptor_type": "gaba", "nodes": Igpe_nodes_ids}, # None means apply to all
+ {"source": "I", "target": "I", # I -> I This is a self-connection for population "Igpi"
+ "synapse_model": synapse_model, "conn_spec": conn_spec_fixed_prob(pGPiGPi), # conn_spec
+ "weight": np.abs(wGPiGPi).item(), "delay": within_node_delay,
+ "receptor_type": "gaba", "nodes": Igpi_nodes_ids}, # None means apply to all
+ {"source": "IdSN", "target": "IdSN", # IdSN -> IdSN This is a self-connection for population "IdSN"
+ "synapse_model": synapse_model, "conn_spec": conn_spec_fixed_prob(pdSNdSN), # conn_spec
+ "weight": np.abs(wdSNdSN).item(), "delay": within_node_delay,
+ "receptor_type": "gaba", "nodes": Istr_nodes_ids},
+ {"source": "IiSN", "target": "IiSN", # IiSN -> IiSN This is a self-connection for population "IiSN"
+ "synapse_model": synapse_model, "conn_spec": conn_spec_fixed_prob(piSNiSN), # conn_spec
+ "weight": np.abs(wiSNiSN).item(), "delay": within_node_delay,
+ "receptor_type": "gaba", "nodes": Istr_nodes_ids},
+ ]
+"""
+
+# Among/Between region-node connections
+# Given that only the AMPA population of one region-node couples to
+# all populations of another region-node,
+# we need only one connection type
+
+# When any of the properties model, conn_spec, weight, delay, receptor_type below
+# depends on regions, set a handle to a function with
+# arguments (source_region_index=None, target_region_index=None)
+
+from tvb_multiscale.core.spiking_models.builders.templates import scale_tvb_weight, tvb_delay
+
+# We set global coupling scaling to 1.0,
+# because we need the Maith et al optimized weights without any scaling:
+spiking_network_builder.global_coupling_scaling = 1.0
+
+# Function that will return the TVB weight with optional scaling:
+class TVBWeightFun(object):
+
+ def __init__(self,
+ global_coupling_scaling=spiking_network_builder.global_coupling_scaling,
+ tvb_weights = spiking_network_builder.tvb_weights):
+ self.global_coupling_scaling = float(global_coupling_scaling)
+ self.tvb_weights = tvb_weights.copy()
+
+ def __call__(self, source_node, target_node):
+ return scale_tvb_weight(source_node, target_node, self.tvb_weights,
+ scale=self.global_coupling_scaling)
+
+# Function that will return the TVB delay unless SPIKING_NODES_DELAYS == False:
+tvb_delay_fun = \
+ lambda source_node, target_node: \
+ np.maximum(spiking_network_builder.tvb_dt,
+ tvb_delay(source_node, target_node, spiking_network_builder.tvb_delays)) \
+ if SPIKING_NODES_DELAYS else within_node_delay
+
+# Creating devices to be able to observe ANNarchy activity:
+
+spiking_network_builder.output_devices = []
+
+period = 1.0
+
+# Creating devices to be able to observe ANNarchy activity:
+params = {} # deepcopy(spiking_network_builder.config.ANNARCHY_OUTPUT_DEVICES_PARAMS_DEF["SpikeMonitor"])
+for pop in spiking_network_builder.populations:
+ connections = OrderedDict({})
+ # label <- target population
+ connections[pop["label"]] = pop["label"]
+ spiking_network_builder.output_devices.append(
+ {"model": "spike_recorder", "params": deepcopy(params),
+ "connections": connections, "nodes": pop["nodes"]}) # None means apply to "all"
+
+# Labels have to be different for every connection to every distinct population
+
+
+# # Create a spike stimulus input device
+spiking_network_builder.input_devices = [] #
+
+# ----------------------------------------------------------------------------------------------------------------
+# ----------------------------------------------------------------------------------------------------------------
+# ----------------------------------------------------------------------------------------------------------------
+
+# This will be transferred to NetPyNE
+config.simulation_length = simulator.simulation_length
+
+spiking_network_builder.configure(netParams=netParams, simConfig=cfg, autoCreateSpikingNodes=False)
+# confirmed: now spiking_network_builder.global_coupling_scaling = 0.00390625
+
+
+def synaptic_weight_scale_func(is_coupling_mode_tvb):
+ # TODO: to be tuned or removed?
+ if is_coupling_mode_tvb: # "TVB"
+ return 1 # 1e-2
+ else: # "spikeNet"
+ return 1 # 5
+
+
+# spiking_network_builder.global_coupling_scaling *= simulator.model.G
+spiking_network_builder.netpyne_synaptic_weight_scale = synaptic_weight_scale_func(is_coupling_mode_tvb=INTERFACE_COUPLING_MODE=="TVB")
+
+
+netpyne_network = spiking_network_builder.build() # set_defaults=False
+# confirmed: 6 spiking pops of 200 neurs
+
+#print("\n\nLOOK HERE\n\n")
+#print("print(spiking_network_builder.netpyne_instance.netParams.connParams.items())")
+#print(spiking_network_builder.netpyne_instance.netParams.connParams.items())
+for i, (connId, conn) in enumerate(spiking_network_builder.netpyne_instance.netParams.connParams.items()):
+ print(f"{i}. {connId}: {conn.get('weight')} {conn.get('probability')}")
+# confirmed: 13 connections between pops (weights and probs confirmed)
+
+# confirmed: 6 spike recorders, as in ANNarchy when "_ts" stuff commented out
+
+# TODO: 6 pops vs 4 in annarchy (but should be okay, because of different way of setting spiking_network_builder.populations)
+populations_sizes = []
+print("Population sizes: ")
+for pop in spiking_network_builder.populations:
+ populations_sizes.append(int(np.round(pop["scale"] * spiking_network_builder.population_order)))
+ print("%s: %d" % (pop["label"], populations_sizes[-1]))
+
+
+from tvb_multiscale.tvb_netpyne.interfaces.builders import TVBNetpyneInterfaceBuilder
+
+# Build a TVB-ANNarchy interface with all the appropriate connections between the
+# TVB and ANNarchy modelled regions
+tvb_spikeNet_model_builder = TVBNetpyneInterfaceBuilder()
+
+tvb_spikeNet_model_builder.config = config
+tvb_spikeNet_model_builder.tvb_cosimulator = simulator
+tvb_spikeNet_model_builder.spiking_network = netpyne_network
+# This can be used to set default tranformer and proxy models:
+tvb_spikeNet_model_builder.model = "RATE" # "RATE" (or "SPIKES", "CURRENT") TVB->ANNarchy interface
+tvb_spikeNet_model_builder.input_flag = True # If True, NetPyNE->TVB update will be implemented
+tvb_spikeNet_model_builder.output_flag = True # If True, TVB->NetPyNE coupling will be implemented
+# If default_coupling_mode = "TVB", large scale coupling towards spiking regions is computed in TVB
+# and then applied with no time delay via a single "TVB proxy node" / ANNarchy device for each spiking region,
+# "1-to-1" TVB->ANNarchy coupling.
+# If any other value, we need 1 "TVB proxy node" / ANNarchy device for each TVB sender region node, and
+# large-scale coupling for spiking regions is computed in ANNarchy,
+# taking into consideration the TVB connectome weights and delays,
+# in this "1-to-many" TVB->ANNarchy coupling.
+tvb_spikeNet_model_builder.default_coupling_mode = INTERFACE_COUPLING_MODE # "spikeNet" # "TVB"
+# Number of neurons per population to be used to compute population mean instantaneous firing rates:
+tvb_spikeNet_model_builder.proxy_inds = np.array(spiking_nodes_ids)
+tvb_spikeNet_model_builder.N_E = spiking_network_builder.population_order
+tvb_spikeNet_model_builder.TC_proxy_inds = np.array(vim_node_id)
+
+# Set exclusive_nodes = True (Default) if the spiking regions substitute for the TVB ones:
+tvb_spikeNet_model_builder.exclusive_nodes = True
+
+tvb_spikeNet_model_builder.output_interfaces = []
+tvb_spikeNet_model_builder.input_interfaces = []
+
+from tvb_multiscale.tvb_netpyne.interfaces.builders import NetpyneInputProxyModels
+
+# TVB applies a global coupling scaling of coupling.a * model.G
+tvb_spikeNet_model_builder.global_coupling_scaling = \
+ tvb_spikeNet_model_builder.tvb_cosimulator.coupling.a[0].item() * simulator.model.G
+print("global_coupling_scaling = %g" % tvb_spikeNet_model_builder.global_coupling_scaling)
+
+"""
+# Total TVB indegree weight to STN:
+wTVBSTNs = simulator.connectivity.weights[Estn_nodes_ids, 5:].squeeze()
+wTVBSTN = wTVBSTNs.sum().item()
+print("wTVBSTN = %g" % wTVBSTN)
+CTXtoSTNinds = 5 + np.where(wTVBSTNs > 0.0)[0] # indices of TVB regions coupling to STN
+
+# Total TVB indegree weight to Striatum:
+wTVBSNs = simulator.connectivity.weights[Istr_nodes_ids, 5:].squeeze()
+wTVBSN = wTVBSNs.sum().item()
+print("wTVBSN = %g" % wTVBSN)
+CTXtoSNinds = 5 + np.where(wTVBSNs > 0.0)[0] # indices of TVB regions coupling to Striatum
+
+# Approximate effective scaling of TVB coupling to STN
+# after normalizing with STN indegree and
+# multiplying with the Maith et al. optimized CTX -> STN weight
+iwCtxSTN = STN_factor*wCtxSTN / wTVBSTN
+print("iwCtxSTN = %g" % iwCtxSTN)
+
+# Approximate effective scaling of TVB coupling to dSN
+# after normalizing with Striatum indegree and
+# multiplying with the Maith et al. optimized CTX -> dSN weight
+iwCtxdSN = dSN_factor*wCtxdSN / wTVBSN
+print("iwCtxdSN = %g" % iwCtxdSN)
+
+# Approximate effective scaling of TVB coupling to iSN
+# after normalizing with Striatum indegree and
+# multiplying with the Maith et al. optimized CTX -> iSN weight
+iwCtxiSN = iSN_factor*wCtxiSN / wTVBSN
+print("iwCtxiSN = %g" % iwCtxiSN)
+"""
+#TVB (cortex) connects to & provides INPUT to spiking region TC_pop (vim)
+# Total TVB indegree weight to TC_pop:
+wTVBTCs = simulator.connectivity.weights[vim_node_id, 5:].squeeze()
+print("simulator.connectivity.weights")
+print(simulator.connectivity.weights[4,5:])
+print("wTVBTCs")
+print(wTVBTCs)
+wTVBTC = wTVBTCs.sum().item()
+print("wTVBTC = %g" % wTVBTC)
+CTXtoTCinds = 5 + np.where(wTVBTC > 0.0)[0] # indices of TVB regions coupling to TC
+
+# Approximate effective scaling of TVB coupling to TC_pop
+# after normalizing with TC indegree and
+# multiplying with the Maith et al. optimized CTX -> TC weight
+# TODO:this isn't in maith. should we still have this line here? if so, how do we set the value for this?
+wCtxTC = 2 * 1e-3 # random temp value I picked in maith value range
+pCtxTC = 0.5 # random temp value I picked (maith value range unknown)
+iwCtxTC = TC_factor*wCtxTC / wTVBTC
+print("iwCtxTC = %g" % iwCtxTC)
+
+# confirmed:
+# global_coupling_scaling = 75.9402
+# wTVBSTN = 0.013739
+# wTVBSN = 0.0503774
+# iwCtxSTN = 0.0530743
+# iwCtxdSN = 0.00666963
+# iwCtxiSN = 0.0045056
+
+# ----------------------------------------------------------------------------------------------------------------
+# ----------------------------------------- BUILDER --------------------------------------------------------------
+# ----------------------------------------------------------------------------------------------------------------
+
+from tvb_multiscale.core.interfaces.base.transformers.models.red_wong_wang import RedWongWangExc
+
+
+# --------For spike transmission from TVB to ANNarchy devices acting as TVB proxy nodes with TVB delays:--------
+
+
+# TVB -> NetPyNE
+
+
+if tvb_spikeNet_model_builder.default_coupling_mode == "spikeNet":
+
+ # If coupling is computing in NetPyNE, we need as many TVB proxies
+ # as TVB regions coupling to Thalamus
+ proxy_inds = np.unique(np.concatenate([CTXtoTCinds]))
+
+ # This is the TVB coupling weight function that will determine the connections' weights
+ # from TVB proxies to the target TC population:
+ class TVBWeightFunInterface(object):
+
+ def __init__(self, scaling):
+ self.scaling = float(scaling)
+
+ def __call__(self, source_node, target_node, tvb_weights):
+ return (scale_tvb_weight(source_node, target_node, tvb_weights, scale=self.scaling))
+
+ # A similar function for TVB coupling delays is also applied in the background
+ # without need to be explicitly defined by the user
+
+# Optionally adjust interface scale factors here
+# to have the same result counter act possible changes to G and coupling.a:
+# STN_factor /= tvb_spikeNet_model_builder.global_coupling_scaling
+# dSN_factor /= tvb_spikeNet_model_builder.global_coupling_scaling
+# iSN_factor /= tvb_spikeNet_model_builder.global_coupling_scaling
+
+#dSN_factor *= 40
+#iSN_factor *= 40
+
+
+
+#
+tvb_spikeNet_model_builder.synaptic_weight_scale_func = synaptic_weight_scale_func
+
+tvb_spikeNet_model_builder.output_interfaces = []
+# Mean spike rates are applied in parallel to all target neurons
+for trg_pop, target_nodes, conn_scaling, this_conn_spec, scale_factor in \
+ zip(["TC_pop"], # NetPyNE target populations
+ # Target region indices in NetPyNE:
+ [tvb_spikeNet_model_builder.TC_proxy_inds],
+ # Maith et al optimized...
+ [wCtxTC], # ...weights
+ # ...and probabilities for CTX -> TC_pop connections
+ [conn_spec_fixed_prob(prob=pCtxTC), # pCtxTC
+ ],
+ # Interface scaling factors scaled by TVB weights' indegree to TC:
+ [TC_factor/wTVBTC]):
+ tvb_spikeNet_model_builder.output_interfaces.append(
+ {"voi": np.array(["R"]), # Source TVB state variable
+ "populations": np.array([trg_pop]), # NetPyNE target population
+ "model": "RATE",
+ "spiking_proxy_inds": target_nodes, # Target region indices in NetPyNE
+ # This spike generator device generates spike trains
+ # with autocorrelation corr at a time scale tau
+ "proxy_model": NetpyneInputProxyModels.RATE, # ANNarchyInputProxyModels.RATE_TO_SPIKES, #
+
+ # TODO: above needed?
+
+ "proxy_params": {"geometry": 600, "record": ["spike"],
+ "corr": 0.3, "tau": 10.0, # comment for RATE_TO_SPIKES
+ "number_of_neurons": tvb_spikeNet_model_builder.N_E, # todo: de-hardcode
+ },
+ 'conn_spec': this_conn_spec, # dictionary of connection properties
+ 'coupling_mode': tvb_spikeNet_model_builder.default_coupling_mode
+ }) # None means all here
+
+ # For both coupling modes, we scale the TVB rate already at the TVB -> ANNarchy transformer
+ # with the interface scale factor (normalized by TVB indegree to STN/Striatum)
+ # and the global coupling scaling.
+ if tvb_spikeNet_model_builder.output_interfaces[-1]["coupling_mode"] == "spikeNet":
+ tvb_spikeNet_model_builder.output_interfaces[-1]["proxy_inds"] = proxy_inds
+ # In this case connections from each TVB proxy to TC
+ # are scaled additionally with the Maith et al. optimized weights
+ tvb_spikeNet_model_builder.output_interfaces[-1]["weights"] = TVBWeightFunInterface(conn_scaling)
+ tvb_spikeNet_model_builder.output_interfaces[-1]["transformer_params"] = \
+ {"scale_factor": scale_factor * tvb_spikeNet_model_builder.global_coupling_scaling}
+ # In total:
+ # From each TVB proxy node we get a rate scaled as (
+ # (coupling.a * G * STN_factor/wTVBSTN) * R_i, (i for all TVB regions)
+ # Then, spikes generated from each TVB proxy are transferred via connections
+ # with weights TVB_w_ji * wCtxSTN or wCtxiSN or wCtxdSN (j for STN or Striatum)
+ # and probabilities pCtxSTN or pCtxiSN or pCtxdSN, respectively
+ else:
+ # In this case connections from each TVB proxy to TC_pop
+ # are equal to the Maith et al. optimized weights
+ tvb_spikeNet_model_builder.output_interfaces[-1]["weights"] = conn_scaling
+ # In this case coupling.a is already applied during computing TVB coupling.
+ # Therefore we scale only with model.G
+ tvb_spikeNet_model_builder.output_interfaces[-1]["transformer_params"] = \
+ {"scale_factor": scale_factor * tvb_spikeNet_model_builder.G}
+ # In total:
+ # From each TVB proxy node we get a total coupling rate scaled
+ # as (coupling.a * G STN_factor/wTVBSTN) * R_j, (j for STN or Striatum)
+ # Then, spikes generated from each TVB proxy are transferred via connections
+ # with weights wCtxSTN or wCtxiSN or wCtxdSN and
+ # probabilities pCtxSTN or pCtxiSN or pCtxdSN, respectively
+
+
+from tvb_multiscale.core.interfaces.base.transformers.models.red_wong_wang import ElephantSpikesRateRedWongWangExc
+
+tvb_spikeNet_model_builder.input_interfaces = []
+# TVB <-- ANNarchy:
+for src_pop, nodes, in zip(
+ # Source populations in ANNarchy:
+ [np.array(["TC_pop"])],
+ # Source regions indices in ANNarchy:
+ [vim_node_id]): # Thalamus
+ # TVB <- NetPyNE
+ tvb_spikeNet_model_builder.input_interfaces.append(
+ {"voi": np.array(["S", "R"]), # Target state variables in TVB
+ "populations": src_pop, # Source populations in NetPyNE
+ # This transformer not only converts spike counts to rates for state variable R,
+ # but also integrates the dS/dt to compute the respective S!:
+ "transformer": ElephantSpikesRateRedWongWangExc,
+ "transformer_params":
+ # Spike counts are converted to rates via:
+ # number_of_spikes / number_of_neurons_per_population / number_of_populations
+ # (mind that there are 2 populations in Striatum)
+ {"scale_factor": np.array([1.0]) / tvb_spikeNet_model_builder.N_E / len(src_pop),
+ # The integrator used to integrate dS/dt
+ "integrator":CONFIGURED.DEFAULT_TRANSFORMER_INTEGRATOR_MODEL(
+ dt=simulator.integrator.dt),
+ "state": np.zeros((2, len(nodes))), # initial condition
+ # Parameters of the dS/dt differential equation:
+ "tau_s": simulator.model.tau_s, # time constant of integrating S
+ "tau_r": np.array([10.0]), # time constant of integrating R to low pass filter it
+ "gamma": simulator.model.gamma},
+ "proxy_inds": np.array(nodes)}) # None means all here
+
+
+# ----------------------------------------------------------------------------------------------------------------
+# ----------------------------------------------------------------------------------------------------------------
+# ----------------------------------------------------------------------------------------------------------------
+
+# Configure and build:
+tvb_spikeNet_model_builder.configure()
+tvb_spikeNet_model_builder.print_summary_info_details(recursive=1)
+
+print("\noutput (TVB->NetPyNE coupling) interfaces' configurations:\n")
+display(tvb_spikeNet_model_builder.output_interfaces)
+
+print("\ninput (NetPyNE->TVB update) interfaces' configurations:\n")
+display(tvb_spikeNet_model_builder.input_interfaces)
+
+simulator = tvb_spikeNet_model_builder.build()
+
+# NetPyNE model is built in two steps. First need to create declarative-style specification for both spiking network itself and TVB-Netpyne proxy devides (interfaces).
+# Once it's done above using builders, network can be instantiated based on the specification:
+
+netpyne.simConfig.duration = simulation_length # TODO: do it properly! (find out why it gets lost)
+netpyne.simConfig.dt = spiking_network_builder.spiking_dt
+netpyne.allowSelfConns = True
+
+netpyne_network.netpyne_instance.instantiateNetwork()
+
+
+simulator.simulate_spiking_simulator = netpyne_network.netpyne_instance.run # set the method to run NetPyNE
+
+# simulator.print_summary_info(recursive=3)
+# simulator.print_summary_info_details(recursive=3)
+
+print("\n\noutput (TVB->NetPyNE coupling) interfaces:\n")
+simulator.output_interfaces.print_summary_info_details(recursive=2)
+
+print("\n\ninput (NetPyNE->TVB update) interfaces:\n")
+simulator.input_interfaces.print_summary_info_details(recursive=2)
+
+# -----------------------------------4. Compile network ---------------------------------------------------------
+# Compile the ANNarchy network...
+# tic_compile = time.time()
+# netpyne_network.configure()
+# print("Compiled! in %f min" % ((time.time() - tic_compile) / 60))
+
+netpyne_network.print_summary_info_details(recursive=2, connectivity=False)
+
+ # -----------------------------------5. Simulate and gather results-------------------------------------------------
+# simulation_length = 1500.0
+# transient = 500.0 # simulation_length/11
+# ...and simulate!
+
+simulator.configure() #
+
+simulation_length = \
+ np.ceil(simulation_length / simulator.synchronization_time) * simulator.synchronization_time
+
+advance_simulation_for_delayed_monitors_output = True
+if simulation_mode == "rs":
+ simulation_length1 = simulation_length
+else:
+ start_stimulus = np.ceil(start_stimulus / simulator.synchronization_time) * simulator.synchronization_time
+ simulation_length1 = start_stimulus
+ advance_simulation_for_delayed_monitors_output = False
+
+t_start = time.time()
+
+results = simulator.run(simulation_length=simulation_length1,
+ advance_simulation_for_delayed_monitors_output=advance_simulation_for_delayed_monitors_output
+ ) # 35.0 with stimulus application
+
+print("\nSimulated in %f secs!" % (time.time() - t_start))
+
+netpyne.finalize()
+
+from netpyne import sim
+popIds = [id for id in sim.net.pops.keys()][4] # only thalamus
+sa = sim.analysis
+#%matplotlib inline
+sim.analysis.plotConn(showFig=True, includePre=popIds, includePost=popIds, feature='weight');
+sim.analysis.plotConn(showFig=True, includePre=popIds, includePost=popIds, feature='prob');
+sim.analysis.plotConn(showFig=True, includePre=popIds, includePost=popIds, feature='probability');
\ No newline at end of file
diff --git a/tvb_multiscale/tvb_netpyne/interfaces/models/wilson_cowan.py b/tvb_multiscale/tvb_netpyne/interfaces/models/wilson_cowan.py
index b72389b3b..1f295eb1d 100644
--- a/tvb_multiscale/tvb_netpyne/interfaces/models/wilson_cowan.py
+++ b/tvb_multiscale/tvb_netpyne/interfaces/models/wilson_cowan.py
@@ -25,4 +25,6 @@ def default_output_config(self):
self.output_interfaces[0]["transformer_params"] = {"scale_factor": tuneable_factor * np.array([1.0])}
def default_input_config(self):
- WilsonCowanTVBSpikeNetInterfaceBuilder.default_input_config(self)
\ No newline at end of file
+ WilsonCowanTVBSpikeNetInterfaceBuilder.default_input_config(self)
+ # TODO: move to dedicated subclass
+ self.input_interfaces[1]['spiking_proxy_inds'] = [8, 361] # M1, cerebellar cortex
\ No newline at end of file
diff --git a/tvb_multiscale/tvb_netpyne/netpyne/module.py b/tvb_multiscale/tvb_netpyne/netpyne/module.py
index 5f437251c..9f348c486 100644
--- a/tvb_multiscale/tvb_netpyne/netpyne/module.py
+++ b/tvb_multiscale/tvb_netpyne/netpyne/module.py
@@ -175,13 +175,10 @@ def run(self, length):
tvbIterationEnd = self.time + length
def _(simTime): pass
if self.nextIntervalFuncCall:
- while (self.nextIntervalFuncCall < tvbIterationEnd):
- if self.time < sim.cfg.duration:
- sim.run.runForInterval(self.nextIntervalFuncCall - self.time, _)
- self.intervalFunc(self.time)
- self.nextIntervalFuncCall = self.time + self.interval
- else:
- break
+ while (self.nextIntervalFuncCall < min(tvbIterationEnd, sim.cfg.duration)):
+ sim.run.runForInterval(self.nextIntervalFuncCall - self.time, _)
+ self.intervalFunc(self.time)
+ self.nextIntervalFuncCall = self.time + self.interval
if tvbIterationEnd > self.time:
if self.time < sim.cfg.duration:
sim.run.runForInterval(tvbIterationEnd - self.time, _)
diff --git a/tvb_multiscale/tvb_netpyne/netpyne_models/models/default_exc_io_inh_i.py b/tvb_multiscale/tvb_netpyne/netpyne_models/models/default_exc_io_inh_i.py
deleted file mode 100644
index 8b6d5f8dc..000000000
--- a/tvb_multiscale/tvb_netpyne/netpyne_models/models/default_exc_io_inh_i.py
+++ /dev/null
@@ -1,203 +0,0 @@
-from tvb_multiscale.tvb_netpyne.netpyne_models.builders.base import NetpyneNetworkBuilder
-from tvb.contrib.scripts.utils.data_structures_utils import ensure_list
-from tvb_multiscale.tvb_netpyne.netpyne_models.builders.netpyne_templates import random_normal_weight, random_normal_tvb_weight, random_uniform_tvb_delay
-from collections import OrderedDict
-import numpy as np
-
-class DefaultExcIOInhIBuilder(NetpyneNetworkBuilder):
-
- def __init__(self, tvb_simulator={}, spiking_nodes_inds=[], netpyne_instance=None, config=None, logger=None):
- super(DefaultExcIOInhIBuilder, self).__init__(tvb_simulator, spiking_nodes_inds, netpyne_instance=netpyne_instance, config=config, logger=logger)
-
- self.scale_e = 1.2
- self.scale_i = 0.4
-
- self.w_ee = 0.02
- self.w_ei = 0.01
- self.w_ie = 0.01
- self.w_ii = 0.01
-
- self.delay = 5
-
- self.conn_spec_all_to_all = {"rule": "all_to_all"}
- self.conn_spec_prob_low = {"rule": {"prob": 0.01}}
- self.conn_spec_prob_high = {"rule": {"prob": 0.03}}
-
- self.params = {} # TODO: pass some info at least about transforming the default delay (because TVB support only delay as single number)
-
- def configure(self):
- from netpyne import specs
- netParams = specs.NetParams()
- cfg = specs.SimConfig()
-
- self.receptor_type_E = 'exc'
- self.receptor_type_I = 'inh'
- netParams.synMechParams[self.receptor_type_E] = {'mod': 'Exp2Syn', 'tau1': 0.8, 'tau2': 5.3, 'e': 0} # NMDA
- netParams.synMechParams[self.receptor_type_I] = {'mod': 'Exp2Syn', 'tau1': 0.6, 'tau2': 8.5, 'e': -75} # GABA
-
- PYRcell = {'secs': {}}
- PYRcell['secs']['soma'] = {'geom': {}, 'mechs': {}} # soma params dict
- PYRcell['secs']['soma']['geom'] = {'diam': 18.8, 'L': 18.8, 'Ra': 123.0} # soma geometry
- PYRcell['secs']['soma']['mechs']['hh'] = {'gnabar': 0.12, 'gkbar': 0.036, 'gl': 0.003, 'el': -70} # soma hh mechanism
- netParams.cellParams['PYR'] = PYRcell
-
- # simConfig.verbose = True
-
- cfg.recordTraces = {'V_soma':{'sec':'soma','loc':0.5,'var':'v'}} # Dict with traces to record
-
- cfg.recordStep = 0.1
- cfg.savePickle = False # Save params, network and sim output to pickle file
- cfg.saveJson = False
-
- super(DefaultExcIOInhIBuilder, self).configure(netParams, cfg, autoCreateSpikingNodes=True)
- self.global_coupling_scaling *= self.tvb_serial_sim.get("model.G", np.array([2.0]))[0].item()
- self.lamda = self.tvb_serial_sim.get("model.lamda", np.array([0.0]))[0].item()
-
- def proxy_node_synaptic_model_funcs(self):
- return {"E": lambda src_node, dst_node: self.receptor_type_E,
- "I": lambda src_node, dst_node: self.receptor_type_E}
-
- def set_defaults(self):
- self.set_populations()
- self.set_populations_connections()
- self.set_nodes_connections()
- self.set_output_devices()
- # self.set_input_devices()
-
-
- def set_populations(self):
- self.populations = [
- {"label": "E", "model": "PYR",
- "nodes": None, # None means "all"
- "params": self.params,
- "scale": self.scale_e},
- {"label": "I", "model": "PYR",
- "nodes": None, # None means "all"
- "params": self.params,
- # "params": netpyne_network_builder.params_I,
- "scale": self.scale_i}]
-
- # connections between E and I populations of same spiking node
-
- def set_populations_connections(self):
- self.populations_connections = [
- self.set_EE_populations_connections(),
- self.set_EI_populations_connections(),
- self.set_IE_populations_connections(),
- self.set_II_populations_connections()
- ]
-
- def set_EE_populations_connections(self):
- connections = \
- {"source": "E", "target": "E",
- # "synapse_model": self.default_populations_connection["synapse_model"], # TODO: here and below, is this needed or `receptor_type` below is just enough?
- "conn_spec": self.conn_spec_prob_low,
- "weight": self.weight_fun_ee,
- "delay": self.delay,
- "receptor_type": self.receptor_type_E, "nodes": None} # None means "all"
- # connections.update(self.pop_conns_EE)
- return connections
-
- def set_EI_populations_connections(self):
- connections = \
- {"source": "E", "target": "I",
- "conn_spec": self.conn_spec_prob_low,
- "weight": self.weight_fun_ei,
- "delay": self.delay,
- "receptor_type": self.receptor_type_E, "nodes": None} # None means "all"
- return connections
-
- def set_IE_populations_connections(self):
- connections = \
- {"source": "I", "target": "E",
- "conn_spec": self.conn_spec_prob_high,
- "weight": self.weight_fun_ie,
- "delay": self.delay,
- "receptor_type": self.receptor_type_I, "nodes": None} # None means "all"
- return connections
-
- def set_II_populations_connections(self):
- connections = \
- {"source": "I", "target": "I",
- "conn_spec": self.conn_spec_prob_high,
- "weight": self.weight_fun_ii,
- "delay": self.delay,
- "receptor_type": self.receptor_type_I, "nodes": None} # None means "all"
- return connections
-
- def weight_fun_ee(self, nodeId):
- return self.within_node_weight_fun(self.w_ee, nodeId)
-
- def weight_fun_ei(self, nodeId):
- return self.within_node_weight_fun(self.w_ei, nodeId)
-
- def weight_fun_ie(self, nodeId):
- return self.within_node_weight_fun(self.w_ie, nodeId)
-
- def weight_fun_ii(self, nodeId):
- return self.within_node_weight_fun(self.w_ii, nodeId)
-
- def within_node_weight_fun(self, weight, nodeId):
- weights = ensure_list(weight)
- weight = weights[nodeId] if (len(weights) > 1) else weights[0]
- return random_normal_weight(weight)
-
- # connections between population of different spiking nodes
-
- def set_nodes_connections(self):
- self.nodes_connections = [
- {"source": "E", "target": "E",
- "conn_spec": self.conn_spec_prob_low,
- "weight": self.tvb_weight_fun,
- "delay": self.tvb_delay_fun,
- # Each region emits spikes in its own port:
- "receptor_type": self.receptor_type_E, "source_nodes": None, "target_nodes": None} # None means "all"
- ]
- if self.lamda > 0:
- self.nodes_connections.append({
- "source": "E", "target": "I",
- "conn_spec": self.conn_spec_prob_low,
- # using lamda to scale connectivity weights (or alternatively, it can be used to downscale connection probability in 'conn_spec' above):
- "weight": lambda source_node, target_node: self.tvb_weight_fun(source_node, target_node, self.lamda),
- "delay": self.tvb_delay_fun,
- # Each region emits spikes in its own port:
- "receptor_type": self.receptor_type_E, "source_nodes": None, "target_nodes": None # None means "all"
- })
-
- # output devices
-
- def set_output_devices(self):
- # Creating devices to be able to observe NetPyNE activity:
- # Labels have to be different
- self.output_devices = [self.set_spike_recorder()]
-
- def set_spike_recorder(self):
- connections = OrderedDict()
- # label <- target population
- connections["E"] = "E"
- connections["I"] = "I"
- params = self.config.NETPYNE_OUTPUT_DEVICES_PARAMS_DEF["spike_recorder"].copy()
- # params["record_to"] = self.output_devices_record_to
- device = {"model": "spike_recorder", "params": params,
- # "neurons_fun": lambda node_id, population: population[:np.minimum(100, len(population))],
- "connections": connections, "nodes": None} # None means all here
- # device.update(self.spike_recorder)
- return device
-
- def tvb_weight_fun(self, source_node, target_node, lamda=None, sigma=0.1):
- scale = self.global_coupling_scaling * self.netpyne_synaptic_weight_scale
- if lamda:
- scale *= lamda
- return random_normal_tvb_weight(source_node, target_node, self.tvb_weights, scale=scale, sigma=sigma)
-
- def tvb_delay_fun(self, source_node, target_node, low=None, high=None, sigma=0.1):
- if low is None:
- low = self.tvb_dt
- if high is None:
- high = 2*low
- return random_uniform_tvb_delay(source_node, target_node, self.tvb_delays, low, high, sigma)
-
- def build(self, set_defaults=True):
- if set_defaults:
- self.set_defaults()
- return super(DefaultExcIOInhIBuilder, self).build()
\ No newline at end of file
diff --git a/tvb_multiscale/tvb_netpyne/netpyne_models/models/thalamic_VIM_exc_io_inh_i.py b/tvb_multiscale/tvb_netpyne/netpyne_models/models/thalamic_VIM_exc_io_inh_i.py
new file mode 100644
index 000000000..beabb6c19
--- /dev/null
+++ b/tvb_multiscale/tvb_netpyne/netpyne_models/models/thalamic_VIM_exc_io_inh_i.py
@@ -0,0 +1,147 @@
+from tvb_multiscale.tvb_netpyne.netpyne_models.builders.base import NetpyneNetworkBuilder
+from tvb_multiscale.tvb_netpyne.netpyne_models.builders.netpyne_templates import random_normal_weight, random_normal_tvb_weight, random_uniform_tvb_delay
+from collections import OrderedDict
+import numpy as np
+
+class ThalamicVIMBuilder(NetpyneNetworkBuilder):
+
+ def __init__(self, tvb_simulator={}, spiking_nodes_inds=[], netpyne_instance=None, config=None, logger=None):
+ super(ThalamicVIMBuilder, self).__init__(tvb_simulator, spiking_nodes_inds, netpyne_instance=netpyne_instance, config=config, logger=logger)
+
+ self.scale_e = 1.2
+ self.scale_i = 0.4
+
+ def configure(self):
+ try:
+ # presuming that network is cloned to ./tvb_multiscale/tvb_netpyne/netpyne_models/models/thalamic_VIM_ET
+ from .thalamic_VIM_ET.src.netParams import netParams
+ from .thalamic_VIM_ET.src.cfg import cfg
+ except ModuleNotFoundError as e:
+ raise Exception(f'Spiking network should be cloned locally and imported here as `netParams` and `cfg` (error: {e})')
+
+ # for external stimuli
+ self.synMechE = 'exc'
+ self.synMechI = 'inh'
+ netParams.synMechParams[self.synMechE] = {'mod': 'Exp2Syn', 'tau1': 0.8, 'tau2': 5.3, 'e': 0} # NMDA
+ netParams.synMechParams[self.synMechI] = {'mod': 'Exp2Syn', 'tau1': 0.6, 'tau2': 8.5, 'e': -75} # GABA
+
+ def intervalFunc(simTime):
+ pass
+ cfg.interval = 1
+ cfg.intervalFunc = intervalFunc
+
+ super(ThalamicVIMBuilder, self).configure(netParams, cfg, autoCreateSpikingNodes=False)
+ self.global_coupling_scaling *= self.tvb_serial_sim.get("model.G", np.array([2.0]))[0].item()
+ self.lamda = self.tvb_serial_sim.get("model.lamda", np.array([0.0]))[0].item()
+
+ def proxy_node_synaptic_model_funcs(self):
+ return {"E": lambda src_node, dst_node: self.synMechE,
+ "I": lambda src_node, dst_node: self.synMechI}
+
+ def set_defaults(self):
+ self.set_populations()
+ # self.set_output_devices()
+
+ def set_populations(self):
+
+ spikingPopsE = {
+ self.primary_motor_cortex_R: 'PY_pop',
+ #self.brainstem: 'ION_pop',
+ #self.thalamus_R: 'TC_pop',
+ #self.cerebellar_cortex_L: 'GrC_pop'
+ }
+
+ spikingPopsI = {
+ self.primary_motor_cortex_R: 'FSI_pop',
+ #self.cerebellar_cortex_L: 'PC_pop'
+ }
+
+ self.populations = [
+ {
+ "label": "E", "model": None, # 'model' not used, as it is needed only in case of autoCreateSpikingNodes is True
+ "nodes": None, # None means "all"
+ "params": lambda node_id: {"global_label": spikingPopsE[node_id]},# {"global_label": "PY_pop"}, # population of spiking network to be interfaced with excitatory population of TVB node
+ "scale": self.scale_e,
+ },
+ {
+ "label": "I", "model": None,
+ "nodes": [self.primary_motor_cortex_R],#, self.cerebellar_cortex_L],
+ "params": lambda node_id: {"global_label": spikingPopsI[node_id]}, # population of spiking network to be interfaced with inhibitory population of TVB node
+ "scale": self.scale_i
+ },
+ ]
+ # self.populations_connections = [] #don't need?
+
+
+ # TODO: seems to be redundant, having input_interfaces in TVBNetpyneInterfaceBuilder?
+ # def set_output_devices(self):
+ # # Creating devices to be able to observe NetPyNE activity:
+ # # Labels have to be different
+ # self.output_devices = [self.set_spike_recorder()]
+
+ # def set_spike_recorder(self):
+ # connections = OrderedDict()
+ # # Keys are arbitrary. Value is either 'E' or 'I' for this model, and will be matched with items in self.populations dict
+ # connections["E"] = "E"
+ # connections["I"] = "I"
+ # params = self.config.NETPYNE_OUTPUT_DEVICES_PARAMS_DEF["spike_recorder"].copy()
+ # # params["record_to"] = self.output_devices_record_to
+ # device = {"model": "spike_recorder", "params": params,
+ # # "neurons_fun": lambda node_id, population: population[:np.minimum(100, len(population))],
+ # "connections": connections, "nodes": None} # None means all here
+ # # device.update(self.spike_recorder)
+ # return device
+
+ def build(self, set_defaults=True):
+ if set_defaults:
+ self.set_defaults()
+ return super(ThalamicVIMBuilder, self).build()
+
+class WilsonCowanThalamicVIMBuilder(ThalamicVIMBuilder):
+
+ def __init__(self, tvb_simulator={}, spiking_nodes_inds=[], netpyne_instance=None,
+ config=None, logger=None):
+ super(WilsonCowanThalamicVIMBuilder, self).__init__(tvb_simulator, spiking_nodes_inds, netpyne_instance, config, logger)
+
+ def set_defaults(self, **kwargs):
+ super(WilsonCowanThalamicVIMBuilder, self).set_defaults()
+
+class WilsonCowanTVBNetpyneInterfaceBuilder(WilsonCowanThalamicVIMBuilder):
+ def __init__(self):
+ super(WilsonCowanTVBNetpyneInterfaceBuilder, self).__init__()
+
+ def a():
+ # Build a TVB-NetPyNE interface with all the appropriate connections between the
+ # TVB and NetPyNE modelled regions
+
+ from tvb_multiscale.tvb_netpyne.interfaces.builders import TVBNetpyneInterfaceBuilder
+ tvb_netpyne_interface_builder = TVBNetpyneInterfaceBuilder() # non opinionated builder
+
+ tvb_netpyne_interface_builder.config = config
+ tvb_netpyne_interface_builder.tvb_cosimulator = simulator
+ tvb_netpyne_interface_builder.spiking_network = netpyne_network
+ # This can be used to set default tranformer and proxy models:
+ tvb_netpyne_interface_builder.model = INTERFACE_MODEL # "RATE"
+ tvb_netpyne_interface_builder.input_flag = True # If True, NetPyNE->TVB update will be implemented
+ tvb_netpyne_interface_builder.output_flag = True # If True, TVB->NetPyNE coupling will be implemented
+ # If default_coupling_mode = "TVB", large scale coupling towards spiking regions is computed in TVB
+ # and then applied with no time delay via a single "TVB proxy node" / NetPyNE device for each spiking region,
+ # "1-to-1" TVB->NetPyNE coupling.
+ # If any other value, we need 1 "TVB proxy node" / NetPyNE device for each TVB sender region node, and
+ # large-scale coupling for spiking regions is computed in NetPyNE,
+ # taking into consideration the TVB connectome weights and delays,
+ # in this "1-to-many" TVB->NetPyNE coupling.
+ tvb_netpyne_interface_builder.default_coupling_mode = INTERFACE_COUPLING_MODE # "spikeNet" # "TVB"
+
+ # Number of neurons per population to be used to compute population mean instantaneous firing rates:
+ tvb_netpyne_interface_builder.N_E = int(netpyne_model_builder.population_order * scale_e)
+ tvb_netpyne_interface_builder.N_I = int(netpyne_model_builder.population_order * scale_i)
+
+ tvb_netpyne_interface_builder.proxy_inds = NETPYNE_NODES_INDS
+ # Set exclusive_nodes = True (Default) if the spiking regions substitute for the TVB ones:
+ tvb_netpyne_interface_builder.exclusive_nodes = True
+
+ tvb_netpyne_interface_builder.synaptic_weight_scale_func = synaptic_weight_scale_func
+
+ tvb_netpyne_interface_builder.output_interfaces = []
+ tvb_netpyne_interface_builder.input_interfaces = []