__author__ = 'Aaron D. Milstein'
# Includes modification of an early version of SWC_neuron.py by Daniele Linaro.
# Includes an extension of BtMorph, created by Ben Torben-Nielsen and modified by Daniele Linaro.
from function_lib import *
from neuron import h # must be found in system $PYTHONPATH
import pprint
import btmorph # must be found in system $PYTHONPATH
# SWC files must use this nonstandard convention to exploit trunk and tuft categorization
swc_types = [soma_type, axon_type, basal_type, apical_type, trunk_type, tuft_type] = [1, 2, 3, 4, 5, 6]
sec_types = ['soma', 'axon_hill', 'ais', 'axon', 'basal', 'trunk', 'apical', 'tuft', 'spine_neck', 'spine_head']
verbose = 0 # Turn on for text reporting during model initialization and simulation
gid_max = 0 # Every new HocCell will receive a global identifier for network simulation, and increment this counter
#-------Wrapper for converting SWC --> BtMorph --> Skeleton morphology in NEURON hoc------------------------
class HocCell(object):
def __init__(self, morph_filename=None, mech_filename=None):
"""
:param morph_filename: str : path to .swc file containing morphology
:param mech_filename: str : path to .pkl file specifying cable parameters and membrane mechanisms
"""
global gid_max
self._gid = gid_max
gid_max += 1
self.tree = btmorph.STree2() # Builds a simple tree to store nodes of type 'SHocNode'
self.index = 0 # Keep track of number of nodes
self._node_dict = {'soma': [], 'axon': [], 'basal': [], 'trunk': [], 'apical': [], 'tuft': [], 'spine': []}
self.mech_dict = self.load_mech_dict(mech_filename) # Refer to function_lib for description of structure of
# mechanism dictionary. loads from .pkl or
# default_mech_dict in function_lib
if not morph_filename is None:
self.load_morphology_from_swc(morph_filename)
self.reinit_mechanisms() # Membrane mechanisms must be reinitialized whenever cable properties (Ra, cm) or
# spatial resolution (nseg) changes.
self.spike_detector = None
if self.axon:
self.init_spike_detector()
self.random = np.random.RandomState()
def load_morphology_from_swc(self, morph_filename):
"""
This method builds an STree2 comprised of SHocNode nodes associated with hoc sections, connects the hoc
sections, and initializes various parameters: Ra, cm, L, diam, nseg
This method implements a standardized soma and axon:
The soma consists of two cylindrical hoc sections of equal length and diameter, connected (0) to (0).
The basal dendritic tree is connected to soma[1](1), and the apical tree is connected to soma[0](1).
The axon is attached to soma[0](0), and consists of three sequential hoc sections:
1) axon[0] : a tapered cylindrical 'axon hillock' section connected to soma[0](0)
2) axon[1] : a tapered cylindrical 'axon initial segment' section connected to axon[0](1)
3) axon[2] : a cylindrical 'axon' section connected to axon[1](1)
"""
raw_tree = btmorph.STree2() # import the full tree from an SWC file
raw_tree.read_SWC_tree_from_file(morph_dir+morph_filename, types=range(10))
soma_length = 14.
soma_diam = 9.
for index in range(2):
node = self.make_section('soma')
node.sec.L = soma_length/2.
node.sec.diam = soma_diam
self._init_cable(node) # consults the mech_dict to initialize Ra, cm, and nseg
self.tree.root = self.soma[0]
self.soma[1].connect(self.soma[0], 0, 0)
for index in range(3):
self.make_section('axon')
self.axon[0].type = 'axon_hill'
self.axon[0].sec.L = 10.
self.axon[0].set_diam_bounds(3., 2.) # stores the diameter boundaries for a tapered cylindrical section
self.axon[1].type = 'ais'
self.axon[1].sec.L = 15.
self.axon[1].set_diam_bounds(2., 0.5) # (2., 1.)
self.axon[2].sec.L = 500.
self.axon[2].sec.diam = 0.5 # 1.
self.axon[0].connect(self.soma[0], 0, 0)
self.axon[1].connect(self.axon[0], 1, 0)
self.axon[2].connect(self.axon[1], 1, 0)
for node in self.axon:
self._init_cable(node)
for child in raw_tree.root.children:
self.make_skeleton(child, self.tree.root)
def make_section(self, sec_type):
"""
Create a new hoc section to associate with this node, and this cell, and store information about it in the
node's content dictionary.
:param sec_type: str
:return node: :class:'SHocNode'
"""
node = SHocNode(self.index)
if self.index == 0:
self.tree.root = node
self.index += 1
node.type = sec_type
if sec_type in ['spine_head', 'spine_neck']:
self._node_dict['spine'].append(node)
else:
self._node_dict[sec_type].append(node)
node.sec = h.Section(name=node.name, cell=self)
return node
def test_sec_properties(self, node=None):
"""
Used for debugging and validating model specification.
:param node:
"""
if node is None:
node = self.tree.root
# node.sec.push()
# h.psection()
# h.pop_section()
print '{} [nseg: {}, Ra: {}]'.format(node.name, node.sec.nseg, node.sec.Ra)
# h('for (x) print (x), diam(x)', sec=node.sec)
for child in node.children:
self.test_sec_properties(child)
def make_skeleton(self, raw_node, parent, length=0, diams=None):
"""
Following construction of soma and axon nodes of type 'SHocNode' in the tree of type 'STree2', this method
recursively converts dendritic 'SNode2' nodes into 'SHocNode' nodes, and connects them to the appropriate
somatic nodes. Skeletonized dendritic nodes have only one hoc section for each stretch of unbranched dendrite,
with length equal to the sum of the lengths of the converted SNode2 nodes.
Nodes that taper more than 1 um remain tapered, otherwise they are converted into untapered cylinders with
diameter equal to the mean diameter of the the converted SNode2 nodes.
Dendrite types that are pre-categorized as basal, apical, trunk, or tuft in the input .swc file are preserved.
:param raw_node: :class:'SNode2'
:param parent: :class:'SHocNode'
:param length: int or float
:param diams: None or (list: float)
"""
global verbose
dend_types=([basal_type, apical_type, trunk_type, tuft_type], ['basal', 'apical', 'trunk', 'tuft'])
swc = raw_node.get_content()['p3d']
swc_type = swc.type
if swc_type in dend_types[0]:
diam = swc.radius*2
length += self.get_node_length_swc(raw_node)
leaves = len(raw_node.children)
# create a new node when encountering 1) branch points, 2) terminal ends, or 3) change in swc_type
if leaves > 1 or leaves == 0 or \
(leaves == 1 and not raw_node.children[0].get_content()['p3d'].type == swc_type):
sec_type = dend_types[1][dend_types[0].index(swc_type)]
new_node = self.make_section(sec_type)
new_node.sec.L = length
if (self.tree.is_root(parent)) and (sec_type == 'basal'):
parent = self.soma[1]
new_node.connect(parent)
if diams is None:
new_node.sec.diam = diam
self._init_cable(new_node)
if verbose:
print '{} [nseg: {}, diam: {}, length: {}, parent: {}]'.format(new_node.name, new_node.sec.nseg,
diam, length, new_node.parent.name)
else:
diams.append(diam)
if len(diams) > 2:
mean = np.mean(diams)
stdev = np.std(diams)
if stdev*2 > 1: # If 95% of the values are within 1 um, don't taper
new_node.set_diam_bounds(mean+stdev, mean-stdev)
self._init_cable(new_node)
if verbose:
print '{} [nseg: {}, diam: ({}:{}), length: {}, parent: {}]'.format(new_node.name,
new_node.sec.nseg, mean+stdev, mean-stdev, length, new_node.parent.name)
else:
new_node.sec.diam = mean
self._init_cable(new_node)
if verbose:
print '{} [nseg: {}, diam: {}, length: {}, parent: {}]'.format(new_node.name,
new_node.sec.nseg, mean, length, new_node.parent.name)
elif abs(diams[0]-diams[1]) > 1:
new_node.set_diam_bounds(diams[0], diams[1])
self._init_cable(new_node)
if verbose:
print '{} [diam: ({}:{}), length: {}, parent: {}]'.format(new_node.name, new_node.sec.nseg,
diams[0], diams[1], length, new_node.parent.name)
else:
mean = np.mean(diams)
new_node.sec.diam = mean
self._init_cable(new_node)
if verbose:
print '{} [nseg: {}, diam: {}, length: {}, parent: {}]'.format(new_node.name,
new_node.sec.nseg, mean, length, new_node.parent.name)
# Follow all branches from this fork
for child in raw_node.children:
self.make_skeleton(child, new_node)
else: # Follow unbranched dendrite
if diams is None:
diams = [diam]
else:
diams.append(diam)
self.make_skeleton(raw_node.children[0], parent, length, diams)
def get_nodes_of_subtype(self, sec_type):
"""
This method searches the node dictionary for nodes of a given sec_type and returns them in a list. Used during
specification of membrane mechanisms.
:param sec_type: str
:return: list of :class:'SHocNode'
"""
if sec_type in ['axon_hill', 'ais', 'axon']:
return [node for node in self.axon if node.type == sec_type]
elif sec_type in ['spine_head', 'spine_neck']:
return [node for node in self.spine if node.type == sec_type]
else:
return self._node_dict[sec_type]
def load_mech_dict(self, mech_filename=None):
"""
This method loads the dictionary specifying membrane mechanism parameters. If a .pkl file is not provided, a
global variable default_mech_dict from function_lib is used.
:param mech_filename: str
"""
if not mech_filename is None:
return read_from_pkl(data_dir+mech_filename+'.pkl')
else:
local_mech_dict = copy.deepcopy(default_mech_dict)
return local_mech_dict
def _init_cable(self, node):
"""
If the mechanism dictionary specifies the cable properties 'Ra' or 'cm', then _modify_mechanism() properly sets
those parameters, and reinitializes the number of segments per section. To avoid redundancy, this
method passes _modify_mechanism() a copy of the dictionary with the spatial_res parameter removed, since this is
consulted in setting nseg. However, if spatial_res is the only parameter being specified, it is passed to
_modify_mechanism()
:param node: :class:'SHocNode'
"""
sec_type = node.type
if sec_type in self.mech_dict and 'cable' in self.mech_dict[sec_type]:
mech_content = copy.deepcopy(self.mech_dict[sec_type]['cable'])
if ('Ra' in mech_content) or ('cm' in mech_content):
if 'spatial_res' in mech_content:
del mech_content['spatial_res']
self._modify_mechanism(node, 'cable', mech_content)
elif 'spatial_res' in mech_content:
self._modify_mechanism(node, 'cable', mech_content)
else:
node.init_nseg()
node.reinit_diam()
def reinit_mechanisms(self, reset_cable=0):
"""
Once a mechanism dictionary has been loaded, and a morphology has been specified, this method traverses through
the tree of SHocNode nodes following order of inheritance and properly sets membrane mechanism parameters,
including gradients and inheritance of parameters from nodes along the path from root. Since cable parameters
are set during specification of morphology, it is not necessary to immediately reinitialize these parameters
again. However, they can be manually reinitialized with the reset_cable flag.
:param reset_cable: boolean
"""
for sec_type in sec_types:
if sec_type in self.mech_dict:
nodes = self.get_nodes_of_subtype(sec_type)
self._reinit_mech(nodes, reset_cable)
def init_synaptic_mechanisms(self):
"""
Spines and synapses are inserted after loading a morphology and specifying membrane mechanisms. This method can
be executed after inserting synapses. It traverses the dendritic tree in order of inheritance and just sets
synaptic mechanism parameters specified in the mechanism dictionary.
"""
for dend_type in ['soma', 'basal', 'trunk', 'apical', 'tuft']:
if dend_type in self.mech_dict and 'synapse' in self.mech_dict[dend_type] and \
self.sec_type_has_synapses(dend_type):
for node in self.get_nodes_of_subtype(dend_type):
self._modify_mechanism(node, 'synapse', self.mech_dict[dend_type]['synapse'])
def node_has_synapses(self, node, syn_type=None):
"""
Checks if a given node contains synapses, or spines with synapses. Can also check for a synaptic point process
of a specific type.
:param node: :class:'SHocNode'
:param syn_type: str
:return: boolean
"""
if [syn for syn in node.synapses if syn_type is None or syn_type in syn._syn]:
return True
else:
for spine in node.spines:
if [syn for syn in spine.synapses if syn_type is None or syn_type in syn._syn]:
return True
return False
def sec_type_has_synapses(self, sec_type, syn_type=None):
"""
Checks if any nodes of a given sec_type contain synapses, or spines with synapses. Can also check for a synaptic
point process of a specific type.
:param sec_type: str
:param syn_type: str
:return: boolean
"""
for node in self.get_nodes_of_subtype(sec_type):
if self.node_has_synapses(node, syn_type):
return True
return False
def _reinit_mech(self, nodes, reset_cable=0):
"""
Given a list of nodes, this method loops through all the mechanisms specified in the mechanism dictionary for
the hoc section type of each node and updates their associated parameters. If the reset_cable flag is set to 1,
cable parameters are modified first, then the parameters for all other mechanisms are reinitialized.
Parameters for synaptic point processes can also be specified in the mechanism dictionary, so one must use the
method init_synaptic_mechanisms() after inserting synapses.
:param nodes: list of :class:'SHocNode'
:param reset_cable: int or boolean
"""
for node in nodes:
sec_type = node.type
if sec_type in self.mech_dict:
if ('cable' in self.mech_dict[sec_type]) and reset_cable: # cable properties must be set first, as they
# can change nseg, which will affect insertion
# of membrane mechanism gradients
self._init_cable(node)
for mech_name in (mech_name for mech_name in self.mech_dict[sec_type]
if not mech_name in ['cable', 'ions']):
self._modify_mechanism(node, mech_name, self.mech_dict[sec_type][mech_name])
# ion-related parameters do not exist until after membrane mechanisms have been inserted
if 'ions' in self.mech_dict[sec_type]:
self._modify_mechanism(node, 'ions', self.mech_dict[sec_type]['ions'])
def reinitialize_subset_mechanisms(self, sec_type, mech_name):
"""
During parameter optimization, it is often convenient to reinitialize all the parameters for a single mechanism
in a subset of compartments. For example, g_pas in basal dendrites that inherit the value from the soma after
modifying the value in the soma compartment.
:param sec_type: str
:param mech_name: str
:return:
"""
if sec_type in self.mech_dict and mech_name in self.mech_dict[sec_type]:
for node in self.get_nodes_of_subtype(sec_type):
self._modify_mechanism(node, mech_name, self.mech_dict[sec_type][mech_name])
def _modify_mechanism(self, node, mech_name, mech_content):
"""
This method loops through all the parameters for a single mechanism specified in the mechanism dictionary and
calls self._parse_mech_content to interpret the rules and set the values for the given node.
:param node: :class:'SHocNode'
:param mech_name: str
:param mech_content: dict
"""
if not mech_content is None:
# only modify synaptic mechanism parameters if synapses have been inserted
if mech_name == 'synapse':
if self.node_has_synapses(node):
for syn_type in mech_content:
for param_name in mech_content[syn_type]:
# accommodate multiple dict entries with different location constraints for a single
# parameter
if type(mech_content[syn_type][param_name]) == dict:
self._parse_mech_content(node, mech_name, param_name,
mech_content[syn_type][param_name], syn_type)
else:
for mech_content_entry in mech_content[syn_type][param_name]:
self._parse_mech_content(node, mech_name, param_name, mech_content_entry, syn_type)
else:
for param_name in mech_content:
# accommodate multiple dict entries with different location constraints for a single parameter
if type(mech_content[param_name]) == dict:
self._parse_mech_content(node, mech_name, param_name, mech_content[param_name])
else:
for mech_content_entry in mech_content[param_name]:
self._parse_mech_content(node, mech_name, param_name, mech_content_entry)
else:
node.sec.insert(mech_name)
def _parse_mech_content(self, node, mech_name, param_name, rules, syn_type=None):
"""
This method loops through all the segments in a node and sets the value(s) for a single mechanism parameter by
interpreting the rules specified in the mechanism dictionary. Properly handles ion channel gradients and
inheritance of values from the closest segment of a specified type of section along the path from root. Also
handles rules with distance boundaries, and rules to set synaptic (point process) parameters. For gradients,
specifying a slope implies a linear gradient, while specifying both a slope and a tau implies an exponential
gradient.
:param node: :class:'SHocNode'
:param mech_name: str
:param param_name: str
:param rules: dict
:param syn_type: str
"""
if mech_name == 'synapse':
if syn_type is None:
raise Exception('Cannot set synaptic mechanism parameter: {} without a specified point process'.
format(param_name))
elif not self.node_has_synapses(node, syn_type):
return None # ignore mechanism dictionary entries for types of synapses that have not been inserted
if 'origin' in rules: # an 'origin' with no 'value' inherits a starting parameter from the origin sec_type
# a 'value' with no 'origin' is independent of other sec_types
# an 'origin' with a 'value' uses the origin sec_type only as a reference point for
# applying a distance-dependent gradient
if rules['origin'] == 'parent':
if node.type == 'spine_head':
donor = node.parent.parent
elif node.type == 'spine_neck':
donor = node.parent
else:
donor = self.get_dendrite_origin(node)
elif rules['origin'] in sec_types:
donor = self._get_node_along_path_to_root(node, rules['origin'])
else:
raise Exception('Mechanism: {} parameter: {} cannot reference unknown sec_type: {}'.format(mech_name,
param_name, rules['origin']))
else:
donor = None
if 'value' in rules:
baseline = rules['value']
elif donor is None:
raise Exception('Cannot set mechanism: {} parameter: {} without a specified origin or value'.format(
mech_name, param_name))
else:
if (mech_name == 'cable') and (param_name == 'spatial_res'):
baseline = self._get_spatial_res(donor)
elif mech_name == 'synapse':
if self.sec_type_has_synapses(donor.type, syn_type):
baseline = self._inherit_mech_param(node, donor, mech_name, param_name, syn_type)
else:
raise Exception('Cannot inherit synaptic mechanism: {} parameter: {} from sec_type: {}'.format(
syn_type, param_name, donor.type))
else:
baseline = self._inherit_mech_param(node, donor, mech_name, param_name)
if mech_name == 'cable': # cable properties can be inherited, but cannot be specified as gradients
if param_name == 'spatial_res':
node.init_nseg(baseline)
else:
setattr(node.sec, param_name, baseline)
node.init_nseg(self._get_spatial_res(node))
node.reinit_diam()
else:
if 'min_loc' in rules or 'max_loc' in rules or 'slope' in rules:
if donor is None:
raise Exception('Cannot follow specifications for mechanism: {} parameter: {} without a provided '
'origin'.format(mech_name, param_name))
if mech_name == 'synapse':
self._specify_synaptic_parameter(node, param_name, baseline, rules, syn_type, donor)
else:
if 'min_loc' in rules:
min_distance = rules['min_loc']
else:
min_distance = None
if 'max_loc' in rules:
max_distance = rules['max_loc']
else:
max_distance = None
min_seg_distance = self.get_distance_to_node(donor, node, 0.5/node.sec.nseg)
max_seg_distance = self.get_distance_to_node(donor, node, (0.5 + node.sec.nseg - 1)/node.sec.nseg)
# if any part of the section is within the location constraints, insert the mechanism, and specify
# the parameter at the segment level
if (min_distance is None or max_seg_distance >= min_distance) and \
(max_distance is None or min_seg_distance <= max_distance):
if not mech_name == 'ions':
node.sec.insert(mech_name)
if min_distance is None:
min_distance = 0.
for seg in node.sec:
seg_loc = self.get_distance_to_node(donor, node, seg.x)
if seg_loc >= min_distance and (max_distance is None or seg_loc <= max_distance):
if 'slope' in rules:
seg_loc -= min_distance
if 'tau' in rules:
if 'xhalf' in rules: # sigmoidal gradient
offset = baseline - rules['slope'] / (1. +
np.exp(rules['xhalf'] / rules['tau']))
value = offset + rules['slope'] / (1. +
np.exp((rules['xhalf'] - seg_loc) / rules['tau']))
else: # exponential gradient
offset = baseline - rules['slope']
value = offset + rules['slope'] * np.exp(seg_loc / rules['tau'])
else: # linear gradient
value = baseline + rules['slope'] * seg_loc
if 'min' in rules and value < rules['min']:
value = rules['min']
#elif value < 0.:
# value = 0.
elif 'max' in rules and value > rules['max']:
value = rules['max']
else:
value = baseline
elif 'outside' in rules: # by default, if only some segments in a section meet the
value = rules['outside'] # location constraints, the parameter inherits the
else: # mechanism's default value. if another value is desired, it
value = None # can be specified via an 'outside' key in the mechanism
if not value is None: # dictionary entry
if mech_name == 'ions':
setattr(seg, param_name, value)
else:
setattr(getattr(seg, mech_name), param_name, value)
elif mech_name == 'ions':
setattr(node.sec, param_name, baseline)
elif mech_name == 'synapse':
self._specify_synaptic_parameter(node, param_name, baseline, rules, syn_type)
else:
node.sec.insert(mech_name)
setattr(node.sec, param_name+"_"+mech_name, baseline)
def _specify_synaptic_parameter(self, node, param_name, baseline, rules, syn_type, donor=None):
"""
This method interprets an entry from the mechanism dictionary to set parameters associated with a synaptic
point_process mechanism that has been inserted either into a spine attached to this node, or inserted directly
in this node. Appropriately implements slopes and inheritances.
:param node: :class:'SHocNode'
:param param_name: str
:param baseline: float
:param rules: dict
:param syn_type: str
:param donor: :class:'SHocNode' or None
"""
syn_list = []
syn_list.extend(node.synapses)
for spine in node.spines:
syn_list.extend(spine.synapses)
if 'min_loc' in rules:
min_distance = rules['min_loc']
else:
min_distance = 0.
if 'max_loc' in rules:
max_distance = rules['max_loc']
else:
max_distance = None
if 'variance' in rules and rules['variance'] == 'normal':
normal = True
else:
normal = False
for syn in syn_list:
if syn_type in syn._syn: # not all synapses contain every synaptic mechanism
target = syn.target(syn_type)
if donor is None:
if normal:
value = self.random.normal(baseline, baseline / 6.)
setattr(target, param_name, value)
else:
setattr(target, param_name, baseline)
else:
distance = self.get_distance_to_node(donor, node, syn.loc)
# note: if only some synapses in a section meet the location constraints, the synaptic parameter
# will maintain its default value in all other locations. values for other locations must be
# specified with an additional entry in the mechanism dictionary
if distance >= min_distance and (max_distance is None or distance <= max_distance):
if 'slope' in rules:
distance -= min_distance
if 'tau' in rules: # exponential gradient
if 'xhalf' in rules: # sigmoidal gradient
offset = baseline - rules['slope'] / (1. + np.exp(rules['xhalf'] / rules['tau']))
value = offset + rules['slope'] / (1. + np.exp((rules['xhalf'] - distance) /
rules['tau']))
else:
offset = baseline - rules['slope']
value = offset + rules['slope'] * np.exp(distance / rules['tau'])
else: # linear gradient
value = baseline + rules['slope'] * distance
if 'min' in rules and value < rules['min']:
value = rules['min']
#elif value < 0.:
# value = 0.
elif 'max' in rules and value > rules['max']:
value = rules['max']
else:
value = baseline
if normal:
value = self.random.normal(value, value / 6.)
setattr(target, param_name, value)
def set_special_mech_param_linear_gradient(self, mech_name, param_name, sec_type_list, criterion, end_val):
"""
This is an admittedly ad-hoc procedure to implement a linearly decreasing gradient of sodium channels in
terminal branches that is not easily accomplished by the general procedures implementing the mechanism
dictionary.
:param mech_name: str
:param param_name: str
:param sec_type_list: list of str
:param criterion: boolean function
:param end_val: float
"""
for sec_type in sec_type_list:
for node in self.get_nodes_of_subtype(sec_type):
if criterion(node):
start_val = getattr(node.sec(0.), param_name+'_'+mech_name)
slope = end_val - start_val
for seg in node.sec:
value = start_val + slope * seg.x
setattr(getattr(seg, mech_name), param_name, value)
def init_spike_detector(self, node=None, loc=1., param='_ref_v', delay=None, weight=None, threshold=None,
target=None):
"""
Converts analog voltage in the specified section to digital spike output. By default, initializes an h.NetCon
object with voltage as a reference parameter and no target. Can later be connected by h.ParallelContext, or
re-initialized with a target on a cell contained within the local processing environment.
:param node: :class:'SHocNode'
:param loc: float
:param param: str
:param delay: float
:param weight: float
:param threshold: float
:return: :class:'h.NetCon'
"""
if node is None:
if self.axon:
node = self.axon[-1]
else:
raise Exception('No source node specified for spike detector.')
if self.spike_detector is not None:
if delay is None:
delay = self.spike_detector.delay
if weight is None:
weight = self.spike_detector.weight[0]
if threshold is None:
threshold = self.spike_detector.threshold
else:
if delay is None:
delay = 0.
if weight is None:
weight = 1.
if threshold is None:
threshold = -30.
self.spike_detector = h.NetCon(getattr(node.sec(loc), param), target, sec=node.sec)
self.spike_detector.delay = delay
self.spike_detector.weight[0] = weight
self.spike_detector.threshold = threshold
def get_dendrite_origin(self, node):
"""
This method determines the section type of the given node, and returns the node representing the primary branch
point for the given section type. Basal and trunk sections originate at the soma, and apical and tuft dendrites
originate at the trunk. For spines, recursively calls with parent node to identify the parent branch first.
:param node: :class:'SHocNode'
:return: :class:'SHocNode'
"""
sec_type = node.type
if sec_type in ['spine_head', 'spine_neck']:
return self.get_dendrite_origin(node.parent)
elif sec_type in ['basal', 'trunk', 'axon_hill', 'ais', 'axon']:
return self._get_node_along_path_to_root(node, 'soma')
elif sec_type in ['apical', 'tuft']:
return self._get_node_along_path_to_root(node, 'trunk')
elif sec_type == 'soma':
return node
def _get_node_along_path_to_root(self, node, sec_type):
"""
This method follows the path from the given node to the root node, and returns the first node with section type
sec_type.
:param node: :class:'SHocNode'
:param sec_type: str
:return: :class:'SHocNode'
"""
parent = node
while not parent is None:
if parent in self.soma and not sec_type == 'soma':
parent = None
elif parent.type == sec_type:
return parent
else:
parent = parent.parent
raise Exception('The path from node: {} to root does not contain sections of type: {}'.format(node.name,
sec_type))
def _get_closest_synapse(self, node, loc, syn_type=None, downstream=True):
"""
This method finds the closest synapse to the specified location within or downstream of the provided node. Used
for inheritance of synaptic mechanism parameters. Can also look upstream instead. Can also find the closest
synapse containing a synaptic point_process of a specific type.
:param node: :class:'SHocNode'
:param loc: float
:param syn_type: str
:return: :class:'Synapse'
"""
syn_list = [syn for syn in node.synapses if syn_type is None or syn_type in syn._syn]
for spine in node.spines:
syn_list.extend([syn for syn in spine.synapses if syn_type is None or syn_type in syn._syn])
if not syn_list:
if downstream:
for child in [child for child in node.children if child.type == node.type]:
target_syn = self._get_closest_synapse(child, 0., syn_type)
if target_syn is not None:
return target_syn
return None
elif node.parent.type == node.type:
return self._get_closest_synapse(node.parent, 1., syn_type, downstream=False)
else:
min_distance = 1.
target_syn = None
for syn in syn_list:
distance = abs(syn.loc - loc)
if distance < min_distance:
min_distance = distance
target_syn = syn
return target_syn
def _inherit_mech_param(self, node, donor, mech_name, param_name, syn_type=None):
"""
When the mechanism dictionary specifies that a node inherit a parameter value from a donor node, this method
returns the value of that parameter found in the section or final segment of the donor node. For synaptic
mechanism parameters, searches for the closest synapse in the donor node. If the donor node does not contain
synapses due to location constraints, this method searches first child branches, then parent nodes of the same
sec_type as the donor node.
:param node: :class:'SHocNode'
:param donor: :class:'SHocNode'
:param mech_name: str
:param param_name: str
:param syn_type: str
:return: float
"""
try:
if mech_name in ['cable', 'ions']:
return getattr(donor.sec, param_name)
elif mech_name == 'synapse':
try: # first look downstream for a nearby synapse, then upstream.
return getattr(self._get_closest_synapse(donor, 1., syn_type).target(syn_type), param_name)
except (AttributeError, KeyError):
return getattr(self._get_closest_synapse(donor, 1., syn_type, downstream=False).target(syn_type),
param_name)
else:
loc = donor.sec.nseg/(donor.sec.nseg + 1.) # accesses the last segment of the section
return getattr(getattr(donor.sec(loc), mech_name), param_name)
except (AttributeError, NameError, KeyError):
if syn_type is None:
print 'Exception: Mechanism: {} parameter: {} cannot be inherited from sec_type: {}'.format(mech_name,
param_name, donor.type)
else:
print 'Exception: Problem inheriting synaptic mechanism: {} parameter {} from sec_type: {}'.format(
syn_type, param_name, donor.type)
raise KeyError
def _get_spatial_res(self, node):
"""
Checks the mechanism dictionary if the section type of this node has a specified spatial resolution factor.
Used to scale the number of segments per section in the hoc model by a factor of an exponent of 3.
:param node: :class:'SHocNode
:return: int
"""
try: # if spatial_res has not been specified for the origin type of section, it defaults to 0
rules = self.mech_dict[node.type]['cable']['spatial_res']
except KeyError:
return 0
if 'value' in rules:
return rules['value']
elif 'origin' in rules:
if rules['origin'] in sec_types: # if this sec_type also inherits the value, continue following the path
return self._get_spatial_res(self._get_node_along_path_to_root(node, rules['origin']))
else:
print 'Exception: Spatial resolution cannot be inherited from sec_type: {}'.format(rules['origin'])
raise KeyError
else:
print 'Exception: Cannot set spatial resolution without a specified origin or value'
raise KeyError
def modify_mech_param(self, sec_type, mech_name, param_name=None, value=None, origin=None, slope=None, tau=None,
xhalf=None, min=None, max=None, min_loc=None, max_loc=None, outside=None, syn_type=None, variance=None,
replace=True):
"""
Modifies or inserts new membrane mechanisms into hoc sections of type sec_type. First updates the mechanism
dictionary, then sets the corresponding hoc parameters. This method is meant to be called manually during
initial model specification, or during parameter optimization. For modifications to persist across simulations,
the mechanism dictionary must be saved to a file using self.export_mech_dict() and re-imported during HocCell
initialization.
:param sec_type: str
:param mech_name: str
:param param_name: str
:param value: float
:param origin: str
:param slope: float
:param tau: float
:param xhalf: float
:param min: float
:param max: float
:param min_loc: float
:param max_loc: float
:param outside: float
:param syn_type: str
:param variance: str
:param replace: bool
"""
global verbose
backup_content = None
mech_content = None
if not sec_type in sec_types:
raise Exception('Cannot specify mechanism: {} parameter: {} for unknown sec_type: {}'.format(mech_name,
param_name, sec_type))
if mech_name in ['cable', 'ions', 'synapse']:
if param_name is None:
raise Exception('No parameter specified for mechanism: {}'.format(mech_name))
if not param_name is None:
if value is None and origin is None:
raise Exception('Cannot set mechanism: {} parameter: {} without a specified origin or value'.format(
mech_name, param_name))
if mech_name == 'synapse':
if syn_type is None:
raise Exception('Cannot set synaptic mechanism parameter: {} without a specified point '
'process'.format(param_name))
else:
self._modify_synaptic_mech_param(sec_type, param_name, value, origin, slope, tau, xhalf, min, max,
min_loc, max_loc, outside, syn_type, variance, replace)
return
rules = {}
if not origin is None:
if not origin in sec_types+['parent']:
raise Exception('Cannot inherit mechanism: {} parameter: {} from unknown sec_type: {}'.format(
mech_name, param_name, origin))
else:
rules['origin'] = origin
if not value is None:
rules['value'] = value
if not slope is None:
rules['slope'] = slope
if not tau is None:
rules['tau'] = tau
if not xhalf is None:
rules['xhalf'] = xhalf
if not min is None:
rules['min'] = min
if not max is None:
rules['max'] = max
if not min_loc is None:
rules['min_loc'] = min_loc
if not max_loc is None:
rules['max_loc'] = max_loc
if not outside is None:
rules['outside'] = outside
# currently only implemented for synaptic parameters
if not variance is None:
rules['variance'] = variance
mech_content = {param_name: rules}
if not sec_type in self.mech_dict: # No mechanisms have been inserted into this type of section yet
self.mech_dict[sec_type] = {mech_name: mech_content}
elif not mech_name in self.mech_dict[sec_type]: # This mechanism has not yet been inserted into
backup_content = copy.deepcopy(self.mech_dict[sec_type]) # this type of section
self.mech_dict[sec_type][mech_name] = mech_content
elif self.mech_dict[sec_type][mech_name] is None: # This mechanism has been inserted, but no
backup_content = copy.deepcopy(self.mech_dict[sec_type]) # parameters have been specified
self.mech_dict[sec_type][mech_name] = mech_content
elif param_name in self.mech_dict[sec_type][mech_name]: # This parameter has already been specified. Now
backup_content = copy.deepcopy(self.mech_dict[sec_type]) # have to determine whether to replace or extend
if replace: # the current dictionary entry.
self.mech_dict[sec_type][mech_name][param_name] = rules
elif type(self.mech_dict[sec_type][mech_name][param_name]) == dict:
self.mech_dict[sec_type][mech_name][param_name] = [self.mech_dict[sec_type][mech_name][param_name],
rules]
elif type(self.mech_dict[sec_type][mech_name][param_name]) == list:
self.mech_dict[sec_type][mech_name][param_name].append(rules)
elif not param_name is None: # This mechanism has been inserted, but this
backup_content = copy.deepcopy(self.mech_dict[sec_type]) # parameter has not yet been specified
self.mech_dict[sec_type][mech_name][param_name] = rules
try:
nodes = self.get_nodes_of_subtype(sec_type)
if mech_name == 'cable': # all membrane mechanisms in sections of type sec_type must be reinitialized after
# changing cable properties
if param_name in ['Ra', 'cm', 'spatial_res']:
self._reinit_mech(nodes, reset_cable=1)
else:
print 'Exception: Unknown cable property: {}'.format(param_name)
raise KeyError
else:
for node in nodes:
try:
self._modify_mechanism(node, mech_name, mech_content)
except (AttributeError, NameError, ValueError, KeyError):
if not param_name is None:
print 'Exception: Problem modifying mechanism: {} parameter: {}'.format(mech_name,
param_name)
else:
print 'Exception: Problem inserting mechanism: {}'.format(mech_name)
raise KeyError
except KeyError:
if backup_content is None:
del self.mech_dict[sec_type]
else:
self.mech_dict[sec_type] = copy.deepcopy(backup_content)
finally:
if verbose:
pprint.pprint(self.mech_dict)
def _modify_synaptic_mech_param(self, sec_type, param_name=None, value=None, origin=None, slope=None, tau=None,
xhalf=None, min=None, max=None, min_loc=None, max_loc=None, outside=None, syn_type=None,
variance=None, replace=True):
"""
Modifies or inserts new synaptic point process mechanisms into hoc sections of type sec_type. First updates the
mechanism dictionary, then sets the corresponding hoc parameters. Internal method handles special nested
dictionary specification for synaptic parameters.
:param sec_type: str
:param mech_name: str
:param param_name: str
:param value: float
:param origin: str
:param slope: float
:param tau: float
:param xhalf: float
:param min: float
:param max: float
:param min_loc: float
:param max_loc: float
:param outside: float
:param syn_type: str
:param variance: str
:param replace: bool
"""
global verbose
backup_content = None
rules = {}
if not origin is None:
if not origin in sec_types + ['parent']:
raise Exception('Cannot inherit synaptic mechanism: {} parameter: {} from unknown sec_type: {}'.format(
syn_type, param_name, origin))
else:
rules['origin'] = origin
if not value is None:
rules['value'] = value
if not slope is None:
rules['slope'] = slope
if not tau is None:
rules['tau'] = tau
if not xhalf is None:
rules['xhalf'] = xhalf
if not min is None:
rules['min'] = min
if not max is None:
rules['max'] = max
if not min_loc is None:
rules['min_loc'] = min_loc
if not max_loc is None:
rules['max_loc'] = max_loc
if not outside is None:
rules['outside'] = outside
if not variance is None:
rules['variance'] = variance
mech_content = {param_name: rules}
if not sec_type in self.mech_dict: # No mechanisms have been inserted into this type of section yet
self.mech_dict[sec_type] = {'synapse': {syn_type: mech_content}}
elif not 'synapse' in self.mech_dict[sec_type]: # No synaptic mechanisms have been specified in this type of
backup_content = copy.deepcopy(self.mech_dict[sec_type]) # section yet
self.mech_dict[sec_type]['synapse'] = {syn_type: mech_content}
elif not syn_type in self.mech_dict[sec_type]['synapse']: # This synaptic mechanism has not yet been inserted
backup_content = copy.deepcopy(self.mech_dict[sec_type]) # into this type of section
self.mech_dict[sec_type]['synapse'][syn_type] = mech_content
elif self.mech_dict[sec_type]['synapse'][syn_type] is None: # This mechanism has been inserted, but no
backup_content = copy.deepcopy(self.mech_dict[sec_type]) # parameters have been specified
self.mech_dict[sec_type]['synapse'][syn_type] = mech_content
elif param_name in self.mech_dict[sec_type]['synapse'][syn_type]: # This parameter has already been specified.
backup_content = copy.deepcopy(self.mech_dict[sec_type]) # Now have to determine whether to replace or
if replace: # extend the current dictionary entry.
self.mech_dict[sec_type]['synapse'][syn_type][param_name] = rules
elif type(self.mech_dict[sec_type]['synapse'][syn_type][param_name]) == dict:
self.mech_dict[sec_type]['synapse'][syn_type][param_name] = \
[self.mech_dict[sec_type]['synapse'][syn_type][param_name], rules]
elif type(self.mech_dict[sec_type]['synapse'][syn_type][param_name]) == list:
self.mech_dict[sec_type]['synapse'][syn_type][param_name].append(rules)
elif not param_name is None: # This mechanism has been inserted, but this
backup_content = copy.deepcopy(self.mech_dict[sec_type]) # parameter has not yet been specified
self.mech_dict[sec_type]['synapse'][syn_type][param_name] = rules
try:
nodes = self.get_nodes_of_subtype(sec_type)
for node in nodes:
try:
self._modify_mechanism(node, 'synapse', {syn_type: mech_content})
except (AttributeError, NameError, ValueError, KeyError):
if not param_name is None:
print 'Exception: Problem modifying synaptic mechanism: {} parameter: {}'.format(syn_type,
param_name)
else:
print 'Exception: Problem inserting synaptic mechanism: {}'.format(syn_type)
raise KeyError
except KeyError:
if backup_content is None:
del self.mech_dict[sec_type]
else:
self.mech_dict[sec_type] = copy.deepcopy(backup_content)
finally:
if verbose:
pprint.pprint(self.mech_dict)
def export_mech_dict(self, mech_filename=None):
"""
Following modifications to the mechanism dictionary either during model specification or parameter optimization,
this method stores the current mech_dict to a pickle file stamped with the date and time. This allows the
current set of mechanism parameters to be recalled later.
"""
if mech_filename is None:
mech_filename = 'mech_dict_'+datetime.datetime.today().strftime('%m%d%Y%H%M')+'.pkl'
write_to_pkl(data_dir+mech_filename+'.pkl', self.mech_dict)
print "Exported mechanism dictionary to "+mech_filename+'.pkl'
def get_node_by_distance_to_soma(self, distance, sec_type):
"""
Gets the first node of the given section type at least the given distance from a soma node.
Not particularly useful, since it will always return the same node.
:param distance: int or float
:param sec_type: str
:return: :class:'SHocNode'
"""
nodes = self._node_dict[sec_type]
for node in nodes:
if self.get_distance_to_node(self.tree.root, node) >= distance:
return node
raise Exception('No node is {} um from a soma node.'.format(distance))
def get_distance_to_node(self, root, node, loc=None):
"""
Returns the distance from the given location on the given node to its connection with a root node.
:param root: :class:'SHocNode'
:param node: :class:'SHocNode'
:param loc: float
:return: int or float
"""
length = 0.
if node in self.soma:
return length
if not loc is None:
length += loc*node.sec.L
if root in self.soma:
while not node.parent in self.soma:
node.sec.push()
loc = h.parent_connection()
h.pop_section()
node = node.parent
length += loc*node.sec.L
elif self.node_in_subtree(root, node):
while not node.parent is root:
node.sec.push()
loc = h.parent_connection()
h.pop_section()
node = node.parent
length += loc*node.sec.L
else:
return None # node is not connected to root
return length
def node_in_subtree(self, root, node):
"""
Checks if a node is contained within a subtree of root.
:param root: 'class':SNode2 or SHocNode
:param node: 'class':SNode2 or SHocNode
:return: boolean
"""
nodelist = []
self.tree._gather_nodes(root, nodelist)
if node in nodelist:
return True
else:
return False
def get_path_length_swc(self, path):
"""
Calculates the distance between nodes given a list of SNode2 nodes connected in a path.
:param path: list : :class:'SNode2'
:return: int or float
"""
distance = 0
for i in range(len(path)-1):
distance += np.sqrt(np.sum((path[i].get_content()['p3d'].xyz - path[i+1].get_content()['p3d'].xyz)**2))
return distance
def get_node_length_swc(self, node):
"""
Calculates the distance between the center points of an SNode2 node and its parent.
:param node: :class:'SNode2'
:return: float
"""
if not self.tree.is_root(node):
return np.sqrt(np.sum((node.get_content()['p3d'].xyz - node.parent.get_content()['p3d'].xyz)**2))
else:
return np.sqrt(np.sum(node.get_content()['p3d'].xyz**2))
def get_branch_order(self, node):
"""
Calculates the branch order of a SHocNode node. The order is defined as 0 for all soma, axon, and apical trunk
dendrite nodes, but defined as 1 for basal dendrites that branch from the soma, and apical and tuft dendrites
that branch from the trunk. Increases by 1 after each additional branch point. Makes sure not to count spines.
:param node: :class:'SHocNode'
:return: int
"""
if node.type in ['soma', 'axon_hill', 'ais', 'axon']:
return 0
elif node.type == 'trunk':
children = [child for child in node.parent.children if not child.type == 'spine_neck']
if len(children) > 1 and children[0].type == 'trunk' and children[1].type == 'trunk':
return 1
else:
return 0
else:
order = 0
path = [branch for branch in self.tree.path_between_nodes(node, self.get_dendrite_origin(node)) if
not branch.type in ['soma', 'trunk']]
for node in path:
if self.is_terminal(node):
order += 1
elif len([child for child in node.parent.children if not child.type == 'spine_neck']) > 1:
order += 1
elif node.parent.type == 'trunk':
order += 1
return order
def is_terminal(self, node):
"""
Calculates if a node is a terminal dendritic branch.
:param node: :class:'SHocNode'
:return: bool
"""
if node.type in ['soma', 'axon_hill', 'ais', 'axon']:
return False
else:
return not bool([child for child in node.children if not child.type == 'spine_neck'])
def is_bifurcation(self, node, child_type):
"""
Calculates if a node bifurcates into at least two children of specified type.
:param node: :class:'SHocNode'
:param child_type: string
:return: bool
"""
return len([child for child in node.children if child.type == child_type]) >= 2
def set_stochastic_synapses(self, value):
"""
This method turns stochastic filtering of presynaptic release on or off for all synapses contained in this cell.
:param value: int in [0, 1]
"""
for nodelist in self._node_dict.itervalues():
for node in nodelist:
for syn in node.synapses:
syn.stochastic = value
@property
def gid(self):
return self._gid
@property
def soma(self):
return self._node_dict['soma']
@property
def axon(self):
return self._node_dict['axon']
@property
def basal(self):
return self._node_dict['basal']
@property
def apical(self):
return self._node_dict['apical']
@property
def trunk(self):
return self._node_dict['trunk']
@property
def tuft(self):
return self._node_dict['tuft']
@property
def spine(self):
return self._node_dict['spine']
#------------------------------Extend SNode2 to interact with NEURON hoc sections------------------------
class SHocNode(btmorph.btstructs2.SNode2):
"""
Extends SNode2 with some methods for storing and retrieving additional information in the node's content
dictionary related to running NEURON models specified in the hoc language.
"""
def __init__(self, index=0):
"""
:param index: int : unique node identifier
"""
btmorph.btstructs2.SNode2.__init__(self, index)
self.content['spines'] = []
self.content['synapses'] = []
def get_sec(self):
"""
Returns the hoc section associated with this node, stored in the node's content dictionary.
:return: :class:'neuron.h.Section'
"""
if 'sec' in self.content:
return self.content['sec']
else:
raise Exception('This node does not yet have an associated hoc section.')
def set_sec(self, sec):
"""
Stores the hoc section associated with this node in the node's content dictionary.
:param sec: :class:'neuron.h.Section'
"""
self.content['sec'] = sec
sec = property(get_sec, set_sec)
def init_nseg(self, spatial_res=0):
"""
Initializes the number of hoc segments in this node's hoc section (nseg) based on the AC length constant.
Must be re-initialized whenever basic cable properties Ra or cm are changed. If the node is a tapered cylinder,
it should contain at least 3 segments. The spatial resolution parameter increases the number of segments per
section by a factor of an exponent of 3.
If a section's nseg has been manually increased beyond the suggestion of the mechanism dictionary, this method
does not decrease it.
:param spatial_res: int
"""
sugg_nseg = d_lambda_nseg(self.sec)
if not self.get_diam_bounds() is None:
sugg_nseg = max(sugg_nseg, 3)
sugg_nseg *= 3**spatial_res
if self.sec.nseg < sugg_nseg:
self.sec.nseg = sugg_nseg
def reinit_diam(self):
"""
For a node associated with a hoc section that is a tapered cylinder, every time the spatial resolution
of the section (nseg) is changed, the section diameters must be reinitialized. This method checks the
node's content dictionary for diameter boundaries and recalibrates the hoc section associated with this node.
"""
if not self.get_diam_bounds() is None:
[diam1, diam2] = self.get_diam_bounds()
h('diam(0:1)={}:{}'.format(diam1, diam2), sec=self.sec)
def get_diam_bounds(self):
"""
If the hoc section associated with this node is a tapered cylinder, this method returns a list containing
the values of the diameters at the 0 and 1 ends of the section, stored in the node's content dictionary.
Otherwise, it returns None (for non-conical cylinders).
:return: (list: int) or None
"""
if 'diam' in self.content:
return self.content['diam']
else:
return None
def set_diam_bounds(self, diam1, diam2):
"""
For a node associated with a hoc section that is a tapered cylinder, this stores a list containing the values
of the diameters at the 0 and 1 ends of the section in the node's content dictionary.
:param diam1: int
:param diam2: int
"""
self.content['diam'] = [diam1, diam2]
self.reinit_diam()
def get_type(self):
"""
NEURON sections are assigned a node type for convenience in order to later specify membrane mechanisms and
properties for each type of compartment.
:return: str
"""
if 'type' in self.content:
return self.content['type']
else:
raise Exception('This node does not yet have a defined type.')
def set_type(self, type):
"""
Checks that type is a string in the list of defined section types, and stores the value in the node's content
dictionary.
:param type: str
"""
if type in sec_types:
self.content['type'] = type
else:
raise Exception('That is not a defined type of section.')
type = property(get_type, set_type)
def connect(self, parent, pindex=1, cindex=0):
"""
Connects this SHocNode node to a parent node, and establishes a connection between their associated
hoc sections.
:param parent: :class:'SHocNode'
:param pindex: int in [0,1] Connect to this end of the parent hoc section.
:param cindex: int in [0,1] Connect this end of the child hoc section
"""
self.parent=parent
parent.add_child(self)
self.sec.connect(parent.sec, pindex, cindex)
@property
def name(self):
"""
Returns a str containing the name of the hoc section associated with this node. Consists of a type descriptor
and an index identifier.
:return: str
"""
if 'type' in self.content:
return '{0.type}{0.index}'.format(self)
else:
raise Exception('This node does not yet have a defined type.')
@property
def spines(self):
"""
Returns a list of the spine head sections attached to the hoc section associated with this node.
:return: list of :class:'SHocNode' of sec_type == 'spine_head'
"""
return self.content['spines']
@property
def synapses(self):
"""
Returns a list of the objects of :class:'Synapse' associated with this node.
:return: list of hoc objects, type depends on .mod file(s) used to implement synapses
"""
return self.content['synapses']
class CA1_Pyr(HocCell):
def __init__(self, morph_filename=None, mech_filename=None, full_spines=True):
HocCell.__init__(self, morph_filename, mech_filename)
self.random.seed(self.gid) # This cell will always have the same spine and GABA_A synapse locations as long as
# they are inserted in the same order
if full_spines:
self.insert_spines_in_subset(['basal', 'trunk', 'apical', 'tuft'])
def insert_spines_in_subset(self, sec_type_list):
"""
This method populates the cell tree with spines following spine density information from Erk Bloss &
Nelson Spruston. Basal dendrites have no spines until the first branch point, and a higher density beyond the
second branch point. Trunk dendrites have no spines until the first branch point, and an increasing density
until the tuft branch point(s). Apical dendrites have a density that varies with the distance from the soma of
their original branch point from the trunk. Terminal tuft branches have a higher density than their parents.
Should standardize the implementation of the rules for each type of dendrite and import the
density dictionary from a file, similar to the implementation of the membrane mechanism dictionary.
:param sec_type_list: list of str
"""
densities = {'trunk': {'min': 0.2418, 'max': 3.8,
'start': min([self.get_distance_to_node(self.tree.root, branch) for branch in
self.apical]),
'end': max([self.get_distance_to_node(self.tree.root, branch) for branch in
self.trunk])},
'basal': {'1': 0., '2': 0.4428, '>2': 1.891},
'apical': {'min': 2.273, 'max': 2.688,
'start': min([self.get_distance_to_node(self.tree.root, branch) for branch in
self.apical]),
'end': max([self.get_distance_to_node(self.tree.root, branch)
for branch in self.apical if self.get_branch_order(branch) == 1])},
'tuft': {'parent': 1.354, 'terminal': 0.7157}
}
if 'basal' in sec_type_list:
for node in self.basal:
order = self.get_branch_order(node)
if order == 2:
self.insert_spines_every(node, densities['basal']['2'])
elif order > 2:
self.insert_spines_every(node, densities['basal']['>2'])
if 'trunk' in sec_type_list:
for node in self.trunk:
distance = self.get_distance_to_node(self.tree.root, node)
if distance >= densities['trunk']['start']:
slope = (densities['trunk']['max'] - densities['trunk']['min']) / \
(densities['trunk']['end'] - densities['trunk']['start'])
density = densities['trunk']['min'] + slope * (distance - densities['trunk']['start'])
self.insert_spines_every(node, density)
if 'apical' in sec_type_list:
for node in self.apical:
distance = self.get_distance_to_node(self.tree.root, self.get_dendrite_origin(node), loc=1.)
slope = (densities['apical']['max'] - densities['apical']['min']) / \
(densities['apical']['end'] - densities['apical']['start'])
density = densities['apical']['min'] + slope * (distance - densities['apical']['start'])
self.insert_spines_every(node, density)
if 'tuft' in sec_type_list:
for node in self.tuft:
if self.is_terminal(node):
self.insert_spines_every(node, densities['tuft']['terminal'])
else:
self.insert_spines_every(node, densities['tuft']['parent'])
self._reinit_mech(self.spine)
def insert_spines_every(self, node, density):
"""
Given a mean spine density in /um, insert spines in the node at the specified density.
:param node: :class:'SHocNode'
:param density: float: mean density in /um
"""
L = node.sec.L
beta = 1./density
interval = self.random.exponential(beta)
while interval < L:
self.insert_spine(node, interval/L)
interval += self.random.exponential(beta)
def insert_spine(self, node, parent_loc, child_loc=0):
"""
Spines consist of two hoc sections: a cylindrical spine head and a cylindrical spine neck.
:param node: :class:'SHocNode'
:param parent_loc: float
:param child_loc: int
"""
neck = self.make_section('spine_neck')
neck.connect(node, parent_loc, child_loc)
neck.sec.L = 1.58
neck.sec.diam = 0.077
self._init_cable(neck)
head = self.make_section('spine_head')
head.connect(neck)
node.spines.append(head)
head.sec.L = 0.5 # open cylinder, matches surface area of sphere with diam = 0.5
head.sec.diam = 0.5
self._init_cable(head)
def insert_inhibitory_synapses_in_subset(self, sec_type_list=None):
"""
:param sec_type_list: str
"""
if sec_type_list is None:
sec_type_list = ['soma', 'ais', 'basal', 'trunk', 'apical', 'tuft']
densities = {'soma': 2.857, # 4.285,
'ais': 0.53,
'trunk': {'min': 0.3022, 'max': 0.0627,
'start': 0.,
'end': max([self.get_distance_to_node(self.tree.root, branch) for branch in
self.trunk])},
'basal': {'primary': 0.3129, 'intermediate': 0.1728, 'terminal': 0.06543},
'apical': {'min': 0.03885, 'max': 0.04512,
'start': min([self.get_distance_to_node(self.tree.root, branch) for branch in
self.apical]),
'end': max([self.get_distance_to_node(self.tree.root, branch)
for branch in self.apical if self.get_branch_order(branch) == 1])},
'tuft': {'parent': 0.2104, 'terminal': 0.1619}
}
if 'soma' in sec_type_list:
for node in self.soma:
self.insert_inhibitory_synapse_every(node, densities['soma'])
if 'ais' in sec_type_list:
for node in self.get_nodes_of_subtype('ais'):
self.insert_inhibitory_synapse_every(node, densities['ais'])
if 'basal' in sec_type_list:
for node in self.basal:
if self.is_terminal(node):
self.insert_inhibitory_synapse_every(node, densities['basal']['terminal'])
else:
order = self.get_branch_order(node)
if order == 1:
self.insert_inhibitory_synapse_every(node, densities['basal']['primary'])
else:
self.insert_inhibitory_synapse_every(node, densities['basal']['intermediate'])
if 'trunk' in sec_type_list:
for node in self.trunk:
distance = self.get_distance_to_node(self.tree.root, node)
if distance >= densities['trunk']['start']:
slope = (densities['trunk']['max'] - densities['trunk']['min']) / \
(densities['trunk']['end'] - densities['trunk']['start'])
density = densities['trunk']['min'] + slope * (distance - densities['trunk']['start'])
self.insert_inhibitory_synapse_every(node, density)
if 'apical' in sec_type_list:
for node in self.apical:
distance = self.get_distance_to_node(self.tree.root, self.get_dendrite_origin(node), loc=1.)
slope = (densities['apical']['max'] - densities['apical']['min']) / \
(densities['apical']['end'] - densities['apical']['start'])
density = densities['apical']['min'] + slope * (distance - densities['apical']['start'])
self.insert_inhibitory_synapse_every(node, density)
if 'tuft' in sec_type_list:
for node in self.tuft:
if self.is_terminal(node):
self.insert_inhibitory_synapse_every(node, densities['tuft']['terminal'])
else:
self.insert_inhibitory_synapse_every(node, densities['tuft']['parent'])
def insert_inhibitory_synapse_every(self, node, density, syn_types=['GABA_A_KIN'], stochastic=0):
"""
:param node: :class:'SHocNode'
:param density: float: mean density in /um
:param syn_types: list of str
:param stochastic: int
"""
L = node.sec.L
beta = 1./density
interval = self.random.exponential(beta)
while interval < L:
syn = Synapse(self, node, type_list=syn_types, stochastic=stochastic, loc=interval/L)
interval += self.random.exponential(beta)
def zero_na(self):
"""
Set na channel conductances to zero in all compartments. Used during parameter optimization.
"""
for sec_type in ['soma', 'axon_hill', 'ais', 'axon', 'basal', 'trunk', 'apical', 'tuft']:
for na_type in (na_type for na_type in ['nas_kin', 'nat_kin', 'nas', 'nax'] if na_type in
self.mech_dict[sec_type]):
self.modify_mech_param(sec_type, na_type, 'gbar', 0.)
def zero_h(self):
"""
Set ih conductances to zero in all compartments. Used during parameter optimization.
"""
self.modify_mech_param('soma', 'h', 'ghbar', 0.)
self.mech_dict['trunk']['h']['ghbar']['value'] = 0.
self.mech_dict['trunk']['h']['ghbar']['slope'] = 0.
for sec_type in ['basal', 'trunk', 'apical', 'tuft']:
self.reinitialize_subset_mechanisms(sec_type, 'h')
def set_terminal_branch_na_gradient(self, gmin=0.):
"""
This is an admittedly ad-hoc procedure to implement a linearly decreasing gradient of sodium channels in
terminal branches that is not easily accomplished by the general procedures implementing the mechanism
dictionary.
"""
na_type = (na_type for na_type in ['nas_kin', 'nat_kin', 'nas', 'nax']
if na_type in self.mech_dict['trunk']).next()
self.set_special_mech_param_linear_gradient(na_type, 'gbar', ['basal', 'apical', 'tuft'],
self.is_terminal, gmin)
class Synapse(object):
"""
The implementation in hoc of synaptic mechanisms that can be triggered is complicated. This container is an attempt
to wrap all the objects required to deliver synaptic events to a section, and have separable synaptic mechanisms
(e.g. GluA-Rs and GluN-Rs) respond with individually specifiable weights and kinetics.
To make model specification and simulation implementation straightforward, synapses are not meant to be moved once
they are initialized.
"""
def __init__(self, cell, node, type_list=None, stochastic=1, loc=0.5, delay=0., weight=1., threshold=-30.,
stochastic_type=None, source=None, source_node=None, source_param=None, source_loc=0.5):
"""
Design goals: A source (like a spike detector in a presynaptic neuron) can be specified. If not, a VecStim
object is used a source, which can be played events at specified times using its .play method. If stochastic,
all spikes are intercepted by a point process with release probability dynamics and its own unique and
independent random variable from a uniform distribution. If not, the specified synaptic mechanisms are connected
directly to the source of spikes.
:param cell: :class:'HocCell'
:param node: :class:'SHoCNode'
:param type_list: list of str
:param stochastic: int in [0, 1]
:param loc: float
:param delay: float
:param weight: float
:param threshold: float
:param stochastic_type: str
:param source: hoc artificial cell or otherwise hoc object not associated with a section
:param source_node: :class:'SHocNode' specifies the section containing a range variable to be used as a source
:param source_param: str
:param source_loc: str
"""
self._cell = cell
self._node = node
self._stochastic = stochastic
self._loc = loc
self._delay = delay
self._weight = weight
self._threshold = threshold
self._syn = {}
self.randObj = None
self._node.synapses.append(self)
if source_param is None:
source_param = '_ref_v'
if stochastic_type is None:
self._stochastic_type = 'Pr'
else:
self._stochastic_type = stochastic_type
if not source is None:
self._source = {'object': source}
elif not source_node is None:
self._source = {'object': getattr(source_node.sec(source_loc), source_param), 'node': source_node}
else:
self._source = {'object': h.VecStim()}
if type_list is None:
type_list = ['AMPA_KIN']
elif type(type_list) is not list:
type_list = [type_list]
if self.stochastic:
self._init_stochastic()
for target in type_list:
syn = getattr(h, target)(self.node.sec(self._loc))
self._syn[target] = {'target': syn}
if self.stochastic:
self._syn[target]['netcon'] = h.NetCon(self.target(self._stochastic_type), syn)
self.netcon(target).delay = delay
self.netcon(target).weight[0] = weight
self.netcon(target).threshold = threshold
else:
self._init_netcon(target)
def _init_stochastic(self):
"""
This method constructs and initializes a stochastic filtering mechanism that intercepts spikes delivered to this
synapse and calculates whether or not to pass a spike to the rest of the specified synaptic mechanisms.
"""
if self.randObj is None: # if this synapse has never been stochastic, it needs a new random number generator
self.randObj = h.Random()
self.randObj.MCellRan4(self.cell.gid*1e4+1, self.node.index*1e4+self.node.synapses.index(self)+1)
# a unique sequence for up to ~10,000 spikes per synapse; ~10,000 synapses per node;
# ~4,290,000 nodes per cell; ~4,290,000 cell in a network
self.randObj.uniform(0, 1)
else: # if this synapse has already been stochastic before, this restarts its random number generator
self.randObj.seq(self.cell.gid*1e4+1)
syn = getattr(h, self._stochastic_type)(self.node.sec(self._loc))
self._syn[self._stochastic_type] = {'target': syn}
self._init_netcon(self._stochastic_type, delay=0.)
self.target(self._stochastic_type).setRandObjRef(self.randObj)
def target(self, target):
"""
Returns the hoc object for the synaptic mechanism of the specified type
:param target: str
:return: :class:'h.HocObject'
"""
if target in self._syn:
return self._syn[target]['target']
else:
raise KeyError('Synapse type: {} not found at a synapse in {}'.format(target, self._node.name))
def _init_netcon(self, target, delay=None, weight=None, threshold=None):
"""
Appropriately initializes new netcon object, depending on whether the current source dictionary specifies a
hocObject without a section, or a reference variable contained within a section.
:param target: str
:param delay = float
:param weight = float
:param threshold = float
"""
if target in self._syn:
source = self._source['object']
syn = self._syn[target]
if weight is None:
weight = self.weight
if threshold is None:
threshold = self.threshold
if delay is None:
delay = self.delay
if 'node' in self._source:
node = self._source['node']
syn['netcon'] = h.NetCon(source, syn['target'], sec=node.sec)
else:
syn['netcon'] = h.NetCon(source, syn['target'])
syn['netcon'].delay = delay
syn['netcon'].weight[0] = weight
syn['netcon'].threshold = threshold
else:
raise KeyError('Synapse type: {} not found at a synapse in {}'.format(target, self._node.name))
def netcon(self, target):
"""
Returns the hoc network connection linking the synaptic mechanism of the specified type to a source of spikes.
:param target: str
:return: :class:'h.NetCon'
"""
if target in self._syn:
return self._syn[target]['netcon']
else:
raise KeyError('Synapse type: {} not found at a synapse in {}'.format(target, self._node.name))
def change_source(self, source=None, node=None, param=None, loc=0.5):
"""
In order to change the source of a synapse from the default VecStim object to a membrane potential or artificial
cell, all the netcon objects must be deleted and replaced with new ones. Preserves previously set values for
delay, weight, and threshold for all synaptic mechanisms associated with this synapse.
:param source: hoc artificial cell or otherwise hoc object not associated with a section
:param node: :class:'SHocNode'
:param param: str corresponding to range variable in section
:param loc: float
"""
netcon_dict = {}
if param is None:
param = '_ref_v'
for target in self._syn:
netcon_dict[target] = {}
netcon_dict[target]['delay'] = self.netcon(target).delay
netcon_dict[target]['weight'] = self.netcon(target).weight[0]
netcon_dict[target]['threshold'] = self.netcon(target).threshold
if source is None:
if node is None:
raise Exception('A source or reference node must be provided to establish a new synaptic connection.')
else:
del self._source
self._source = {'object': getattr(node.sec(loc), param), 'node': node}
else:
del self._source
self._source = {'object': source}
if self._stochastic:
del self._syn[self._stochastic_type]['netcon']
delay = netcon_dict[self._stochastic_type]['delay']
weight = netcon_dict[self._stochastic_type]['weight']
threshold = netcon_dict[self._stochastic_type]['threshold']
self._init_netcon(self._stochastic_type, delay=delay, weight=weight, threshold=threshold)
else:
for target in (target for target in self._syn if not target == self._stochastic_type):
del self._syn[target]['netcon']
delay = netcon_dict[target]['delay']
weight = netcon_dict[target]['weight']
threshold = netcon_dict[target]['threshold']
self._init_netcon(target, delay=delay, weight=weight, threshold=threshold)
def get_stochastic(self):
"""
Returns the value of an internal variable indicating if this synapse has a stochastic filter for spikes.
:return: int in [0, 1]
"""
return self._stochastic
def set_stochastic(self, value):
"""
Turns on or off stochastic filtering of spikes, preserving delay, weight, and threshold for all synaptic
mechanisms associated with this synapse.
:param value: int in [0, 1]
"""
if not (value == self._stochastic):
self._stochastic = value
if value:
self._init_stochastic()
for target in (target for target in self._syn if not target == self._stochastic_type):
delay = self.netcon(target).delay
weight = self.netcon(target).weight[0]
threshold = self.netcon(target).threshold
del self._syn[target]['netcon']
self._syn[target]['netcon'] = h.NetCon(self.target(self._stochastic_type), self.target(target))
self.netcon(target).delay = delay
self.netcon(target).weight[0] = weight
self.netcon(target).threshold = threshold
else:
for target in (target for target in self._syn if not target == self._stochastic_type):
delay = self.netcon(target).delay
weight = self.netcon(target).weight[0]
threshold = self.netcon(target).threshold
del self._syn[target]['netcon']
self._init_netcon(target, delay=delay, weight=weight, threshold=threshold)
del self._syn[self._stochastic_type]
stochastic = property(get_stochastic, set_stochastic)
def get_delay(self):
"""
Returns the default value of the time delay (ms) between spike and activation for this synapse.
:return: int or float
"""
return self._delay
def set_delay(self, value):
"""
Changes the value of the time delay (ms) between spike and activation for all synaptic mechanisms associated
with this synapse, except self._stochastic_type, which retains its current value until set manually.
:param value: int or float
"""
self._delay = value
for target in (target for target in self._syn if not target == self._stochastic_type):
self.netcon(target).delay = value
delay = property(get_delay, set_delay)
def get_weight(self):
"""
Returns the default value of the activation weight for this synapse.
:return: float
"""
return self._weight
def set_weight(self, value):
"""
Changes the value of the activation weight for all synaptic mechanisms associated with this synapse.
:param value: float
"""
self._weight = value
for target in (target for target in self._syn):
self.netcon(target).weight[0] = value
weight = property(get_weight, set_weight)
def get_threshold(self):
"""
Returns the value of the activation threshold for this synapse.
:return: float
"""
return self._threshold
def set_threshold(self, value):
"""
Changes the value of the activation threshold for this synapse.
:param value: float
"""
self._threshold = value
for target in (target for target in self._syn):
self.netcon(target).threshold = value
threshold = property(get_threshold, set_threshold)
@property
def source(self):
"""
Returns the hocObject currently being used as a source.
:return: :class:'hocObject'
"""
return self._source['object']
@property
def cell(self):
"""
Returns the cell containing this synapse.
:return: :class:'HocCell'
"""
return self._cell
@property
def node(self):
"""
Returns the node containing this synapse.
:return: :class:'SHocNode'
"""
return self._node
@property
def loc(self):
"""
Returns the location along the hoc section containing this synapse. For convenience, if the synapse is
contained in a spine_head, this property method returns the location along the branch section where the
spine_neck is connected.
:return: int or float
"""
if self.node.type == 'spine_head':
self.node.parent.sec.push()
loc = h.parent_connection()
h.pop_section()
return loc
else:
return self._loc