Chapter 3: Multi-cell, single population network (with BioNet)

In this tutorial, we will create a more complex network that contains multiple biophysical cells, but all of them having the same cell-type (we will cover hetergenous networks in the next tutorial). The network will contain recurrent connections, as well as external input that provide the next with stimululation.

Note - scripts and files for running this tutorial can be found in the directory sources/chapter03/

requirements: * bmtk * NEURON 7.4+

1. Building the Network

First we will build our internal network, which consists of 100 different cells. All the cells are of the same type (we’ll show how to build a heterogeneous network in the next tutorial), however they all have a different location and y-axis rotation.

nodes

import numpy as np
from bmtk.builder.networks import NetworkBuilder
from bmtk.builder.auxi.node_params import positions_columinar, xiter_random

cortex = NetworkBuilder('mcortex')
cortex.add_nodes(
    N=100,
    pop_name='Scnn1a',
    positions=positions_columinar(N=100, center=[0, 50.0, 0], max_radius=30.0, height=100.0),
    rotation_angle_yaxis=xiter_random(N=100, min_x=0.0, max_x=2*np.pi),
    rotation_angle_zaxis=3.646878266,
    potental='exc',
    model_type='biophysical',
    model_template='ctdb:Biophys1.hoc',
    model_processing='aibs_perisomatic',
    dynamics_params='472363762_fit.json',
    morphology='Scnn1a_473845048_m.swc'
)

The parameter N is used to indicate the number of cells in our population. The positions of each cell is defined by the columinar built-in method, which will random place our cells in a column (users can define their own positions as shown here). The rotation_angel_yaxis is similarl defined by a built-in function that will randomly assign each cell a given y angle.

One thing to note is that while yaxis is defined by a function which returns a lists of values, the zaxis is defined by a single value. This means that all cells will share the zaxis. we could alteratively give all cells the same y-axis rotation:

rotation_angle_yaxis=rotation_value

or give all cells a unique z-rotation angle

rotation_angle_zaxis=xiter_random(N=100, min_x=0.0, max_x=2*np.pi)

and in general, it is at the discretion of the modeler to choose what parameters are unqiue to each cell, and what parameters are global to the cell-type.

edges

Next we want to add recurrent edges. To create the connections we will use the built-in distance_connector function, which will assign the number of connections between two cells randomly (between range nsyn_min and nsysn_max) but weighted by distance. The other parameters, including the synaptic model (AMPA_ExcToExc) will be shared by all connections.

To use this, or even customized, connection functions, we must pass in the name of our connection function using the “connection_rule” parameter, and the function parameters through “connection_params” as a dictionary, which will looks something like:

connection_rule=<name_of_function>
connection_params={'param_arg1': val1, 'param_arg2': val2, ...}

The connection_rule method isn’t explicitly called by the script. Rather when the build() method is called, the connection_rule will iterate through every source/target node pair, and use the rule and build a connection matrix.

After building the connections based on our connection function, we will save the nodes and edges files into the network/ directory.

from bmtk.builder.auxi.edge_connectors import distance_connector

cortex.add_edges(
    source={'pop_name': 'Scnn1a'}, target={'pop_name': 'Scnn1a'},
    connection_rule=distance_connector,
    connection_params={'d_weight_min': 0.0, 'd_weight_max': 0.34, 'd_max': 50.0, 'nsyn_min': 0, 'nsyn_max': 10},
    syn_weight=2.0e-04,
    distance_range=[30.0, 150.0],
    target_sections=['basal', 'apical', 'soma'],
    delay=2.0,
    dynamics_params='AMPA_ExcToExc.json',
    model_template='exp2syn'
)

<bmtk.builder.connection_map.ConnectionMap at 0x7f89003e26d0>

cortex.build()
cortex.save_nodes(output_dir='sim_ch03/network')
cortex.save_edges(output_dir='sim_ch03/network')

External network

After building our internal network, we will build the external thalamic network which will provide input (see previous tutorial for more detail). Our thalamic network will consist of 100 “filter” cells, which aren’t actual cells by just place holders for spike-trains.

thalamus = NetworkBuilder('mthalamus')
thalamus.add_nodes(
    N=100,
    pop_name='tON',
    potential='exc',
    model_type='virtual'
)

The external network doesn’t have recurrent connections. Rather all the cells are feedforward onto the internal network. To do this is in a separate script which must reload the saved mcortex cell files using the import function. Then we create an edge with the thalamus nodes as the sources and the cortext nodes as the targets. This time we use the built-in connect_random connection rule, which will randomly assign each thalamus –> cortex connection between 0 and 12 synaptic connections.

from bmtk.builder.auxi.edge_connectors import connect_random

thalamus.add_edges(
    source=thalamus.nodes(), target=cortex.nodes(),
    connection_rule=connect_random,
    connection_params={'nsyn_min': 0, 'nsyn_max': 12},
    syn_weight=5.0e-05,
    distance_range=[0.0, 150.0],
    target_sections=['basal', 'apical'],
    delay=2.0,
    dynamics_params='AMPA_ExcToExc.json',
    model_template='exp2syn'
)

thalamus.build()
thalamus.save_nodes(output_dir='sim_ch03/network')
thalamus.save_edges(output_dir='sim_ch03/network')

Spike Trains

We next need to create the individual spike trains for our thalamic filter cells. We will use a Poission distrubition to create a random distribution of spikes for our 300 hundred cells each firing at ~ 15 Hz over a 3 second window. Then we can save our spike trains as a SONATA file under sim_ch03/inputs directory.

from bmtk.utils.reports.spike_trains import PoissonSpikeGenerator

psg = PoissonSpikeGenerator(population='mthalamus')
psg.add(node_ids=range(100),  # Have 10 nodes to match mthalamus
        firing_rate=15.0,    # 15 Hz, we can also pass in a nonhomogeneous function/array
        times=(0.0, 3.0))    # Firing starts at 0 s up to 3 s
psg.to_sonata('sim_ch03/inputs/mthalamus_spikes.h5')

# Let's do a quick check that we have reasonable results. Should see somewhere on the order of 15*3*100 = 4500
# spikes
psg.to_dataframe().head()

	node_ids	timestamps	population
0	0	87.943425	mthalamus
1	0	110.660282	mthalamus
2	0	140.897319	mthalamus
3	0	170.834006	mthalamus
4	0	174.629060	mthalamus

2. Setting up BioNet

file structure.

Before running a simulation, we will need to create the runtime environment, including parameter files, run-script and configuration files. You’ve already completed Chapter 02 tutorial you can just copy the files to sim_ch03 (just make sure not to overwrite the network and inputs directory).

Or create them from scracth by either running the command:

$ python -m bmtk.utils.sim_setup  \
   --report-vars v,cai            \
   --network sim_ch03/network     \
   --spikes-inputs mthalamus:sim_ch03/inputs/mthalamus_spikes.h5 \
   --dt 0.1             \
   --tstop 3000.0       \
   --include-examples   \
   --compile-mechanisms \
   bionet sim_ch03

from bmtk.utils.sim_setup import build_env_bionet

build_env_bionet(
    base_dir='sim_ch03',
    config_file='config.json',
    network_dir='sim_ch03/network',
    tstop=3000.0, dt=0.1,
    report_vars=['v', 'cai'],     # Record membrane potential and calcium (default soma)
    spikes_inputs=[('mthalamus',   # Name of population which spikes will be generated for
                    'sim_ch03/inputs/mthalamus_spikes.h5')],
    include_examples=True,    # Copies components files
    compile_mechanisms=True   # Will try to compile NEURON mechanisms
)

WARNING:root:Configuration file /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/config.json already exists, skipping. Please delete existing file, use a different name, or use overwrite_config=True.

/home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms
Mod files: "modfiles/CaDynamics.mod" "modfiles/Ca_HVA.mod" "modfiles/Ca_LVA.mod" "modfiles/Ih.mod" "modfiles/Im.mod" "modfiles/Im_v2.mod" "modfiles/Kd.mod" "modfiles/K_P.mod" "modfiles/K_T.mod" "modfiles/Kv2like.mod" "modfiles/Kv3_1.mod" "modfiles/Nap.mod" "modfiles/NaTa.mod" "modfiles/NaTs.mod" "modfiles/NaV.mod" "modfiles/SK.mod" "modfiles/vecevent.mod"

COBJS=''
 -> Compiling mod_func.cpp
 -> NMODL ../modfiles/CaDynamics.mod
 -> NMODL ../modfiles/Ca_HVA.mod
 -> NMODL ../modfiles/Ca_LVA.mod
 -> NMODL ../modfiles/Ih.mod
 -> NMODL ../modfiles/Im.mod
 -> NMODL ../modfiles/Im_v2.mod
 -> NMODL ../modfiles/Kd.mod
 -> NMODL ../modfiles/K_P.mod
 -> NMODL ../modfiles/K_T.mod
 -> NMODL ../modfiles/Kv2like.mod
 -> NMODL ../modfiles/Kv3_1.mod
 -> NMODL ../modfiles/Nap.mod
 -> NMODL ../modfiles/NaTs.mod
 -> NMODL ../modfiles/NaTa.mod
 -> NMODL ../modfiles/NaV.mod
 -> NMODL ../modfiles/SK.mod
 -> NMODL ../modfiles/vecevent.mod
 -> Compiling CaDynamics.c
 -> Compiling Ca_HVA.c
 -> Compiling Ca_LVA.c
 -> Compiling Ih.c
 -> Compiling Im.c
 -> Compiling Im_v2.c
 -> Compiling Kd.c
 -> Compiling K_P.c

Translating Ca_HVA.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/Ca_HVA.c
Translating CaDynamics.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/CaDynamics.c
Translating Ca_LVA.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/Ca_LVA.c
Thread Safe
Thread Safe
Thread Safe
Translating Ih.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/Ih.c
Translating Im_v2.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/Im_v2.c
Translating Im.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/Im.c
Thread Safe
Thread Safe
Thread Safe
Translating Kd.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/Kd.c
Translating K_P.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/K_P.c
Translating K_T.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/K_T.c
Thread Safe
Thread Safe
Thread Safe
Translating Kv2like.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/Kv2like.c
Translating Kv3_1.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/Kv3_1.c
Translating Nap.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/Nap.c
Thread Safe
Thread Safe
Thread Safe
Translating NaTs.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/NaTs.c
Translating NaTa.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/NaTa.c
Thread Safe
Thread Safe
Translating NaV.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/NaV.c
NEURON's CVode method ignores conservation
Notice: LINEAR is not thread safe.
Translating SK.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/SK.c
Translating vecevent.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch03/components/mechanisms/x86_64/vecevent.c
Notice: VERBATIM blocks are not thread safe
Thread Safe

 -> Compiling K_T.c
 -> Compiling Kv2like.c
 -> Compiling Kv3_1.c
 -> Compiling Nap.c
 -> Compiling NaTa.c
 -> Compiling NaTs.c
 -> Compiling NaV.c
 -> Compiling SK.c
 -> Compiling vecevent.c
 => LINKING shared library ./libnrnmech.so
Successfully created x86_64/special

It’s a good idea to check the configuration files sim_ch03/circuit_config.json and sim_ch03/simulation_config.json, especially to make sure that bmtk will know to use our generated spikes file (if you don’t see the below section in the simulation_config.json file go ahead and add it).

{
  "inputs": {
    "tc_spikes": {
      "input_type": "spikes",
      "module": "csv",
      "input_file": "${BASE_DIR}/mthalamus_spikes.csv",
      "node_set": "mthalamus"
    }
  }
}

3. Running the simulation

Once our config file is setup we can run a simulation either through the command line:

$ cd sim_ch03
$ python run_bionet.py config.json

or through the script

from bmtk.simulator import bionet


conf = bionet.Config.from_json('sim_ch03/config.json')
conf.build_env()
net = bionet.BioNetwork.from_config(conf)
sim = bionet.BioSimulator.from_config(conf, network=net)
sim.run()

2022-08-09 21:44:46,934 [INFO] Created log file

INFO:NEURONIOUtils:Created log file

numprocs=1
2022-08-09 21:44:47,012 [INFO] Building cells.

INFO:NEURONIOUtils:Building cells.

2022-08-09 21:44:54,838 [INFO] Building recurrent connections

INFO:NEURONIOUtils:Building recurrent connections

2022-08-09 21:44:55,189 [INFO] Building virtual cell stimulations for mthalamus_spikes

INFO:NEURONIOUtils:Building virtual cell stimulations for mthalamus_spikes

2022-08-09 21:44:58,773 [INFO] Running simulation for 3000.000 ms with the time step 0.100 ms

INFO:NEURONIOUtils:Running simulation for 3000.000 ms with the time step 0.100 ms

2022-08-09 21:44:58,774 [INFO] Starting timestep: 0 at t_sim: 0.000 ms

INFO:NEURONIOUtils:Starting timestep: 0 at t_sim: 0.000 ms

2022-08-09 21:44:58,775 [INFO] Block save every 5000 steps

INFO:NEURONIOUtils:Block save every 5000 steps

2022-08-09 21:45:26,516 [INFO]     step:5000 t_sim:500.00 ms

INFO:NEURONIOUtils:    step:5000 t_sim:500.00 ms

2022-08-09 21:45:54,497 [INFO]     step:10000 t_sim:1000.00 ms

INFO:NEURONIOUtils:    step:10000 t_sim:1000.00 ms

2022-08-09 21:46:26,947 [INFO]     step:15000 t_sim:1500.00 ms

INFO:NEURONIOUtils:    step:15000 t_sim:1500.00 ms

2022-08-09 21:46:56,036 [INFO]     step:20000 t_sim:2000.00 ms

INFO:NEURONIOUtils:    step:20000 t_sim:2000.00 ms

2022-08-09 21:47:24,824 [INFO]     step:25000 t_sim:2500.00 ms

INFO:NEURONIOUtils:    step:25000 t_sim:2500.00 ms

2022-08-09 21:47:55,361 [INFO]     step:30000 t_sim:3000.00 ms

INFO:NEURONIOUtils:    step:30000 t_sim:3000.00 ms

2022-08-09 21:47:55,393 [INFO] Simulation completed in 2.0 minutes, 56.62 seconds

INFO:NEURONIOUtils:Simulation completed in 2.0 minutes, 56.62 seconds

4. Analyzing the run.

If successful, we should have our results in the output directory. We can use the analyzer to plot a raster of the spikes over time:

from bmtk.analyzer.spike_trains import plot_raster


_ = plot_raster(config_file='sim_ch03/config.json')

In our config file we used the cell_vars and node_id_selections parameters to save the calcium influx and membrane potential of selected cells. We can also use the analyzer to display these traces:

from bmtk.analyzer.compartment import plot_traces

_ = plot_traces(config_file='sim_ch03/config.json', report_name='v_report')
_ = plot_traces(config_file='sim_ch03/config.json', report_name='v_report', node_ids=[50])
_ = plot_traces(config_file='sim_ch03/config.json', report_name='cai_report')

5. Modifying the network

Customized node params

When building our cortex nodes, we used some built-in functions to set certain parameters like positions and y-axis rotations:

cortex.add_nodes(N=100,
                 pop_name='Scnn1a',
                 positions=positions_columinar(N=100, center=[0, 50.0, 0], max_radius=30.0, height=100.0),
                 rotation_angle_yaxis=xiter_random(N=100, min_x=0.0, max_x=2*np.pi),
                 ...

These functions will assign every cell a unique value in the positions and rotation_angle_yaxis parameters, unlike the pop_name parameter which will be the same for all 100 cells. We can verify by the following code:

cortex_nodes = list(cortex.nodes())
n0 = cortex_nodes[0]
n1 = cortex_nodes[1]
print('cell 0: pop_name: {}, positions: {}, angle_yaxis: {}'.format(n0['pop_name'], n0['positions'], n0['rotation_angle_yaxis']))
print('cell 1: pop_name: {}, positions: {}, angle_yaxis: {}'.format(n1['pop_name'], n1['positions'], n1['rotation_angle_yaxis']))

cell 0: pop_name: Scnn1a, positions: [-12.67936906  16.99325185  -8.20277884], angle_yaxis: 5.127788634883025
cell 1: pop_name: Scnn1a, positions: [ -1.22804411  91.62129377 -11.58852979], angle_yaxis: 1.0965697825257734

The Network Builder contains a growing number of built-in functions. However for advanced networks a modeler will probably want to assign parameters using their own functions. To do so, a modeler only needs to passes in, or alternatively create a function that returns, a list of size N. When saving the network, each individual position will be saved in the nodes.h5 file assigned to each cell by gid.

def cortex_positions(N):
    # codex to create a list/numpy array of N (x, y, z) positions.
    return [...]

cortex.add_nodes(N=100,
                 positions=cortex_positions(100),
                 ...

or if we wanted we could give all cells the same position (The builder has no restrictions on this, however this may cause issues if you’re trying to create connections based on distance). When saving the network, the same position is assigned as a global cell-type property, and thus saved in the node_types.csv file.

cortex.add_nodes(N=100,
                 positions=np.ndarray([100.23, -50.67, 89.01]),
                 ...

We can use the same logic not just for positions and rotation_angle, but for any parameter we choose.

Customized connector functions

When creating edges, we used the built-in distance_connector function to help create the connection matrix. There are a number of built-in connection functions, but we also allow modelers to create their own. To do so, the modeler must create a function that takes in a source, target, and a variable number of parameters, and pass back a natural number representing the number of connections.

The Builder will iterate over that function passing in every source/target node pair (filtered by the source and target parameters in add_edges()). The source and target parameters are essentially dictionaries that can be used to fetch properties of the nodes. A typical example would look like:

def customized_connector(source, target, param1, param2, param3):
    if source.node_id == target.node_id:
        # necessary if we don't want autapses
        return 0
    source_pot = source['potential']
    target_pot = target['potential']
    # some code to determine number of connections
    return n_synapses

...
cortex.add_edges(source=<source_nodes>, target=<target_nodes>,
                 connection_rule=customized_connector,
                 connection_params={'param1': <p1>, 'param2': <p2>, 'param3': <p3>},
                 ...