Source code for cspaceSteadyState

#########################################################################
# crossComptOscillator.py ---
#
# Filename:  crossComptOscillator.py
# Author: Upinder S. Bhalla
# Maintainer:
# Created: Oct  12 16:26:05 2014 (+0530)
# Version:
# Last-Updated: May 15 2017
#           By:
#     Update #:
# URL:
# Keywords:
# Compatibility:
#
#
# Commentary:
#
#
#
#
# Change log: Indentation clean up
#
## This program is part of 'MOOSE', the
## Messaging Object Oriented Simulation Environment.
##           Copyright (C) 2013 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 pylab
import numpy
import moose

def displayPlots():
    for x in moose.wildcardFind( '/model/graphs/#' ):
        t = numpy.arange( 0, x.vector.size, 1 ) #sec
        pylab.plot( t, x.vector, label=x.name )
    pylab.legend()
    pylab.show()

def getState( ksolve, state ):
    state.randomInit()
    moose.start( 0.1 ) # Run the model for 2 seconds.
    state.settle()
    '''
    scale = 1.0 / ( 1e-15 * 6.022e23 )
    for x in ksolve.nVec[0]:
            print x * scale,
    # print ksolve.nVec[0]
    print state.nIter, state.status, state.stateType, state.nNegEigenvalues, state.nPosEigenvalues, state.solutionStatus
     '''
    moose.start( 20.0 ) # Run model for 10 seconds, just for display


[docs]def main(): """ This example sets up the kinetic solver and steady-state finder, on a bistable model. It looks for the fixed points 100 times, as follows: - Set up the random initial condition that fits the conservation laws - Run for 2 seconds. This should not be mathematically necessary, but for obscure numerical reasons it is much more likely that the steady state solver will succeed in finding a state. - Find the fixed point - Print out the fixed point vector and various diagnostics. - Run for 10 seconds. This is completely unnecessary, and is done here just so that the resultant graph will show what kind of state has been found. After it does all this, the program runs for 100 more seconds on the last found fixed point (which turns out to be a saddle node), then is hard-switched in the script to the first attractor basin from which it runs for another 100 seconds till it settles there, and then is hard-switched yet again to the second attractor and runs for 100 seconds. Looking at the output you will see many features of note: - the first attractor (stable point) and the saddle point (unstable fixed point) are both found quite often. But the second attractor is found just once. Has a very small basin of attraction. - The values found for each of the fixed points match well with the values found by running the system to steady-state at the end. - There are a large number of failures to find a fixed point. These are found and reported in the diagnostics. They show up on the plot as cases where the 10-second runs are not flat. If you wanted to find fixed points in a production model, you would not need to do the 10-second runs, and you would need to eliminate the cases where the state-finder failed. Then you could identify the good points and keep track of how many of each were found. There is no way to guarantee that all fixed points have been found using this algorithm! You may wish to sample concentration space logarithmically rather than linearly. """ # The wildcard uses # for single level, and ## for recursive. #compartment = makeModel() moose.loadModel( '../genesis/M1719.cspace', '/model', 'ee' ) compartment = moose.element( 'model/kinetics' ) compartment.name = 'compartment' ksolve = moose.Ksolve( '/model/compartment/ksolve' ) stoich = moose.Stoich( '/model/compartment/stoich' ) stoich.compartment = compartment stoich.ksolve = ksolve #ksolve.stoich = stoich stoich.path = "/model/compartment/##" state = moose.SteadyState( '/model/compartment/state' ) moose.reinit() state.stoich = stoich #state.showMatrices() state.convergenceCriterion = 1e-7 moose.le( '/model/graphs' ) a = moose.element( '/model/compartment/a' ) b = moose.element( '/model/compartment/b' ) c = moose.element( '/model/compartment/c' ) for i in range( 0, 100 ): getState( ksolve, state ) moose.start( 100.0 ) # Run the model for 100 seconds. b = moose.element( '/model/compartment/b' ) c = moose.element( '/model/compartment/c' ) # move most molecules over to b b.conc = b.conc + c.conc * 0.95 c.conc = c.conc * 0.05 moose.start( 100.0 ) # Run the model for 100 seconds. # move most molecules back to a c.conc = c.conc + b.conc * 0.95 b.conc = b.conc * 0.05 moose.start( 100.0 ) # Run the model for 100 seconds. # Iterate through all plots, dump their contents to data.plot. displayPlots() quit()
# Run the 'main' if this script is executed standalone. if __name__ == '__main__': main()