Hi,
We are extending the reaction-diffusion module (rxd) to support extracellular space, it’s not currently implement in the distributed version of NEURON, but the source code is available on github
https://github.com/ramcdougal/nrn.git
Here is a simple example of two cells in the extracellular space defined by rxd.Extracellular, the volume fraction is the relative proportion of free space (typically 20%) and the tortuosity reduces extracellular diffusion due to obstacles (typically ~1.6) in brain tissue (See
http://www.sciencedirect.com/science/ar ... 3698012612)
If your problem doesn’t require extracellular diffusion, you could make dx equal to the length of the grid, so it is all in one voxel.
The initial concentrations are taken from NEURON defaults (h.nai0_na_ion, h.nao0_na_ion, etc.) and the concentration can be accessed via the NEURON sections or the rxd Species objects.
Code: Select all
from neuron import rxd, h
import numpy
from matplotlib import pyplot
h.load_file('stdrun.hoc')
soma1 = h.Section(name='soma1')
soma2 = h.Section(name='soma2')
# Geometry
soma1.nseg = 5
soma1.L = soma1.diam = 20
soma2.nseg = 5
soma2.L = soma2.diam = 20
# insert mechanisms into the cells
soma1.insert('hh')
soma2.insert('hh')
# stimulation
stim1 = h.IClamp(soma1(0.5))
stim1.delay = 0
stim1.dur = 50
stim1.amp = 0.015
stim2 = h.IClamp(soma2(0.5))
stim2.delay = 50
stim2.dur = 50
stim2.amp = 0.03
# enable the extracellular space
rxd.options.enable.extracellular = True
# define the grid with the volume fraction (the proportion of free volume)
# and the tortuosity (the increase in path length due to obstacles)
extracellular = rxd.Extracellular(xlo=-25, ylo=-25, zlo=-100, xhi = 25, yhi = 25, zhi = 100, dx = 10, volume_fraction=0.2, tortuosity=1.6)
# the intracellular region
cyt = rxd.Region(h.allsec(), nrn_region='i')
# the ion species
rxd_na = rxd.Species([extracellular,cyt], name = 'na', charge = 1, d = 1)
rxd_k = rxd.Species([extracellular, cyt], name = 'k', charge = 1, d = 1)
#record into vectors
#intracellular (i) and extracellular (o) concentrations
nai_vec1 = h.Vector()
nai_vec1.record(soma1(0.5)._ref_nai)
nao_vec1 = h.Vector()
nao_vec1.record(soma1(0.5)._ref_nao)
ki_vec1 = h.Vector()
ki_vec1.record(soma1(0.5)._ref_ki)
ko_vec1 = h.Vector()
ko_vec1.record(soma1(0.5)._ref_ko)
nai_vec2 = h.Vector()
nai_vec2.record(soma2(0.5)._ref_nai)
nao_vec2 = h.Vector()
nao_vec2.record(soma2(0.5)._ref_nao)
ki_vec2 = h.Vector()
ki_vec2.record(soma2(0.5)._ref_ki)
ko_vec2 = h.Vector()
ko_vec2.record(soma2(0.5)._ref_ko)
# membrane potential
v_vec1 = h.Vector()
v_vec2 = h.Vector()
v_vec1.record(soma1(0.5)._ref_v)
v_vec2.record(soma2(0.5)._ref_v)
t_vec = h.Vector()
t_vec.record(h._ref_t)
h.finitialize()
# access the intracellular and extracellular species
na_ics = rxd_na[cyt]
k_ics = rxd_k[cyt]
na_ecs = rxd_na[extracellular]
k_ecs = rxd_k[extracellular]
print "Initial concentrations\t\tintracellular\t\textracellular"
print "Na+\t\t\t\t%g\t\t\t%g" % (numpy.mean(na_ics.nodes.concentration), na_ecs.states3d.mean())
print "K+\t\t\t\t%g\t\t\t%g" % (numpy.mean(k_ics.nodes.concentration), k_ecs.states3d.mean())
h.tstop=100
h.run()
print "Final concentrations\t\tintracellular\t\textracellular"
print "Na+\t\t\t\t%g\t\t\t%g" % (numpy.mean(na_ics.nodes.concentration), na_ecs.states3d.mean())
print "K+\t\t\t\t%g\t\t\t%g" % (numpy.mean(k_ics.nodes.concentration), k_ecs.states3d.mean())
pyplot.figure(figsize=(11,6))
pyplot.subplot(3,2,(1,2))
pyplot.plot(t_vec, v_vec1)
pyplot.plot(t_vec, v_vec2)
pyplot.xlabel('time (ms)')
pyplot.ylabel('mV')
pyplot.xlim(0,100)
pyplot.subplot(3,2,3)
pyplot.plot(t_vec, ki_vec1)
pyplot.plot(t_vec, ki_vec2)
pyplot.xlabel('time (ms)')
pyplot.ylabel('intracellular K+ mM')
pyplot.xlim(0,100)
pyplot.subplot(3,2,4)
pyplot.plot(t_vec, nai_vec1)
pyplot.plot(t_vec, nai_vec2)
pyplot.xlabel('time (ms)')
pyplot.ylabel('intracellular Na+ mM')
pyplot.xlim(0,100)
pyplot.subplot(3,2,5)
pyplot.plot(t_vec, ko_vec1)
pyplot.plot(t_vec, ko_vec2)
pyplot.xlabel('time (ms)')
pyplot.ylabel('extracellular K+ mM')
pyplot.xlim(0,100)
pyplot.subplot(3,2,6)
pyplot.plot(t_vec, nao_vec1)
pyplot.plot(t_vec, nao_vec2)
pyplot.xlabel('time (ms)')
pyplot.ylabel('extracellular Na+ mM')
pyplot.xlim(0,100)
pyplot.show()
https://slack-files.com/T0DM7EUDV-F5XJXSNGH-2aa6646349
Best wishes,
Adam.