Complex Geometry

Example notebook to load a complex geometry from a swc file, and add synapses on it.

import sinaps as sn

import numpy as np
import holoviews as hv
import random

Creating Neuron structure from swc file

We choose a neuron from the neuromorpho database:

#filename = "http://neuromorpho.org/dableFiles/chang/CNG%20version/V247fs-MT-Untreated-56.CNG.swc"
filename = "https://raw.githubusercontent.com/ngltr/sinaps/master/docs/sample_data/V247fs-MT-Untreated-56.CNG.swc"

Use the function read_swc from io module to create a neuron from this file:

nrn = sn.io.read_swc(filename)
nrn.plot()

Accessing data

The neuron structure is stored using the networkx library, and you can access this underlying structure using:

nrn.graph

<networkx.classes.graph.Graph at 0x7f9761cbc4f0>

For example :

nrn.graph.edges

EdgeView([(1, 2), (1, 3), (1, 4), (1, 29), (1, 40), (1, 50), (1, 64), (1, 89), (1, 96), (4, 5), (5, 6), (5, 10), (6, 7), (7, 8), (7, 9), (10, 11), (10, 12), (12, 13), (12, 14), (14, 15), (15, 16), (15, 17), (17, 18), (18, 19), (19, 20), (19, 21), (21, 22), (22, 23), (23, 24), (24, 25), (25, 26), (26, 27), (27, 28), (29, 30), (30, 31), (31, 32), (31, 36), (32, 33), (32, 34), (34, 35), (36, 37), (37, 38), (37, 39), (40, 41), (41, 42), (41, 43), (43, 44), (44, 45), (45, 46), (45, 47), (47, 48), (48, 49), (50, 51), (51, 52), (52, 53), (52, 54), (54, 55), (54, 58), (55, 56), (56, 57), (58, 59), (59, 60), (59, 61), (61, 62), (62, 63), (64, 65), (65, 66), (66, 67), (66, 75), (67, 68), (68, 69), (68, 70), (70, 71), (70, 72), (72, 73), (73, 74), (75, 76), (75, 80), (76, 77), (77, 78), (78, 79), (80, 81), (81, 82), (81, 83), (83, 84), (83, 85), (85, 86), (86, 87), (86, 88), (89, 90), (90, 91), (91, 92), (92, 93), (93, 94), (94, 95), (96, 97), (97, 98), (97, 104), (98, 99), (98, 100), (100, 101), (100, 102), (102, 103), (104, 105), (105, 106), (106, 107), (107, 108), (108, 109), (109, 110), (110, 111), (111, 112), (111, 114), (112, 113), (114, 115), (115, 116), (116, 117), (116, 122), (117, 118), (118, 119), (119, 120), (120, 121), (122, 123), (122, 124), (124, 125), (124, 126), (126, 127)])

nrn.graph.nodes

NodeView((1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127))

You can access the section corresponding to an edge in the graph structure using the section attribute of edges

nrn.graph[1][2]['section']

Section soma_2

• L: 4.2914 um

• a: 4.316 um

• C_m: 1 uF/cm² (c_m=271.18 fF/μm)

• R_l: 150 Ω.cm (r_l=25.632 aΩ/μm)

• channels:

• point_channels:

Setting custom radius

The swc file contains the radius for each node of the dendrite. In the following, we redefine the radius of each section, assuming that this radius decreases at each bifurcation. We define for that a recursive function.

def set_radius(G,dep,a,factor):
    if G.degree(dep)>1:
        a=max(a*(G.degree(dep)-1)**(factor),0.1)
        for node in G[dep]:
            s = G[dep][node]['section']
            if s.a == 1000:# 1000 = marker
                s.a=a
                set_radius(G,node,a,factor)

nrn['dendrite'].a = 1000 #initialisation

# Setting radius of each node
set_radius(nrn.graph,1,5,-0.5)

nrn.plot()

The neuron is plotted using the coordinates of the swc file, but we can force the recalculation of the layout to optimize the graph visualization:

nrn.plot.layout(force=True)
nrn.plot().opts(node_size=0)

Calculating layout...[OK]

By defaut, the layout is calculated using the Kamada-Kawai algorithm but you can use any other node positioning algorithm setting the layoutparameter

Modeling synapses

Selecting leaves

We now randomly choose 10 leaves (excluding the soma) where we put synpases: we first remove the Hodgkin-Huxley channels, and add AMPA receptors.

We first find the node indices of the soma, to exclude it:

s = nrn['soma']
[nrn.sections[s] for s in s]

[(1, 2), (1, 3)]

Then we use the nrn.leaves function to find leaves of the neuron structure

leaves = nrn.leaves() # Finding the leaves
leaves.remove(2) # Removing the soma
leaves.remove(3) # Removing the soma
stim_leaves = random.sample(leaves,10) # getting randomly 10 leaves
stim_sec = nrn[stim_leaves] # list of the stimulated sections

stim_sec.a

[0.6250000000000002,
 0.31250000000000017,
 0.6250000000000002,
 1.2500000000000002,
 0.8838834764831847,
 0.31250000000000017,
 0.8838834764831847,
 0.31250000000000017,
 0.8838834764831847,
 0.8838834764831847]

stim_sec is a list of sections, of type section_list, on which one can operate. For example, stim_sec.a will give the radius of all sections in stim_sec. It also allows the operations on all the series in list, such as adding channels, or changing the radius.

Setting up the channels

We add Hodgkin Huxley channel everywhere in the neuron and AMPAR channel on previously chosen synapses.

The nrn[:] command select all sections and allws to see or set parameters in one command

nrn[:].clear_channels() #reset all channels
nrn[:].add_channel(sn.channels.Hodgkin_Huxley()) # add HH channel everywhere
stim_sec.clear_channels() #remove all channels on stim_sec sections
stim_sec.add_channel(sn.channels.AMPAR(0.5,gampa=0.5),1) # add NMDAR channel on stim_sec sections

Plotting stimulated leaves

nrn.add_node_data(type=0) #default node
nrn.add_node_data(stim_leaves,type=2) # synapses are put in a different color, different size.
plt = nrn.plot()
plt.opts(
    node_size = hv.dim('type')*4,
    node_color = 'red',

)

Running the simulation

# Initialisation of the simulation
sim=sn.Simulation(nrn,dx=50)

/home/docs/checkouts/readthedocs.org/user_builds/sinaps/envs/stable/lib/python3.8/site-packages/sinaps/core/model.py:402: RuntimeWarning: invalid value encountered in true_divide
  k = k / k.sum(axis=1, keepdims=True)

# Runing the simulation
sim.run((0,100))

100%|██████████| 100.0/100 [00:09<00:00, 10.32ms/s]

Plotting

# Plotting the potential
sim.plot()

sim.plot.I(sn.channels.Hodgkin_Huxley)

sim.plot.V_field()

# Extracting some sections
sim['soma'].plot()

sim.plot.I_field(sn.channels.AMPAR)