Adding 300 NEST example

examples/300_pointneurons/build_network.py

@@ -0,0 +1,157 @@
+import os
+import numpy as np
+
+from bmtk.builder import NetworkBuilder
+from bmtk.builder.bionet import SWCReader
+from bmtk.utils.io.spike_trains import PoissonSpikesGenerator
+from bmtk.builder.aux.node_params import positions_columinar, xiter_random
+
+build_recurrent_edges = True
+
+print('Building internal network')
+# List of non-virtual cell models
+cell_models = [
+    {
+        'model_name': 'Scnn1a', 'ei': 'e',
+        'model_template': 'nest:iaf_psc_alpha',
+        'dynamics_params': '472363762_point.json'
+    },
+    {
+        'model_name': 'Rorb', 'ei': 'e',
+        'model_template': 'nest:iaf_psc_alpha',
+        'dynamics_params': '473863510_point.json'
+    },
+    {
+        'model_name': 'Nr5a1', 'ei': 'e',
+        'model_template': 'nest:iaf_psc_alpha',
+        'dynamics_params': '473863035_point.json'
+    },
+    {
+        'model_name': 'PV1', 'ei': 'i',
+        'model_template': 'nest:iaf_psc_alpha',
+        'dynamics_params': '472912177_point.json'
+    },
+    {
+        'model_name': 'PV2', 'ei': 'i',
+        'model_template': 'nest:iaf_psc_alpha',
+        'dynamics_params': '473862421_point.json'
+    }
+]
+
+'''
+morphologies = {p['model_name']: SWCReader(os.path.join('../shared_components/morphologies',
+                                                        '{}.swc'.format(p['morphology'])))
+                for p in cell_models}
+def build_edges(src, trg, sections=['basal', 'apical'], dist_range=[50.0, 150.0]):
+    """Function used to randomly assign a synaptic location based on the section (soma, basal, apical) and an
+    arc-length dist_range from the soma. This function should be passed into the network and called during the build
+    process.
+
+    :param src: source cell (dict)
+    :param trg: target cell (dict)
+    :param sections: list of target cell sections to synapse onto
+    :param dist_range: range (distance from soma center) to place
+    :return:
+    """
+    # Get morphology and soma center for the target cell
+    swc_reader = morphologies[trg['model_name']]
+    target_coords = [trg['x'], trg['y'], trg['z']]
+
+    sec_ids, sec_xs = swc_reader.choose_sections(sections, dist_range)  # randomly choose sec_ids
+    coords = swc_reader.get_coord(sec_ids, sec_xs, soma_center=target_coords)  # get coords of sec_ids
+    dist = swc_reader.get_dist(sec_ids)
+    swctype = swc_reader.get_type(sec_ids)
+    return sec_ids, sec_xs, coords[0][0], coords[0][1], coords[0][2], dist[0], swctype[0]
+'''
+
+# Build a network of 300 biophysical cells to simulate
+internal = NetworkBuilder("internal")
+for i, model_props in enumerate(cell_models):
+    n_cells = 80 if model_props['ei'] == 'e' else 30  # 80% excitatory, 20% inhib
+
+    # Randomly get positions uniformly distributed in a column
+    positions = positions_columinar(N=n_cells, center=[0, 10.0, 0], max_radius=50.0, height=200.0)
+
+    internal.add_nodes(N=n_cells,
+                       x=positions[:, 0], y=positions[:, 1], z=positions[:, 2],
+                       rotation_angle_yaxis=xiter_random(N=n_cells, min_x=0.0, max_x=2 * np.pi),  # randomly rotate y axis
+                       model_type='point_process',
+                       **model_props)
+
+
+
+if build_recurrent_edges:
+    def n_connections(src, trg, prob=0.5, min_syns=2, max_syns=7):
+        return 0 if np.random.uniform() > prob else np.random.randint(min_syns, max_syns)
+
+    # exc --> exc connections
+    cm = internal.add_edges(source={'ei': 'e'}, target={'ei': 'e'},
+                            connection_rule=lambda s, t: np.random.binomial(1, 0.2),
+                            dynamics_params='ExcToExc.json',
+                            model_template='static_synapse',
+                            # syn_weight=2.5,
+                            delay=2.0)
+    cm.add_properties('syn_weight', rule=2.5, dtypes=np.float)
+
+    # exc --> inh connections
+    cm = internal.add_edges(source={'ei': 'e'}, target={'ei': 'i'},
+                            connection_rule=lambda s, t: np.random.binomial(1, 0.5),
+                            dynamics_params='ExcToInh.json',
+                            model_template='static_synapse',
+                            # syn_weight=7.0,
+                            delay=2.0)
+    cm.add_properties('syn_weight', rule=7.0, dtypes=np.float)
+
+    # inh --> exc connections
+    cm = internal.add_edges(source={'ei': 'i'}, target={'ei': 'e'},
+                            connection_rule=lambda s, t: np.random.binomial(1, 0.5),
+                            #connection_rule=lambda *_: np.random.randint(0, 4),
+                            dynamics_params='InhToExc.json',
+                            model_template='static_synapse',
+                            # syn_weight=-7.5,
+                            delay=2.0)
+    cm.add_properties('syn_weight', rule=-7.5, dtypes=np.float)
+
+    # inh --> inh connections
+    cm = internal.add_edges(source={'ei': 'i'}, target={'ei': 'i'},
+                            connection_rule=lambda s, t: np.random.binomial(1, 0.5),
+                            dynamics_params='InhToInh.json',
+                            model_template='static_synapse',
+                            # syn_weight=-3.0,
+                            delay=2.0)
+    cm.add_properties('syn_weight', rule=-3.0, dtypes=np.float)
+
+
+internal.build()
+
+print('Saving internal')
+internal.save(output_dir='network')
+
+
+print('Building external connections')
+external = NetworkBuilder("external")
+external.add_nodes(N=100, model_type='virtual', ei='e')
+cm = external.add_edges(target=internal.nodes(ei='e'), source=external.nodes(),
+                        connection_rule=lambda *_: np.random.binomial(1, .7),
+                        dynamics_params='ExcToExc.json',
+                        model_template='static_synapse',
+                        # syn_weight=50.0,
+                        delay=2.0)
+cm.add_properties('syn_weight', rule=50.0, dtypes=np.float)
+
+
+cm = external.add_edges(target=internal.nodes(ei='i'), source=external.nodes(),
+                        connection_rule=lambda *_: np.random.binomial(1, .7),
+                        dynamics_params='ExcToInh.json',
+                        model_template='static_synapse',
+                        # syn_weight=75.0,
+                        delay=2.0)
+cm.add_properties('syn_weight', rule=65.0, dtypes=np.float)
+
+
+external.build()
+
+print('Saving external')
+external.save(output_dir='network')
+
+

examples/300_pointneurons/circuit_config.json

@@ -0,0 +1,37 @@
+{
+  "manifest": {
+    "$NETWORK_DIR": "./network",
+    "$COMPONENT_DIR": "../shared_components/nest_models"
+  },
+
+  "components": {
+    "point_neuron_models_dir": "$COMPONENT_DIR/cell_models",
+    "synaptic_models_dir": "$COMPONENT_DIR/synaptic_models"
+  },
+
+  "networks": {
+    "nodes": [
+      {
+        "nodes_file": "$NETWORK_DIR/internal_nodes.h5",
+        "node_types_file": "$NETWORK_DIR/internal_node_types.csv"
+      },
+      {
+        "nodes_file": "$NETWORK_DIR/external_nodes.h5",
+        "node_types_file": "$NETWORK_DIR/external_node_types.csv"
+      }
+    ],
+
+    "edges": [
+      {
+        "edges_file": "$NETWORK_DIR/internal_internal_edges.h5",
+        "edge_types_file": "$NETWORK_DIR/internal_internal_edge_types.csv",
+        "enabled": true
+      },
+      {
+        "edges_file": "$NETWORK_DIR/external_internal_edges.h5",
+        "edge_types_file": "$NETWORK_DIR/external_internal_edge_types.csv"
+      }
+    ]
+
+  }
+}

examples/300_pointneurons/config.json

@@ -0,0 +1,4 @@
+{
+  "network": "./circuit_config.json",
+  "simulation": "./simulation_config.json"
+}

examples/300_pointneurons/generate_spikes.py

@@ -0,0 +1,7 @@
+import os
+from bmtk.utils.io.spike_trains import PoissonSpikesGenerator
+
+if not os.path.exists('inputs'):
+    os.mkdir('inputs')
+psg = PoissonSpikesGenerator(range(100), 14.0, tstop=3000.0, precision=3)
+psg.to_hdf5('inputs/external_spike_trains.h5')

examples/300_pointneurons/network/external_internal_edge_types.csv

@@ -0,0 +1,3 @@
+edge_type_id target_query source_query delay dynamics_params model_template
+100 ei=='e' * 2.0 ExcToExc.json static_synapse
+101 ei=='i' * 2.0 ExcToInh.json static_synapse

examples/300_pointneurons/network/external_node_types.csv

@@ -0,0 +1,2 @@
+node_type_id model_type ei
+100 virtual e

examples/300_pointneurons/network/internal_internal_edge_types.csv

@@ -0,0 +1,5 @@
+edge_type_id target_query source_query delay dynamics_params model_template
+100 ei=='e' ei=='e' 2.0 ExcToExc.json static_synapse
+101 ei=='i' ei=='e' 2.0 ExcToInh.json static_synapse
+102 ei=='e' ei=='i' 2.0 InhToExc.json static_synapse
+103 ei=='i' ei=='i' 2.0 InhToInh.json static_synapse

examples/300_pointneurons/network/internal_node_types.csv

@@ -0,0 +1,6 @@
+node_type_id ei model_template model_type dynamics_params model_name
+104 i nest:iaf_psc_alpha point_process 473862421_point.json PV2
+100 e nest:iaf_psc_alpha point_process 472363762_point.json Scnn1a
+101 e nest:iaf_psc_alpha point_process 473863510_point.json Rorb
+102 e nest:iaf_psc_alpha point_process 473863035_point.json Nr5a1
+103 i nest:iaf_psc_alpha point_process 472912177_point.json PV1

examples/300_pointneurons/node_sets.json

@@ -0,0 +1,7 @@
+{
+  "external": {
+    "population": "external"
+  },
+
+  "recorded_cells": [0, 80, 160, 240, 270]
+}

examples/300_pointneurons/plot_potential.py

@@ -0,0 +1,76 @@
+"""
+A script for plotting cell variable (eg voltage, calcium) traces from a SONATA hdf5 file.
+
+To plot a single ouput:
+ $ python plot_potential.py path/to/cell_var.h5
+
+If you want to perform a side-by-side comparison of multiple output files:
+ $ python plot_potential.py cell_var_0.h5 cell_var_1.h5 ... cell_var_n.h5
+cell_var_0.h5 is used as the base-line for plotting number of gid/subplots and time-traces.
+
+"""
+
+import os
+import matplotlib.pyplot as plt
+from optparse import OptionParser
+from bmtk.utils.cell_vars import CellVarsFile
+
+
+def plot_vars(file_names, cell_var='v', gid_list=[]):
+    """Plots variable traces for a SONATA h5 file. If multiple spike files are specified will do a side-by-side
+    comparsion for each gid.
+
+    :param file_names: list of cell_var file names
+    :param cell_var: cell variable to plot trace
+    :param gid_list: used to set what gid/subplots to show (if empty list just plot all possible gids)
+    """
+    # convert to list if single spike file passed in
+    file_names = [file_names] if not isinstance(file_names, (tuple, list)) else file_names
+    assert(len(file_names) > 0)
+
+    # Use bmtk to parse the cell-var files
+    cell_var_files = [CellVarsFile(fn) for fn in file_names]
+
+    # get first spike file and properties
+    cvf_base = cell_var_files[0]
+    xlim = [cvf_base.time_trace[0], cvf_base.time_trace[-1]]  # Use the same x-axis across all subplots
+    gid_list = cvf_base.gids if not gid_list else gid_list  # if gid_list is None just get all gids in first file
+    n_cells = len(cvf_base.gids)
+
+    fig, ax = plt.subplots(n_cells, 1, figsize=(10, 10))
+    for subplot, gid in enumerate(gid_list):
+        for i, cvf in enumerate(cell_var_files):
+            # plot all traces
+            ax[subplot].plot(cvf.time_trace, cvf.data(gid, cell_var), label=file_names[i])
+
+        ax[subplot].yaxis.set_label_position("right")
+        ax[subplot].set_ylabel('gid {}'.format(gid), fontsize='xx-small')
+        ax[subplot].set_xlim(xlim)
+        if subplot + 1 < n_cells:
+            # remove x-axis labels on all but the last plot
+            ax[subplot].set_xticklabels([])
+        else:
+            # Use the last plot to get the legend
+            handles, labels = ax[subplot].get_legend_handles_labels()
+            fig.legend(handles, labels, loc='upper right')
+
+    plt.show()
+
+
+if __name__ == '__main__':
+    parser = OptionParser(usage="Usage: python %prog [options] <cell_vars>.h5 [<cell_vars2>.h5 ...]")
+    parser.add_option('--variable', type=str, dest='cell_var', default='V_m', help='Cell variable to compare (v, cai, etc)')
+    parser.add_option('--gids', type='string', dest='gids', default=[],
+                      action='callback',
+                      callback=lambda option, opt, value, parser: setattr(parser.values, option.dest, [int(v) for v in value.split(',')]),
+                      help='comma seperated list of gids to plot')
+    options, args = parser.parse_args()
+
+    if len(args) == 0:
+        # If no file is specified see if there is a output/membrane_potential.h5 file to plot membrane voltage
+        if not os.path.exists('output/membrane_potential.h5'):
+            raise Exception('Please specifiy hdf5 file to read in arguments. Exiting!')
+        else:
+            plot_vars('output/membrane_potential.h5', cell_var=options.cell_var)
+    else:
+        plot_vars(file_names=args, cell_var=options.cell_var, gid_list=options.gids)

examples/300_pointneurons/plot_spikes.py

@@ -0,0 +1,3 @@
+from bmtk.analyzer.visualization.spikes import plot_spikes
+
+plot_spikes('network/internal_nodes.h5', 'network/internal_node_types.csv', 'output/spikes.h5', group_key='model_name')

examples/300_pointneurons/run_pointnet.py

@@ -0,0 +1,18 @@
+import os, sys
+from bmtk.simulator import pointnet
+from bmtk.analyzer.visualization.spikes import plot_spikes
+
+
+def main(config_file):
+    configure = pointnet.Config.from_json(config_file)
+    configure.build_env()
+
+    graph = pointnet.PointNetwork.from_config(configure)
+    sim = pointnet.PointSimulator.from_config(configure, graph)
+    sim.run()
+
+    # plot_spikes('network/internal_nodes.h5', 'network/internal_node_types.csv', 'output/spikes.h5', group_key='model_name')
+
+
+if __name__ == '__main__':
+    main('config.json')

examples/300_pointneurons/simulation_config.json

@@ -0,0 +1,49 @@
+{
+  "manifest": {
+    "$BASE_DIR": ".",
+    "$OUTPUT_DIR": "$BASE_DIR/output",
+    "$INPUT_DIR": "$BASE_DIR/inputs"
+  },
+
+  "run": {
+    "tstop": 1500.0,
+    "dt": 0.001,
+    "nsteps_block": 5000
+  },
+
+  "target_simulator":"NEST",
+
+  "network": "$BASE_DIR/circuit_config.json",
+
+  "conditions": {
+    "celsius": 34.0,
+    "v_init": -80
+  },
+
+  "node_sets_file": "node_sets.json",
+
+  "inputs": {
+    "external_spike_trains": {
+      "input_type": "spikes",
+      "module": "h5",
+      "input_file": "$INPUT_DIR/external_spike_trains.h5",
+      "node_set": "external"
+    }
+  },
+
+  "output": {
+    "log_file": "log.txt",
+    "output_dir": "$OUTPUT_DIR",
+    "spikes_file": "spikes.h5",
+    "spikes_sort_order": "time"
+  },
+
+  "reports": {
+    "membrane_potential": {
+      "cells": "recorded_cells",
+      "variable_name": "V_m",
+      "module": "membrane_report",
+      "sections": "multimeter_report"
+    }
+  }
+}

examples/shared_components/nest_models/cell_models/472363762_point.json

@@ -0,0 +1,9 @@
+{
+  "I_e": 0.0,
+  "tau_m": 44.9,
+  "C_m": 239.0,
+  "t_ref": 3.0,
+  "E_L": -78.0,
+  "V_th": -43.0,
+  "V_reset": -55.0
+}