Source code for cylinderDiffusion

#########################################################################
## This program is part of 'MOOSE', the
## Messaging Object Oriented Simulation Environment.
##           Copyright (C) 2014 Upinder S. Bhalla. and NCBS
## It is made available under the terms of the
## GNU Lesser General Public License version 2.1
## See the file COPYING.LIB for the full notice.
#########################################################################


import sys
sys.path.append('../../python')
import math
import pylab
import numpy
import matplotlib.pyplot as plt
import moose

import os
import signal
PID = os.getpid()

def doNothing( *args ):
    pass

signal.signal( signal.SIGUSR1, doNothing )

def makeModel():
    # create container for model
    r0 = 2e-6        # m
    r1 = 1e-6        # m
    num = 100
    diffLength = 1e-6 # m
    len = num * diffLength        # m
    diffConst = 10e-12 # m^2/sec
    motorRate = 10e-6 # m/sec
    concA = 1 # millimolar

    model = moose.Neutral( 'model' )
    compartment = moose.CylMesh( '/model/compartment' )
    compartment.r0 = r0
    compartment.r1 = r1
    compartment.x0 = 0
    compartment.x1 = len
    compartment.diffLength = diffLength

    assert( compartment.numDiffCompts == num )

    # create molecules and reactions
    a = moose.Pool( '/model/compartment/a' )
    b = moose.Pool( '/model/compartment/b' )
    c = moose.Pool( '/model/compartment/c' )
    d = moose.BufPool( '/model/compartment/d' )
    r1 = moose.Reac( '/model/compartment/r1' )
    moose.connect( r1, 'sub', b, 'reac' )
    moose.connect( r1, 'sub', d, 'reac' )
    moose.connect( r1, 'prd', c, 'reac' )
    r1.Kf = 1.0 # 1/(mM.sec)
    r1.Kb = 1.0 # 1/sec

    # Assign parameters
    a.diffConst = diffConst
    b.diffConst = diffConst / 2.0
    b.motorConst = motorRate
    c.diffConst = 0
    d.diffConst = diffConst


    # Make solvers
    ksolve = moose.Ksolve( '/model/compartment/ksolve' )
    dsolve = moose.Dsolve( '/model/compartment/dsolve' )
    stoich = moose.Stoich( '/model/compartment/stoich' )
    stoich.compartment = compartment
    stoich.ksolve = ksolve
    stoich.dsolve = dsolve
    os.kill( PID, signal.SIGUSR1 )
    stoich.path = "/model/compartment/##"

    print((dsolve.numPools))
    assert( dsolve.numPools == 3 )
    a.vec[0].concInit = concA
    b.vec[0].concInit = concA
    c.vec[0].concInit = concA
    d.vec.concInit = concA / 5.0
    d.vec[num-1].concInit = concA

def makePlots():
    plt.ion()
    fig = plt.figure( figsize=(12,6) )
    dynamic = fig.add_subplot( 111 )

    a = moose.vec( '/model/compartment/a' )
    b = moose.vec( '/model/compartment/b' )
    c = moose.vec( '/model/compartment/c' )
    d = moose.vec( '/model/compartment/d' )
    pos = numpy.arange( 0, a.conc.size, 1 )
    aline, = dynamic.plot( pos, a.conc, label='a' )
    bline, = dynamic.plot( pos, b.conc, label='b' )
    cline, = dynamic.plot( pos, c.conc, label='c' )
    dline, = dynamic.plot( pos, d.conc, label='d' )

    plt.ylabel( 'Conc (mM)' )
    plt.xlabel( 'Cylinder voxel #' )
    plt.legend()
    timelabel = plt.text( 10, 0.8, 'time = 0.0' )
    fig.canvas.draw()
    return( fig, dynamic, timelabel, aline, bline, cline, dline )

def updatePlots( plotlist, time ):
    a = moose.vec( '/model/compartment/a' )
    b = moose.vec( '/model/compartment/b' )
    c = moose.vec( '/model/compartment/c' )
    d = moose.vec( '/model/compartment/d' )

    plotlist[2].set_text( "time = %g" % time )
    plotlist[3].set_ydata( a.conc )
    plotlist[4].set_ydata( b.conc )
    plotlist[5].set_ydata( c.conc )
    plotlist[6].set_ydata( d.conc )
    plotlist[0].canvas.draw()


[docs]def main(): """ This example illustrates how to set up a diffusion/transport model with a simple reaction-diffusion system in a tapering cylinder: | Molecule **a** diffuses with diffConst of 10e-12 m^2/s. | Molecule **b** diffuses with diffConst of 5e-12 m^2/s. | Molecule **b** also undergoes motor transport with a rate of 10e-6 m/s, | Thus it 'piles up' at the end of the cylinder. | Molecule **c** does not move: diffConst = 0.0 | Molecule **d** does not move: diffConst = 10.0e-12 but it is buffered. | Because it is buffered, it is treated as non-diffusing. All molecules other than **d** start out only in the leftmost (first) voxel, with a concentration of 1 mM. **d** is present throughout at 0.2 mM, except in the last voxel, where it is at 1.0 mM. The cylinder has a starting radius of 2 microns, and end radius of 1 micron. So when the molecule undergoing motor transport gets to the narrower end, its concentration goes up. There is a little reaction in all compartments: ``b + d <===> c`` As there is a high concentration of **d** in the last compartment, when the molecule **b** reaches the end of the cylinder, the reaction produces lots of **c**. Note that molecule **a** does not participate in this reaction. The concentrations of all molecules are displayed in an animation. """ runtime = 20.0 diffdt = 0.005 plotdt = 0.1 makeModel() # Set up clocks. The dsolver to know before assigning stoich moose.setClock( 10, diffdt ) # 10 is the standard clock for Dsolve. moose.setClock( 16, plotdt ) # 16 is the standard clock for Ksolve. a = moose.element( '/model/compartment/a' ) b = moose.element( '/model/compartment/b' ) c = moose.element( '/model/compartment/c' ) d = moose.element( '/model/compartment/d' ) moose.reinit() atot = sum( a.vec.n ) btot = sum( b.vec.n ) ctot = sum( c.vec.n ) dtot = sum( d.vec.n ) plotlist = makePlots() for t in numpy.arange( 0, runtime, plotdt ): moose.start( plotdt ) updatePlots( plotlist, t ) # moose.start( runtime ) # Run the model atot2 = sum( a.vec.n ) btot2 = sum( b.vec.n ) ctot2 = sum( c.vec.n ) dtot2 = sum( d.vec.n ) print('Ratio of initial to final total numbers of of a, b, c, d = ') print((atot2/atot, btot2/btot, ctot2/ctot, dtot2/dtot)) print(('Initial to final (b+c)=', (btot2 + ctot2) / (btot + ctot ))) print("\nHit '0' to exit") try: raw_input( ) except NameError as e: # python3 input( ) quit()
# Run the 'main' if this script is executed standalone. if __name__ == '__main__': main()