Source code for chemDoseResponse

#########################################################################
## 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

[docs]def main(): """ This example builds a dose-response of a bistable model of a chemical system. It uses the kinetic solver *Ksolve* and the steady-state finder *SteadyState*. The model is set up within the script. The basic approach is to increment the control variable, **a** in this case, while monitoring **b**. The algorithm marches through a series of values of the buffered pool **a** and measures resultant values of pool **b**. At each cycle the algorithm calls the steady-state finder. Since **a** is incremented only a small amount on each step, each new steady state is (usually) quite close to the previous one. The exception is when there is a state transition. Here we plot three dose-response curves to illustrate the bistable nature of the system. On the upward going curve in blue, **a** starts low. Here, **b** follows the low arm of the curve and then jumps up to the high value at roughly *log( [a] ) = -0.55*. On the downward going curve in green, **b** follows the high arm of the curve forming a nice hysteretic loop. Eventually **b** has to fall to the low state at about *log( [a] ) = -0.83* Through nasty concentration manipulations, we find the third arm of the curve, which tracks the unstable fixed point. This is in red. We find this arm by setting an initial point close to the unstable fixed point, which the steady-state finder duly locates. We then follow a dose-response curve as with the other arms of the curve. Note that the steady-state solver doesn't always succeed in finding a good solution, despite moving only in small steps. Nevertheless the resultant curves are smooth because it gives up pretty close to the correct value, simply because the successive points are close together. Overall, the system is pretty robust despite the core root-finder computations in GSL being temperamental. In doing a production dose-response series you may wish to sample concentration space logarithmically rather than linearly. """ compartment = makeModel() ksolve = moose.Ksolve( '/model/compartment/ksolve' ) stoich = moose.Stoich( '/model/compartment/stoich' ) stoich.compartment = compartment stoich.ksolve = ksolve stoich.path = "/model/compartment/##" state = moose.SteadyState( '/model/compartment/state' ) moose.reinit() state.stoich = stoich state.convergenceCriterion = 1e-6 moose.seed( 111 ) # Used when generating the samples in state space b = moose.element( '/model/compartment/b' ) a = moose.element( '/model/compartment/a' ) c = moose.element( '/model/compartment/c' ) a.concInit = 0.1 deltaA = 0.002 num = 150 avec = [] bvec = [] moose.reinit() # Now go up. for i in range( 0, num ): moose.start( 1.0 ) # Run the model for 1 seconds. state.settle() # This function finds the steady states. avec.append( a.conc ) bvec.append( b.conc ) a.concInit += deltaA #print i, a.conc, b.conc pylab.plot( numpy.log10( avec ), numpy.log10( bvec ), label='b vs a up' ) # Now go down. avec = [] bvec = [] for i in range( 0, num ): moose.start( 1.0 ) # Run the model for 1 seconds. state.settle() # This function finds the steady states. avec.append( a.conc ) bvec.append( b.conc ) a.concInit -= deltaA #print i, a.conc, b.conc pylab.plot( numpy.log10( avec ), numpy.log10( bvec ), label='b vs a down' ) # Now aim for the middle. We do this by judiciously choosing a # start point that should be closer to the unstable fixed point. avec = [] bvec = [] a.concInit = 0.28 b.conc = 0.15 for i in range( 0, 65 ): moose.start( 1.0 ) # Run the model for 1 seconds. state.settle() # This function finds the steady states. avec.append( a.conc ) bvec.append( b.conc ) a.concInit -= deltaA #print i, a.conc, b.conc pylab.plot( numpy.log10( avec ), numpy.log10( bvec ), label='b vs a mid' ) pylab.ylim( [-1.7, 1.2] ) pylab.legend() pylab.show() quit()
[docs]def makeModel(): """ This function creates a bistable reaction system using explicit MOOSE calls rather than load from a file. The reaction is:: a ---b---> 2b # b catalyzes a to form more of b. 2b ---c---> a # c catalyzes b to form a. a <======> 2b # a interconverts to b. """ # create container for model model = moose.Neutral( 'model' ) compartment = moose.CubeMesh( '/model/compartment' ) compartment.volume = 1e-15 # the mesh is created automatically by the compartment mesh = moose.element( '/model/compartment/mesh' ) # create molecules and reactions a = moose.BufPool( '/model/compartment/a' ) b = moose.Pool( '/model/compartment/b' ) c = moose.Pool( '/model/compartment/c' ) enz1 = moose.Enz( '/model/compartment/b/enz1' ) enz2 = moose.Enz( '/model/compartment/c/enz2' ) cplx1 = moose.Pool( '/model/compartment/b/enz1/cplx' ) cplx2 = moose.Pool( '/model/compartment/c/enz2/cplx' ) reac = moose.Reac( '/model/compartment/reac' ) # connect them up for reactions moose.connect( enz1, 'sub', a, 'reac' ) moose.connect( enz1, 'prd', b, 'reac' ) moose.connect( enz1, 'prd', b, 'reac' ) # Note 2 molecules of b. moose.connect( enz1, 'enz', b, 'reac' ) moose.connect( enz1, 'cplx', cplx1, 'reac' ) moose.connect( enz2, 'sub', b, 'reac' ) moose.connect( enz2, 'sub', b, 'reac' ) # Note 2 molecules of b. moose.connect( enz2, 'prd', a, 'reac' ) moose.connect( enz2, 'enz', c, 'reac' ) moose.connect( enz2, 'cplx', cplx2, 'reac' ) moose.connect( reac, 'sub', a, 'reac' ) moose.connect( reac, 'prd', b, 'reac' ) moose.connect( reac, 'prd', b, 'reac' ) # Note 2 order in b. # Assign parameters a.concInit = 1 b.concInit = 0 c.concInit = 0.01 enz1.kcat = 0.4 enz1.Km = 4 enz2.kcat = 0.6 enz2.Km = 0.01 reac.Kf = 0.001 reac.Kb = 0.01 return compartment
# Run the 'main' if this script is executed standalone. if __name__ == '__main__': main()