Source code for funcReacLotkaVolterra

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()