# -*- coding: utf-8 -*-
# crossComptSimpleReacGSSA.py
#
# Filename:crossComptSimpleReacGSSA.py
# Author: Upinder S. Bhalla
# Maintainer:
# Created: Oct 12 16:26:05 2014 (+0530)
# Version:
# Last-Updated: May 16 2017
# By: Upinder S. Bhalla
# Update #:
# URL:
# Keywords:
# Compatibility:
#
# Commentary:
# __BROKEN__
#
# Change log:
# 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 makeModel():
# create container for model
model = moose.Neutral( 'model' )
compt0 = moose.CubeMesh( '/model/compt0' )
compt0.volume = 1e-18
compt1 = moose.CubeMesh( '/model/compt1' )
compt1.volume = 1e-19
compt2 = moose.CubeMesh( '/model/compt2' )
compt2.volume = 1e-20
# Position containers so that they abut each other, with
# compt1 in the middle.
side = compt1.dy
compt0.y1 += side
compt0.y0 += side
compt2.x1 += side
compt2.x0 += side
print(('Volumes = ', compt0.volume, compt1.volume, compt2.volume))
# create molecules and reactions
a = moose.Pool( '/model/compt0/a' )
b = moose.Pool( '/model/compt1/b' )
c = moose.Pool( '/model/compt2/c' )
reac0 = moose.Reac( '/model/compt1/reac0' )
reac1 = moose.Reac( '/model/compt1/reac1' )
# connect them up for reactions
moose.connect( reac0, 'sub', a, 'reac' )
moose.connect( reac0, 'prd', b, 'reac' )
moose.connect( reac1, 'sub', b, 'reac' )
moose.connect( reac1, 'prd', c, 'reac' )
# Assign parameters
a.concInit = 1
b.concInit = 12.1
c.concInit = 1
reac0.Kf = 0.1
reac0.Kb = 0.1
reac1.Kf = 0.1
reac1.Kb = 0.1
# Create the output tables
graphs = moose.Neutral( '/model/graphs' )
outputA = moose.Table2 ( '/model/graphs/concA' )
outputB = moose.Table2 ( '/model/graphs/concB' )
outputC = moose.Table2 ( '/model/graphs/concC' )
# connect up the tables
moose.connect( outputA, 'requestOut', a, 'getConc' );
moose.connect( outputB, 'requestOut', b, 'getConc' );
moose.connect( outputC, 'requestOut', c, 'getConc' );
# Build the solvers. No need for diffusion in this version.
ksolve0 = moose.Gsolve( '/model/compt0/ksolve0' )
ksolve1 = moose.Gsolve( '/model/compt1/ksolve1' )
ksolve2 = moose.Gsolve( '/model/compt2/ksolve2' )
stoich0 = moose.Stoich( '/model/compt0/stoich0' )
stoich1 = moose.Stoich( '/model/compt1/stoich1' )
stoich2 = moose.Stoich( '/model/compt2/stoich2' )
# Configure solvers
stoich0.compartment = compt0
stoich1.compartment = compt1
stoich2.compartment = compt2
stoich0.ksolve = ksolve0
stoich1.ksolve = ksolve1
stoich2.ksolve = ksolve2
stoich0.path = '/model/compt0/#'
stoich1.path = '/model/compt1/#'
stoich2.path = '/model/compt2/#'
stoich1.buildXreacs( stoich0 )
stoich1.buildXreacs( stoich2 )
stoich0.filterXreacs()
stoich1.filterXreacs()
stoich2.filterXreacs()
[docs]def main():
"""
This example illustrates a simple cross compartment reaction::
a <===> b <===> c
Here each molecule is in a different compartment.
The initial conditions are such that the end conc on all compartments
should be 2.0.
The time course depends on which compartment the Reac object is
embedded in.
The cleanest thing numerically and also conceptually is to have both
reactions in the same compartment, in this case the middle one
(**compt1**).
The initial conditions have a lot of **B**. The equilibrium with
**C** is fast and so **C** shoots up and passes **B**, peaking at
about (2.5,9). This is also just
about the crossover point.
**A** starts low and slowly climbs up to equilibrate.
If we put **reac0** in **compt0** and **reac1** in **compt1**,
it behaves the same
qualitiatively but now the peak is at around (1, 5.2)
This configuration of reactions makes sense from the viewpoint of
having the
reactions always in the compartment with the smaller volume, which is
important if we need to have junctions where many small voxels talk to
one big voxel in another compartment.
Note that putting the reacs in other compartments doesn't work and in
some combinations (e.g., **reac0** in **compt0** and **reac1** in
**compt2**) give
numerical instability.
"""
simdt = 0.1
plotdt = 0.1
runtime = 100.0
makeModel()
# MOOSE autoschedules everything.
moose.reinit()
moose.start( runtime ) # Run the model for 100 seconds.
print("All concs should converge to 2.0 even though vols differ:")
for x in moose.wildcardFind( '/model/compt#/#[ISA=PoolBase]' ):
print((x.name, x.conc))
# Iterate through all plots, dump their contents to data.plot.
for x in moose.wildcardFind( '/model/graphs/conc#' ):
t = numpy.linspace( 0, runtime, x.vector.size ) # sec
pylab.plot( t, x.vector, label=x.name )
pylab.legend()
pylab.show()
quit()
# Run the 'main' if this script is executed standalone.
if __name__ == '__main__':
main()