Chapter 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 0x7f7188ed8910>

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.

psg.to_dataframe()

	node_id	population	timestamps
0	0	mthalamus	0.035761
1	0	mthalamus	0.040679
2	0	mthalamus	0.096620
3	0	mthalamus	0.244019
4	0	mthalamus	0.331973
5	0	mthalamus	0.386464
6	0	mthalamus	0.474339
7	0	mthalamus	0.532093
8	0	mthalamus	0.641981
9	0	mthalamus	0.661213
10	0	mthalamus	0.703497
11	0	mthalamus	0.773827
12	0	mthalamus	0.924455
13	0	mthalamus	1.018524
14	0	mthalamus	1.100684
15	0	mthalamus	1.202221
16	0	mthalamus	1.236319
17	0	mthalamus	1.362414
18	0	mthalamus	1.468461
19	0	mthalamus	1.518727
20	0	mthalamus	1.832532
21	0	mthalamus	1.924180
22	0	mthalamus	1.963368
23	0	mthalamus	2.007959
24	0	mthalamus	2.041782
25	0	mthalamus	2.049386
26	0	mthalamus	2.051970
27	0	mthalamus	2.064528
28	0	mthalamus	2.546912
29	0	mthalamus	2.590909
...	...	...	...
280	9	mthalamus	0.196973
281	9	mthalamus	0.201703
282	9	mthalamus	0.226206
283	9	mthalamus	0.264307
284	9	mthalamus	0.267023
285	9	mthalamus	0.662817
286	9	mthalamus	0.696093
287	9	mthalamus	0.727064
288	9	mthalamus	0.856552
289	9	mthalamus	0.869575
290	9	mthalamus	0.929920
291	9	mthalamus	0.985744
292	9	mthalamus	1.152974
293	9	mthalamus	1.266891
294	9	mthalamus	1.390516
295	9	mthalamus	1.416569
296	9	mthalamus	1.427497
297	9	mthalamus	1.563913
298	9	mthalamus	1.711437
299	9	mthalamus	1.730856
300	9	mthalamus	1.796984
301	9	mthalamus	1.821040
302	9	mthalamus	1.830675
303	9	mthalamus	1.987896
304	9	mthalamus	2.151594
305	9	mthalamus	2.313487
306	9	mthalamus	2.366667
307	9	mthalamus	2.792077
308	9	mthalamus	2.825268
309	9	mthalamus	2.873350

310 rows × 3 columns

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',
                 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
                )

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:

$ python run_bionet.py simulation_config.json

or through the script

from bmtk.simulator import bionet


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

2019-07-21 13:13:44,864 [INFO] Created log file
Mechanisms already loaded from path: /home/kael/Workspace/bmtk/docs/tutorial/sim_ch02/components/mechanisms.  Aborting.
2019-07-21 13:13:44,899 [INFO] Building cells.
2019-07-21 13:13:45,075 [INFO] Building recurrent connections
2019-07-21 13:13:45,076 [INFO] Building virtual cell stimulations for mthalamus_spikes
2019-07-21 13:13:45,151 [INFO] Running simulation for 3000.000 ms with the time step 0.100 ms
2019-07-21 13:13:45,152 [INFO] Starting timestep: 0 at t_sim: 0.000 ms
2019-07-21 13:13:45,154 [INFO] Block save every 5000 steps
2019-07-21 13:13:45,342 [INFO]     step:5000 t_sim:500.00 ms
2019-07-21 13:13:45,525 [INFO]     step:10000 t_sim:1000.00 ms
2019-07-21 13:13:45,707 [INFO]     step:15000 t_sim:1500.00 ms
2019-07-21 13:13:45,886 [INFO]     step:20000 t_sim:2000.00 ms
2019-07-21 13:13:46,068 [INFO]     step:25000 t_sim:2500.00 ms
2019-07-21 13:13:46,253 [INFO]     step:30000 t_sim:3000.00 ms
2019-07-21 13:13:46,272 [INFO] Simulation completed in 1.121 seconds

4. Analyzing the run

from bmtk.analyzer.spike_trains import to_dataframe
to_dataframe(config_file='sim_ch02/simulation_config.json')

	node_ids	population	timestamps
0	0	mcortex	45.8
1	0	mcortex	86.7
2	0	mcortex	121.9
3	0	mcortex	205.9
4	0	mcortex	252.3
5	0	mcortex	287.3
6	0	mcortex	312.7
7	0	mcortex	406.3
8	0	mcortex	483.2
9	0	mcortex	592.4
10	0	mcortex	700.0
11	0	mcortex	734.3
12	0	mcortex	874.4
13	0	mcortex	934.9
14	0	mcortex	1001.2
15	0	mcortex	1023.7
16	0	mcortex	1109.1
17	0	mcortex	1126.3
18	0	mcortex	1194.6
19	0	mcortex	1206.0
20	0	mcortex	1227.3
21	0	mcortex	1277.4
22	0	mcortex	1308.4
23	0	mcortex	1348.6
24	0	mcortex	1369.2
25	0	mcortex	1401.5
26	0	mcortex	1421.0
27	0	mcortex	1450.0
28	0	mcortex	1472.8
29	0	mcortex	1657.1
30	0	mcortex	1672.0
31	0	mcortex	1754.7
32	0	mcortex	1795.6
33	0	mcortex	1837.3
34	0	mcortex	1913.4
35	0	mcortex	1930.0
36	0	mcortex	1956.9
37	0	mcortex	1999.6
38	0	mcortex	2051.5
39	0	mcortex	2130.1
40	0	mcortex	2172.8
41	0	mcortex	2303.4
42	0	mcortex	2455.7
43	0	mcortex	2508.4
44	0	mcortex	2566.1
45	0	mcortex	2616.4
46	0	mcortex	2646.7
47	0	mcortex	2675.6
48	0	mcortex	2793.9
49	0	mcortex	2836.1
50	0	mcortex	2970.5

from bmtk.analyzer.cell_vars import plot_report

plot_report(config_file='sim_ch02/simulation_config.json')

5. Additional things to do:

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('network/mthalamus_mcortex_edge_types.csv', sep=' ')

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

(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"
    },
}