Source code for rxdReacDiffusion

import numpy
import pylab
import moose
import time
[docs]def main(): """ This example implements a reaction-diffusion like system which is bistable and propagates losslessly. It is based on the NEURON example rxdrun.py, but incorporates more compartments and runs for a longer time. The system is implemented as a hybrid of a reaction and a function which sets its rates. Please see rxdFuncDiffusion.py for a variant that uses just a function object to set up the system. """ dt = 0.1 # define the geometry compt = moose.CylMesh( '/cylinder' ) compt.r0 = compt.r1 = 1 compt.x1 = 100 compt.diffLength = 0.2 assert( compt.numDiffCompts == compt.x1/compt.diffLength ) #define the molecule. Its geometry is defined by its parent volume, cylinder c = moose.Pool( '/cylinder/pool' ) c.diffConst = 1 # define diffusion constant # There is an implicit reaction substrate/product. MOOSE makes it explicit. buf = moose.BufPool( '/cylinder/buf' ) buf.nInit = 1 # The reaction is something entirely peculiar, not a chemical thing. reaction = moose.Reac( '/cylinder/reac' ) reaction.Kb = 0 # so here we set up a function calculation to do the same thing. func = moose.Function( '/cylinder/reac/func' ) func.expr = "(1 - x0) * (0.3 - x0)" func.x.num = 1 #specify number of input variables. #Connect the reaction to the pools moose.connect( reaction, 'sub', c, 'reac' ) moose.connect( reaction, 'prd', buf, 'reac' ) #Connect the function to the reaction moose.connect( func, 'valueOut', reaction, 'setNumKf' ) #Connect the molecules to the func moose.connect( c, 'nOut', func.x[0], 'input' ) #Set up solvers ksolve = moose.Ksolve( '/cylinder/ksolve' ) dsolve = moose.Dsolve( '/cylinder/dsolve' ) stoich = moose.Stoich( '/cylinder/stoich' ) stoich.compartment = compt stoich.ksolve = ksolve stoich.dsolve = dsolve stoich.path = '/cylinder/##' for i in range( 10, 18 ): moose.setClock( i, dt ) #initialize x = numpy.arange( 0, compt.x1, compt.diffLength ) c.vec.nInit = [ (q < 0.2 * compt.x1) for q in x ] # Run and plot it. moose.reinit() updateDt = 50 runtime = updateDt * 4 plt = pylab.plot( x, c.vec.n, label='t = 0 ') t1 = time.time() for t in range( 0, runtime-1, updateDt ): moose.start( updateDt ) plt = pylab.plot( x, c.vec.n, label='t = '+str(t + updateDt) ) print(("Time = ", time.time() - t1)) pylab.ylim( 0, 1.05 ) pylab.legend() pylab.show()
if __name__ == '__main__': main()