Note: requires CRNSimulatorSSA.m to be loaded.

Note that the input format to SimulateRxnsysSSA is consistent with SimulateRxnsys in that a reaction system is created with rxn and conc statements. Arbitrary order (zero-order, unimolecular, bimolecular, etc) reactions are supported. The output of SimulateRxnsysSSA is a TimeSeries object (list of discrete times and corresponding molecular counts).

See documentation for SimulateRxnsysSSA:

?SimulateRxnsysSSA

Approximate majority

n=1000;

rsys={
rxn[x+y,x+b,1/n],
rxn[x+y,y+b,1/n],
rxn[x+b,2x,1/n],
rxn[y+b,2y,1/n],

conc[{x,y},n/2]};

(* Simulate until no further reactions are possible *)
trajectory=SimulateRxnsysSSA[rsys,∞];

ListLinePlot[
trajectory,
PlotRange->{0,All},
PlotLegends->SpeciesInRxnsys[rsys],
MaxPlotPoints->500 (* Limit the number of points to plot faster *)]

Predator-prey

n=300;

tmax=50;

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

conc[{x1,x2},n/2]};

(* Simulate for tmax time *)
trajectory=SimulateRxnsysSSA[rsys,tmax];

ListLinePlot[
trajectory,
PlotRange->{0,All},
PlotLegends->SpeciesInRxnsys[rsys],
MaxPlotPoints->500 (* Limit the number of points to plot faster *)]

The TimeSeries object “trajectory” can be used as either a list of times and counts, or counts as a function of t (a continuous variable). Thus we can also pass it to Plot as follows.

(* trajectory[t] represents the counts of species as a function of t *)
Plot[trajectory[t][[1]], {t,0,10}, PlotRange->All]