Skip to content

Instantly share code, notes, and snippets.

@dilawar

dilawar/M1719.cspace

Forked from nikalras/M1719.cspace
Created August 23, 2014 04:56
Show Gist options
  • Save dilawar/3bbe9879fe35eabcfe23 to your computer and use it in GitHub Desktop.
Save dilawar/3bbe9879fe35eabcfe23 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
cspaceSteadyState.py
#########################################################################
## 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.
#########################################################################
# 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 makes it 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.
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
def main():
# Schedule the whole lot
moose.setClock( 4, 0.1 ) # for the computational objects
moose.setClock( 5, 0.1 ) # clock for the solver
moose.setClock( 8, 1.0 ) # for the plots
# The wildcard uses # for single level, and ## for recursive.
#compartment = makeModel()
moose.loadModel( '../Genesis_files/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.useClock( 5, '/model/compartment/ksolve', 'process' )
moose.useClock( 8, '/model/graphs/#', 'process' )
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()
Raw
M1719.cspace
M1719: |Labc|Ldce|Ebze| 0.20 0.30 0.32 0.32 0.32 0.32 0.263296394963 0.111387109786 0.470177600114 0.029771700191 0.173209325064 3.28129207574
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment