import math
import pylab
import numpy
import moose
runtime = 50.0
def makeModel():
# create container for model
model = moose.Neutral( 'model' )
harmonic = moose.CubeMesh( '/model/harmonic' )
harmonic.volume = 1e-15
lotka = moose.CubeMesh( '/model/lotka' )
lotka.volume = 1e-15
# create molecules and reactions
x = moose.Pool( '/model/lotka/x' )
y = moose.Pool( '/model/lotka/y' )
z = moose.BufPool( '/model/lotka/z' ) # Dummy molecule.
xreac = moose.Reac( '/model/lotka/xreac' )
yreac = moose.Reac( '/model/lotka/yreac' )
xrate = moose.Function( '/model/lotka/xreac/func' )
yrate = moose.Function( '/model/lotka/yreac/func' )
# Parameters
alpha = 1.0
beta = 1.0
gamma = 1.0
delta = 1.0
k = 1.0
x.nInit = 2.0
y.nInit = 1.0
z.nInit = 0.0
xrate.x.num = 1
yrate.x.num = 1
xrate.expr = "x0 * " + str( beta ) + " - " + str( alpha )
yrate.expr = str( gamma ) + " - x0 * " + str( delta )
xreac.Kf = k
yreac.Kf = k
xreac.Kb = 0
yreac.Kb = 0
# connect them up for reactions
moose.connect( y, 'nOut', xrate.x[0], 'input' )
moose.connect( x, 'nOut', yrate.x[0], 'input' )
moose.connect( xrate, 'valueOut', xreac, 'setNumKf' )
moose.connect( yrate, 'valueOut', yreac, 'setNumKf' )
moose.connect( xreac, 'sub', x, 'reac' )
moose.connect( xreac, 'prd', z, 'reac' )
moose.connect( yreac, 'sub', y, 'reac' )
moose.connect( yreac, 'prd', z, 'reac' )
# Create the output tables
graphs = moose.Neutral( '/model/graphs' )
xplot = moose.Table2 ( '/model/graphs/x' )
yplot = moose.Table2 ( '/model/graphs/y' )
# connect up the tables
moose.connect( xplot, 'requestOut', x, 'getN' );
moose.connect( yplot, 'requestOut', y, 'getN' );
[docs]def main():
"""
The funcReacLotkaVolterra example shows how to use function objects
as part of differential equation systems in the framework of the MOOSE
kinetic solvers. Here the system is set up explicitly using the
scripting, in normal use one would expect to use SBML.
In this example we set up a Lotka-Volterra system. The equations
are readily expressed as a pair of reactions each of whose rate is
governed by a function::
x' = x( alpha - beta.y )
y' = -y( gamma - delta.x )
This translates into two reactions::
x ---> z Kf = beta.y - alpha
y ---> z Kf = gamma - delta.x
Here z is a dummy molecule whose concentration is buffered to zero.
The model first runs using default Exponential Euler integration.
This is not particularly accurate even with a small timestep.
The model is then converted to use the deterministic Kinetic solver
Ksolve. This is accurate and faster.\n
Note that we cannot use the stochastic GSSA solver for this system, it
cannot handle a reaction term whose rate keeps changing.
"""
makeModel()
for i in range( 11, 18 ):
moose.setClock( i, 0.001 )
moose.setClock( 18, 0.1 )
moose.reinit()
moose.start( runtime ) # Run the model
# Iterate through all plots, dump their contents to data.plot.
for x in moose.wildcardFind( '/model/graphs/#' ):
#x.xplot( 'scriptKineticModel.plot', x.name )
t = numpy.arange( 0, x.vector.size, 1 ) * x.dt # sec
pylab.plot( t, x.vector, label=x.name )
pylab.ylim( 0, 2.5 )
pylab.title( "Exponential Euler solution. Note slight error buildup" )
pylab.legend()
pylab.figure()
compt = moose.element( '/model/lotka' )
ksolve = moose.Ksolve( '/model/lotka/ksolve' )
stoich = moose.Stoich( '/model/lotka/stoich' )
stoich.compartment = compt
stoich.ksolve = ksolve
stoich.path = '/model/lotka/##'
moose.reinit()
moose.start( runtime ) # Run the model
for i in range( 11, 18 ):
moose.setClock( i, 0.1 )
for x in moose.wildcardFind( '/model/graphs/#' ):
t = numpy.arange( 0, x.vector.size, 1 ) * x.dt # sec
pylab.plot( t, x.vector, label=x.name )
pylab.ylim( 0, 2.5 )
pylab.title( "Runge-Kutta solution." )
pylab.legend()
pylab.show()
quit()
# Run the 'main' if this script is executed standalone.
if __name__ == '__main__':
main()