Skip to content

Commit 2219c40

Browse files
adamjhnGoran Jelic-Cizmek
adamjhn
authored and
Goran Jelic-Cizmek
committed
Fix for dynamic ECS diffusion characteristics. (#2229)
* Fix for dynamic ECS diffusion characteristics. * Update extracellular.ipynb
1 parent 1d87d5d commit 2219c40

5 files changed

Lines changed: 136 additions & 22 deletions

File tree

‎docs/rxd-tutorials/extracellular.ipynb‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -428,6 +428,31 @@
428428
"metadata": {},
429429
"source": [
430430
"### Inhomogeneities\n",
431+
"#### Initial conditions\n",
432+
"Initial condition in the ECS need not be a scalar concentration. Suppose we only want to apply the buffer shown above to a sphere in the middle of the ECS, this can be achieved by passing function that takes a node as an argument."
433+
]
434+
},
435+
{
436+
"cell_type": "code",
437+
"execution_count": null,
438+
"metadata": {},
439+
"outputs": [],
440+
"source": [
441+
"del A, buffering\n",
442+
"A = rxd.Species(\n",
443+
" ecs,\n",
444+
" name=\"buffer\",\n",
445+
" charge=1,\n",
446+
" d=0,\n",
447+
" initial=lambda nd: 10 if nd.x3d**2 + nd.y3d**2 + nd.z3d**2 < 25**2 else 0,\n",
448+
")\n",
449+
"buffering = rxd.Reaction(k + A, AK, kf, kb)"
450+
]
451+
},
452+
{
453+
"cell_type": "markdown",
454+
"metadata": {},
455+
"source": [
431456
"#### Anisotropy\n",
432457
"The diffusion coefficient for a species can be different in each direction. e.g. to limit diffusion to the x,y-plane;"
433458
]

‎share/lib/python/neuron/rxd/region.py‎

Lines changed: 23 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -384,25 +384,32 @@ def _parse_tortuosity(self, value, is_permeability=False):
384384
ecs_permeability = h.Vector(parsed_value.flatten())
385385
else:
386386
ecs_permeability = h.Vector(parsed_value.flatten()).pow(-2)
387-
elif hasattr(value, "_extracellular_instances") or hasattr(
388-
self, "_extracellular"
387+
elif hasattr(value, "_extracellular_regions") or hasattr(
388+
value, "_extracellular"
389389
):
390390
# Check Species or SpeciesOnExtracellular is defined on this region
391391
if (
392-
hasattr(value, "_extracellular_instances")
393-
and self not in value._extracellular_instances
392+
hasattr(value, "_extracellular_regions")
393+
and self not in value._extracellular_regions
394394
) or (
395395
(
396396
hasattr(value, "_extracellular")
397397
and value._extracellular()
398-
and value._extracellular() != self
398+
and value._extracellular()._region != self
399399
)
400400
):
401401
raise RxDException(
402402
"permeability can be set to a State or Parameter, but it must be defined on this Extracellular region %r"
403403
% self
404404
)
405405
else:
406+
# make sure permeability has been initialized
407+
if hasattr(value, "_extracellular_regions") and not hasattr(
408+
value, "_extracellular_instances"
409+
):
410+
initializer._do_init()
411+
value._finitialize()
412+
406413
if is_permeability:
407414
parsed_value = value
408415
ecs_permeability = None
@@ -455,25 +462,31 @@ def _parse_volume_fraction(self, volume_fraction):
455462
warnings.warn(
456463
"Dynamic changes to the volume fraction does not change the concentrations of species already in the region."
457464
)
458-
elif hasattr(volume_fraction, "_extracellular_instances") or hasattr(
459-
self, "_extracellular"
465+
elif hasattr(volume_fraction, "_extracellular_regions") or hasattr(
466+
volume_fraction, "_extracellular"
460467
):
461468
# Check Species or SpeciesOnExtracellular is defined on this region
462469
if (
463-
hasattr(volume_fraction, "_extracellular_instances")
464-
and self not in volume_fraction._extracellular_instances
470+
hasattr(volume_fraction, "_extracellular_regions")
471+
and self not in volume_fraction._extracellular_regions
465472
) or (
466473
(
467474
hasattr(volume_fraction, "_extracellular")
468475
and volume_fraction._extracellular()
469-
and volume_fraction._extracellular() != self
476+
and volume_fraction._extracellular()._region != self
470477
)
471478
):
472479
raise RxDException(
473480
"volume fraction can be set to a State or Parameter, but it must be defined on this Extracellular region %r"
474481
% self
475482
)
476483
else:
484+
# make sure volume_fraction has been initialized
485+
if hasattr(volume_fraction, "_extracellular_regions") and not hasattr(
486+
volume_fraction, "_extracellular_instances"
487+
):
488+
initializer._do_init()
489+
volume_fraction._finitialize()
477490
alpha = volume_fraction
478491
ecs_alpha = None
479492
warnings.warn(

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

Lines changed: 27 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -133,6 +133,10 @@
133133
_set_tortuosity = nrn_dll_sym("set_tortuosity")
134134
_set_tortuosity.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.py_object]
135135

136+
# function to change extracellular volume fraction
137+
_set_volume_fraction = nrn_dll_sym("set_volume_fraction")
138+
_set_volume_fraction.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.py_object]
139+
136140

137141
# The difference here is that defined species only exists after rxd initialization
138142
_all_species = []
@@ -160,11 +164,14 @@ def _update_tortuosity(region):
160164
for s, r in _extracellular_diffusion_objects.items():
161165
if (
162166
r == region
163-
and hasattr(s, "_id")
167+
and hasattr(s, "_grid_id")
164168
and not hasattr(s, "_deleted")
165169
and not s._diffusion_characteristic
166170
):
167-
_set_tortuosity(s._id, region._permeability_vector)
171+
if hasattr(region._permeability_vector, "_ref_x"):
172+
_set_tortuosity(0, s._grid_id, region._permeability_vector._ref_x[0])
173+
else:
174+
_set_tortuosity(0, s._grid_id, region._permeability_vector)
168175

169176

170177
def _update_volume_fraction(region):
@@ -173,11 +180,16 @@ def _update_volume_fraction(region):
173180
for s, r in _extracellular_diffusion_objects.items():
174181
if (
175182
r == region
176-
and hasattr(s, "_id")
183+
and hasattr(s, "_grid_id")
177184
and not hasattr(s, "_deleted")
178185
and not s._diffusion_characteristic
179186
):
180-
_set_volume_fraction(s._id, region._volume_faction_vector)
187+
if hasattr(region._volume_fraction_vector, "_ref_x"):
188+
_set_volume_fraction(
189+
0, s._grid_id, region._volume_fraction_vector._ref_x[0]
190+
)
191+
else:
192+
_set_volume_fraction(0, s._grid_id, region._volume_fraction_vector)
181193

182194

183195
def _1d_submatrix_n():
@@ -399,12 +411,15 @@ def alpha_by_location(self, locs):
399411
e = self._extracellular()._region
400412
if numpy.isscalar(e.alpha):
401413
return e.alpha
402-
alphas = []
403-
for loc in locs:
404-
i = int((loc[0] - e._xlo) / e._dx[0])
405-
j = int((loc[1] - e._ylo) / e._dx[1])
406-
k = int((loc[2] - e._zlo) / e._dx[2])
407-
alphas.append(e.alpha[i, j, k])
414+
elif hasattr(e.alpha, "nodes"):
415+
alphas = [e.alpha.nodes((x, y, z)).value[0] for (x, y, z) in locs]
416+
else:
417+
alphas = []
418+
for loc in locs:
419+
i = int((loc[0] - e._xlo) / e._dx[0])
420+
j = int((loc[1] - e._ylo) / e._dx[1])
421+
k = int((loc[2] - e._zlo) / e._dx[2])
422+
alphas.append(e.alpha[i, j, k])
408423
return numpy.array(alphas)
409424

410425
def node_by_ijk(self, i, j, k):
@@ -1176,7 +1191,7 @@ def __init__(
11761191
11771192
11781193
"""
1179-
_extracellular_diffusion_objects[self] = None
1194+
_extracellular_diffusion_objects[self] = region
11801195

11811196
# ensure 3D points exist
11821197
h.define_shape()
@@ -1632,6 +1647,7 @@ def __init__(
16321647
# initialize self if the rest of rxd is already initialized
16331648
if initializer.is_initialized():
16341649
self._do_init()
1650+
self._finitialize()
16351651
rxd._update_node_data(True, True)
16361652

16371653
def _do_init(self):

‎test/rxd/ecs/test_ecs_example.py‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -196,3 +196,63 @@ def test_ecs_nodelists(ecs_example):
196196
# test accessing specific node by location
197197
nd = x[ecs].nodes((-38, 27, 27))[0]
198198
assert (nd.x3d, nd.y3d, nd.z3d) == (-38.5, 27.5, 27.5)
199+
200+
201+
def test_ecs_example_dynamic_tort(ecs_example):
202+
"""Test ecs_example with dynamic tortuosity"""
203+
204+
(h, rxd, data, save_path), make_model = ecs_example
205+
model = make_model(lambda x, y, z: 1.0, 1)
206+
h.finitialize(-65)
207+
(
208+
cell1,
209+
cell2,
210+
cyt,
211+
org,
212+
cyt_org_membrane,
213+
ecs,
214+
x,
215+
Xcyt,
216+
Xorg,
217+
createX,
218+
cell1_param,
219+
createX,
220+
cyt_org_leak,
221+
) = model
222+
perm = rxd.Species(ecs, name="perm", initial=1.0 / 1.6**2)
223+
ecs.permeability = perm
224+
h.continuerun(1000)
225+
226+
if not save_path:
227+
max_err = compare_data(data)
228+
assert max_err < tol
229+
230+
231+
def test_ecs_example_dynamic_alpha(ecs_example):
232+
"""Test ecs_example with fixed step and inhomogeneous tortuosity methods"""
233+
234+
(h, rxd, data, save_path), make_model = ecs_example
235+
model = make_model(lambda x, y, z: 1.0, 1.6)
236+
h.finitialize(-65)
237+
(
238+
cell1,
239+
cell2,
240+
cyt,
241+
org,
242+
cyt_org_membrane,
243+
ecs,
244+
x,
245+
Xcyt,
246+
Xorg,
247+
createX,
248+
cell1_param,
249+
createX,
250+
cyt_org_leak,
251+
) = model
252+
alpha = rxd.Species(ecs, name="alpha", initial=0.2)
253+
ecs.alpha = alpha
254+
h.continuerun(1000)
255+
256+
if not save_path:
257+
max_err = compare_data(data)
258+
assert max_err < tol

‎test/rxd/testdata‎

0 commit comments

Comments
 (0)