Fix rxd initialization of extracellular species. (#1326)

share/lib/python/neuron/rxd/species.py

@@ -1118,7 +1118,7 @@ def _locate_segments(self):
     def _update_pointers(self):
         # TODO: call this anytime the _grid_id changes and anytime the structure_change_count change
         self._isalive()
-        if self._species is None or self._charge == 0:
+        if self._species is None:
             return
         self._seg_indices = self._locate_segments()
         from .geometry import _surface_areas1d
@@ -1133,7 +1133,7 @@ def _update_pointers(self):
                     grid_indices.append(i)
                     neuron_pointers.append(getattr(seg, stateo))
         _set_grid_concentrations(grid_list, self._grid_id, grid_indices, neuron_pointers)
-        if isinstance(_defined_species[self._species][self._region](), Parameter):
+        if self._charge == 0 or isinstance(_defined_species[self._species][self._region](), Parameter):
             _set_grid_currents(grid_list, self._grid_id, [], [], [])
         else:
             tenthousand_over_charge_faraday = 10000. / (self._charge * h.FARADAY)

share/lib/python/neuron/rxdtests/tests/hh_param.py

@@ -75,6 +75,9 @@ def init(ics,ecs):
 
 x = rxd.Species([cyt, mem, ecs], name='x', charge=1, initial=1e9)
 
+# a parameter without a charge to test initialization 
+dump = rxd.Parameter([cyt, ecs], name='dump')
+
 ki, ko, nai, nao, xi, xo = k[cyt], kecs[ecs], na[cyt], naecs[ecs], x[cyt], x[ecs]
 
 # gates

test/rxd/ecs/test_ecs_reinit.py

@@ -13,30 +13,35 @@ def simple_model(neuron_instance):
     cyt = rxd.Region(dend, name='cyt', nrn_region='i')
     k = rxd.Species([cyt, ecs], name='k', d=1, charge=1,
                     initial=lambda nd: 140 if nd.region == cyt else 3)
+    # test initialization for parameters with and without names
+    paramA = rxd.Parameter([ecs], name='paramA', initial=1)
+    paramB = rxd.Parameter([ecs], initial=0)
     decay = rxd.Rate(k, -0.1*k)
-    model = (dend, cyt, ecs, k, decay)
+    model = (dend, cyt, ecs, k, paramA, paramB,  decay)
     yield (neuron_instance, model)
 
 def test_ecs_reinit(simple_model):
     """Test rxd.re_init updates extracellular node values from NEURON values"""
 
     neuron_instance, model = simple_model
     h, rxd, data, save_path = neuron_instance
-    dend, cyt, ecs, k, decay = model
+    dend, cyt, ecs, k, paramA, paramB, decay = model
     h.finitialize(-65)
     dend(0.2).ko = 0
     rxd.re_init()
     assert(k[ecs].nodes((0,0,0)).value == [0])
+    assert(dend(0.2).paramAo == 1)
 
 def test_ecs_reinit_cvode(simple_model):
     """Test rxd.re_init updates extracellular node values from NEURON segments
        with CVode """
 
     neuron_instance, model = simple_model
     h, rxd, data, save_path = neuron_instance
-    dend, cyt, ecs, k, decay = model
+    dend, cyt, ecs, k, paramA, paramB, decay = model
     h.CVode().active(True)
     h.finitialize(-65)
     dend(0.2).ko = 0
     rxd.re_init()
     assert(k[ecs].nodes((0,0,0)).value == [0])
+    assert(dend(0.2).paramAo == 1)

test/rxd/test_reinit.py

@@ -12,16 +12,18 @@ def simple_model(neuron_instance):
     cyt = rxd.Region(dend, name='cyt', nrn_region='i')
     k = rxd.Species([cyt], name='k', d=1, charge=1,
                     initial=140)
+    paramA = rxd.Parameter([cyt], name='paramA', initial=1)
+    paramB = rxd.Parameter([cyt], initial=0)
     decay = rxd.Rate(k, -0.1*k)
-    model = (dend, cyt, k, decay)
+    model = (dend, cyt, k, paramA, paramB, decay)
     yield (neuron_instance, model)
 
 def test_reinit(simple_model):
     """Test rxd.re_init updates node values from NEURON values"""
 
     neuron_instance, model = simple_model
     h, rxd, data, save_path = neuron_instance
-    dend, cyt, k, decay = model
+    dend, cyt, k, paramA, paramB, decay = model
     h.finitialize(-65)
     dend(0.5).ki = 0
     rxd.re_init()
@@ -32,7 +34,7 @@ def test_reinit_cvode(simple_model):
 
     neuron_instance, model = simple_model
     h, rxd, data, save_path = neuron_instance
-    dend, cyt, k, decay = model
+    dend, cyt, k, paramA, paramB, decay = model
     h.finitialize(-65)
     h.CVode().active(True)
     dend(0.5).ki = 0
@@ -44,7 +46,7 @@ def test_reinit_3d(simple_model):
 
     neuron_instance, model = simple_model
     h, rxd, data, save_path = neuron_instance
-    dend, cyt, k, decay = model
+    dend, cyt, k, paramA, paramB, decay = model
     rxd.set_solve_type(dimension=3)
     # check changing the units after initialization
     h.finitialize(-65)
@@ -60,7 +62,7 @@ def test_reinit_3d_cvode(simple_model):
 
     neuron_instance, model = simple_model
     h, rxd, data, save_path = neuron_instance
-    dend, cyt, k, decay = model
+    dend, cyt, k, paramA, paramB, decay = model
     rxd.set_solve_type(dimension=3)
     h.CVode().active(True)
     # check changing the units after initialization