Tutorial 2: Single cell simulation with external feedfoward input (with BioNet)

In the previous tutorial we built a single cell and stimulated it with a current injection. In this example we will keep our single-cell network, but instead of stimulation by a step current, we’ll set-up an external network that synapses onto our cell.

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

Requirements: * Python 2.7 or 3.6+ * bmtk * NEURON 7.4+

Step 1: Building the network.

Similar to the previous tutorial, we want to build and save a network consisting of a single biophysically detailed cell.

from bmtk.builder.networks import NetworkBuilder


cortex = NetworkBuilder('mcortex')
cortex.add_nodes(
    cell_name='Scnn1a_473845048',
    potental='exc',
    model_type='biophysical',
    model_template='ctdb:Biophys1.hoc',
    model_processing='aibs_perisomatic',
    dynamics_params='472363762_fit.json',
    morphology='Scnn1a_473845048_m.swc'
)

cortex.build()
cortex.save_nodes(output_dir='sim_ch02/network')

But we will also want a collection of external spike-generating cells that will synapse onto our cell. To do this we create a second network which can represent thalamic input. We will call our network “mthalamus”, and it will consist of 10 cells. These cells are not biophysical but instead “virtual” cells. Virtual cells don’t have a morphology or the normal properties of a neuron, but rather act as spike generators.

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

Now that we built our nodes, we want to create a connection between our 10 thalamic cells onto our cortex cell. To do so we use the add_edges function like so:

thalamus.add_edges(
    source={'pop_name': 'tON'}, target=cortex.nodes(),
    connection_rule=5,
    syn_weight=0.001,
    delay=2.0,
    weight_function=None,
    target_sections=['basal', 'apical'],
    distance_range=[0.0, 150.0],
    dynamics_params='AMPA_ExcToExc.json',
    model_template='exp2syn'
)

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

Let us break down how this method call works:

thalamus.add_edges(source={'pop_name': 'tON'}, target=cortex.nodes(),

• Here we specify which set of nodes to use as sources and targets. Our source/pre-synaptic cells are all thamalus cells with the property “pop_name=tON”, which in this case is every thalmus cell (We could also use source=thalamus.nodes(), or source={‘level_of_detail’: ‘filter’}). The target/post-synaptic is all cell(s) of the “cortex” network.

connection_rule=5,

• The connection_rule parameter determines how many synapses exists between every source/target pair. In this very trivial case we are indicating that between every thamalic –> cortical cell connection, there are 5 synapatic connections. In future tutorials we will show how we can create more complex customized rules.

syn_weight=0.001,
delay=2.0,
weight_function=None,

• Here we are specifying the connection weight. For every connection in this edge-type, there is a connection strength of 5e-05 (units) and a connection delay of 2 ms. The weight function is used to adjust the weights before runtime. Later we will show how to create customized weight functions.

target_sections=['basal', 'apical'],
distance_range=[0.0, 150.0],

• This is used by BioNet to determine where on the post-synaptic cell to place the synapse. By default placement is random within the given section and range.

  dynamics_params='AMPA_ExcToExc.json',
  model_template='exp2syn')
  

• The params_file give the parameters of the synpase, including the time constant and potential. Here we are using an AMPA type synaptic model with an Excitatory connection. The set_params_function is used by BioNet to convert the model into a valid NEURON synaptic object.

Finally we are ready to build the model and save the thalamic nodes and edges.

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

The network/ directory will contain multiple nodes and edges files. It should have nodes (and node-types) files for both the thalamus and cortex network. And edges (and edge-types) files for the thalamus –> cortex connections. Nodes and edges for different networks and their connections are spread out across different files which allows us in the future to rebuild, edit or replace part of setup in a piecemeal and efficent manner.

Spike Trains

We need to give our 10 thalamic cells spike trains. There are multiple ways to do this, but an easy way to use a sonata hdf5 file. The following function will create a file to provide the spikes for our 10 cells.

from bmtk.utils.reports.spike_trains import PoissonSpikeGenerator

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

The spikes files are save as a Sonata formated spikes train. BMTK also allow us to use csv files and NWB files to set the spike train - we just have to tell the configuration file the appropiate format to read.

print('Number of spikes: {}'.format(psg.n_spikes()))
print('Units: {}'.format(psg.units()))

psg.to_dataframe().head()

Number of spikes: 315
Units: ms

	node_ids	timestamps	population
0	0	45.327109	mthalamus
1	0	77.644564	mthalamus
2	0	163.482462	mthalamus
3	0	517.296866	mthalamus
4	0	535.154362	mthalamus

Step 2: Setting up BioNet environment.

file 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, however you should run the following anyway in a command-line:

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

Or call the function directly in python

from bmtk.utils.sim_setup import build_env_bionet

build_env_bionet(
    base_dir='sim_ch02',
    config_file='config.json',
    network_dir='sim_ch02/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_ch02/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_ch02/config.json already exists, skipping. Please delete existing file, use a different name, or use overwrite_config=True.
WARNING:root:Configuration file /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/circuit_config.json already exists, skipping. Please delete existing file, use a different name, or use overwrite_config=True.
WARNING:root:Configuration file /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/simulation_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_ch02/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/NaTa.mod
 -> NMODL ../modfiles/NaTs.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_LVA.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/Ca_LVA.c
Translating Ca_HVA.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/Ca_HVA.c
Translating CaDynamics.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/CaDynamics.c
Thread Safe
Thread Safe
Thread Safe
Translating Ih.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/Ih.c
Translating Im.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/Im.c
Translating Im_v2.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/Im_v2.c
Thread Safe
Thread Safe
Thread Safe
Translating K_P.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/K_P.c
Translating Kd.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/Kd.c
Translating K_T.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/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_ch02/components/mechanisms/x86_64/Kv2like.c
Translating Kv3_1.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/Kv3_1.c
Translating Nap.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/Nap.c
Thread Safe
Thread Safe
Thread Safe
Translating NaTs.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/NaTs.c
Translating NaTa.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/NaTa.c
Translating NaV.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/NaV.c
Thread Safe
Thread Safe
NEURON's CVode method ignores conservation
Notice: LINEAR is not thread safe.
Translating SK.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/components/mechanisms/x86_64/SK.c
Translating vecevent.mod into /home/ping/bmtk_change/bmtk/docs/tutorial/sim_ch02/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

The last thing that we need to do is to update the configuration file to read “thalamus_spikes.csv”. To do so we open simulation_config.json in a text editor and add the following to the input section.

{
  "inputs": {
    "mthalamus_spikes": {
      "input_type": "spikes",
      "module": "sonata",
      "input_file": "${BASE_DIR}/input/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_ch02/
$ python run_bionet.py config.json

or through the script

from bmtk.simulator import bionet


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

2022-08-08 16:06:48,667 [INFO] Created log file

INFO:NEURONIOUtils:Created log file

numprocs=1
2022-08-08 16:06:48,750 [INFO] Building cells.

INFO:NEURONIOUtils:Building cells.

2022-08-08 16:06:48,919 [INFO] Building recurrent connections

INFO:NEURONIOUtils:Building recurrent connections

2022-08-08 16:06:48,922 [INFO] Building virtual cell stimulations for mthalamus_spikes

INFO:NEURONIOUtils:Building virtual cell stimulations for mthalamus_spikes

2022-08-08 16:06:48,955 [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-08 16:06:48,956 [INFO] Starting timestep: 0 at t_sim: 0.000 ms

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

2022-08-08 16:06:48,957 [INFO] Block save every 5000 steps

INFO:NEURONIOUtils:Block save every 5000 steps

2022-08-08 16:06:49,289 [INFO]     step:5000 t_sim:500.00 ms

INFO:NEURONIOUtils:    step:5000 t_sim:500.00 ms

2022-08-08 16:06:49,582 [INFO]     step:10000 t_sim:1000.00 ms

INFO:NEURONIOUtils:    step:10000 t_sim:1000.00 ms

2022-08-08 16:06:49,881 [INFO]     step:15000 t_sim:1500.00 ms

INFO:NEURONIOUtils:    step:15000 t_sim:1500.00 ms

2022-08-08 16:06:50,176 [INFO]     step:20000 t_sim:2000.00 ms

INFO:NEURONIOUtils:    step:20000 t_sim:2000.00 ms

2022-08-08 16:06:50,495 [INFO]     step:25000 t_sim:2500.00 ms

INFO:NEURONIOUtils:    step:25000 t_sim:2500.00 ms

2022-08-08 16:06:50,800 [INFO]     step:30000 t_sim:3000.00 ms

INFO:NEURONIOUtils:    step:30000 t_sim:3000.00 ms

2022-08-08 16:06:50,808 [INFO] Simulation completed in 1.853 seconds

INFO:NEURONIOUtils:Simulation completed in 1.853 seconds

4. Analyzing the run

from bmtk.analyzer.spike_trains import to_dataframe

results_df = to_dataframe(config_file='sim_ch02/config.json')
print('Number of Spikes: {}'.format(len(results_df)))
results_df.head()

Number of Spikes: 57

	timestamps	node_ids	population
0	45.0	0	mcortex
1	1375.1	0	mcortex
2	1532.4	0	mcortex
3	1686.9	0	mcortex
4	1780.4	0	mcortex

from bmtk.analyzer.compartment import plot_traces

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

5. Additional things to do:

Changing edge properties

When using the Network Builder add_edges method, we gave all the edges the same parameter values (delay, weight, target_section, etc.). All connection created by this method call constitute an single edge-type that share the same parameters, and are specified in the mthalamic_edge_type.csv file

import pandas as pd
pd.read_csv('sim_ch02/network/mthalamus_mcortex_edge_types.csv', sep=' ')

	edge_type_id	target_query	source_query	model_template	weight_function	target_sections	delay	dynamics_params	syn_weight	distance_range
0	100	*	pop_name=='tON'	exp2syn	NaN	['basal', 'apical']	2.0	AMPA_ExcToExc.json	0.001	[0.0, 150.0]

(if in the build script we called add_edges multiple times, we’d have multiple edge-types).

Using a simple text-editor we can modify this file directly, change parameters before a simulation run without having to rebuild the entire network (although for a network this small it may not be beneficial).

weight_function

By default BioNet uses the value in syn_weight to set a synaptic weight, which is a constant stored in the network files. Often we will want to adjust the synaptic weight between simulations, but don’t want to have to regenerate the network. BioNet allows us to specify custom synaptic weight functions that will calculate synaptic weight before each simulation.

To do so first we must set the value of ‘weight_function’ column. Either we can open up the file mthalamus_mcortex_edge_types.csv with a text-editor and change the column.

edge_type_id	target_query	source_query	…	weight_function
100	*	pop_name==’tON’	…	adjusted_weight

or we can rebuild the edges

thalamus.add_edges(source={'pop_name': 'tON'}, target=cortex.nodes(),
                   connection_rule=5,
                   syn_weight=0.001,
                   weight_function=adjusted_weight,
                   delay=2.0,
                   target_sections=['basal', 'apical'],
                   distance_range=[0.0, 150.0],
                   dynamics_params='AMPA_ExcToExc.json',
                   model_template='exp2syn')

Then we write a custom weight function. The weight functions will be called during the simulation when building each synapse, and requires three parameters - target_cell, source_cell, and edge_props. These three parameters are dictionaries which can be used to access properties of the source node, target node, and edge, respectively. The function must return a floating point number which will be used to set the synaptic weight

def adjusted_weights(target_cell, source_cell, edge_props):
    if target_cell['cell_name'] == 'Scnn1a':
        return edge_prop["weight_max"]*0.5
    elif target_cell['cell_name'] == 'Rorb'
        return edge_prop["weight_max"]*1.5
    else:
        ...

Finally we must tell BioNet where to access the function which we can do by using the add_weight_function.

from bmtk.simulator import bionet

bionet.nrn.add_weight_function(adjusted_weights)

conf = bionet.Config.from_json('config.json')
...

Using NWB for spike trains

Instead of using csv files to set the spike trains of our external network, we can also use nwb files. The typical setup would look like the following in the config file:

{
  "inputs": {
    "LGN_spikes": {
      "input_type": "spikes",
      "module": "nwb",
      "input_file": "$INPUT_DIR/lgn_spikes.nwb",
      "node_set": "lgn",
      "trial": "trial_0"
    },
  }
}