diff --git a/python/moose/neuroml2/reader.py b/python/moose/neuroml2/reader.py index b68d4d473f..44aa1948f2 100644 --- a/python/moose/neuroml2/reader.py +++ b/python/moose/neuroml2/reader.py @@ -5,7 +5,7 @@ # Implementation of reader for NeuroML 2 models. # TODO: handle morphologies of more than one segment... # Author: Subhasis Ray, Padraig Gleeson -# Maintainer: Dilawar Singh , Padraig Gleeson +# Maintainer: Padraig Gleeson, Subhasis Ray # Created: Wed Jul 24 15:55:54 2013 (+0530) # Notes: # For update/log, please see git-blame documentation or browse the github @@ -15,25 +15,29 @@ import math import logging import numpy as np +from collections import defaultdict + import moose -logger_ = logging.getLogger('moose.nml2') +logger_ = logging.getLogger("moose.nml2") -nml_not_available_msg = '' +nml_not_available_msg = "" try: import neuroml as nml import pyneuroml.pynml as pynml except ImportError as error: - raise ImportError(f'Could not import neuroml/pyneuroml. Please make sure you have pyneuroml installed (`pip install pyneuroml`)') from error - - + raise ImportError( + f"Could not import neuroml/pyneuroml. Please make sure you have pyneuroml installed (`pip install pyneuroml`)" + ) from error + + from moose.neuroml2.units import SI -def _write_flattened_nml( doc, outfile ): +def _write_flattened_nml(doc, outfile): """_write_flattened_nml Concat all NML2 read by moose and generate one flattened NML file. Only useful when debugging. @@ -42,10 +46,12 @@ def _write_flattened_nml( doc, outfile ): :param outfile: Name of the output file. """ import neuroml.writers - neuroml.writers.NeuroMLWriter.write( doc, outfile ) - logger_.debug( "Wrote flattened NML model to %s" % outfile ) -def _gates_sorted( all_gates ): + neuroml.writers.NeuroMLWriter.write(doc, outfile) + logger_.debug("Wrote flattened NML model to %s" % outfile) + + +def _gates_sorted(all_gates): """_gates_sorted Parameters @@ -59,23 +65,25 @@ def _gates_sorted( all_gates ): X, Y or Z gate respectively. Otherwise do not touch them i.e. first gate will be loaded into X, second into Y and so on. """ - allMooseGates = ['x', 'y', 'z'] - allGatesDict = { g.id : g for g in all_gates } - gateNames = [ g.id.lower() for g in all_gates ] + allMooseGates = ["x", "y", "z"] + allGatesDict = {g.id: g for g in all_gates} + gateNames = [g.id.lower() for g in all_gates] if set(gateNames).issubset(set(allMooseGates)): sortedGates = [] for gid in allMooseGates: - sortedGates.append( allGatesDict.get(gid) ) + sortedGates.append(allGatesDict.get(gid)) return sortedGates return all_gates -def _unique( ls ): - res = [ ] + +def _unique(ls): + res = [] for l in ls: if l not in res: - res.append( l ) + res.append(l) return res + def _isConcDep(ct): """_isConcDep Check if componet is dependent on concentration. Most HHGates are @@ -87,23 +95,27 @@ def _isConcDep(ct): :return: True if Component is depenant on conc, False otherwise. """ # logger_.debug(f"{'#' * 10} EXTENDS {ct.extends}") - if 'ConcDep' in ct.extends: + if "ConcDep" in ct.extends: return True return False + def _findCaConcVariableName(): """_findCaConcVariableName Find a suitable CaConc for computing HHGate tables. This is a hack, though it is likely to work in most cases. """ - caConcs = moose.wildcardFind( '/library/##[TYPE=CaConc]' ) - assert len(caConcs) >= 1, "No moose.CaConc found. Currently moose \ + caConcs = moose.wildcardFind("/library/##[TYPE=CaConc]") + assert ( + len(caConcs) >= 1 + ), "No moose.CaConc found. Currently moose \ supports HHChannel which depends only on moose.CaConc ." return caConcs[0].name + def sarea(comp): """sarea - Return the surface area (2ϖrL) of compartment from length and diameter. + Return the surface area (2 * pi * r * L) of compartment from length and diameter. :param comp: Compartment instance. :type comp: str @@ -115,9 +127,10 @@ def sarea(comp): else: return math.pi * comp.diameter * comp.diameter + def xarea(compt): """xarea - Return the cross sectional area (𝜋r²) from the diameter of the compartment. + Return the cross sectional area (pi * r^2) from the diameter of the compartment. Note: ---- @@ -128,7 +141,8 @@ def xarea(compt): :return: cross sectional area. :rtype: float """ - return math.pi * (compt.diameter/2.0)**2.0 + return math.pi * (compt.diameter / 2.0) ** 2.0 + def setRa(comp, resistivity): """Calculate total raxial from specific value `resistivity`""" @@ -137,20 +151,23 @@ def setRa(comp, resistivity): else: comp.Ra = resistivity * 4.0 / (comp.diameter * np.pi) + def setRm(comp, condDensity): """Set membrane resistance""" - comp.Rm = 1/(condDensity * sarea(comp)) + comp.Rm = 1 / (condDensity * sarea(comp)) + def setEk(comp, erev): """Set reversal potential""" comp.setEm(erev) + def getSegments(nmlcell, component, sg_to_segments): """Get the list of segments the `component` is applied to""" sg = component.segment_groups if sg is None: segments = nmlcell.morphology.segments - elif sg == 'all': + elif sg == "all": segments = [seg for seglist in sg_to_segments.values() for seg in seglist] else: segments = sg_to_segments[sg] @@ -171,27 +188,33 @@ class NML2Reader(object): creates a passive neuronal morphology `/library/Purk2M9s`. """ + def __init__(self, verbose=False): global logger_ - self.lunit = 1e-6 # micron is the default length unit + self.lunit = 1e-6 # micron is the default length unit self.verbose = verbose if self.verbose: - logger_.setLevel( logging.DEBUG ) + logger_.setLevel(logging.DEBUG) self.doc = None self.filename = None - self.nml_cells_to_moose = {} # NeuroML object to MOOSE object - self.nml_segs_to_moose = {} # NeuroML object to MOOSE object - self.nml_chans_to_moose = {} # NeuroML object to MOOSE object - self.nml_conc_to_moose = {} # NeuroML object to MOOSE object - self.moose_to_nml = {} # Moose object to NeuroML object + self.nml_cells_to_moose = {} # NeuroML object to MOOSE object + self.nml_segs_to_moose = {} # NeuroML object to MOOSE object + self.nml_chans_to_moose = {} # NeuroML object to MOOSE object + self.nml_conc_to_moose = {} # NeuroML object to MOOSE object + self.moose_to_nml = {} # Moose object to NeuroML object self.proto_cells = {} # map id to prototype cell in moose self.proto_chans = {} # map id to prototype channels in moose self.proto_pools = {} # map id to prototype pools (Ca2+, Mg2+) - self.includes = {} # Included files mapped to other readers - self.lib = moose.element('/library') if moose.exists('/library') \ - else moose.Neutral( '/library' ) + self.includes = {} # Included files mapped to other readers + self.lib = ( + moose.element("/library") + if moose.exists("/library") + else moose.Neutral("/library") + ) self.id_to_ionChannel = {} - self._cell_to_sg = {} # nml cell to dict - the dict maps segment groups to segments + self._cell_to_sg = ( + {} + ) # nml cell to dict - the dict maps segment groups to segments self._variables = {} self.cells_in_populations = {} @@ -201,17 +224,18 @@ def __init__(self, verbose=False): self.network = None def read(self, filename, symmetric=True): - filename = os.path.realpath( filename ) + filename = os.path.realpath(filename) self.doc = nml.loaders.read_neuroml2_file( - filename, include_includes=True, verbose=self.verbose) + filename, include_includes=True, verbose=self.verbose + ) self.filename = filename - logger_.info('Parsed the NeuroML2 file: %s'% filename) + logger_.info("Parsed the NeuroML2 file: %s" % filename) if self.verbose: - _write_flattened_nml( self.doc, '%s__flattened.xml' % self.filename ) + _write_flattened_nml(self.doc, "%s__flattened.xml" % self.filename) - if len(self.doc.networks)>=1: + if len(self.doc.networks) >= 1: self.network = self.doc.networks[0] moose.celsius = self._getTemperature() @@ -223,64 +247,75 @@ def read(self, filename, symmetric=True): # logger_.debug(f"{'%' * 10} Creating cell prototype {cell}") self.createCellPrototype(cell, symmetric=symmetric) - if len(self.doc.networks)>=1: + if len(self.doc.networks) >= 1: self.createPopulations() self.createInputs() def _getTemperature(self): if self.network is not None: - if self.network.type=="networkWithTemperature": + if self.network.type == "networkWithTemperature": return SI(self.network.temperature) else: # Why not, if there's no temp dependence in nml..? return 0 - return SI('25') + return SI("25") def getCellInPopulation(self, pop_id, index): return self.cells_in_populations[pop_id][index] def getComp(self, pop_id, cellIndex, segId): - return moose.element('%s/%s/%s/%s' % ( self.lib.path, pop_id, cellIndex - , self.seg_id_to_comp_name[self.pop_to_cell_type[pop_id]][segId]) + return moose.element( + "%s/%s/%s/%s" + % ( + self.lib.path, + pop_id, + cellIndex, + self.seg_id_to_comp_name[self.pop_to_cell_type[pop_id]][segId], ) + ) def createPopulations(self): for pop in self.network.populations: - ep = '%s/%s' % (self.lib.path, pop.id) + ep = "%s/%s" % (self.lib.path, pop.id) mpop = moose.element(ep) if moose.exists(ep) else moose.Neutral(ep) - self.cells_in_populations[pop.id] ={} + self.cells_in_populations[pop.id] = {} for i in range(pop.size): - logger_.info("Creating %s/%s instances of cell: %s under %s"%(i,pop.size,pop.component, mpop)) - self.pop_to_cell_type[pop.id]=pop.component - chid = moose.copy(self.proto_cells[pop.component], mpop, '%s'%(i)) - self.cells_in_populations[pop.id][i]=chid - + logger_.info( + "Creating %s/%s instances of cell: %s under %s" + % (i, pop.size, pop.component, mpop) + ) + self.pop_to_cell_type[pop.id] = pop.component + chid = moose.copy(self.proto_cells[pop.component], mpop, "%s" % (i)) + self.cells_in_populations[pop.id][i] = chid def getInput(self, input_id): - return moose.element('%s/inputs/%s'%(self.lib.path,input_id)) - + return moose.element("%s/inputs/%s" % (self.lib.path, input_id)) def createInputs(self): for el in self.network.explicit_inputs: - pop_id = el.target.split('[')[0] - i = el.target.split('[')[1].split(']')[0] + pop_id = el.target.split("[")[0] + i = el.target.split("[")[1].split("]")[0] seg_id = 0 - if '/' in el.target: - seg_id = el.target.split('/')[1] + if "/" in el.target: + seg_id = el.target.split("/")[1] input = self.getInput(el.input) - moose.connect(input, 'output', self.getComp(pop_id,i,seg_id), 'injectMsg') + moose.connect(input, "output", self.getComp(pop_id, i, seg_id), "injectMsg") for il in self.network.input_lists: for ii in il.input: input = self.getInput(il.component) - moose.connect(input, 'output' - , self.getComp(il.populations,ii.get_target_cell_id(),ii.get_segment_id()) - , 'injectMsg') - + moose.connect( + input, + "output", + self.getComp( + il.populations, ii.get_target_cell_id(), ii.get_segment_id() + ), + "injectMsg", + ) def createCellPrototype(self, cell, symmetric=True): """To be completed - create the morphology, channels in prototype""" - ep = '%s/%s' % (self.lib.path, cell.id) + ep = "%s/%s" % (self.lib.path, cell.id) nrn = moose.element(ep) if moose.exists(ep) else moose.Neuron(ep) self.proto_cells[cell.id] = nrn self.nml_cells_to_moose[cell.id] = nrn @@ -289,7 +324,6 @@ def createCellPrototype(self, cell, symmetric=True): self.importBiophysics(cell, nrn) return cell, nrn - def createMorphology(self, nmlcell, moosecell, symmetric=True): """Create the MOOSE compartmental morphology in `moosecell` using the segments in NeuroML2 cell `nmlcell`. Create symmetric @@ -298,31 +332,26 @@ def createMorphology(self, nmlcell, moosecell, symmetric=True): """ morphology = nmlcell.morphology segments = morphology.segments - id_to_segment = dict([(seg.id, seg) for seg in segments]) + id_to_segment = {seg.id: seg for seg in segments} if symmetric: compclass = moose.SymCompartment + src, dst = "proximal", "distal" else: compclass = moose.Compartment + src, dst = "axial", "raxial" # segment names are used as compartment names - assuming # naming convention does not clash with that in MOOSE cellpath = moosecell.path id_to_comp = {} - self.seg_id_to_comp_name[nmlcell.id]={} + self.seg_id_to_comp_name[nmlcell.id] = {} + # create or access compartments by segment name (or id if name + # is not set) and keep two-way mapping for seg in segments: - if seg.name is not None: - ep = '%s/%s' % (cellpath, seg.name) - id_to_comp[seg.id] = moose.element(ep) if moose.exists(ep) else compclass(ep) - self.seg_id_to_comp_name[nmlcell.id][seg.id] = seg.name - else: - name = 'comp_%s'%seg.id - ep = '%s/%s' % (cellpath, name) - id_to_comp[seg.id] = moose.element(ep) if moose.exists(ep) else compclass(ep) - self.seg_id_to_comp_name[nmlcell.id][seg.id] = name + name = seg.name if seg.name is not None else f"comp_{seg.id}" + comp = compclass(f"{cellpath}/{name}") + id_to_comp[seg.id] = comp + self.seg_id_to_comp_name[nmlcell.id][seg.id] = name # Now assign the positions and connect up axial resistance - if not symmetric: - src, dst = 'axial', 'raxial' - else: - src, dst = 'proximal', 'distal' for segid, comp in id_to_comp.items(): segment = id_to_segment[segid] try: @@ -337,33 +366,37 @@ def createMorphology(self, nmlcell, moosecell, symmetric=True): p0 = parent.distal else: raise Exception( - 'No proximal point and no parent segment for segment:'+\ - 'name=%s, id=%s' % (segment.name, segment.id) - ) - comp.x0, comp.y0, comp.z0 = (x*self.lunit for x in map(float, (p0.x, p0.y, p0.z))) + "No proximal point and no parent segment for segment: " + f"name={segment.name}, id={segment.id}" + ) + comp.x0, comp.y0, comp.z0 = ( + float(x) * self.lunit for x in (p0.x, p0.y, p0.z) + ) p1 = segment.distal - comp.x, comp.y, comp.z = (x*self.lunit for x in map(float, (p1.x, p1.y, p1.z))) - comp.length = np.sqrt((comp.x-comp.x0)**2+(comp.y-comp.y0)**2+(comp.z-comp.z0)**2) + comp.x, comp.y, comp.z = (float(x) * self.lunit for x in (p1.x, p1.y, p1.z)) + comp.length = np.sqrt( + (comp.x - comp.x0) ** 2 + + (comp.y - comp.y0) ** 2 + + (comp.z - comp.z0) ** 2 + ) - # This can pose problem with moose where both ends of - # compartment have same diameter. We are averaging the two - # - may be splitting the compartment into two is better? - comp.diameter = (float(p0.diameter)+float(p1.diameter)) * self.lunit / 2 + # NOTE: moose compartments are cylindrical (both ends of a + # compartment have same diameter). So taking the average + # of the two ends in case of truncated-cone. + comp.diameter = (float(p0.diameter) + float(p1.diameter)) * self.lunit / 2 if parent: pcomp = id_to_comp[parent.id] moose.connect(comp, src, pcomp, dst) - sg_to_segments = {} + # map segment-group to segments + sg_to_segments = defaultdict(list) for sg in morphology.segment_groups: - sg_to_segments[sg.id] = [id_to_segment[m.segments] for m in sg.members] - for sg in morphology.segment_groups: - if not sg.id in sg_to_segments: - sg_to_segments[sg.id] = [] + for m in sg.members: + sg_to_segments[sg.id].append(id_to_segment[m.segments]) for inc in sg.includes: for cseg in sg_to_segments[inc.segment_groups]: sg_to_segments[sg.id].append(cseg) - - if not 'all' in sg_to_segments: - sg_to_segments['all'] = [ s for s in segments ] + if "all" not in sg_to_segments: + sg_to_segments["all"] = [s for s in segments] self._cell_to_sg[nmlcell.id] = sg_to_segments return id_to_comp, id_to_segment, sg_to_segments @@ -373,15 +406,20 @@ def importBiophysics(self, nmlcell, moosecell): according to NeuroML2 cell `nmlcell`.""" bp = nmlcell.biophysical_properties if bp is None: - logger_.info('Warning: %s in %s has no biophysical properties' % (nmlcell.id, self.filename)) + logger_.info( + "Warning: %s in %s has no biophysical properties" + % (nmlcell.id, self.filename) + ) return self.importMembraneProperties(nmlcell, moosecell, bp.membrane_properties) - self.importIntracellularProperties(nmlcell, moosecell, bp.intracellular_properties) + self.importIntracellularProperties( + nmlcell, moosecell, bp.intracellular_properties + ) def importMembraneProperties(self, nmlcell, moosecell, mp): """Create the membrane properties from nmlcell in moosecell""" if self.verbose: - logger_.info('Importing membrane properties') + logger_.info("Importing membrane properties") self.importCapacitances(nmlcell, moosecell, mp.specific_capacitances) self.importChannelsToCell(nmlcell, moosecell, mp) self.importInitMembPotential(nmlcell, moosecell, mp) @@ -414,7 +452,7 @@ def importSpecies(self, nmlcell, properties): # nml file here. concModel = species.concentration_model if (concModel is not None) and (concModel not in self.proto_pools): - logger_.warn("No concentrationModel '%s' found."%concModel) + logger_.warn("No concentrationModel '%s' found." % concModel) continue segments = getSegments(nmlcell, species, sg_to_segments) for seg in segments: @@ -434,18 +472,28 @@ def copySpecies(self, species, compartment): proto_pool = innerReader.proto_pools[concModel] break if not proto_pool: - msg = 'No prototype pool for %s referred to by %s' % (concModel, species.id) + msg = "No prototype pool for %s referred to by %s" % (concModel, species.id) logger_.error(msg) raise RuntimeError(msg) pool_id = moose.copy(proto_pool, compartment, species.id) pool = moose.element(pool_id) # print('&' * 10, compartment.path, compartment.length, compartment.diameter, pool.thick) if compartment.length <= 0: - vol = 4 * np.pi * (0.5 * compartment.diameter**3 - (0.5 * compartment.diameter - pool.thick)**3) / 3 + vol = ( + 4 + * np.pi + * ( + 0.5 * compartment.diameter**3 + - (0.5 * compartment.diameter - pool.thick) ** 3 + ) + / 3 + ) else: - vol = (np.pi * compartment.length * ( - 0.5 * compartment.diameter + pool.thick) * - (0.5 * compartment.diameter - pool.thick) + vol = ( + np.pi + * compartment.length + * (0.5 * compartment.diameter + pool.thick) + * (0.5 * compartment.diameter - pool.thick) ) pool.B = pool.B / vol return pool @@ -458,18 +506,18 @@ def importAxialResistance(self, nmlcell, intracellularProperties): comp = self.nml_segs_to_moose[seg.id] setRa(comp, SI(r.value)) - def isPassiveChan(self,chan): - if chan.type == 'ionChannelPassive': + def isPassiveChan(self, chan): + if chan.type == "ionChannelPassive": return True - if hasattr(chan,'gates'): - return len(chan.gate_hh_rates)+len(chan.gates)==0 + if hasattr(chan, "gates"): + return len(chan.gate_hh_rates) + len(chan.gates) == 0 return False def evaluate_moose_component(self, ct, variables): - print( "[INFO ] Not implemented." ) + print("[INFO ] Not implemented.") return False - def calculateRateFn(self, ratefn, vmin, vmax, tablen=3000, vShift='0mV'): + def calculateRateFn(self, ratefn, vmin, vmax, tablen=3000, vShift="0mV"): """Returns A / B table from ngate.""" from moose.neuroml2.hhfit import exponential2 from moose.neuroml2.hhfit import sigmoid2 @@ -477,22 +525,24 @@ def calculateRateFn(self, ratefn, vmin, vmax, tablen=3000, vShift='0mV'): from moose.neuroml2.units import SI rate_fn_map = { - 'HHExpRate': exponential2, - 'HHSigmoidRate': sigmoid2, - 'HHSigmoidVariable': sigmoid2, - 'HHExpLinearRate': linoid2 - } + "HHExpRate": exponential2, + "HHSigmoidRate": sigmoid2, + "HHSigmoidVariable": sigmoid2, + "HHExpLinearRate": linoid2, + } tab = np.linspace(vmin, vmax, tablen) if self._is_standard_nml_rate(ratefn): - midpoint, rate, scale = map(SI, (ratefn.midpoint, ratefn.rate, ratefn.scale)) + midpoint, rate, scale = map( + SI, (ratefn.midpoint, ratefn.rate, ratefn.scale) + ) return rate_fn_map[ratefn.type](tab, rate, scale, midpoint) for ct in self.doc.ComponentType: if ratefn.type != ct.name: continue - logger_.info(f"Using %s to evaluate rate"%ct.name) + logger_.info(f"Using %s to evaluate rate" % ct.name) rate = [] for v in tab: # Note: MOOSE HHGate are voltage and/or concentration @@ -504,37 +554,56 @@ def calculateRateFn(self, ratefn, vmin, vmax, tablen=3000, vShift='0mV'): # Find a suitable CaConc from the /library. Currently on Ca # dependent channels are allowed. caConcName = _findCaConcVariableName() - req_vars = {'v': '0.0V', 'caConc':f'{max(1e-11,v):g}', caConcName:f'{max(1e-11,v):g}','vShift':vShift,'temperature':self._getTemperature()} + req_vars = { + "v": "0.0V", + "caConc": f"{max(1e-11,v):g}", + caConcName: f"{max(1e-11,v):g}", + "vShift": vShift, + "temperature": self._getTemperature(), + } # logger_.debug(f"{'A' * 30} {req_vars}") else: - req_vars = {'v':'%sV'%v,'vShift':vShift,'temperature':self._getTemperature()} - req_vars.update( self._variables ) + req_vars = { + "v": "%sV" % v, + "vShift": vShift, + "temperature": self._getTemperature(), + } + req_vars.update(self._variables) vals = pynml.evaluate_component(ct, req_variables=req_vars) - v = vals.get('x', vals.get('t', vals.get('r', None))) + v = vals.get("x", vals.get("t", vals.get("r", None))) if v is not None: rate.append(v) return np.array(rate) - print( "[WARN ] Could not determine rate: %s %s %s" %(ratefn.type,vmin,vmax)) + print("[WARN ] Could not determine rate: %s %s %s" % (ratefn.type, vmin, vmax)) return np.array([]) - def importChannelsToCell(self, nmlcell, moosecell, membrane_properties): sg_to_segments = self._cell_to_sg[nmlcell.id] - for chdens in membrane_properties.channel_densities + membrane_properties.channel_density_v_shifts: + for chdens in ( + membrane_properties.channel_densities + + membrane_properties.channel_density_v_shifts + ): segments = getSegments(nmlcell, chdens, sg_to_segments) condDensity = SI(chdens.cond_density) erev = SI(chdens.erev) try: ionChannel = self.id_to_ionChannel[chdens.ion_channel] except KeyError: - logger_.info('No channel with id: %s' % chdens.ion_channel) + logger_.info("No channel with id: %s" % chdens.ion_channel) continue if self.verbose: - logger_.info('Setting density of channel %s in %s to %s; erev=%s (passive: %s)'%( - chdens.id, segments, condDensity,erev,self.isPassiveChan(ionChannel)) + logger_.info( + "Setting density of channel %s in %s to %s; erev=%s (passive: %s)" + % ( + chdens.id, + segments, + condDensity, + erev, + self.isPassiveChan(ionChannel), ) + ) if self.isPassiveChan(ionChannel): for seg in segments: @@ -543,9 +612,11 @@ def importChannelsToCell(self, nmlcell, moosecell, membrane_properties): setEk(comp, erev) else: for seg in segments: - self.copyChannel(chdens, self.nml_segs_to_moose[seg.id], condDensity, erev) - '''moose.le(self.nml_to_moose[seg]) - moose.showfield(self.nml_to_moose[seg], field="*", showtype=True)''' + self.copyChannel( + chdens, self.nml_segs_to_moose[seg.id], condDensity, erev + ) + """moose.le(self.nml_to_moose[seg]) + moose.showfield(self.nml_to_moose[seg], field="*", showtype=True)""" def copyChannel(self, chdens, comp, condDensity, erev): """Copy moose prototype for `chdens` condutcance density to `comp` @@ -561,49 +632,62 @@ def copyChannel(self, chdens, comp, condDensity, erev): proto_chan = innerReader.proto_chans[chdens.ion_channel] break if not proto_chan: - raise Exception('No prototype channel for %s referred to by %s' % (chdens.ion_channel, chdens.id)) + raise Exception( + "No prototype channel for %s referred to by %s" + % (chdens.ion_channel, chdens.id) + ) if self.verbose: - logger_.info('Copying the %s to %s, %s; erev=%s'%(chdens.id, comp, condDensity, erev)) + logger_.info( + "Copying the %s to %s, %s; erev=%s" + % (chdens.id, comp, condDensity, erev) + ) orig = chdens.id chid = moose.copy(proto_chan, comp, chdens.id) chan = moose.element(chid) els = list(self.paths_to_chan_elements.keys()) for p in els: - pp = p.replace('%s/'%chdens.ion_channel,'%s/'%orig) - self.paths_to_chan_elements[pp] = self.paths_to_chan_elements[p].replace('%s/'%chdens.ion_channel,'%s/'%orig) - #logger_.info(self.paths_to_chan_elements) + pp = p.replace("%s/" % chdens.ion_channel, "%s/" % orig) + self.paths_to_chan_elements[pp] = self.paths_to_chan_elements[p].replace( + "%s/" % chdens.ion_channel, "%s/" % orig + ) + # logger_.info(self.paths_to_chan_elements) chan.Gbar = sarea(comp) * condDensity chan.Ek = erev - moose.connect(chan, 'channel', comp, 'channel') + moose.connect(chan, "channel", comp, "channel") return chan def _is_standard_nml_rate(self, rate): - return rate.type=='HHExpLinearRate' \ - or rate.type=='HHExpRate' or \ - rate.type=='HHSigmoidRate' or \ - rate.type=='HHSigmoidVariable' + return ( + rate.type == "HHExpLinearRate" + or rate.type == "HHExpRate" + or rate.type == "HHSigmoidRate" + or rate.type == "HHSigmoidVariable" + ) def createHHChannel(self, chan, vmin=-150e-3, vmax=100e-3, vdivs=5000): - mchan = moose.HHChannel('%s/%s' % (self.lib.path, chan.id)) + mchan = moose.HHChannel("%s/%s" % (self.lib.path, chan.id)) mgates = [moose.element(x) for x in [mchan.gateX, mchan.gateY, mchan.gateZ]] - assert len(chan.gate_hh_rates)<=3, "We handle only up to 3 gates in HHCHannel" + assert len(chan.gate_hh_rates) <= 3, "We handle only up to 3 gates in HHCHannel" if self.verbose: - logger_.info('== Creating channel: %s (%s) -> %s (%s)'%(chan.id, chan.gate_hh_rates, mchan, mgates)) + logger_.info( + "== Creating channel: %s (%s) -> %s (%s)" + % (chan.id, chan.gate_hh_rates, mchan, mgates) + ) all_gates = chan.gates + chan.gate_hh_rates # Sort all_gates such that they come in x, y, z order. - all_gates = _gates_sorted( all_gates ) + all_gates = _gates_sorted(all_gates) for ngate, mgate in zip(all_gates, mgates): if ngate is None: continue - if mgate.name.endswith('X'): + if mgate.name.endswith("X"): mchan.Xpower = ngate.instances - elif mgate.name.endswith('Y'): + elif mgate.name.endswith("Y"): mchan.Ypower = ngate.instances - elif mgate.name.endswith('Z'): + elif mgate.name.endswith("Z"): mchan.Zpower = ngate.instances mgate.min = vmin mgate.max = vmax @@ -613,28 +697,38 @@ def createHHChannel(self, chan, vmin=-150e-3, vmax=100e-3, vdivs=5000): # HH-channels, the meaning of forwardRate and # reverseRate and steadyState are not clear in the # classes GateHHRatesInf, GateHHRatesTau and in - # FateHHTauInf the meaning of timeCourse and + # GateHHTauInf the meaning of timeCourse and # steady state is not obvious. Is the last one # refering to tau_inf and m_inf?? fwd = ngate.forward_rate rev = ngate.reverse_rate - self.paths_to_chan_elements['%s/%s'%(chan.id,ngate.id)] = '%s/%s'%(chan.id,mgate.name) + self.paths_to_chan_elements["%s/%s" % (chan.id, ngate.id)] = "%s/%s" % ( + chan.id, + mgate.name, + ) q10_scale = 1 if ngate.q10_settings: - if ngate.q10_settings.type == 'q10Fixed': - q10_scale= float(ngate.q10_settings.fixed_q10) - elif ngate.q10_settings.type == 'q10ExpTemp': - q10_scale = math.pow( float(ngate.q10_settings.q10_factor) - , (self._getTemperature()- SI(ngate.q10_settings.experimental_temp))/10 - ) - else: - raise Exception('Unknown Q10 scaling type %s: %s'%( - ngate.q10_settings.type,ngate.q10_settings) + if ngate.q10_settings.type == "q10Fixed": + q10_scale = float(ngate.q10_settings.fixed_q10) + elif ngate.q10_settings.type == "q10ExpTemp": + q10_scale = math.pow( + float(ngate.q10_settings.q10_factor), + ( + self._getTemperature() + - SI(ngate.q10_settings.experimental_temp) ) + / 10, + ) + else: + raise Exception( + "Unknown Q10 scaling type %s: %s" + % (ngate.q10_settings.type, ngate.q10_settings) + ) - logger_.debug('+ Gate: %s; %s; %s; %s; %s; scale=%s'% ( - ngate.id, mgate.path, mchan.Xpower, fwd, rev, q10_scale) - ) + logger_.debug( + "+ Gate: %s; %s; %s; %s; %s; scale=%s" + % (ngate.id, mgate.path, mchan.Xpower, fwd, rev, q10_scale) + ) if (fwd is not None) and (rev is not None): alpha = self.calculateRateFn(fwd, vmin, vmax, vdivs) @@ -644,8 +738,12 @@ def createHHChannel(self, chan, vmin=-150e-3, vmax=100e-3, vdivs=5000): mgate.tableB = q10_scale * (alpha + beta) # Assuming the meaning of the elements in GateHHTauInf ... - if hasattr(ngate,'time_course') and hasattr(ngate,'steady_state') \ - and (ngate.time_course is not None) and (ngate.steady_state is not None): + if ( + hasattr(ngate, "time_course") + and hasattr(ngate, "steady_state") + and (ngate.time_course is not None) + and (ngate.steady_state is not None) + ): tau = ngate.time_course inf = ngate.steady_state tau = self.calculateRateFn(tau, vmin, vmax, vdivs) @@ -653,7 +751,11 @@ def createHHChannel(self, chan, vmin=-150e-3, vmax=100e-3, vdivs=5000): mgate.tableA = q10_scale * (inf / tau) mgate.tableB = q10_scale * (1 / tau) - if hasattr(ngate,'steady_state') and (ngate.time_course is None) and (ngate.steady_state is not None): + if ( + hasattr(ngate, "steady_state") + and (ngate.time_course is None) + and (ngate.steady_state is not None) + ): inf = ngate.steady_state tau = 1 / (alpha + beta) if inf is not None: @@ -662,27 +764,27 @@ def createHHChannel(self, chan, vmin=-150e-3, vmax=100e-3, vdivs=5000): mgate.tableA = q10_scale * (inf / tau) mgate.tableB = q10_scale * (1 / tau) - logger_.info('%s: Created %s for %s'%(self.filename,mchan.path,chan.id)) + logger_.info("%s: Created %s for %s" % (self.filename, mchan.path, chan.id)) return mchan def createPassiveChannel(self, chan): - epath = '%s/%s' % (self.lib.path, chan.id) - if moose.exists( epath ): + epath = "%s/%s" % (self.lib.path, chan.id) + if moose.exists(epath): mchan = moose.element(epath) else: mchan = moose.Leakage(epath) - logger_.info('%s: Created %s for %s'%(self.filename,mchan.path,chan.id)) + logger_.info("%s: Created %s for %s" % (self.filename, mchan.path, chan.id)) return mchan def importInputs(self, doc): - epath = '%s/inputs' % (self.lib.path) - if moose.exists( epath ): - minputs = moose.element( epath ) + epath = "%s/inputs" % (self.lib.path) + if moose.exists(epath): + minputs = moose.element(epath) else: - minputs = moose.Neutral( epath ) + minputs = moose.Neutral(epath) for pg_nml in doc.pulse_generators: - epath = '%s/%s' % (minputs.path, pg_nml.id) + epath = "%s/%s" % (minputs.path, pg_nml.id) pg = moose.element(epath) if moose.exists(epath) else moose.PulseGen(epath) pg.firstDelay = SI(pg_nml.delay) pg.firstWidth = SI(pg_nml.duration) @@ -690,12 +792,11 @@ def importInputs(self, doc): pg.secondDelay = 1e9 logger_.debug(f'{"$" * 10} Created input {epath}') - def importIonChannels(self, doc, vmin=-150e-3, vmax=100e-3, vdivs=5000): - logger_.info('%s: Importing the ion channels' % self.filename ) + logger_.info("%s: Importing the ion channels" % self.filename) - for chan in doc.ion_channel+doc.ion_channel_hhs: - if chan.type == 'ionChannelHH': + for chan in doc.ion_channel + doc.ion_channel_hhs: + if chan.type == "ionChannelHH": mchan = self.createHHChannel(chan) elif self.isPassiveChan(chan): mchan = self.createPassiveChannel(chan) @@ -705,8 +806,11 @@ def importIonChannels(self, doc, vmin=-150e-3, vmax=100e-3, vdivs=5000): self.id_to_ionChannel[chan.id] = chan self.nml_chans_to_moose[chan.id] = mchan self.proto_chans[chan.id] = mchan - logger_.info(self.filename + ': Created ion channel %s for %s %s'%( - mchan.path, chan.type, chan.id)) + logger_.info( + self.filename + + ": Created ion channel %s for %s %s" + % (mchan.path, chan.type, chan.id) + ) def importConcentrationModels(self, doc): for concModel in doc.decaying_pool_concentration_models: @@ -714,19 +818,21 @@ def importConcentrationModels(self, doc): def createDecayingPoolConcentrationModel(self, concModel): """Create prototype for concentration model""" - if hasattr(concModel, 'name') and concModel.name is not None: + if hasattr(concModel, "name") and concModel.name is not None: name = concModel.name else: name = concModel.id - ca = moose.CaConc('%s/%s' % (self.lib.path, name)) + ca = moose.CaConc("%s/%s" % (self.lib.path, name)) ca.CaBasal = SI(concModel.resting_conc) ca.tau = SI(concModel.decay_constant) ca.thick = SI(concModel.shell_thickness) - ca.B = 5.2e-6 # B = 5.2e-6/(Ad) where A is the area of the - # shell and d is thickness - must divide by - # shell volume when copying + ca.B = 5.2e-6 # B = 5.2e-6/(Ad) where A is the area of the + # shell and d is thickness - must divide by + # shell volume when copying self.proto_pools[concModel.id] = ca self.nml_conc_to_moose[concModel.id] = ca self.moose_to_nml[ca] = concModel - logger_.debug('Created moose element: %s for nml conc %s' % (ca.path, concModel.id)) + logger_.debug( + "Created moose element: %s for nml conc %s" % (ca.path, concModel.id) + ) diff --git a/python/moose/neuroml2/test_granule98.py b/python/moose/neuroml2/test_granule98.py index c300edd103..85662f2e11 100644 --- a/python/moose/neuroml2/test_granule98.py +++ b/python/moose/neuroml2/test_granule98.py @@ -4,7 +4,7 @@ # Description: # Author: Subhasis Ray # Created: Mon Apr 8 21:41:22 2024 (+0530) -# Last-Updated: Wed Apr 10 20:16:04 2024 (+0530) +# Last-Updated: Wed Jul 17 15:51:22 2024 (+0530) # By: Subhasis Ray # @@ -40,15 +40,15 @@ def run(modeldir, nogui=True): moose.connect(vm, 'requestOut', soma, 'getVm') print('A' * 10, soma) - simtime = 300e-3 + simtime = 700e-3 moose.reinit() moose.start(simtime) t = np.arange(len(vm.vector)) * vm.dt print('%' * 10, len(vm.vector), len(inj.vector)) - results = np.array([t, vm.vector, inj.vector], dtype=[('time', float), ('Vm', float), ('Im', float)]) - fname = 'Granule_98.npy' - np.save(fname, results) + results = np.block([t[:, np.newaxis], vm.vector[:, np.newaxis], inj.vector[:, np.newaxis]]) + fname = 'Granule_98.dat' + np.savetxt(fname, X=results, header='time Vm Im', delimiter=' ') print(f'Saved results in {fname}') if not nogui: import matplotlib.pyplot as plt