Note: requires CRNSimulator.m to be loaded.

See basic documentation for the key functions:

?SimulateRxnsys

?RxnsysToOdesys

?SpeciesInRxnsys

See more details for the different statements allowed in rxnsys: 
rxn[], revrxn[], rxnl[], term[], conc[]
by querying them with “?” like this:

?revrxn

Predator-Prey

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

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

tmax=45;

sol=SimulateRxnsys[rsys,tmax];

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

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

Showing: reversible reactions, nanomolar concentrations

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: defining modules

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]

]]

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]
};

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.

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

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 non-mass action kinetics: Repressilator

A reaction system can include term[] statements (eg. [x, y[t]/(1+y[t])]), which are additively combined
together to generate the ODEs. You can intermix rxn[] and term[] statements.

tmax=20;

n=4;γ=1;β=10;k=5;

rsys=
{term[x2,((β (k^n))/((k^n)+((x1[t])^n)))],term[x3,((β (k^n))/((k^n)+((x2[t])^n)))],term[x1,((β (k^n))/((k^n)+((x3[t])^n)))],
term[x1,-γ x1[t]],term[x2,-γ x2[t]],term[x3,-γ x3[t]],

conc[x1,1],conc[x2,0.9],conc[x3,0.1]};

sol=SimulateRxnsys[rsys,tmax];

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

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

Example of providing an external input

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]