I pushed a version that supports extracellular diffusion to the development branch at
http://github.com/ramcdougal/nrn.
To try it, clone from there, compile, and run the code below. You should get two windows, one that shows a red square on a blue background (the initial condition), and one where diffusion over time has turned the square into a circle. (It's actually a 3D simulation, but the square and circle are the cross-sections.)
Code: Select all
from neuron import h, rxd
from matplotlib import pyplot
import time
# extracellular disabled by default until fully tested, so must override
rxd.options.enable.extracellular = True
# Need to have at least one Section or CVode.statistics will crash
s = h.Section()
h.load_file('stdrun.hoc')
extracellular = rxd.Extracellular(-100, -100, -50, 100, 100, 50, dx=4)
ca = rxd.Species(extracellular, name='ca', charge=2, d=1, initial=0)
h.CVode().active(1)
h.finitialize()
ca_ext = ca[extracellular]
ca_ext.states3d[5:15, 5:15, :] = 1
# NEURON fact-of-life: always have to call this when manually changing values in a CVode simulation
h.CVode().re_init()
pyplot.figure()
pyplot.imshow(ca_ext.states3d[:, :, 0].copy(), interpolation='nearest', vmin=0, vmax=1, extent=ca_ext.extent('xy'), origin='lower')
pyplot.title('Initial, plane 0')
print 'initial sum: %g' % ca_ext.states3d.sum()
print 'advancing...'
start = time.time()
h.continuerun(100)
print '... advance took %gs' % (time.time() - start)
print 'final sum: %g' % ca_ext.states3d.sum()
print 'final max: %g' % ca_ext.states3d.max()
h.CVode().statistics()
pyplot.figure()
pyplot.imshow(ca_ext.states3d[:, :, 0], interpolation='nearest', vmin=0, vmax=1, extent=ca_ext.extent('xy'), origin='lower')
pyplot.title('t = %g, plane 0' % h.t)
pyplot.show()
The test assumes that the extracellular space occupies a volume_fraction=1 (i.e. 100% of the space; that there is nothing else out there) and that the tortuosity=1 (again, basically, that there's nothing else out there), but you can change this by sending in the appropriate options (scalars or arrays) to rxd.Extracellular.
At the beginning, there's a line that specifically enables extracellular support. Why? Because it's an experimental, incomplete version and that's there to remind you of that... in particular note that reactions on the extracellular space are not yet implemented in this version, and that non-unity volume fractions and tortuosity currently are only respected in fixed step simulations. (That's coming in a future commit though.)