Tutorial 5: Point-Neuron Network Models (with PointNet)

In this chapter we will create a heterogeneous network of point-model neurons and use the PointNet simulator which will run the network using the NEST simulator. As with the previous BioNet examples will create both a internal recurrently-connected network of different node types, and an external network of “virtual” neurons that will drive the firing of the internal neurons. And we’ll show how to drive network activity by using a current clamp.

PointNet, like BioNet and the other simulators, use the SONATA data format for representing networks, setting up simulations and saving results. Thus the tools used to build and display biophysically detailed networks in the previous chapters will work just the same.

Requirements: * bmtk * NEST 2.11+

1. Building the network

There are two ways of generating a network of point-neurons. Either we can take the existing biophysical network created in the previous chapters and make some minor adjustments to the neuron models being used. Or we can build a new network from scratch using the BMTK Builder.

Converting networks

We want to take the BioNet V1 network and change parameters so that the individual neurons are using point models. Luckily there parameters are stored in the node and edge “types” csv files, thus we can easily change them with a simple text editor (emacs, vi, sublime-text, etc). Here is an example of the old V1_node_types.csv:

[1]:
import pandas as pd

pd.read_csv('sources/chapter05/converted_network/V1_node_types_bionet.csv', sep=' ')
[1]:
node_type_id ei morphology_file model_processing pop_name location model_template model_type dynamics_params
0 100 e Scnn1a.swc aibs_perisomatic Scnn1a L4 ctdb:Biophys1.hoc biophysical 472363762_fit.json
1 101 i Pvalb.swc aibs_perisomatic PV L4 ctdb:Biophys1.hoc biophysical 472912177_fit.json
2 102 e NaN NaN LIF_exc VisL4 nrn:IntFire1 point_process IntFire1_exc_1.json
3 103 i NaN NaN LIF_inh VisL4 nrn:IntFire1 point_process IntFire1_inh_1.json

and here is the V1_node_types.csv used for PointNet:

[2]:
pd.read_csv('sources/chapter05/converted_network/V1_node_types.csv', sep=' ')
[2]:
node_type_id model_type model_template ei pop_name location model_processing dynamics_params
0 100 point_process nest:iaf_psc_alpha e Scnn1a VisL4 NaN 472363762_point.json
1 101 point_process nest:iaf_psc_alpha i PV VisL4 NaN 472912177_point.json
2 102 point_process nest:iaf_psc_alpha e LIF_exc VisL4 NaN IntFire1_exc_point.json
3 103 point_process nest:iaf_psc_alpha i LIF_inh VisL4 NaN IntFire1_inh_point.json

Changes: * model_type - PointNet will not support the “biophysical” model_type and only support “point_process” neuron models. * model_template - nrn:IntFire1 and ctdb:Biophys1.hoc are special directives for running NEURON based models. Instead we replaced them with the “nest:<nest-model>” directive (note we can replace iaf_psc_alpha with any valid NEST model). * dynamics_params - We have new json parameters files for the new NEST based models. * model_processing - “aibs_perisomatic” is a special command for adjusting the morphology of biophysical models, and since our NEST-based models do not have a morphology we set it to none which tells the simulator to use the models as-is (note: you can implement custom model_processing functions for PointNet that will be explained later).

We must also adjust the edges_types.csv files:

[3]:
pd.read_csv('sources/chapter05/converted_network/V1_V1_edge_types.csv', sep=' ')
[3]:
edge_type_id target_query source_query syn_weight model_template delay weight_sigma weight_function dynamics_params
0 100 pop_name=='Scnn1a' ei=='e' 2 static_synapse 2 50.0 wmax ExcToExc.json
1 101 pop_name=='LIF_exc' ei=='e' 2 static_synapse 2 50.0 wmax instantaneousExc.json
2 102 ei=='i'&model_type=='point_process' ei=='i' 2 static_synapse 2 NaN wmax InhToInh.json
3 103 ei=='i'&model_type=='point_process' ei=='i' 2 static_synapse 2 NaN wmax instantaneousInh.json
4 104 ei=='e'&model_type=='point_process' ei=='i' 2 static_synapse 2 NaN wmax InhToExc.json
5 105 ei=='e'&model_type=='point_process' ei=='i' 2 static_synapse 2 NaN wmax instantaneousInh.json
6 106 pop_name=='PV' ei=='e' 2 static_synapse 2 NaN wmax ExcToInh.json
7 107 pop_name=='LIF_inh' ei=='e' 2 static_synapse 2 NaN wmax instantaneousExc.json
  • model_template has been changed to use a NEST based model type (static_synapse)

  • Use different dynamics_parameter files

  • It’s important to readjust syn_weight as values appropiate for NEURON based models are oftern wrong for NEST based models.

Notice we don’t have to change any of the hdf5 files. The network topology remains the same making it a powerful tool for comparing networks of different levels of resolution.

Building a model from scratch.

We can use the BMTK Network Builder to create new network files just for point-based modeling

V1 Network

First lets build a “V1” network of 400 cells, split into 4 different populations

[4]:
from bmtk.builder.networks import NetworkBuilder
from bmtk.builder.auxi.node_params import positions_columinar

net = NetworkBuilder("V1")
net.add_nodes(N=80,  # Create a population of 80 neurons
              positions=positions_columinar(N=80, center=[0, 50.0, 0], max_radius=30.0, height=100.0),
              pop_name='Scnn1a', location='VisL4', ei='e',  # optional parameters
              model_type='point_process',  # Tells the simulator to use point-based neurons
              model_template='nest:iaf_psc_alpha',  # tells the simulator to use NEST iaf_psc_alpha models
              dynamics_params='472363762_point.json'  # File containing iaf_psc_alpha mdoel parameters
             )

net.add_nodes(N=20, pop_name='PV', location='VisL4', ei='i',
              positions=positions_columinar(N=20, center=[0, 50.0, 0], max_radius=30.0, height=100.0),
              model_type='point_process',
              model_template='nest:iaf_psc_alpha',
              dynamics_params='472912177_point.json')

net.add_nodes(N=200, pop_name='LIF_exc', location='L4', ei='e',
              positions=positions_columinar(N=200, center=[0, 50.0, 0], min_radius=30.0, max_radius=60.0, height=100.0),
              model_type='point_process',
              model_template='nest:iaf_psc_alpha',
              dynamics_params='IntFire1_exc_point.json')

net.add_nodes(N=100, pop_name='LIF_inh', location='L4', ei='i',
              positions=positions_columinar(N=100, center=[0, 50.0, 0], min_radius=30.0, max_radius=60.0, height=100.0),
              model_type='point_process',
              model_template='nest:iaf_psc_alpha',
              dynamics_params='IntFire1_inh_point.json')

We can now go ahead and created synaptic connections then build and save our network.

[5]:
from bmtk.builder.auxi.edge_connectors import distance_connector

## E-to-E connections
net.add_edges(source={'ei': 'e'}, target={'pop_name': 'Scnn1a'},
              connection_rule=distance_connector,
              connection_params={'d_weight_min': 0.0, 'd_weight_max': 0.34, 'd_max': 300.0, 'nsyn_min': 3, 'nsyn_max': 7},
              syn_weight=3.0,
              delay=2.0,
              dynamics_params='ExcToExc.json',
              model_template='static_synapse')

net.add_edges(source={'ei': 'e'}, target={'pop_name': 'LIF_exc'},
              connection_rule=distance_connector,
              connection_params={'d_weight_min': 0.0, 'd_weight_max': 0.34, 'd_max': 300.0, 'nsyn_min': 3, 'nsyn_max': 7},
              syn_weight=3.0,
              delay=2.0,
              dynamics_params='instantaneousExc.json',
              model_template='static_synapse')


### Generating I-to-I connections
net.add_edges(source={'ei': 'i'}, target={'pop_name': 'PV'},
              connection_rule=distance_connector,
              connection_params={'d_weight_min': 0.0, 'd_weight_max': 1.0, 'd_max': 160.0, 'nsyn_min': 3, 'nsyn_max': 7},
              syn_weight=-3.0,
              delay=2.0,
              dynamics_params='InhToInh.json',
              model_template='static_synapse')

net.add_edges(source={'ei': 'i'}, target={'ei': 'i', 'pop_name': 'LIF_inh'},
              connection_rule=distance_connector,
              connection_params={'d_weight_min': 0.0, 'd_weight_max': 1.0, 'd_max': 160.0, 'nsyn_min': 3, 'nsyn_max': 7},
              syn_weight=-3.0,
              delay=2.0,
              dynamics_params='instantaneousInh.json',
              model_template='static_synapse')

### Generating I-to-E connections
net.add_edges(source={'ei': 'i'}, target={'ei': 'e', 'pop_name': 'Scnn1a'},
              connection_rule=distance_connector,
              connection_params={'d_weight_min': 0.0, 'd_weight_max': 1.0, 'd_max': 160.0, 'nsyn_min': 3, 'nsyn_max': 7},
              syn_weight=-3.0,
              delay=2.0,
              dynamics_params='InhToExc.json',
              model_template='static_synapse')

net.add_edges(source={'ei': 'i'}, target={'ei': 'e', 'pop_name': 'LIF_exc'},
              connection_rule=distance_connector,
              connection_params={'d_weight_min': 0.0, 'd_weight_max': 1.0, 'd_max': 160.0, 'nsyn_min': 3, 'nsyn_max': 7},
              syn_weight=-3.0,
              delay=2.0,
              dynamics_params='instantaneousInh.json',
              model_template='static_synapse')

### Generating E-to-I connections
net.add_edges(source={'ei': 'e'}, target={'pop_name': 'PV'},
              connection_rule=distance_connector,
              connection_params={'d_weight_min': 0.0, 'd_weight_max': 1.0, 'd_max': 300.0, 'nsyn_min': 3, 'nsyn_max': 25},
              syn_weight=6.0,
              delay=2.0,
              dynamics_params='ExcToInh.json',
              model_template='static_synapse')


net.add_edges(source={'ei': 'e'}, target={'pop_name': 'LIF_inh'},
              connection_rule=distance_connector,
              connection_params={'d_weight_min': 0.0, 'd_weight_max': 0.26, 'd_max': 300.0, 'nsyn_min': 3, 'nsyn_max': 7},
              syn_weight=3.0,
              delay=2.0,
              dynamics_params='instantaneousExc.json',
              model_template='static_synapse')
[5]:
<bmtk.builder.connection_map.ConnectionMap at 0x7fa91747ee80>
[6]:
net.build()
net.save_nodes(output_dir='sim_ch05/network')
net.save_edges(output_dir='sim_ch05/network')

Building external network

Next we want to create an external network of “virtual cells” with spike-trains that will synapse onto our V1 cells and drive activity. We will call this external network “LGN” and contains 500 excitatory cells.

[7]:
lgn = NetworkBuilder('LGN')
lgn.add_nodes(N=500, pop_name='tON', potential='exc', model_type='virtual')

We will use a special function for setting the number of synapses between the LGN –> V1 cells. The select_source_cells function will be called during the build process.

[8]:
import numpy as np

def select_source_cells(sources, target, nsources_min=10, nsources_max=30, nsyns_min=3, nsyns_max=12):
    total_sources = len(sources)
    nsources = np.random.randint(nsources_min, nsources_max)
    selected_sources = np.random.choice(total_sources, nsources, replace=False)
    syns = np.zeros(total_sources)
    syns[selected_sources] = np.random.randint(nsyns_min, nsyns_max, size=nsources)
    return syns

lgn.add_edges(source=lgn.nodes(), target=net.nodes(pop_name='Scnn1a'),
              iterator='all_to_one',
              connection_rule=select_source_cells,
              connection_params={'nsources_min': 10, 'nsources_max': 25},
              syn_weight=15.0,
              delay=2.0,
              dynamics_params='ExcToExc.json',
              model_template='static_synapse')

lgn.add_edges(source=lgn.nodes(), target=net.nodes(pop_name='PV'),
              connection_rule=select_source_cells,
              connection_params={'nsources_min': 15, 'nsources_max': 35},
              iterator='all_to_one',
              syn_weight=6.0,
              delay=2.0,
              dynamics_params='ExcToInh.json',
              model_template='static_synapse')

lgn.add_edges(source=lgn.nodes(),  target=net.nodes(pop_name='LIF_exc'),
              connection_rule=select_source_cells,
              connection_params={'nsources_min': 10, 'nsources_max': 25},
              iterator='all_to_one',
              syn_weight= 10.0,
              delay=2.0,
              dynamics_params='instantaneousExc.json',
              model_template='static_synapse')

lgn.add_edges(source=lgn.nodes(),  target=net.nodes(pop_name='LIF_inh'),
              connection_rule=select_source_cells,
              connection_params={'nsources_min': 15, 'nsources_max': 30},
              iterator='all_to_one',
              syn_weight=10.0,
              delay=2.0,
              dynamics_params='instantaneousExc.json',
              model_template='static_synapse')
[8]:
<bmtk.builder.connection_map.ConnectionMap at 0x7fa9173cfd30>

Finally we build and save our lgn network.

[9]:
lgn.build()
lgn.save_nodes(output_dir='sim_ch05/network')
lgn.save_edges(output_dir='sim_ch05/network')

2. Setting up PointNet Environment

Directory Structure

Before running a simulation, we will need to create the runtime environment, including parameter files, run-script and configuration files. If using the tutorial these files will already be in place. Otherwise we can use a command-line:

$ python -m bmtk.utils.sim_setup        \
    --network sim_ch05/network/         \
    --include-examples                  \
    --tstop 3000.0                      \
    pointnet sim_ch05/

or

[10]:
from bmtk.utils.sim_setup import build_env_pointnet

build_env_pointnet(base_dir='sim_ch05',
                   network_dir='sim_ch05/network',
                   tstop=3000.0,
                   dt=0.01,
                   include_examples=True,         # Copies components files
                  )

The network files are written to circuit_config.json and the simulation parameters are set in simulation_config, which are both linked to the main config.json file. The simulation time is set to run for 3000.0 ms (tstop). In general, all the parameters needed to setup and start a simulation are found in the config files, and adjusting network/simulation conditions can be done by editing these json files in a text editor.

lgn input

We need to provide our LGN external network cells with spike-trains so they can activate our recurrent network. Previously we showed how to do this by generating csv files. We can also use SONATA spike files.

First download the file:

$ cd sim_ch05
$ wget https://github.com/AllenInstitute/bmtk/raw/develop/examples/spikes_inputs/lgn_spikes.h5

Then we must edit the simulation_config.json file to tell the simulator to find the file and which network to associate it with.

{
  "inputs": {
    "LGN_spikes": {
      "input_type": "spikes",
      "module": "h5",
      "input_file": "$BASE_DIR/lgn_spikes.h5",
      "node_set": "LGN"
    }
  }
}

3. Running the simulation

The call to sim_setup created a file run_pointnet.py which we can run directly in a command line:

$ python run_pointnet.py config.json

or if you have mpi setup:

$ mpirun -np $NCORES python run_pointnet.py config.json

Or we can run it directly

[16]:
from bmtk.simulator import pointnet

configure = pointnet.Config.from_json('sim_ch05/config.json')
configure.build_env()
network = pointnet.PointNetwork.from_config(configure)
sim = pointnet.PointSimulator.from_config(configure, network)
sim.run()
2022-10-25 08:57:58,631 [INFO] Created log file
INFO:NestIOUtils:Created log file
2022-10-25 08:57:58,651 [INFO] Batch processing nodes for V1/0.
INFO:NestIOUtils:Batch processing nodes for V1/0.
2022-10-25 08:57:58,660 [INFO] Batch processing nodes for LGN/0.
INFO:NestIOUtils:Batch processing nodes for LGN/0.
2022-10-25 08:57:58,693 [INFO] Setting up output directory
INFO:NestIOUtils:Setting up output directory
2022-10-25 08:57:58,694 [INFO] Building cells.
INFO:NestIOUtils:Building cells.
2022-10-25 08:57:58,710 [INFO] Building recurrent connections
INFO:NestIOUtils:Building recurrent connections
2022-10-25 08:57:58,776 [INFO] Build virtual cell stimulations for LGN_spikes
INFO:NestIOUtils:Build virtual cell stimulations for LGN_spikes
2022-10-25 08:57:59,970 [INFO] Network created.
INFO:NestIOUtils:Network created.
2022-10-25 08:58:00,001 [INFO] Starting Simulation
INFO:NestIOUtils:Starting Simulation
2022-10-25 08:58:03,008 [INFO] Simulation finished, finalizing results.
INFO:NestIOUtils:Simulation finished, finalizing results.
2022-10-25 08:58:03,284 [INFO] Done.
INFO:NestIOUtils:Done.

4. Analyzing results

Results of the simulation, as specified in the config, are saved into the output directory. Using the analyzer functions, we can do things like plot the raster plot

[17]:
from bmtk.analyzer.spike_trains import plot_raster, plot_rates

_= plot_raster(config_file='sim_ch05/config.json', group_by='pop_name', plt_style='seaborn-muted')
_images/tutorial_pointnet_modeling_27_0.png

Or we can plot the rates of the different populations

[18]:
_= plot_rates(config_file='sim_ch05/config.json', group_by='pop_name', plt_style='seaborn-muted')
_images/tutorial_pointnet_modeling_29_0.png

In our simulation_config.json in the reports section, we can see we also record the V_m (i.e membrane potential) of a select sample of cells. By default these files are written to an hdf5 file with the same name as the report (membrane_potential.h5), and we can use the analyzer to show the time course of some of these cells.

5. Additional Information

Recording membrane potential (and other cell properties).

By default bmtk will automatically record and save the spike-trains of our (non-virtual) cells. Sometimes it can be useful to record other properties besides spikes, like the membrane voltage traces, of all or a subsection of the cells. To do so we only need to open the simulation_config.json with our prefered text editor and insert the following section:

{
  "reports": {
    "membrane_potential": {
      "cells": {
        "population": "V1",
        "node_id": [0, 80, 100]
      },
      "variable_name": "V_m",
      "module": "multimeter_report",
      "sections": "soma"
    }
  }
}

This will tell bmtk to record the “V_m” variable for selected cells with ids 0, 80, 100 - and save the results in a file output/membrane_potential.h5. After re-running the simulation the file will be generated, according to the SONATA format, which we can read using hdf5 or bmtk like below:

[19]:
from bmtk.analyzer.compartment import plot_traces

_ = plot_traces(config_file='sim_ch05/config.json', group_by='pop_name', plt_style='seaborn-muted')
_images/tutorial_pointnet_modeling_33_0.png

Perturbation simulations

A common use case is to mimick a network where the firing rate of select cells are being depressed or stimulated using current clamps or optogenetic methods. The best way to do this is by applying a current clamp in a similar manner used in bionet.

The point_120cells network in the bmtk examples consists of 120 recurrently connected excitatory and inhibitory cells with randomized inputs. Normally the results without any type of perturbations will look as follows: 30213ab18c0745f699e02999de8bb24d

To artifically excite some of the cells we can open the config.json file and add the following section to the “inputs”:

{
  "exc_perturbation": {
    "input_type": "current_clamp",
    "module": "IClamp",
    "node_set": {
      "population": "cortex",
      "node_ids": [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]
    },
    "amp": 230.0,
    "delay": 1.0,
    "duration": 3000.0
  }
}

This will apply a large step current for the duration of the simulation for cells 20 through 39. In this case we knew the exact node_ids of the cells we want to excite. However it is also possible to get a subset using a specific property, for example if we wanted to only select inhibitory L2/3 cells:

{
  "node_set": {
    "population": "cortex",
    "ei": "i",
    "location": "L2/3"
  }
}

To simulate depression of cells 40 to 49 we will add another current_clamp but with a large negative current:

{
  "inh_perturbation": {
    "input_type": "current_clamp",
    "module": "IClamp",
    "node_set": {
      "population": "cortex",
      "node_ids": [40, 41, 42, 43, 44, 45, 46, 47, 48, 49]
    },
    "amp": -230.0,
    "delay": 1.0,
    "duration": 3000.0
  }
}

After we rerun the simulation: b638f017c4a34b869b8d3f961909b702