Skip to content

Reaction

Tim Vaughan edited this page Mar 16, 2015 · 18 revisions

The <reaction> element is defined in the following way:

    <reaction spec='Reaction' reactionName='STRING' rate='RATESPEC'>
      <!-- Reaction String -->
      <!-- Optional Predicate -->
      <!-- Optional RateMultiplier -->
    </reaction>

The optional reactionName attribute can be used to specify a unique identifying name of the reaction.

The rate attribute specifies one or more rates at which the reaction will occur, per individual reactant configuration. (If a reaction requiring two individual members of a population of size 10 has a rate of 0.1, the overall propensity of that reaction will be 0.1*10*(10-1)=9 per unit time.) The value of this attribute is a RATESPEC which may either be a single double-precision number, or a comma-delimited list of RATE:TIME pairs:

RATE1:TIME1, RATE2:TIME2, RATE3:TIME3, ...

where RATE1 is the rate between times TIME1 and TIME2, RATE2 applies between TIME2 and TIME3 and so-on. The rate is assumed to be zero for the period between time zero (the start of the simulation) and TIME1. Additionally, one can omit :TIME1 in which case RATE1 is assumed to apply at time zero.

The value of the element contains the "reaction string" which is a human-readable description of the reactant and product individuals involved in the reaction. It is an ASCII representation of a simple chemical reaction, and takes the form

REACTANT [ + REACTANT [ + ... ] ] -> PRODUCT [ + PRODUCT [ + ... ] ]

where REACTANT and PRODUCT specify individuals involved in the reaction. For the sake of the kinds of birth-death processes we deal with here, only the population that these individuals belongs to matters. This is specified using the following format:

TYPE[loc]

where TYPE is the unique name given to the population type and [loc] is the specific location of the population within that type, given as a list of integers. For example, in the case of a population type named "X" defined using the attribute dim="2 2" a sensible REACTANT or PRODUCT would be "X[1,0]".

In the case that the population type has only a single population, the square brackets can be omitted. Also, if two or more members of a population exists on the same side of the ->, they can be collapsed together by prefixing the population type name with an integer representing the number of times that the individual occurs. Finally, if no individuals are present on one side of the ->, the place-holder 0 must be used there instead. For example, a linear decay process could be represented by the reaction X -> 0.

The optional <predicate> and <rateMultiplier> elements are discussed below in the section "Specifying multiple reactions".

Inheritance-specific considerations

For InheritanceTrajectory or InheritanceEnsemble simulations, there is some additional complexity due to the fact that inheritance relationships between reactants and products need to be specified. To do this, we use the following augmented form of the reaction specification string:

REACTANT:ID [ + REACTANT:ID [ + ... ] ] -> PRODUCT:ID [ + PRODUCT:ID [ + ... ] ]

Each ID must be a positive integer and has the following effect: all products associated with each unique ID on the right-hand-side of the reaction are considered direct descendents of all of the reactants having that same ID on the left-hand side of the reaction.

The following two reaction strings, which correspond to the same process 2X -> X under population-only simulation types, generate very distinct modifications to an inheritance tree. Firstly

X:1 + X:2 -> X:1

results in a continuation of the lineage belonging to the first reactant individual and a termination of the lineage belonging to the second. Alternatively,

X:1 + X:1 -> X:1

results in a joining of the two lineages belonging to the reactants.

If the :IDs are omitted from the individual specifiers, the following "greedy" algorithm is used to assign children to parents:

  1. Select the first reactant individual, naming it R.
  2. Assign all product individuals having the same population type as R to be children of R.
  3. If unprocessed reactants remain, set R to be the next reactant and go back to 2.
  4. End.

If included, the :IDs have no effect on non-inheritance simulations.

Specifying multiple reactions

Occasionally it is necessary to specify a large number of reactions between identical reagent PopulationTypes but which differ with respect to the actual populations involved. A simple example is a model which describes migration between all pairs of 3 locations (islands). The reaction-based description MASTER employs requires enumerating 6 distinct reactions:

X[0] -> X[1]
X[0] -> X[2]
X[1] -> X[0]
X[1] -> X[2]
X[2] -> X[0]
X[2] -> X[1]

If instead we had 5 islands, we would need to specify 20 reactions. At 10 islands, we'd need 90. Specifying these reactions would quickly become tedious, especially as each would need its own <reaction> element!

To avoid this tedium, MASTER allows you to replace location vector indices with named variables such as "i" or "var". When encountered in a reaction specification, these variables are allowed to range over all of the values allowed by the population type definitions. Thus, the following single line:

X[i] -> X[j]

produces reactions for every combination of i and j allowed for the population type X.

Predicates

While use of these index variables dramatically shortens the description of some models, it is insufficient on its own as it often produces reactions that are not wanted. In the migration example above, for instance, the use of variables i and j produces meaningless reactions such as

X[0] -> X[0]

which do not occur in the desired migration model.

To address this, one or more <predicate> elements can be included in the <reaction> element body. These elements allow users to specify an expression that must evaluate to true for the reaction to be included in the model. The reaction is included only if all predicate expressions evaluate to true.

For example, the following reaction element uses a predicate to avoid the meaningless "X[i] -> X[i]" migrations:

    <reaction spec="Reaction" reactionName="migration" rate="1.0">
        X[i] -> X[j]
        <predicate spec="Predicate" value="i != j"/>
    </reaction>

Refer to the MASTER expression syntax page for details on how to construct valid expressions.

Reaction rate multipliers

When specifying multiple reactions in this way, it is often desirable to allow the reaction rate to be a function of the exact locations involved. This can be accomplished through the use of a <rateMultiplier> element, which allows the user to provide an expression involving index variables. This expression is evaulated for each reaction produced by the reaction string and the result is multiplied by the value of the reactionRate attribute to produce a rate specific to that reaction.

For instance, suppose we want the migrations to vary inversely with the absolute difference between i and j. The following reaction element uses a <rateMultiplier> element to accomplish this:

    <reaction spec="Reaction" reactionName="migration" rate="1.0">
        X[i] -> X[j]
        <predicate spec="Predicate" value="i != j"/>
        <rateMultiplier spec="RateMultiplier" value="1/abs(i-j)"/>
    </reaction>

User-defined functions

Certain expressions may appear multiple times in predicates or reaction rate multipliers. To further simplify the specification of the model, MASTER allows users to define named functions of zero or more variables by adding function elements to the model.