mcell/mdl/testsuite/reactions/01-volume_highconc.mdl

/****************************************************************************
 * Test 01: All volume rxn types, simple test of equilibrium and reaction
 *          rates.
 *
 *    A quantitative test of the basic volume-volume reaction types (uni, bi,
 *    tri).  Note that the bi-directional pathways are in here twice.  This is
 *    so that they can be tested via both the non-trimolecular and trimolecular
 *    code pathways.  The latter is accomplished by adding a bogus trimolecular
 *    reaction involving the relevant molecules.
 *
 * Author: Jed Wing <[email protected]>
 * Date:   2008-11-07
 ****************************************************************************/

basename = "01-volume_highconc"
countdir = "dat/" & basename & "/"
vizdir   = "./viz_dat/" & basename & "/"

iterations = 10000
timestep = 1e-6
boxdims = 0.25
voldiffconst = 2.0e-6

num_u = 1000
num_b = 1000
num_t = 1000

// Set reverse rates to expect approx 1 rxn every 10 timesteps
rr_u = 2.0e2
rr_b = 2.0e2
rr_t = 2.0e2

// Compute forward rates to give equilibrium of 50% reacted
boxvol = 8*boxdims*boxdims*boxdims
c_b = (num_b / 2. / boxvol) * 1e15 / 6.022e23
c_t = (num_t / 2. / boxvol) * 1e15 / 6.022e23
rf_u = rr_u
rf_b = rr_b / c_b
rf_t = rr_t / c_t / c_t

PARTITION_X = [ [ -0.3 TO 0.3 STEP 0.025 ] ]
PARTITION_Y = [ [ -0.3 TO 0.3 STEP 0.025 ] ]
PARTITION_Z = [ [ -0.3 TO 0.3 STEP 0.025 ] ]

ITERATIONS = iterations
TIME_STEP = timestep
TIME_STEP_MAX = 1.0

INCLUDE_FILE = "volmols.mdl"

DEFINE_REACTIONS {
  V_u                     -> I_u   [rf_u] : fwd_u    /* unimol */
  V_b1  + V_b2            -> I_b   [rf_b] : fwd_b    /* bimol  */
  V_bt1 + V_bt2           -> I_bt  [rf_b] : fwd_bt   /* bimol (via trimol code) */
  V_b3  + V_b4to          -> I_b2  [rf_b] : fwd_b2   /* bimol (1 target only) */
  V_bt3 + V_bt4to         -> I_bt2 [rf_b] : fwd_bt2  /* bimol (1 target only, via trimol code) */
  V_t1  + V_t2 + V_t3     -> I_t   [rf_t] : fwd_t    /* trimol */
  V_t4  + V_t5to + V_t6   -> I_t2  [rf_t] : fwd_t2   /* trimol (1 target only) */
  V_t7  + V_t8to + V_t9to -> I_t3  [rf_t/3] : fwd_t3   /* trimol (2 target only) */

  I_u   -> V_u                    [rr_u] : rev_u
  I_b   -> V_b1 + V_b2            [rr_b] : rev_b
  I_bt  -> V_bt1 + V_bt2          [rr_b] : rev_bt
  I_b2  -> V_b3 + V_b4to          [rr_b] : rev_b2
  I_bt2 -> V_bt3 + V_bt4to        [rr_b] : rev_bt2
  I_t   -> V_t1 + V_t2 + V_t3     [rr_t] : rev_t
  I_t2  -> V_t4 + V_t5to + V_t6   [rr_t] : rev_t2
  I_t3  -> V_t7 + V_t8to + V_t9to [rr_t/3] : rev_t3

  /* Dummy reactions to trick mcell into using trimol. pathway for bimol rxn. */
  V_bt1 + V_a + V_a -> V_a   [1e2]
  V_bt2 + V_a + V_a -> V_a   [1e2]
  V_bt3 + V_a + V_a -> V_a   [1e2]
}

main_box BOX
{
    CORNERS = [-boxdims, -boxdims, -boxdims], [boxdims, boxdims, boxdims]
}

INSTANTIATE world OBJECT
{
    main_box OBJECT main_box {}

    release_site_1r RELEASE_SITE     { SHAPE = world.main_box   MOLECULE = V_u     NUMBER_TO_RELEASE = num_u / 2 }
    release_site_1p RELEASE_SITE     { SHAPE = world.main_box   MOLECULE = I_u     NUMBER_TO_RELEASE = num_u / 2 }

    /* V_b1 + V_b2 <-> I_b */
    release_site_2r1 RELEASE_SITE    { SHAPE = world.main_box   MOLECULE = V_b1    NUMBER_TO_RELEASE = num_b / 2 }
    release_site_2r2 RELEASE_SITE    { SHAPE = world.main_box   MOLECULE = V_b2    NUMBER_TO_RELEASE = num_b / 2 }
    release_site_2p RELEASE_SITE     { SHAPE = world.main_box   MOLECULE = I_b     NUMBER_TO_RELEASE = num_b / 2 }

    /* V_bt1 + V_bt2 <-> I_bt */
    release_site_2tr1 RELEASE_SITE   { SHAPE = world.main_box   MOLECULE = V_bt1   NUMBER_TO_RELEASE = num_b / 2 }
    release_site_2tr2 RELEASE_SITE   { SHAPE = world.main_box   MOLECULE = V_bt2   NUMBER_TO_RELEASE = num_b / 2 }
    release_site_2tp RELEASE_SITE    { SHAPE = world.main_box   MOLECULE = I_bt    NUMBER_TO_RELEASE = num_b / 2 }

    /* V_b3 + V_b4 <-> I_b2 */
    release_site_2r1t RELEASE_SITE   { SHAPE = world.main_box   MOLECULE = V_b3    NUMBER_TO_RELEASE = num_b / 2 }
    release_site_2r2t RELEASE_SITE   { SHAPE = world.main_box   MOLECULE = V_b4to  NUMBER_TO_RELEASE = num_b / 2 }
    release_site_2pt RELEASE_SITE    { SHAPE = world.main_box   MOLECULE = I_b2    NUMBER_TO_RELEASE = num_b / 2 }

    /* V_bt3 + V_bt4 <-> I_bt2 */
    release_site_2tr1t RELEASE_SITE  { SHAPE = world.main_box   MOLECULE = V_bt3   NUMBER_TO_RELEASE = num_b / 2 }
    release_site_2tr2t RELEASE_SITE  { SHAPE = world.main_box   MOLECULE = V_bt4to NUMBER_TO_RELEASE = num_b / 2 }
    release_site_2tpt RELEASE_SITE   { SHAPE = world.main_box   MOLECULE = I_bt2   NUMBER_TO_RELEASE = num_b / 2 }

    /* V_t1 + V_t2 + V_t3 <-> I_t */
    release_site_3r1 RELEASE_SITE    { SHAPE = world.main_box   MOLECULE = V_t1    NUMBER_TO_RELEASE = num_t / 2 }
    release_site_3r2 RELEASE_SITE    { SHAPE = world.main_box   MOLECULE = V_t2    NUMBER_TO_RELEASE = num_t / 2 }
    release_site_3r3 RELEASE_SITE    { SHAPE = world.main_box   MOLECULE = V_t3    NUMBER_TO_RELEASE = num_t / 2 }
    release_site_3p RELEASE_SITE     { SHAPE = world.main_box   MOLECULE = I_t     NUMBER_TO_RELEASE = num_t / 2 }

    /* V_t4 + V_t5to + V_t6 <-> I_t2 */
    release_site_3r1t RELEASE_SITE   { SHAPE = world.main_box   MOLECULE = V_t4    NUMBER_TO_RELEASE = num_t / 2 }
    release_site_3r2t RELEASE_SITE   { SHAPE = world.main_box   MOLECULE = V_t5to  NUMBER_TO_RELEASE = num_t / 2 }
    release_site_3r3t RELEASE_SITE   { SHAPE = world.main_box   MOLECULE = V_t6    NUMBER_TO_RELEASE = num_t / 2 }
    release_site_3pt RELEASE_SITE    { SHAPE = world.main_box   MOLECULE = I_t2    NUMBER_TO_RELEASE = num_t / 2 }

    /* V_t7 + V_t8to + V_t9to <-> I_t3 */
    release_site_3r1tt RELEASE_SITE  { SHAPE = world.main_box   MOLECULE = V_t7    NUMBER_TO_RELEASE = num_t / 2 }
    release_site_3r2tt RELEASE_SITE  { SHAPE = world.main_box   MOLECULE = V_t8to  NUMBER_TO_RELEASE = num_t / 2 }
    release_site_3r3tt RELEASE_SITE  { SHAPE = world.main_box   MOLECULE = V_t9to  NUMBER_TO_RELEASE = num_t / 2 }
    release_site_3ptt RELEASE_SITE   { SHAPE = world.main_box   MOLECULE = I_t3    NUMBER_TO_RELEASE = num_t / 2 }
}

REACTION_DATA_OUTPUT {
  OUTPUT_BUFFER_SIZE = 50
  STEP = timestep*10
  HEADER = "# "
  {
    COUNT[V_u,world.main_box]     : "V_u",
    COUNT[V_b1,world.main_box]    : "V_b1",
    COUNT[V_b2,world.main_box]    : "V_b2",
    COUNT[V_bt1,world.main_box]   : "V_bt1",
    COUNT[V_bt2,world.main_box]   : "V_bt2",
    COUNT[V_b3,world.main_box]    : "V_b3",
    COUNT[V_b4to,world.main_box]  : "V_b4",
    COUNT[V_bt3,world.main_box]   : "V_bt3",
    COUNT[V_bt4to,world.main_box] : "V_bt4",
    COUNT[V_t1,world.main_box]    : "V_t1",
    COUNT[V_t2,world.main_box]    : "V_t2",
    COUNT[V_t3,world.main_box]    : "V_t3",
    COUNT[V_t4,world.main_box]    : "V_t4",
    COUNT[V_t5to,world.main_box]  : "V_t5",
    COUNT[V_t6,world.main_box]    : "V_t6",
    COUNT[V_t7,world.main_box]    : "V_t7",
    COUNT[V_t8to,world.main_box]  : "V_t8",
    COUNT[V_t9to,world.main_box]  : "V_t9",
    COUNT[I_u,world.main_box]     : "I_u",
    COUNT[I_b,world.main_box]     : "I_b",
    COUNT[I_bt,world.main_box]    : "I_bt",
    COUNT[I_b2,world.main_box]    : "I_b2",
    COUNT[I_bt2,world.main_box]   : "I_bt2",
    COUNT[I_t,world.main_box]     : "I_t",
    COUNT[I_t2,world.main_box]    : "I_t2",
    COUNT[I_t3,world.main_box]    : "I_t3"
  }     => countdir & "V_out.dat"
  {
    COUNT[fwd_u,WORLD]         : "fwd_u",
    COUNT[fwd_b,WORLD]         : "fwd_b",
    COUNT[fwd_bt,WORLD]        : "fwd_bt",
    COUNT[fwd_b2,WORLD]        : "fwd_b2",
    COUNT[fwd_bt2,WORLD]       : "fwd_bt2",
    COUNT[fwd_t,WORLD]         : "fwd_t",
    COUNT[fwd_t2,WORLD]        : "fwd_t2",
    COUNT[fwd_t3,WORLD]        : "fwd_t3",
    COUNT[rev_u,WORLD]         : "rev_u",
    COUNT[rev_b,WORLD]         : "rev_b",
    COUNT[rev_bt,WORLD]        : "rev_bt",
    COUNT[rev_b2,WORLD]        : "rev_b2",
    COUNT[rev_bt2,WORLD]       : "rev_bt2",
    COUNT[rev_t,WORLD]         : "rev_t",
    COUNT[rev_t2,WORLD]        : "rev_t2",
    COUNT[rev_t3,WORLD]        : "rev_t3"
  }     => countdir & "rxn_out.dat"
}