Reaction Schema Simulator Mathematica Package

Version 1.0 notebook: ReactionSchemaSimulator.nb
Last modified: July 14, 2020 by Erik Winfree.
Background: A few years ago, David Soloveichik
wrote CRN
Simulator package for Mathematica, which I find to be
exceptionally convenient for quickly exploring CRN behaviors. It can
simulate both with bulk deterministic (ODE) semantics and with
finitevolume stochastic (CTMC) semantics. Chemical species are just
Mathematica expressions, so it's particularly easy to, e.g. index a
species as X[5], or build a CRN using Mathemtica code. Data analysis
is instantly available, as is documenting one's explorations in a
notebook. For stochastic simulations, the specific CRN is seamlessly
compiled down to C code, so it's pretty fast (although the Gillespie
direct method is used, so CRNs with lots of reactions could be
simulated more efficiently with more advanced algorithms, such as the
one I used
here.
Anyway, I use it a lot, and I teach with it.
Motivation: For simulators such as the above, CRNs must
have a finite number of species and a finite number of reactions 
the representation of the CRN is a list of the reactions. There are
also many cases where it makes sense to consider an infinite number of species,
for example polymers with arbitrary length, or heteropolymers with arbitrary
sequence. In that case, there will also be an infinite number of reactions.
To provide a finite specification of such a system, one needs a finite set of
rules decribing what reactions are possible. In the context of
species being represented as Mathematica expressions, it is natural
for things like polymers to be represented as expressions of arbitrary
length (e.g. DNA[G,A,T,T,A,C,A]) and for reaction to be obtained using
pattern matching rules (e.g. SmaI + DNA[x___,C,C,C,G,G,G,y___] > SmaI
+ DNA[x,C,C,C] + DNA[G,G,G,y]). As illustrated by this example,
Mathematica has a powerful builtin language for pattern matching,
which allows us to concisely express a wide class of such rules 
which we call reaction schemata  in a natural and intuitive way.
Therefore, I prototyped a package for simulating generalized reaction
schemata that can model the activities of polymers, polymermodifying
enzymes, and more abstract processes.
Related work: DNA strand displacement simulators (Visual
DSD, RuleDSD, Peppercorn, etc.), tile selfassembly simulators (xgrow,
ISU TAS, etc.), systems biology simulators that handle the
combinatorics of protein modification and subunit assembly (BioNetGen,
kappa, etc.) all deal with the issue of generalized reaction
mechanisms and combinatorial or infinite numbers of potential species.
Often, optimizations are possible when considering a specifical
application domain. Our simulator is not meant to compete with these;
indeed, it is not particularly suitable for branched polymers and 2D
crystals. Rather, it is meant more as a quickanddirty approach that
enables easy exploration, in Mathematica, of systems that naturally
fall into its domain. This includes the linear Polymer Reaction
Networks
studied here.
Status: Currently, the Reaction Schema Simulator is in the form of a notebook that combines the core code with a variety of example systems. When I have time, I will dissect it into a proper package and an examples notebook.
A product of the The DNA and Natural Algorithms Group.