Source code for TuringOneDim

#########################################################################
## 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 math
import numpy
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import moose

[docs]def makeModel(): """ This example illustrates how to set up a oscillatory Turing pattern in 1-D using reaction diffusion calculations. Reaction system is:: s ---a---> a // s goes to a, catalyzed by a. s ---a---> b // s goes to b, catalyzed by a. a ---b---> s // a goes to s, catalyzed by b. b -------> s // b is degraded irreversibly to s. in sum, **a** has a positive feedback onto itself and also forms **b**. **b** has a negative feedback onto **a**. Finally, the diffusion constant for **a** is 1/10 that of **b**. This chemical system is present in a 1-dimensional (cylindrical) compartment. The entire reaction-diffusion system is set up within the script. """ # create container for model r0 = 1e-6 # m r1 = 1e-6 # m num = 100 diffLength = 1e-6 # m len = num * diffLength # m diffConst = 5e-12 # m^2/sec motorRate = 1e-6 # m/sec concA = 1 # millimolar dt4 = 0.02 # for the diffusion dt5 = 0.2 # for the reaction 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' ) s = moose.Pool( '/model/compartment/s' ) e1 = moose.MMenz( '/model/compartment/e1' ) e2 = moose.MMenz( '/model/compartment/e2' ) e3 = moose.MMenz( '/model/compartment/e3' ) r1 = moose.Reac( '/model/compartment/r1' ) moose.connect( e1, 'sub', s, 'reac' ) moose.connect( e1, 'prd', a, 'reac' ) moose.connect( a, 'nOut', e1, 'enzDest' ) e1.Km = 1 e1.kcat = 1 moose.connect( e2, 'sub', s, 'reac' ) moose.connect( e2, 'prd', b, 'reac' ) moose.connect( a, 'nOut', e2, 'enzDest' ) e2.Km = 1 e2.kcat = 0.5 moose.connect( e3, 'sub', a, 'reac' ) moose.connect( e3, 'prd', s, 'reac' ) moose.connect( b, 'nOut', e3, 'enzDest' ) e3.Km = 0.1 e3.kcat = 1 moose.connect( r1, 'sub', b, 'reac' ) moose.connect( r1, 'prd', s, 'reac' ) r1.Kf = 0.3 # 1/sec r1.Kb = 0 # 1/sec # Assign parameters a.diffConst = diffConst/10 b.diffConst = diffConst s.diffConst = 0 # Make solvers ksolve = moose.Ksolve( '/model/compartment/ksolve' ) dsolve = moose.Dsolve( '/model/dsolve' ) # Set up clocks. The dsolver to know before assigning stoich moose.setClock( 4, dt4 ) moose.setClock( 5, dt5 ) moose.useClock( 4, '/model/dsolve', 'process' ) # Ksolve must be scheduled after dsolve. moose.useClock( 5, '/model/compartment/ksolve', 'process' ) stoich = moose.Stoich( '/model/compartment/stoich' ) stoich.compartment = compartment stoich.ksolve = ksolve stoich.dsolve = dsolve stoich.path = "/model/compartment/##" assert( dsolve.numPools == 3 ) a.vec.concInit = [0.1]*num a.vec[0].concInit *= 1.2 # slight perturbation at one end. b.vec.concInit = [0.1]*num s.vec.concInit = [1]*num
def displayPlots(): a = moose.element( '/model/compartment/a' ) b = moose.element( '/model/compartment/b' ) pos = numpy.arange( 0, a.vec.conc.size, 1 ) pylab.plot( pos, a.vec.conc, label='a' ) pylab.plot( pos, b.vec.conc, label='b' ) pylab.legend() pylab.show() def main(): runtime = 400 displayInterval = 2 makeModel() dsolve = moose.element( '/model/dsolve' ) moose.reinit() #moose.start( runtime ) # Run the model for 10 seconds. a = moose.element( '/model/compartment/a' ) b = moose.element( '/model/compartment/b' ) s = moose.element( '/model/compartment/s' ) img = mpimg.imread( 'turingPatternTut.png' ) #imgplot = plt.imshow( img ) #plt.show() plt.ion() fig = plt.figure( figsize=(12,10) ) png = fig.add_subplot(211) imgplot = plt.imshow( img ) ax = fig.add_subplot(212) ax.set_ylim( 0, 0.5 ) plt.ylabel( 'Conc (mM)' ) plt.xlabel( 'Position along cylinder (microns)' ) pos = numpy.arange( 0, a.vec.conc.size, 1 ) line1, = ax.plot( pos, a.vec.conc, label='a' ) line2, = ax.plot( pos, b.vec.conc, label='b' ) timeLabel = plt.text(60, 0.4, 'time = 0') plt.legend() fig.canvas.draw() for t in range( displayInterval, runtime, displayInterval ): moose.start( displayInterval ) line1.set_ydata( a.vec.conc ) line2.set_ydata( b.vec.conc ) timeLabel.set_text( "time = %d" % t ) fig.canvas.draw() print( "Hit 'enter' to exit" ) try: raw_input( ) except NameError as e: # python3 input( ) # Run the 'main' if this script is executed standalone. if __name__ == '__main__': main()