Note: requires CRNSimulator.m to be loaded.
Note that when you close Mathematica, it may ask you whether you want to save CRNSimulator.m because it has changed.  There is no need to save it;  answer “no”.

See basic documentation for the key functions:

?SimulateRxnsys

?RxnsysToOdesys

?SpeciesInRxnsys

For this package, a model reaction system consists of an unordered set of reactions and concentration statements.  Species names are simply Mathematica symbols (i.e. variables) that have no values or rules assigned to them.  Thus, for example, the symbols “C” and “D” and “N” and “I” should not be used as species names, because Mathematica has built-in meanings for these terms.  Similarly, if your code assigns e.g. x=17, then “x” cannot be used as a species name unless you first execute Clear[x].  Once a reaction system has been defined, SimulateRxnsys[ ] builds the corresponding ODE equations for continuous deterministic mass-action chemical kinetics and performs the simulation.  The numerical integration of ODEs in Mathematica is an interpolating function that can be queried to obtain the value at a particular point.  These functions can be passed to Plot[ ] to draw the results, or they can be used in others ways to interrogate the simulation’s results.

What follows are a few annotated examples.

Predator-Prey

The classical Lotka-Volterra model for ecological dynamics can be represented as a CRN with three reactions.  The rxn[ ] statement indicates the reactants, products, and rate constant.  Note that when a rxn[ ] statement is displayed, it gets printed such that it looks more like you would expect.

Note that when a reaction has no products, that is indicated by placing a constant where products would be listed -- in this case 0, but any other constant would give the same result.  Similarly, for zero-order reactions that consume no reactants, a constant is placed in the reactant slot.

rsys={
rxn[x1+x2,2x2,1.5],
rxn[x1,2x1,1],
rxn[x2,0,1],

conc[x1,2],
conc[x2,1]
}

tmax=45;

sol=SimulateRxnsys[rsys,tmax]

Mathematica functions often return “rewriting rules” to express the solutions to problems.  In this case, it means that the variable “x1”, representing a species concentration, can be rewritten as a numerically-calculated function of time.  Below, we use “/.” to apply the rewriting rule and thus obtain an expression encapsulating the concentrations as a function of time, with the time t being explicitly written as needed by the Plot[ ] function.

plotter={x1[t],x2[t]}/.sol

Plot[plotter,{t,0,tmax},PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue}}]

You can set up the same predator-prey ODEs by hand if you like; this is what SimulateRxnsys is doing under the hood automatically.

tmax=45;
sol=NDSolve[{
x1'[t]== x1[t]-1.5*x1[t]*x2[t],
x2'[t]==1.5*x1[t]*x2[t]-x2[t],
x1[0]==2,x2[0]==1},
{x1,x2},{t,0,tmax}]

Plot[Evaluate[{x1[t],x2[t]}/.sol],{t,0,tmax},PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue}}]

ParametricPlot[Evaluate[{x1[t],x2[t]}/.sol],{t,0,tmax},AxesLabel->{x1,x2}]

Showing: reversible reactions, nanomolar concentrations

It is important to understand that the user is responsible for keeping track of units.  In the example below, the forward rate constants for bimolecular reactions are 1.5*10^6 /M/s and 10^5 /M/s respectively, while the reverse reaction is 10 /s.  The concentration is 10^-9 M, i.e. 1 nM.  Time is expressed in seconds.  

We also introduce here the revrxn[ ] statement, which is merely a shorthand for both the forward reaction and the reverse reaction simultaneously.

rsys={
revrxn[a+b,c,1.5 10^6,10],
rxn[2a,b,10^5],

conc[{a,b,c},10^-9]
}

tmax=10*60*60;

sol=SimulateRxnsys[rsys,tmax];

plotter={a[t],b[t],c[t]}/.sol;

Plot[plotter,{t,0,tmax},PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue}}]

Showing: analog computational circuit

Here we show the CRN implementation of the function discussed in the lecture:  
Z_infinity = max(min(A_0,B_0), A_0 / 2 + C_0).
Note that the species C is represented by “c” because Mathematica has a built-in meaning for “C”.
This simulation time is such that the CRN hasn’t quite gone to completion yet, but is close.

rsys={
rxn[A,X+Y,1],
rxn[2X,P,1],
rxn[c,P,1],
rxn[Y+B,Q,1],
rxn[P,Z+R,1],
rxn[Q,Z+S,1],
rxn[R+S,T,1],
rxn[T+Z,0,1],
conc[A,5],conc[B,2],conc[c,3]
}

tmax=10;

sol=SimulateRxnsys[rsys,tmax];

plotter={A[t],B[t],c[t],Z[t]}/.sol;

Plot[plotter,{t,0,tmax},PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue},{Black}}]

{Z[tmax],Max[Min[A[0],B[0]],A[0]/2+c[0]]}/.sol

Showing: defining modules

An important feature of the CRNSimulator package is that it is possible to define CRN “subroutines” or “modules”, and systematically built larger CRNs based on those generic components.

Here, bibi[ ] produces a CRN representing a molecular implementation of “a + b --> c + d” via the consumption of a set of fuels, “gate”, “backward”, “agate”, “out”, “atrans”, and “helper”.  These fuels are specific to this particular reaction.  In Mathematica, local variables in a Module[ ] statement are given unique instance names, e.g. “gate$4120”, so the reactions returned by different bibi[ ] function calls contain different oddly-named species that won’t interfere with each other.

Also note that “trans” is an indexed species name, again representing a fuel, but in this case one that we want to conveniently refer to later as we look at the results of the simulation.

In the module definition, Seq[ ] is used to wrap the reactions and concentrations (rather than { }). This is so that when multiple modules are combined in rsys we won’t get a nested list. Alternatively, Flatten[ ] could be used to flatten the resulting nested list.

keff=300000;

tmax=5*24*60*60;

bibi[a_,b_,c_,d_]:=
Module[{gate,backward,agate,out,atrans,helper},
Seq[
revrxn[a+gate,backward+agate,keff/10,keff],
rxn[b+agate,out,keff],
revrxn[out+trans[c], c+atrans,keff/10,keff],
rxn[helper+atrans,d,keff],

conc[{gate,trans[c]}, 200 10^-9],
conc[backward,100 10^-9],
conc[helper, 200 10^-9]
]]

The system that we will simulate corresponds to an implementation of { b+a --> b+b, c+b --> c+c, a+c --> a+a }, which goes by the cyclic predator-prey oscillator, the rock-paper-scissors CRN, and other names.  In our case, the fuel is in roughly 10-fold excess of the signal species, and after about 6 oscillations the fuel is used up.

rsys=
{bibi[b,a,b,b],
bibi[c,b,c,c],
bibi[a,c,a,a],

conc[b,20 10^-9],
conc[{a,c},5 10^-9]
}

sol=SimulateRxnsys[rsys,tmax];

plotSpecies=
Join[
{a,b,c},
{trans[a],trans[b],trans[c]},
SpeciesInRxnsysStringPattern[rsys,"gate$*"]];

plotter=#[t]/.sol&/@plotSpecies;
Plot[plotter,{t,0,tmax},AxesOrigin->{0,0},PlotRange->{0,100 10^-9},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Red,Thick,Dashed},{Green,Thick,Dashed},{Blue,Thick,Dashed}}]

Example of adding/removing species during experiment

Note: Also requires CRNSimulatorExtensions.m to be loaded

rsys={
revrxn[a+b,c,0.1,30],
rxn[2a,b,0.01],

conc[c,1]
};

Our interest here is to see how a CRN responds to input that is given at various times, not just as the initial conditions.

We don’t specify the initial concentrations of a, b because these are controlled by a schedule. Positive numbers means add this amount, negative numbers means remove this amount. The first column is the time to wait until next addition.  The remaining columns indicate the amounts to add (or remove) for the respective species listed in the SimulateRxnsysWithSchedule[ ] call.

schedule={
{50,4,0},
{20,0,10},
{20,10,0},
{20,0,-9},
{20,-2,0}
};

SpeciesInRxnsys[rsys]

tmax=Total[schedule[[All,1]]]

sol=SimulateRxnsysWithSchedule[rsys,schedule,{a,b}];

plotter={a[t],b[t],c[t]}/.sol;

Plot[plotter,{t,0,tmax},PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue}}]

Example of providing an external input

You might not want to provide discrete pulse inputs, but rather set a time-varying concentration.

You can intermix definitions of the form x[t]==..., x’[t]==..., with rxn, term, and conc statements.

tmax=20;

rsys=
{u[t]==0.1 Sin[t]+1,
revrxn[x1+u,x2,1,1],
conc[x1,1]
};

sol=SimulateRxnsys[rsys,tmax];

plotter={x1[t],x2[t],u[t]}/.sol;

Plot[plotter,{t,0,tmax},PlotRange->All]