If nml user has set id of the gate to be 'x', 'y' or 'z' then load the

dynamics into HHChannels's gateX, gateY, or gateZ dyamics.

Author: Dilawar Singh

python/moose/neuroml2/reader.py

@@ -43,15 +43,42 @@
         "********************************************************************")
 
 
+# these gates are available in moose. These are prefixed by 'gate'
+_validMooseHHGateIds = ['X', 'Y', 'Z'] 
+
 def _unique(ls):
-    res = [ ]
+    res = []
     for l in ls:
         if l not in res:
             res.append(l)
     return res
 
+
 def _whichGate(chan):
-    return chan.name[-1]
+    global _validMooseHHGateIds
+    c = chan.name[-1]
+    assert c in _validMooseHHGateIds
+    return c
+
+def _pairNmlGateWithMooseGates(mGates, nmlGates):
+    """Return moose gate id from nml.HHGate
+    """
+    global _validMooseHHGateIds
+    # deep copy
+    mooseGatesMap = {_whichGate(x) : x for x in mGates}
+    availableMooseGates = _validMooseHHGateIds[:]
+    mapping = {}
+    for nmlGate in nmlGates: 
+        if nmlGate is None:
+            continue
+        if hasattr(nmlGate, 'id') and nmlGate.id:
+            mapping[nmlGate.id.upper()] = nmlGate
+            availableMooseGates.remove(nmlGate.id.upper())
+        else:
+            mapping[availableMooseGates.pop(0)] = nmlGate
+
+    # Now replace 'X', 'Y', 'Z' with moose gates.
+    return [(mooseGatesMap[x], mapping[x]) for x in mapping]
 
 def _isConcDep(ct):
     """_isConcDep
@@ -75,9 +102,10 @@ def _findCaConc():
     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]' )
+    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]
 
 
@@ -200,7 +228,6 @@ def __init__(self, verbose=False):
         # Just in case.
         self._variables = {}
 
-
     def read(self, filename, symmetric=True):
         self.doc = loaders.read_neuroml2_file(filename,
                                               include_includes=True,
@@ -464,8 +491,9 @@ def isPassiveChan(self, chan):
         'HHExpLinearRate': linoid2
     }
 
-    def calculateRateFn(self, ratefn, vmin, vmax, tablen=3000, vShift='0mV'):
+    def calculateRateFn(self, ratefn, mgate, vmin, vmax, tablen=3000, vShift='0mV'):
         """Returns A / B table from ngate."""
+
         tab = np.linspace(vmin, vmax, tablen)
         if self._is_standard_nml_rate(ratefn):
             midpoint, rate, scale = map(
@@ -476,24 +504,62 @@ def calculateRateFn(self, ratefn, vmin, vmax, tablen=3000, vShift='0mV'):
             if ratefn.type != ct.name:
                 continue
             logger_.info("Using %s to evaluate rate" % ct.name)
-            rate = []
-            for v in tab:
-                req_vars = {'v':'%sV'%v,'vShift':vShift,'temperature':self._getTemperature()}
-                req_vars.update(self._variables)
-                vals = pynml.evaluate_component(ct, req_variables=req_vars)
-                '''print(vals)'''
-                if 'x' in vals:
-                    rate.append(vals['x'])
-                if 't' in vals:
-                    rate.append(vals['t'])
-                if 'r' in vals:
-                    rate.append(vals['r'])
-            return np.array(rate)
-
-    def calculateRateFnCaDep(self, ratefn, ca, camin, camax, tablen=3000, vShift='0mV'):
+            if not _isConcDep(ct):
+                return self._computeRateFn(ct, tab)
+            else:
+                ca = _findCaConc()
+                assert _whichGate(mgate) == 'Z', mgate
+                return self._computeRateFnCa(ca, ct, tab, vShift=vShift)
+
+    def _computeRateFnCa(self, ca, ct, tab, vShift):
+        rate = []
+        for v in tab:
+            req_vars = {
+                ca.name: '%sV' % v,
+                'vShift': vShift,
+                'temperature': self._getTemperature()
+            }
+            req_vars.update(self._variables)
+            vals = pynml.evaluate_component(ct, req_variables=req_vars)
+            '''print(vals)'''
+            if 'x' in vals:
+                rate.append(vals['x'])
+            if 't' in vals:
+                rate.append(vals['t'])
+            if 'r' in vals:
+                rate.append(vals['r'])
+        return np.array(rate)
+
+    def _computeRateFn(self, ct, tab, vShift):
+        rate = []
+        for v in tab:
+            req_vars = {
+                'v': '%sV' % v,
+                'vShift': vShift,
+                'temperature': self._getTemperature()
+            }
+            req_vars.update(self._variables)
+            vals = pynml.evaluate_component(ct, req_variables=req_vars)
+            '''print(vals)'''
+            if 'x' in vals:
+                rate.append(vals['x'])
+            if 't' in vals:
+                rate.append(vals['t'])
+            if 'r' in vals:
+                rate.append(vals['r'])
+        return np.array(rate)
+
+    def calculateRateFnCaDep(self,
+                             ratefn,
+                             ca,
+                             camin,
+                             camax,
+                             tablen=3000,
+                             vShift='0mV'):
         """Returns A / B table from ngate.
 
-        FIXME: Merge with calculateRateFn
+        This function compute rate which depends on ca. 
+        This must go into Z gate.
         """
         tab = np.linspace(camin, camax, tablen)
         if self._is_standard_nml_rate(ratefn):
@@ -504,10 +570,15 @@ def calculateRateFnCaDep(self, ratefn, ca, camin, camax, tablen=3000, vShift='0m
         for ct in self.doc.ComponentType:
             if ratefn.type != ct.name:
                 continue
-            logger_.info("Using %s to evaluate rate (caConc dependant)" % ct.name)
+            logger_.info("Using %s to evaluate rate (caConc dependant)" %
+                         ct.name)
             rate = []
             for v in tab:
-                req_vars = {ca.name:'%sV'%v,'vShift':vShift,'temperature':self._getTemperature()}
+                req_vars = {
+                    ca.name: '%sV' % v,
+                    'vShift': vShift,
+                    'temperature': self._getTemperature()
+                }
                 vals = pynml.evaluate_component(ct, req_variables=req_vars)
                 '''print(vals)'''
                 if 'x' in vals:
@@ -615,107 +686,109 @@ def _is_standard_nml_rate(self, rate):
                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))
-        mgates = map(moose.element, (mchan.gateX, mchan.gateY, mchan.gateZ))
+        path = '%s/%s' % (self.lib.path, chan.id)
+        if moose.exists(path):
+            mchan = moose.element(path)
+        else:
+            mchan =  moose.HHChannel(path)
+        mgates = [
+            moose.element(g) for g in [mchan.gateX, mchan.gateY, mchan.gateZ]
+        ]
 
         # We handle only up to 3 gates in HHCHannel
-        assert (len(chan.gate_hh_rates) <= 3)  
+        assert len(chan.gate_hh_rates) <= 3, "No more than 3 gates"
 
         if self.verbose:
             print('== Creating channel: %s (%s) -> %s (%s)' %
                   (chan.id, chan.gate_hh_rates, mchan, mgates))
+
         all_gates = chan.gates + chan.gate_hh_rates
-        for ngate, mgate in zip(all_gates, mgates):
-            if ngate is None:
-                continue
-            logger_.info('whichGate(mgate) %s %s' %(mgate, _whichGate(mgate)))
-            if _whichGate(mgate) == 'X':
-                mchan.Xpower = ngate.instances
-            elif _whichGate(mgate) == 'Y':
-                mchan.Ypower = ngate.instances
-            elif _whichGate(mgate) == 'Z':
-                mchan.Zpower = ngate.instance
-            mgate.min = vmin
-            mgate.max = vmax
-            mgate.divs = vdivs
-
-            # I saw only examples of GateHHRates in
-            # 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
-            # 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)
-
-            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)
-                    #print('Q10: %s; %s; %s; %s'%(ngate.q10_settings.q10_factor, self._getTemperature(),SI(ngate.q10_settings.experimental_temp),q10_scale))
-                else:
-                    raise Exception(
-                        'Unknown Q10 scaling type %s: %s' %
-                        (ngate.q10_settings.type, ngate.q10_settings))
 
-            logger_.info(' === Gate: %s; %s; %s; %s; %s; scale=%s' %
-                    (ngate.id, mgate.path, mchan.Xpower, fwd, rev, q10_scale))
+        # If user set bnml channels' id to 'x', 'y' or 'z' then pair this gate
+        # with moose.HHChannel's gateX, gateY, gateZ respectively. Else pair
+        # them with gateX, gateY, gateZ acording to list order.
+        for mgate, ngate in _pairNmlGateWithMooseGates(mgates, all_gates):
+            self._addGateToHHChannel(chan, mchan, mgate, ngate, vmin, vmax,
+                                     vdivs)
+        logger_.debug('== Created', mchan.path, 'for', chan.id)
+        return mchan
 
-            print(fwd, rev)
-            quit()
+    def _addGateToHHChannel(self, chan, mchan, mgate, ngate, vmin, vmax,
+                            vdivs):
+        """Add gateX, gateY, gateZ etc to moose.HHChannel (mchan). 
 
-            if (fwd is not None) and (rev is not None):
-                # Note: MOOSE HHGate are either voltage of concentration
-                # dependant. Here we figure out if nml description of gate is
-                # concentration dependant or not. 
-                if not _isConcDep(fwd):
-                    alpha = self.calculateRateFn(fwd, vmin, vmax, vdivs)
-                else:
-                    ca = _findCaConc()
-                    alpha = self.calculateRateFnCaDep(fwd, ca, vmin, vmax, vdivs)
+        Each gate can be voltage dependant and/or concentration dependant.
+        Only caConc dependant channels are supported.
+        """
 
-                if not _isConcDep(rev):
-                    beta = self.calculateRateFnCaDep(rev, vmin, vmax, vdivs)
-                else:
-                    ca = _findCaConc()
-                    beta = self.calculateRateFn(rev, ca, vmin, vmax, vdivs)
-
-                mgate.tableA = q10_scale * (alpha)
-                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):
-                tau = ngate.time_course
-                inf = ngate.steady_state
-                tau = self.calculateRateFn(tau, vmin, vmax, vdivs)
+        # set mgate.Xpower, .Ypower etc.
+        setattr(mchan, _whichGate(mgate) + 'power', ngate.instances)
+
+        mgate.min = vmin
+        mgate.max = vmax
+        mgate.divs = vdivs
+
+        # Note by Padraig: 
+        # ---------------
+        # I saw only examples of GateHHRates in 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 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)
+
+        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)
+                logger_.debug(
+                    'Q10: %s; %s; %s; %s' %
+                    (ngate.q10_settings.q10_factor, self._getTemperature(),
+                     SI(ngate.q10_settings.experimental_temp), q10_scale))
+            else:
+                raise Exception('Unknown Q10 scaling type %s: %s' %
+                                (ngate.q10_settings.type, ngate.q10_settings))
+        logger_.info(' === 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):
+            # Note: MOOSE HHGate are either voltage of concentration
+            # dependant. Here we figure out if nml description of gate is
+            # concentration dependant or not.
+            alpha = self.calculateRateFn(fwd, mgate, vmin, vmax, vdivs)
+            beta = self.calculateRateFn(rev, mgate, vmin, vmax, vdivs)
+
+            mgate.tableA = q10_scale * (alpha)
+            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):
+            tau = ngate.time_course
+            inf = ngate.steady_state
+            tau = self.calculateRateFn(tau, mgate, vmin, vmax, vdivs)
+            inf = self.calculateRateFn(inf, mgate, vmin, vmax, vdivs)
+            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):
+            inf = ngate.steady_state
+            tau = 1 / (alpha + beta)
+            if (inf is not None):
                 inf = self.calculateRateFn(inf, vmin, vmax, vdivs)
                 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):
-                inf = ngate.steady_state
-                tau = 1 / (alpha + beta)
-                if (inf is not None):
-                    inf = self.calculateRateFn(inf, vmin, vmax, vdivs)
-                    mgate.tableA = q10_scale * (inf / tau)
-                    mgate.tableB = q10_scale * (1 / tau)
-
-        if self.verbose:
-            print(self.filename, '== Created', mchan.path, 'for', chan.id)
-        return mchan
-
     def createPassiveChannel(self, chan):
         mchan = moose.Leakage('%s/%s' % (self.lib.path, chan.id))
         if self.verbose: