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¶
[1]:
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.
[2]:
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')
[2]:
<bmtk.builder.connection_map.ConnectionMap at 0x7f62f3911090>
[3]:
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.
[4]:
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.
[5]:
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=1.0e-04,
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.
[6]:
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 nonhomoegenous 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()
[6]:
| node_id | population | timestamps | |
|---|---|---|---|
| 0 | 0 | mthalamus | 0.001397 |
| 1 | 0 | mthalamus | 0.099114 |
| 2 | 0 | mthalamus | 0.145707 |
| 3 | 0 | mthalamus | 0.147085 |
| 4 | 0 | mthalamus | 0.162860 |
| 5 | 0 | mthalamus | 0.208985 |
| 6 | 0 | mthalamus | 0.209418 |
| 7 | 0 | mthalamus | 0.225538 |
| 8 | 0 | mthalamus | 0.232885 |
| 9 | 0 | mthalamus | 0.320728 |
| 10 | 0 | mthalamus | 0.435443 |
| 11 | 0 | mthalamus | 0.456890 |
| 12 | 0 | mthalamus | 0.484684 |
| 13 | 0 | mthalamus | 0.523754 |
| 14 | 0 | mthalamus | 0.542172 |
| 15 | 0 | mthalamus | 0.542389 |
| 16 | 0 | mthalamus | 0.549574 |
| 17 | 0 | mthalamus | 0.686661 |
| 18 | 0 | mthalamus | 0.744817 |
| 19 | 0 | mthalamus | 0.745667 |
| 20 | 0 | mthalamus | 0.814251 |
| 21 | 0 | mthalamus | 0.891552 |
| 22 | 0 | mthalamus | 0.911168 |
| 23 | 0 | mthalamus | 0.923863 |
| 24 | 0 | mthalamus | 1.206938 |
| 25 | 0 | mthalamus | 1.234080 |
| 26 | 0 | mthalamus | 1.284708 |
| 27 | 0 | mthalamus | 1.361369 |
| 28 | 0 | mthalamus | 1.376251 |
| 29 | 0 | mthalamus | 1.441294 |
| ... | ... | ... | ... |
| 4461 | 99 | mthalamus | 1.637971 |
| 4462 | 99 | mthalamus | 1.686438 |
| 4463 | 99 | mthalamus | 1.834531 |
| 4464 | 99 | mthalamus | 1.892833 |
| 4465 | 99 | mthalamus | 1.896758 |
| 4466 | 99 | mthalamus | 1.929408 |
| 4467 | 99 | mthalamus | 1.954290 |
| 4468 | 99 | mthalamus | 1.977157 |
| 4469 | 99 | mthalamus | 1.978290 |
| 4470 | 99 | mthalamus | 1.984737 |
| 4471 | 99 | mthalamus | 2.104143 |
| 4472 | 99 | mthalamus | 2.111367 |
| 4473 | 99 | mthalamus | 2.112607 |
| 4474 | 99 | mthalamus | 2.192592 |
| 4475 | 99 | mthalamus | 2.276194 |
| 4476 | 99 | mthalamus | 2.284277 |
| 4477 | 99 | mthalamus | 2.306376 |
| 4478 | 99 | mthalamus | 2.306636 |
| 4479 | 99 | mthalamus | 2.307304 |
| 4480 | 99 | mthalamus | 2.334148 |
| 4481 | 99 | mthalamus | 2.370546 |
| 4482 | 99 | mthalamus | 2.648120 |
| 4483 | 99 | mthalamus | 2.661539 |
| 4484 | 99 | mthalamus | 2.748451 |
| 4485 | 99 | mthalamus | 2.803204 |
| 4486 | 99 | mthalamus | 2.845400 |
| 4487 | 99 | mthalamus | 2.874962 |
| 4488 | 99 | mthalamus | 2.923604 |
| 4489 | 99 | mthalamus | 2.943302 |
| 4490 | 99 | mthalamus | 2.987504 |
4491 rows × 3 columns
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
[1]:
from bmtk.utils.sim_setup import build_env_bionet
build_env_bionet(base_dir='sim_ch03',
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
)
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:
$ python run_bionet.py simulation_config.json
or through the script
[1]:
from bmtk.simulator import bionet
conf = bionet.Config.from_json('sim_ch03/simulation_config.json')
conf.build_env()
net = bionet.BioNetwork.from_config(conf)
sim = bionet.BioSimulator.from_config(conf, network=net)
sim.run()
2019-07-22 09:45:15,914 [INFO] Created log file
2019-07-22 09:45:16,009 [INFO] Building cells.
2019-07-22 09:45:23,968 [INFO] Building recurrent connections
2019-07-22 09:45:24,609 [INFO] Building virtual cell stimulations for mthalamus_spikes
2019-07-22 09:45:33,506 [INFO] Running simulation for 3000.000 ms with the time step 0.100 ms
2019-07-22 09:45:33,507 [INFO] Starting timestep: 0 at t_sim: 0.000 ms
2019-07-22 09:45:33,508 [INFO] Block save every 5000 steps
2019-07-22 09:46:03,198 [INFO] step:5000 t_sim:500.00 ms
2019-07-22 09:46:33,324 [INFO] step:10000 t_sim:1000.00 ms
2019-07-22 09:47:03,400 [INFO] step:15000 t_sim:1500.00 ms
2019-07-22 09:47:33,230 [INFO] step:20000 t_sim:2000.00 ms
2019-07-22 09:48:03,147 [INFO] step:25000 t_sim:2500.00 ms
2019-07-22 09:48:33,519 [INFO] step:30000 t_sim:3000.00 ms
2019-07-22 09:48:36,802 [INFO] Simulation completed in 3.0 minutes, 3.297 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:
[2]:
from bmtk.analyzer.spike_trains import plot_raster
plot_raster(config_file='sim_ch03/simulation_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:
[5]:
from bmtk.analyzer.cell_vars import plot_report
plot_report(config_file='sim_ch03/simulation_config.json', node_ids=[0, 49, 99])
5. Modifying the network¶
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:
[20]:
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: [ -1.99484437 41.49527042 22.33923077], angle_yaxis: 5.29759513272
cell 1: pop_name: Scnn1a, positions: [-25.72073426 36.01835631 2.43526216], angle_yaxis: 2.94311607964
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.
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>},
...
[ ]: