# Reaction Schema Simulator Mathematica Package

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 finite-volume stochastic (CTMC) semantics. Chemical species are just Mathematica expressions, so it's particularly easy to, e.g. index a species as X, 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 built-in 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, polymer-modifying enzymes, and more abstract processes.

Related work: DNA strand displacement simulators (Visual DSD, RuleDSD, Peppercorn, etc.), tile self-assembly 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 quick-and-dirty 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.

• Version 1.0 notebook: ReactionSchemaSimulator.nb