diff --git a/doc/htmldoc/examples/index.rst b/doc/htmldoc/examples/index.rst index 88a93aed4a..5ef22afc0a 100644 --- a/doc/htmldoc/examples/index.rst +++ b/doc/htmldoc/examples/index.rst @@ -80,6 +80,17 @@ PyNEST examples * :doc:`../auto_examples/Potjans_2014/index` + + + .. grid-item-card:: EI clustered network (Rostami et al) + :img-top: ../static/img/pynest/EI_clustered_network_schematic.png + + :doc:`../auto_examples/EI_clustered_network/index` + + + +.. grid:: 1 1 2 3 + .. grid-item-card:: GLIF (from Allen institute) :img-top: ../static/img/pynest/glif_cond.png @@ -87,10 +98,6 @@ PyNEST examples * :doc:`../auto_examples/glif_psc_neuron` * :doc:`../auto_examples/glif_psc_double_alpha_neuron` - - -.. grid:: 1 1 2 3 - .. grid-item-card:: Compartmental neurons :img-top: ../static/img/pynest/dendritic_synapse_conductances.png @@ -98,13 +105,6 @@ PyNEST examples * :doc:`../auto_examples/compartmental_model/two_comps` - .. grid-item-card:: Rate neurons - :img-top: ../static/img/pynest/rate_neuron.png - - * :doc:`../auto_examples/lin_rate_ipn_network` - * :doc:`../auto_examples/rate_neuron_dm` - - .. grid-item-card:: GIF (from Gerstner lab) :img-top: ../static/img/pynest/gif_pop.png @@ -116,6 +116,13 @@ PyNEST examples .. grid:: 1 1 2 3 + .. grid-item-card:: Rate neurons + :img-top: ../static/img/pynest/rate_neuron.png + + * :doc:`../auto_examples/lin_rate_ipn_network` + * :doc:`../auto_examples/rate_neuron_dm` + + .. grid-item-card:: Hodgkin-Huxley :img-top: ../static/img/pynest/hh_phase.png @@ -127,14 +134,14 @@ PyNEST examples * :doc:`../auto_examples/BrodyHopfield` +.. grid:: 1 1 2 3 + .. grid-item-card:: Brette and Gerstner :img-top: ../static/img/pynest/brette_gerstner2c.png * :doc:`../auto_examples/brette_gerstner_fig_2c` * :doc:`../auto_examples/brette_gerstner_fig_3d` -.. grid:: 1 1 2 3 - .. grid-item-card:: Precise spiking :img-top: ../static/img/pynest/precisespiking.png @@ -146,11 +153,6 @@ PyNEST examples * :doc:`../auto_examples/CampbellSiegert` - .. grid-item-card:: SONATA network - :img-top: ../static/img/300_pointneurons.png - - * :doc:`../auto_examples/sonata_example/sonata_network` - .. grid:: 1 1 2 3 @@ -242,6 +244,11 @@ PyNEST examples .. grid:: 1 1 2 3 + .. grid-item-card:: SONATA network + :img-top: ../static/img/300_pointneurons.png + + * :doc:`../auto_examples/sonata_example/sonata_network` + .. grid-item-card:: HPC benchmark :img-top: ../static/img/nest_logo-faded.png @@ -342,6 +349,7 @@ PyNEST examples ../auto_examples/astrocytes/astrocyte_interaction ../auto_examples/astrocytes/astrocyte_small_network ../auto_examples/astrocytes/astrocyte_brunel + ../auto_examples/EI_clustered_network/index ../auto_examples/eprop_plasticity/index .. toctree:: diff --git a/doc/htmldoc/static/img/pynest/EI_clustered_network_schematic.png b/doc/htmldoc/static/img/pynest/EI_clustered_network_schematic.png new file mode 100644 index 0000000000..26e0fe5527 Binary files /dev/null and b/doc/htmldoc/static/img/pynest/EI_clustered_network_schematic.png differ diff --git a/pynest/examples/EI_clustered_network/README.rst b/pynest/examples/EI_clustered_network/README.rst new file mode 100644 index 0000000000..4f7d5a43d7 --- /dev/null +++ b/pynest/examples/EI_clustered_network/README.rst @@ -0,0 +1,69 @@ +EI-clustered circuit model +========================== + +This is PyNEST implementation of the EI-clustered circuit model described by Rostami et al. [1]_. + +.. figure:: /static/img/pynest/EI_clustered_network_schematic.png + :alt: EI-clustered circuit model. + + Schematic of the EI-clustered circuit model. The network consists of `n_clusters` with one excitatory and one inhibitory population each. + +Citing this code +---------------- + +If you use this code, we ask you to cite the paper by Rostami et al. [1]_ and the NEST release on Zenodo. + +File structure +-------------- + +* :doc:`run_simulation.py `: an example script to try out the EI-clustered circuit model +* :doc:`network.py `: the main ``Network`` class with functions to build and simulate the network +* :doc:`helper.py `: helper functions for calculation of synaptic weights and currents and plot function for raster plots +* :doc:`network_params.py `: network and neuron parameters +* :doc:`stimulus_params.py `: parameters for optional external stimulation +* :doc:`sim_params.py `: simulation parameters + +Running the simulation +---------------------- + +.. code-block:: bash + + python run_simulation.py + +A raster plot of the network activity is saved as ``clustered_ei_raster.png``. + +The code can be parallelized by using multiple threads during the NEST simulation. +This can be done by setting the parameter ``n_vp`` in the ``run_simulation.py`` script. + +Contributions to this PyNEST model implementation +------------------------------------------------- + +2023: initial version of code and documentation by Felix J. Schmitt, Vahid Rostami and Martin Nawrot. + +Acknowledgments +--------------- + +Funding for the study by Rostami et al. [1]_: This work was supported by the German Research Foundation (DFG), +in parts through the Collaborative Research Center ’Motor Control in Health and Disease’ +(DFG-SFB 1451, Project-ID 431549029) and under the Institutional Strategy of the University of Cologne within the +German Excellence Initiative (DFG-ZUK 81/1) and in parts through the DFG graduate school +’Neural Circuit Analysis’ (DFG-RTG 1960, ID 365082554) and through the European Union’s Horizon 2020 Framework +Programme for Research and Innovation under grant agreement number 945539 (Human Brain Project SGA3). +The figure is created with BioRender.com. + +Other implementations of the EI-clustered model +----------------------------------------------- + +A `GeNN version `__ by Felix J. Schmitt, Vahid Rostami and Martin Nawrot [2]_. + +References +---------- + +.. [1] Rostami V, Rost T, Riehle A, van Albada SJ and Nawrot MP. 2020. + Excitatory and inhibitory motor cortical clusters account for balance, variability, and task performance. + bioRxiv 2020.02.27.968339. DOI: `10.1101/2020.02.27.968339 `__. + + +.. [2] Schmitt FJ, Rostami V and Nawrot MP. 2023. + Efficient parameter calibration and real-time simulation of large-scale spiking neural networks with GeNN + and NEST. Front. Neuroinform. 17:941696. DOI: `10.3389/fninf.2023.941696 `__. diff --git a/pynest/examples/EI_clustered_network/helper.py b/pynest/examples/EI_clustered_network/helper.py new file mode 100644 index 0000000000..902535849e --- /dev/null +++ b/pynest/examples/EI_clustered_network/helper.py @@ -0,0 +1,190 @@ +# -*- coding: utf-8 -*- +# +# helper.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +"""PyNEST EI-clustered network: Helper Functions +------------------------------------------------ + +Helper functions to calculate synaptic weights to construct +random balanced networks and plot raster plot with color groups. +""" + +import matplotlib.pyplot as plt +import numpy as np + + +def postsynaptic_current_to_potential(tau_m, tau_syn, c_m=1.0, e_l=0.0): + """Maximum post-synaptic potential amplitude + for exponential synapses and a synaptic efficacy J of 1 pA. + + Parameters + ---------- + tau_m: float + Membrane time constant [ms] + tau_syn: float + Synapse time constant [ms] + c_m: float (optional) + Membrane capacity [pF] (default: 1.0) + e_l: float (optional) + Resting potential [mV] (default: 0.0) + + Returns + ------- + float + maximum psp amplitude (for a 1 pA current spike) [mV] + """ + # time of maximum deflection of the psp [ms] + tmax = np.log(tau_syn / tau_m) / (1 / tau_m - 1 / tau_syn) + # we assume here the current spike is 1 pA, otherwise [mV/pA] + pre = tau_m * tau_syn / c_m / (tau_syn - tau_m) + return (e_l - pre) * np.exp(-tmax / tau_m) + pre * np.exp(-tmax / tau_syn) + + +def calculate_RBN_weights(params): + """Calculate synaptic weights for a random balanced network. + + The synaptic weights are calculated according to the method + described in Eqs. 7-10 [1]_. + + References + ---------- + .. [1] Rostami V, Rost T, Riehle A, van Albada SJ and Nawrot MP. 2020. + Excitatory and inhibitory motor cortical clusters account for balance, + variability, and task performance. bioRxiv 2020.02.27.968339. + DOI: `10.1101/2020.02.27.968339 + `__. + + Parameters + ---------- + params: dict + Dictionary of network parameters + + Returns + ------- + ndarray + synaptic weights 2x2 matrix [[EE, EI], [IE, II]] + """ + + N_E = params.get("N_E") # excitatory units + N_I = params.get("N_I") # inhibitory units + N = N_E + N_I # total units + + E_L = params.get("E_L") + V_th_E = params.get("V_th_E") # threshold voltage + V_th_I = params.get("V_th_I") + + tau_E = params.get("tau_E") + tau_I = params.get("tau_I") + + tau_syn_ex = params.get("tau_syn_ex") + tau_syn_in = params.get("tau_syn_in") + + gei = params.get("gei") + gii = params.get("gii") + gie = params.get("gie") + + amp_EE = postsynaptic_current_to_potential(tau_E, tau_syn_ex) + amp_EI = postsynaptic_current_to_potential(tau_E, tau_syn_in) + amp_IE = postsynaptic_current_to_potential(tau_I, tau_syn_ex) + amp_II = postsynaptic_current_to_potential(tau_I, tau_syn_in) + + baseline_conn_prob = params.get("baseline_conn_prob") # connection probs + + js = np.zeros((2, 2)) + K_EE = N_E * baseline_conn_prob[0, 0] + js[0, 0] = (V_th_E - E_L) * (K_EE**-0.5) * N**0.5 / amp_EE + js[0, 1] = -gei * js[0, 0] * baseline_conn_prob[0, 0] * N_E * amp_EE / (baseline_conn_prob[0, 1] * N_I * amp_EI) + K_IE = N_E * baseline_conn_prob[1, 0] + js[1, 0] = gie * (V_th_I - E_L) * (K_IE**-0.5) * N**0.5 / amp_IE + js[1, 1] = -gii * js[1, 0] * baseline_conn_prob[1, 0] * N_E * amp_IE / (baseline_conn_prob[1, 1] * N_I * amp_II) + return js + + +def rheobase_current(tau_m, e_l, v_th, c_m): + """Rheobase current for membrane time constant and resting potential. + + Parameters + ---------- + tau_m: float + Membrane time constant [ms] + e_l: float + Resting potential [mV] + v_th: float + Threshold potential [mV] + c_m: float + Membrane capacity [pF] + + Returns + ------- + float + rheobase current [pA] + """ + return (v_th - e_l) * c_m / tau_m + + +def raster_plot(spiketimes, tlim=None, colorgroups=None, ax=None, markersize=0.5): + """Raster plot of spiketimes. + + Plots raster plot of spiketimes withing given time limits and + colors neurons according to colorgroups. + + Parameters + ---------- + spiketimes: ndarray + 2D array [2xN_Spikes] + of spiketimes with spiketimes in row 0 and neuron IDs in row 1. + tlim: list of floats (optional) + Time limits of plot: [tmin, tmax], + if None: [min(spiketimes), max(spiketimes)] + colorgroups: list of tuples (optional) + Each element is a tuple in the format + (color, start_neuron_ID, stop_neuron_ID]) for coloring neurons in + given range, if None is given all neurons are black. + ax: axis object (optional) + If None a new figure is created, + else the plot is added to the given axis. + markersize: float (optional) + Size of markers. (default: 0.5) + + Returns + ------- + axis object + Axis object with raster plot. + """ + if ax is None: + fig, ax = plt.subplots() + if tlim is None: + tlim = [min(spiketimes[0]), max(spiketimes[0])] + if colorgroups is None: + colorgroups = [("k", 0, max(spiketimes[1]))] + for color, start, stop in colorgroups: + ax.plot( + spiketimes[0][np.logical_and(spiketimes[1] >= start, spiketimes[1] < stop)], + spiketimes[1][np.logical_and(spiketimes[1] >= start, spiketimes[1] < stop)], + color=color, + marker=".", + linestyle="None", + markersize=markersize, + ) + ax.set_xlim(tlim) + ax.set_ylim([0, max(spiketimes[1])]) + ax.set_xlabel("Time [ms]") + ax.set_ylabel("Neuron ID") + return ax diff --git a/pynest/examples/EI_clustered_network/network.py b/pynest/examples/EI_clustered_network/network.py new file mode 100644 index 0000000000..306014db17 --- /dev/null +++ b/pynest/examples/EI_clustered_network/network.py @@ -0,0 +1,824 @@ +# -*- coding: utf-8 -*- +# +# network.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +"""PyNEST EI-clustered network: Network Class +--------------------------------------------- + +``ClusteredNetwork`` class with functions to build and simulate +the EI-clustered network. +""" + +import pickle + +import helper +import nest +import numpy as np + + +class ClusteredNetwork: + """EI-clustered network objeect to build and simulate the network. + + Provides functions to create neuron populations, + stimulation devices and recording devices for an + EI-clustered network and setups the simulation in + NEST (v3.x). + + Attributes + ---------- + _params: dict + Dictionary with parameters used to construct network. + _populations: list + List of neuron population groups. + _recording_devices: list + List of recording devices. + _currentsources: list + List of current sources. + _model_build_pipeline: list + List of functions to build the network. + """ + + def __init__(self, sim_dict, net_dict, stim_dict): + """Initialize the ClusteredNetwork object. + + Parameters are given and explained in the files network_params.py, + sim_params.py and stimulus_params.py. + + Parameters + ---------- + sim_dict: dict + Dictionary with simulation parameters. + net_dict: dict + Dictionary with network parameters. + stim_dict: dict + Dictionary with stimulus parameters. + """ + + # merge dictionaries of simulation, network and stimulus parameters + self._params = {**sim_dict, **net_dict, **stim_dict} + # list of neuron population groups [E_pops, I_pops] + self._populations = [] + self._recording_devices = [] + self._currentsources = [] + self._model_build_pipeline = [ + self.setup_nest, + self.create_populations, + self.create_stimulation, + self.create_recording_devices, + self.connect, + ] + + if self._params["clustering"] == "weight": + jep = self._params["rep"] + jip = 1.0 + (jep - 1) * self._params["rj"] + self._params["jplus"] = np.array([[jep, jip], [jip, jip]]) + elif self._params["clustering"] == "probabilities": + pep = self._params["rep"] + pip = 1.0 + (pep - 1) ** self._params["rj"] + self._params["pplus"] = np.array([[pep, pip], [pip, pip]]) + else: + raise ValueError("Clustering type not recognized") + + def setup_nest(self): + """Initializes the NEST kernel. + + Reset the NEST kernel and pass parameters to it. + Updates randseed of parameters to the actual + used one if none is supplied. + """ + + nest.ResetKernel() + nest.set_verbosity("M_WARNING") + nest.local_num_threads = self._params.get("n_vp", 4) + nest.resolution = self._params.get("dt") + self._params["randseed"] = self._params.get("randseed") + nest.rng_seed = self._params.get("randseed") + + def create_populations(self): + """Create all neuron populations. + + n_clusters excitatory and inhibitory neuron populations + with the parameters of the network are created. + """ + + # make sure number of clusters and units are compatible + if self._params["N_E"] % self._params["n_clusters"] != 0: + raise ValueError("N_E must be a multiple of Q") + if self._params["N_I"] % self._params["n_clusters"] != 0: + raise ValueError("N_E must be a multiple of Q") + if self._params["neuron_type"] != "iaf_psc_exp": + raise ValueError("Model only implemented for iaf_psc_exp neuron model") + + if self._params["I_th_E"] is None: + I_xE = self._params["I_xE"] # I_xE is the feed forward excitatory input in pA + else: + I_xE = self._params["I_th_E"] * helper.rheobase_current( + self._params["tau_E"], self._params["E_L"], self._params["V_th_E"], self._params["C_m"] + ) + + if self._params["I_th_I"] is None: + I_xI = self._params["I_xI"] + else: + I_xI = self._params["I_th_I"] * helper.rheobase_current( + self._params["tau_I"], self._params["E_L"], self._params["V_th_I"], self._params["C_m"] + ) + + E_neuron_params = { + "E_L": self._params["E_L"], + "C_m": self._params["C_m"], + "tau_m": self._params["tau_E"], + "t_ref": self._params["t_ref"], + "V_th": self._params["V_th_E"], + "V_reset": self._params["V_r"], + "I_e": I_xE + if self._params["delta_I_xE"] == 0 + else I_xE * nest.random.uniform(1 - self._params["delta_I_xE"] / 2, 1 + self._params["delta_I_xE"] / 2), + "tau_syn_ex": self._params["tau_syn_ex"], + "tau_syn_in": self._params["tau_syn_in"], + "V_m": self._params["V_m"] + if not self._params["V_m"] == "rand" + else self._params["V_th_E"] - 20 * nest.random.lognormal(0, 1), + } + I_neuron_params = { + "E_L": self._params["E_L"], + "C_m": self._params["C_m"], + "tau_m": self._params["tau_I"], + "t_ref": self._params["t_ref"], + "V_th": self._params["V_th_I"], + "V_reset": self._params["V_r"], + "I_e": I_xI + if self._params["delta_I_xE"] == 0 + else I_xI * nest.random.uniform(1 - self._params["delta_I_xE"] / 2, 1 + self._params["delta_I_xE"] / 2), + "tau_syn_ex": self._params["tau_syn_ex"], + "tau_syn_in": self._params["tau_syn_in"], + "V_m": self._params["V_m"] + if not self._params["V_m"] == "rand" + else self._params["V_th_I"] - 20 * nest.random.lognormal(0, 1), + } + + # iaf_psc_exp allows stochasticity, if not used - don't supply the parameters and use + # iaf_psc_exp as deterministic model + if (self._params.get("delta") is not None) and (self._params.get("rho") is not None): + E_neuron_params["delta"] = self._params["delta"] + I_neuron_params["delta"] = self._params["delta"] + E_neuron_params["rho"] = self._params["rho"] + I_neuron_params["rho"] = self._params["rho"] + + # create the neuron populations + pop_size_E = self._params["N_E"] // self._params["n_clusters"] + pop_size_I = self._params["N_I"] // self._params["n_clusters"] + E_pops = [ + nest.Create(self._params["neuron_type"], n=pop_size_E, params=E_neuron_params) + for _ in range(self._params["n_clusters"]) + ] + I_pops = [ + nest.Create(self._params["neuron_type"], n=pop_size_I, params=I_neuron_params) + for _ in range(self._params["n_clusters"]) + ] + + self._populations = [E_pops, I_pops] + + def connect(self): + """Connect the excitatory and inhibitory populations with each other + in the EI-clustered scheme + + Raises + ------ + ValueError + If the clustering method is not recognized + """ + + if "clustering" not in self._params or self._params["clustering"] == "weight": + self.connect_weight() + elif self._params["clustering"] == "probabilities": + self.connect_probabilities() + else: + raise ValueError("Clustering method %s not implemented" % self._params["clustering"]) + + def connect_probabilities(self): + """Connect the clusters with a probability EI-cluster scheme + + Connects the excitatory and inhibitory populations with each other + in the EI-clustered scheme by increasing the probabilities of the + connections within the clusters and decreasing the probabilities of the + connections between the clusters. The weights are calculated so that + the total input to a neuron is balanced. + """ + + # self._populations[0] -> Excitatory population + # self._populations[1] -> Inhibitory population + + N = self._params["N_E"] + self._params["N_I"] # total units + # if js are not given compute them so that sqrt(K) spikes equal v_thr-E_L and rows are balanced + # if any of the js is nan or not given + if self._params.get("js") is None or np.isnan(self._params.get("js")).any(): + js = helper.calculate_RBN_weights(self._params) + js *= self._params["s"] + + if self._params["n_clusters"] > 1: + pminus = (self._params["n_clusters"] - self._params["pplus"]) / float(self._params["n_clusters"] - 1) + else: + self._params["pplus"] = np.ones((2, 2)) + pminus = np.ones((2, 2)) + + p_plus = self._params["pplus"] * self._params["baseline_conn_prob"] + p_minus = pminus * self._params["baseline_conn_prob"] + + # Connection probabilities within clusters can exceed 1. In this case, we iteratively split + # the connections in multiple synapse populations with probabilities < 1. + iterations = np.ones((2, 2), dtype=int) + # test if any of the probabilities is larger than 1 + if np.any(p_plus > 1): + print("The probability of some connections is larger than 1.") + print("Pre-splitting the connections in multiple synapse populations:") + printoptions = np.get_printoptions() + np.set_printoptions(precision=2, floatmode="fixed") + print("p_plus:\n", p_plus) + print("p_minus:\n", p_minus) + for i in range(2): + for j in range(2): + if p_plus[i, j] > 1: + iterations[i, j] = int(np.ceil(p_plus[i, j])) + p_plus[i, j] /= iterations[i, j] + print("\nPost-splitting the connections in multiple synapse populations:") + print("p_plus:\n", p_plus) + print("Number of synapse populations:\n", iterations) + np.set_printoptions(**printoptions) + + # define the synapses and connect the populations + + # Excitatory to excitatory neuron connections + j_ee = js[0, 0] / np.sqrt(N) + nest.CopyModel("static_synapse", "EE", {"weight": j_ee, "delay": self._params["delay"]}) + + if self._params["fixed_indegree"]: + K_EE_plus = int(p_plus[0, 0] * self._params["N_E"] / self._params["n_clusters"]) + print("K_EE+: ", K_EE_plus) + K_EE_minus = int(p_minus[0, 0] * self._params["N_E"] / self._params["n_clusters"]) + print("K_EE-: ", K_EE_minus) + conn_params_EE_plus = { + "rule": "fixed_indegree", + "indegree": K_EE_plus, + "allow_autapses": False, + "allow_multapses": True, + } + conn_params_EE_minus = { + "rule": "fixed_indegree", + "indegree": K_EE_minus, + "allow_autapses": False, + "allow_multapses": True, + } + + else: + conn_params_EE_plus = { + "rule": "pairwise_bernoulli", + "p": p_plus[0, 0], + "allow_autapses": False, + "allow_multapses": True, + } + conn_params_EE_minus = { + "rule": "pairwise_bernoulli", + "p": p_minus[0, 0], + "allow_autapses": False, + "allow_multapses": True, + } + for i, pre in enumerate(self._populations[0]): + for j, post in enumerate(self._populations[0]): + if i == j: + # same cluster + for n in range(iterations[0, 0]): + nest.Connect(pre, post, conn_params_EE_plus, "EE") + else: + nest.Connect(pre, post, conn_params_EE_minus, "EE") + + # Inhibitory to excitatory neuron connections + j_ei = js[0, 1] / np.sqrt(N) + nest.CopyModel("static_synapse", "EI", {"weight": j_ei, "delay": self._params["delay"]}) + + if self._params["fixed_indegree"]: + K_EI_plus = int(p_plus[0, 1] * self._params["N_I"] / self._params["n_clusters"]) + print("K_EI+: ", K_EI_plus) + K_EI_minus = int(p_minus[0, 1] * self._params["N_I"] / self._params["n_clusters"]) + print("K_EI-: ", K_EI_minus) + conn_params_EI_plus = { + "rule": "fixed_indegree", + "indegree": K_EI_plus, + "allow_autapses": False, + "allow_multapses": True, + } + conn_params_EI_minus = { + "rule": "fixed_indegree", + "indegree": K_EI_minus, + "allow_autapses": False, + "allow_multapses": True, + } + + else: + conn_params_EI_plus = { + "rule": "pairwise_bernoulli", + "p": p_plus[0, 1], + "allow_autapses": False, + "allow_multapses": True, + } + conn_params_EI_minus = { + "rule": "pairwise_bernoulli", + "p": p_minus[0, 1], + "allow_autapses": False, + "allow_multapses": True, + } + for i, pre in enumerate(self._populations[1]): + for j, post in enumerate(self._populations[0]): + if i == j: + # same cluster + for n in range(iterations[0, 1]): + nest.Connect(pre, post, conn_params_EI_plus, "EI") + else: + nest.Connect(pre, post, conn_params_EI_minus, "EI") + + # Excitatory to inhibitory neuron connections + j_ie = js[1, 0] / np.sqrt(N) + nest.CopyModel("static_synapse", "IE", {"weight": j_ie, "delay": self._params["delay"]}) + + if self._params["fixed_indegree"]: + K_IE_plus = int(p_plus[1, 0] * self._params["N_E"] / self._params["n_clusters"]) + print("K_IE+: ", K_IE_plus) + K_IE_minus = int(p_minus[1, 0] * self._params["N_E"] / self._params["n_clusters"]) + print("K_IE-: ", K_IE_minus) + conn_params_IE_plus = { + "rule": "fixed_indegree", + "indegree": K_IE_plus, + "allow_autapses": False, + "allow_multapses": True, + } + conn_params_IE_minus = { + "rule": "fixed_indegree", + "indegree": K_IE_minus, + "allow_autapses": False, + "allow_multapses": True, + } + + else: + conn_params_IE_plus = { + "rule": "pairwise_bernoulli", + "p": p_plus[1, 0], + "allow_autapses": False, + "allow_multapses": True, + } + conn_params_IE_minus = { + "rule": "pairwise_bernoulli", + "p": p_minus[1, 0], + "allow_autapses": False, + "allow_multapses": True, + } + for i, pre in enumerate(self._populations[0]): + for j, post in enumerate(self._populations[1]): + if i == j: + # same cluster + for n in range(iterations[1, 0]): + nest.Connect(pre, post, conn_params_IE_plus, "IE") + else: + nest.Connect(pre, post, conn_params_IE_minus, "IE") + + # Inhibitory to inhibitory neuron connections + j_ii = js[1, 1] / np.sqrt(N) + nest.CopyModel("static_synapse", "II", {"weight": j_ii, "delay": self._params["delay"]}) + + if self._params["fixed_indegree"]: + K_II_plus = int(p_plus[1, 1] * self._params["N_I"] / self._params["n_clusters"]) + print("K_II+: ", K_II_plus) + K_II_minus = int(p_minus[1, 1] * self._params["N_I"] / self._params["n_clusters"]) + print("K_II-: ", K_II_minus) + conn_params_II_plus = { + "rule": "fixed_indegree", + "indegree": K_II_plus, + "allow_autapses": False, + "allow_multapses": True, + } + conn_params_II_minus = { + "rule": "fixed_indegree", + "indegree": K_II_minus, + "allow_autapses": False, + "allow_multapses": True, + } + + else: + conn_params_II_plus = { + "rule": "pairwise_bernoulli", + "p": p_plus[1, 1], + "allow_autapses": False, + "allow_multapses": True, + } + conn_params_II_minus = { + "rule": "pairwise_bernoulli", + "p": p_minus[1, 1], + "allow_autapses": False, + "allow_multapses": True, + } + for i, pre in enumerate(self._populations[1]): + for j, post in enumerate(self._populations[1]): + if i == j: + # same cluster + for n in range(iterations[1, 1]): + nest.Connect(pre, post, conn_params_II_plus, "II") + else: + nest.Connect(pre, post, conn_params_II_minus, "II") + + def connect_weight(self): + """Connect the clusters with a weight EI-cluster scheme + + Connects the excitatory and inhibitory populations with + each other in the EI-clustered scheme by increasing the weights + of the connections within the clusters and decreasing the weights + of the connections between the clusters. The weights are calculated + so that the total input to a neuron is balanced. + """ + + # self._populations[0] -> Excitatory population + # self._populations[1] -> Inhibitory population + + N = self._params["N_E"] + self._params["N_I"] # total units + + # if js are not given compute them so that sqrt(K) spikes equal v_thr-E_L and rows are balanced + # if any of the js is nan or not given + if self._params.get("js") is None or np.isnan(self._params.get("js")).any(): + js = helper.calculate_RBN_weights(self._params) + js *= self._params["s"] + + # jminus is calculated so that row sums remain constant + if self._params["n_clusters"] > 1: + jminus = (self._params["n_clusters"] - self._params["jplus"]) / float(self._params["n_clusters"] - 1) + else: + self._params["jplus"] = np.ones((2, 2)) + jminus = np.ones((2, 2)) + + # define the synapses and connect the populations + + # Excitatory to excitatory neuron connections + j_ee = js[0, 0] / np.sqrt(N) + nest.CopyModel( + "static_synapse", + "EE_plus", + { + "weight": self._params["jplus"][0, 0] * j_ee, + "delay": self._params["delay"], + }, + ) + nest.CopyModel( + "static_synapse", + "EE_minus", + {"weight": jminus[0, 0] * j_ee, "delay": self._params["delay"]}, + ) + if self._params["fixed_indegree"]: + K_EE = int(self._params["baseline_conn_prob"][0, 0] * self._params["N_E"] / self._params["n_clusters"]) + print("K_EE: ", K_EE) + conn_params_EE = { + "rule": "fixed_indegree", + "indegree": K_EE, + "allow_autapses": False, + "allow_multapses": False, + } + + else: + conn_params_EE = { + "rule": "pairwise_bernoulli", + "p": self._params["baseline_conn_prob"][0, 0], + "allow_autapses": False, + "allow_multapses": False, + } + for i, pre in enumerate(self._populations[0]): + for j, post in enumerate(self._populations[0]): + if i == j: + # same cluster + nest.Connect(pre, post, conn_params_EE, "EE_plus") + else: + nest.Connect(pre, post, conn_params_EE, "EE_minus") + + # Inhibitory to excitatory neuron connections + j_ei = js[0, 1] / np.sqrt(N) + nest.CopyModel( + "static_synapse", + "EI_plus", + { + "weight": j_ei * self._params["jplus"][0, 1], + "delay": self._params["delay"], + }, + ) + nest.CopyModel( + "static_synapse", + "EI_minus", + {"weight": j_ei * jminus[0, 1], "delay": self._params["delay"]}, + ) + if self._params["fixed_indegree"]: + K_EI = int(self._params["baseline_conn_prob"][0, 1] * self._params["N_I"] / self._params["n_clusters"]) + print("K_EI: ", K_EI) + conn_params_EI = { + "rule": "fixed_indegree", + "indegree": K_EI, + "allow_autapses": False, + "allow_multapses": False, + } + else: + conn_params_EI = { + "rule": "pairwise_bernoulli", + "p": self._params["baseline_conn_prob"][0, 1], + "allow_autapses": False, + "allow_multapses": False, + } + for i, pre in enumerate(self._populations[1]): + for j, post in enumerate(self._populations[0]): + if i == j: + # same cluster + nest.Connect(pre, post, conn_params_EI, "EI_plus") + else: + nest.Connect(pre, post, conn_params_EI, "EI_minus") + + # Excitatory to inhibitory neuron connections + j_ie = js[1, 0] / np.sqrt(N) + nest.CopyModel( + "static_synapse", + "IE_plus", + { + "weight": j_ie * self._params["jplus"][1, 0], + "delay": self._params["delay"], + }, + ) + nest.CopyModel( + "static_synapse", + "IE_minus", + {"weight": j_ie * jminus[1, 0], "delay": self._params["delay"]}, + ) + + if self._params["fixed_indegree"]: + K_IE = int(self._params["baseline_conn_prob"][1, 0] * self._params["N_E"] / self._params["n_clusters"]) + print("K_IE: ", K_IE) + conn_params_IE = { + "rule": "fixed_indegree", + "indegree": K_IE, + "allow_autapses": False, + "allow_multapses": False, + } + else: + conn_params_IE = { + "rule": "pairwise_bernoulli", + "p": self._params["baseline_conn_prob"][1, 0], + "allow_autapses": False, + "allow_multapses": False, + } + for i, pre in enumerate(self._populations[0]): + for j, post in enumerate(self._populations[1]): + if i == j: + # same cluster + nest.Connect(pre, post, conn_params_IE, "IE_plus") + else: + nest.Connect(pre, post, conn_params_IE, "IE_minus") + + # Inhibitory to inhibitory neuron connections + j_ii = js[1, 1] / np.sqrt(N) + nest.CopyModel( + "static_synapse", + "II_plus", + { + "weight": j_ii * self._params["jplus"][1, 1], + "delay": self._params["delay"], + }, + ) + nest.CopyModel( + "static_synapse", + "II_minus", + {"weight": j_ii * jminus[1, 1], "delay": self._params["delay"]}, + ) + if self._params["fixed_indegree"]: + K_II = int(self._params["baseline_conn_prob"][1, 1] * self._params["N_I"] / self._params["n_clusters"]) + print("K_II: ", K_II) + conn_params_II = { + "rule": "fixed_indegree", + "indegree": K_II, + "allow_autapses": False, + "allow_multapses": False, + } + else: + conn_params_II = { + "rule": "pairwise_bernoulli", + "p": self._params["baseline_conn_prob"][1, 1], + "allow_autapses": False, + "allow_multapses": False, + } + for i, pre in enumerate(self._populations[1]): + for j, post in enumerate(self._populations[1]): + if i == j: + # same cluster + nest.Connect(pre, post, conn_params_II, "II_plus") + else: + nest.Connect(pre, post, conn_params_II, "II_minus") + + def create_stimulation(self): + """Create a current source and connect it to clusters.""" + if self._params["stim_clusters"] is not None: + stim_amp = self._params["stim_amp"] # amplitude of the stimulation current in pA + stim_starts = self._params["stim_starts"] # list of stimulation start times + stim_ends = self._params["stim_ends"] # list of stimulation end times + amplitude_values = [] + amplitude_times = [] + for start, end in zip(stim_starts, stim_ends): + amplitude_times.append(start + self._params["warmup"]) + amplitude_values.append(stim_amp) + amplitude_times.append(end + self._params["warmup"]) + amplitude_values.append(0.0) + self._currentsources = [nest.Create("step_current_generator")] + for stim_cluster in self._params["stim_clusters"]: + nest.Connect(self._currentsources[0], self._populations[0][stim_cluster]) + nest.SetStatus( + self._currentsources[0], + { + "amplitude_times": amplitude_times, + "amplitude_values": amplitude_values, + }, + ) + + def create_recording_devices(self): + """Creates a spike recorder + + Create and connect a spike recorder to all neuron populations + in self._populations. + """ + self._recording_devices = [nest.Create("spike_recorder")] + self._recording_devices[0].record_to = "memory" + + all_units = self._populations[0][0] + for E_pop in self._populations[0][1:]: + all_units += E_pop + for I_pop in self._populations[1]: + all_units += I_pop + nest.Connect(all_units, self._recording_devices[0], "all_to_all") # Spikerecorder + + def set_model_build_pipeline(self, pipeline): + """Set _model_build_pipeline + + Parameters + ---------- + pipeline: list + ordered list of functions executed to build the network model + """ + self._model_build_pipeline = pipeline + + def setup_network(self): + """Setup network in NEST + + Initializes NEST and creates + the network in NEST, ready to be simulated. + Functions saved in _model_build_pipeline are executed. + """ + for func in self._model_build_pipeline: + func() + + def simulate(self): + """Simulates network for a period of warmup+simtime""" + nest.Simulate(self._params["warmup"] + self._params["simtime"]) + + def get_recordings(self): + """Extract spikes from Spikerecorder + + Extract spikes form the Spikerecorder connected + to all populations created in create_populations. + Cuts the warmup period away and sets time relative to end of warmup. + Ids 1:N_E correspond to excitatory neurons, + N_E+1:N_E+N_I correspond to inhibitory neurons. + + Returns + ------- + spiketimes: ndarray + 2D array [2xN_Spikes] + of spiketimes with spiketimes in row 0 and neuron IDs in row 1. + """ + events = nest.GetStatus(self._recording_devices[0], "events")[0] + # convert them to the format accepted by spiketools + spiketimes = np.append(events["times"][None, :], events["senders"][None, :], axis=0) + spiketimes[1] -= 1 + # remove the pre warmup spikes + spiketimes = spiketimes[:, spiketimes[0] >= self._params["warmup"]] + spiketimes[0] -= self._params["warmup"] + return spiketimes + + def get_parameter(self): + """Get all parameters used to create the network. + Returns + ------- + dict + Dictionary with all parameters of the network and the simulation. + """ + return self._params + + def create_and_simulate(self): + """Create and simulate the EI-clustered network. + + Returns + ------- + spiketimes: ndarray + 2D array [2xN_Spikes] + of spiketimes with spiketimes in row 0 and neuron IDs in row 1. + """ + self.setup_network() + self.simulate() + return self.get_recordings() + + def get_firing_rates(self, spiketimes=None): + """Calculates the average firing rates of + all excitatory and inhibitory neurons. + + Calculates the firing rates of all excitatory neurons + and the firing rates of all inhibitory neurons + created by self.create_populations. + If spiketimes are not supplied, they get extracted. + + Parameters + ---------- + spiketimes: ndarray + 2D array [2xN_Spikes] of spiketimes + with spiketimes in row 0 and neuron IDs in row 1. + + Returns + ------- + tuple[float, float] + average firing rates of excitatory (0) + and inhibitory (1) neurons (spikes/s) + """ + if spiketimes is None: + spiketimes = self.get_recordings() + e_count = spiketimes[:, spiketimes[1] < self._params["N_E"]].shape[1] + i_count = spiketimes[:, spiketimes[1] >= self._params["N_E"]].shape[1] + e_rate = e_count / float(self._params["N_E"]) / float(self._params["simtime"]) * 1000.0 + i_rate = i_count / float(self._params["N_I"]) / float(self._params["simtime"]) * 1000.0 + return e_rate, i_rate + + def set_I_x(self, I_XE, I_XI): + """Set DC currents for excitatory and inhibitory neurons + Adds DC currents for the excitatory and inhibitory neurons. + The DC currents are added to the currents already + present in the populations. + + Parameters + ---------- + I_XE: float + extra DC current for excitatory neurons [pA] + I_XI: float + extra DC current for inhibitory neurons [pA] + """ + for E_pop in self._populations[0]: + I_e_loc = E_pop.get("I_e") + E_pop.set({"I_e": I_e_loc + I_XE}) + for I_pop in self._populations[1]: + I_e_loc = I_pop.get("I_e") + I_pop.set({"I_e": I_e_loc + I_XI}) + + def get_simulation(self, PathSpikes=None): + """Create network, simulate and return results + + Creates the EI-clustered network and simulates it with + the parameters supplied in the object creation. + Returns a dictionary with firing rates, + timing information (dict) and parameters (dict). + If PathSpikes is supplied the spikes get saved to a pickle file. + + Parameters + ---------- + PathSpikes: str (optional) + Path of file for spiketimes, if None, no file is saved + + Returns + ------- + dict + Dictionary with firing rates, + spiketimes (ndarray) and parameters (dict) + """ + + self.setup_network() + self.simulate() + spiketimes = self.get_recordings() + e_rate, i_rate = self.get_firing_rates(spiketimes) + + if PathSpikes is not None: + with open(PathSpikes, "wb") as outfile: + pickle.dump(spiketimes, outfile) + return { + "e_rate": e_rate, + "i_rate": i_rate, + "_params": self.get_parameter(), + "spiketimes": spiketimes, + } diff --git a/pynest/examples/EI_clustered_network/network_params.py b/pynest/examples/EI_clustered_network/network_params.py new file mode 100644 index 0000000000..c40eb9cb48 --- /dev/null +++ b/pynest/examples/EI_clustered_network/network_params.py @@ -0,0 +1,102 @@ +# -*- coding: utf-8 -*- +# +# network_params.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +"""PyNEST EI-clustered network: Network Parameters +------------------------------------------------ + +A dictionary with parameters defining the network and neuron parameters. + +""" + +import numpy as np + +net_dict = { + ############################################ + # neuron parameters + ############################################ + # neuron model + "neuron_type": "iaf_psc_exp", + # Resting potential [mV] + "E_L": 0.0, + # Membrane capacitance [pF] + "C_m": 1.0, + # Membrane time constant for excitatory neurons [ms] + "tau_E": 20.0, + # Membrane time constant for inhibitory neurons [ms] + "tau_I": 10.0, + # Refractory period [ms] + "t_ref": 5.0, + # Threshold for excitatory neurons [mV] + "V_th_E": 20.0, + # Threshold for inhibitory neurons [mV] + "V_th_I": 20.0, + # Reset potential [mV] + "V_r": 0.0, + # synaptic time constant for excitatory synapses [ms] + "tau_syn_ex": 5.0, + # synaptic time constant for inhibitory synapses [ms] + "tau_syn_in": 5.0, + # synaptic delay [ms] + "delay": 0.1, + # Feed forward excitatory input [rheobase current] + "I_th_E": 1.25, + # Feed forward inhibitory input [rheobase current] + "I_th_I": 0.78, + # distribution of feed forward input, + # I_th*[1-delta_I_../2, 1+delta_I_../2] + "delta_I_xE": 0.0, # excitatory + "delta_I_xI": 0.0, # inhibitory + # initial membrane potential: either a float (in mV) to initialize all neurons to a fixed value + # or "rand" for randomized values: "V_th_{E,I}" - 20 * nest.random.lognormal(0, 1) + "V_m": "rand", + ############################################ + # network parameters + ############################################ + # number of excitatory neurons in the network + # Neurons per cluster N_E/n_clusters + "N_E": 4000, + # number of inhibitory neurons in the network + "N_I": 1000, + # Number of clusters + "n_clusters": 20, + # connection probabilities + # baseline_conn_prob[0, 0] E to E, baseline_conn_prob[0, 1] I to E, + # baseline_conn_prob[1, 0] E to I, baseline_conn_prob[1, 1] I to I + "baseline_conn_prob": np.array([[0.2, 0.5], [0.5, 0.5]]), + # inhibitory weight ratios - scaling like random balanced network + "gei": 1.2, # I to E + "gie": 1.0, # E to I + "gii": 1.0, # I to I + # additional scaling factor for all weights + # - can be used to scale weights with network size + "s": 1.0, + # fixed indegree - otherwise established with probability ps + "fixed_indegree": False, + # cluster network by "weight" or "probabilities" + "clustering": "weight", + # ratio excitatory to inhibitory clustering, + # rj = 0 means no clustering, which resembles a clustered network + # with a blanket of inhibition + "rj": 0.82, + # excitatory clustering factor, + # rep = 1 means no clustering, reselmbles a balanced random network + "rep": 6.0, +} diff --git a/pynest/examples/EI_clustered_network/run_simulation.py b/pynest/examples/EI_clustered_network/run_simulation.py new file mode 100644 index 0000000000..ecafa112df --- /dev/null +++ b/pynest/examples/EI_clustered_network/run_simulation.py @@ -0,0 +1,56 @@ +# -*- coding: utf-8 -*- +# +# run_simulation.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +"""PyNEST EI-clustered network: Run Simulation +----------------------------------------------- + +This is an example script for running the EI-clustered model +with two stimulations and generating a raster plot. +""" + +import matplotlib.pyplot as plt +import network +from helper import raster_plot +from network_params import net_dict +from sim_params import sim_dict +from stimulus_params import stim_dict + +if __name__ == "__main__": + # Creates object which creates the EI clustered network in NEST + ei_network = network.ClusteredNetwork(sim_dict, net_dict, stim_dict) + + # Runs the simulation and returns the spiketimes + # get simulation initializes the network in NEST + # and runs the simulation + # it returns a dict with the average rates, + # the spiketimes and the used parameters + result = ei_network.get_simulation() + ax = raster_plot( + result["spiketimes"], + tlim=(0, sim_dict["simtime"]), + colorgroups=[ + ("k", 0, net_dict["N_E"]), + ("darkred", net_dict["N_E"], net_dict["N_E"] + net_dict["N_I"]), + ], + ) + plt.savefig("clustered_ei_raster.png") + print(f"Firing rate of excitatory neurons: {result['e_rate']:6.2f} spikes/s") + print(f"Firing rate of inhibitory neurons: {result['i_rate']:6.2f} spikes/s") diff --git a/pynest/examples/EI_clustered_network/sim_params.py b/pynest/examples/EI_clustered_network/sim_params.py new file mode 100644 index 0000000000..30ac2cd330 --- /dev/null +++ b/pynest/examples/EI_clustered_network/sim_params.py @@ -0,0 +1,41 @@ +# -*- coding: utf-8 -*- +# +# sim_params.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +"""PyNEST EI-clustered network: Simulation Parameters +------------------------------------------------ + +A dictionary with parameters defining the simulation. + +""" + +sim_dict = { + # The full simulation time is the sum of a presimulation time and the main + # simulation time. + # presimulation time (in ms) + "warmup": 1000.0, + # simulation time (in ms) + "simtime": 10000.0, + # resolution of the simulation (in ms) + "dt": 0.1, + "randseed": 55, + # Number of virtual processes + "n_vp": 4, +} diff --git a/pynest/examples/EI_clustered_network/stimulus_params.py b/pynest/examples/EI_clustered_network/stimulus_params.py new file mode 100644 index 0000000000..f5110a11e5 --- /dev/null +++ b/pynest/examples/EI_clustered_network/stimulus_params.py @@ -0,0 +1,38 @@ +# -*- coding: utf-8 -*- +# +# stimulus_params.py +# +# This file is part of NEST. +# +# Copyright (C) 2004 The NEST Initiative +# +# NEST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# NEST is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with NEST. If not, see . + +""" PyNEST EI-clustered network: Stimulus Parameters +----------------------------------------------- + +A dictionary with parameters for an optinal stimulation of clusters. + +""" + +stim_dict = { + # list of clusters to be stimulated (None: no stimulation, 0-n_clusters-1) + "stim_clusters": [2, 3, 4], + # stimulus amplitude (in pA) + "stim_amp": 0.15, + # stimulus start times in ms: list (warmup time is added automatically) + "stim_starts": [2000, 6000], + # list of stimulus end times in ms (warmup time is added automatically) + "stim_ends": [3500, 7500], +}