Note: requires CRNSimulator.m and 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

Brusselator, deterministic (for comparison)

alpha=1.0;  (* units nM/s *)
beta = 2.5;  (* units /s *)
rsys={
rxn[X+X+Y,X+X+X,1], (* units 1 /nM/nM/s *)
rxn[X,Y,beta], (* units /s *)
rxn[X,0,1],  (* units /s *)
rxn[0,X,alpha],  (* units nM/s *)
conc[X,5], (* units nM *)
conc[Y,5]     (* units nM *)
}

tmax=45;
sol=SimulateRxnsys[rsys,tmax];
plotter={X[t],Y[t]}/.sol;Plot[plotter,{t,0,tmax},PlotRange->{0,All},PlotStyle->{{Red},{Green}}]

Just for completeness, it is nice to compare the CRN ODE simulations to direct simulations of the ODEs without the CRNSimulator package.  We’ll first write down the ODEs by hand, then add initial conditions and use numerical integration.

Brusselator={x'[t]==1-3.5x[t]+x[t]^2y[t],y'[t]==2.5x[t]-x[t]^2y[t]}

tmax=45;
sol=NDSolve[Join[Brusselator,{x[0]==5,y[0]==5}],{x,y},{t,0,tmax}];
Plot[Evaluate[{x[t],y[t]}/.sol],{t,0,tmax},PlotRange->{{0,tmax},{0,10}},PlotStyle->{{Red},{Green}}]

Brusselator, stochastic

Concentration c in volume V has c*A*V molecules, where A is Avogadro’s number.
E.g. c = 1 nM, A = 6.02*10^23 molecules/mole, V = 1 fL, 
       c*A*V = 10^-9 moles/L * 6.02 * 10^23 molecules/mole * 10^-15 L = 0.602 molecules.
To put it another way, 1 molecule in a volume of 1 fL has concentration of 1.66 nM,
or, 1 molecule in a volume of 1.66 fL has a concentration of 1 nM.
       
The propensity of a reaction indicates its probability per unit time, e.g. 5.3 per second.
The stochastic simulator calculates the propensity of a reaction by multiplying the counts of each reactant species times the given rate constant.  E.g.
      reaction X --> Y with rate constant 5.3 and #X = 3 has propensity 15.9 /s;
      reaction X+Y --> Z+A with rate constant 1.2 and #X=3, #Y=2 has propensity 7.2 /s;
      reaction 0 --> A+B with rate constant 2.7 has propensity 2.7 /s.
If a reactant appears in higher-than-unity stoichiometry, decremented counts are used:
      reaction X+X --> Z+A with rate constant 1.2 and #X=4 has propensity 1.2*4*3 = 14.4.
      
If the original mass action rate constants had units chosen with respect to nM and s, then this choice corresponds to a default volume of 1.66 fL.   

You can verify this by examining how many molecules per second are produced by a given reaction, according to the deterministic or stochastic formulation.  Consider the reaction 0 --> A+B with rate constant 2.7 nM/s.  In volume 1.66 fL, that’s 2.7 molecules/s, following the above calculation.

To run a stochastic simulation in a larger or smaller volume, with state variables being molecular counts rather than concentrations, the given rate constants must be modified accordingly.  Say the volume is Vm-fold larger than the default volume, then:
      for zero-order reactions, multiply the rate constant by Vm;
      for first-order reactions, don’t change the rate constant;
      for second-order reactions, decrease the rate constant by Vm;
      for third-order reactions, decrease the rate constant by Vm^2;
      etc.
      
    The example below is interesting for Vm=1 and Vm=100 as well...

Vm = 10; (* how many units of 1.66 fL *)
alpha=1.0;  (* units counts/s *)
beta = 2.5;  (* units /s *)
rsys={
rxn[X+X+Y,X+X+X,1/Vm/Vm], (* units 1 /count/count/s *)
rxn[X,Y,beta], (* units /s *)
rxn[X,0,1],  (* units /s *)
rxn[0,X,alpha*Vm],  (* units counts/s *)
conc[X,5*Vm], (* units molecule counts *)
conc[Y,5*Vm]     (* units molecule counts *)
}

tmax=45;
trajectory=SimulateRxnsysSSA[rsys,tmax];
ListLinePlot[trajectory,PlotRange->{0,All},
PlotLegends->SpeciesInRxnsys[rsys]]

Instead of ListPlot or ListLinePlot, which just use the {time, counts} points and optionally draws lines between then, it is also possible to use trajectory[t] as a function that returns the current counts for any given time.  Thus, rather than a saw-tooth plot, you will see a step plot, which more accurately captures the model.  It’s easier to see with smaller Vm and smaller range of times.

Plot[{trajectory[t][[1]],trajectory[t][[2]]}, {t,0,45}, PlotRange->All]

The trajectory[t] notation also makes it easy to superimpose the continuous deterministic ODEs with the stochastic simulation results, after scaling the counts back to nM concentrations.  With larger Vm, the agreement is closer -- but you won’t see indistinguishable behavior until Vm is so large that the stochastic simulations are too slow.

Show[
Plot[trajectory[t][[1]]/Vm,{t,0,tmax},PlotRange->{0,All},PlotStyle->{{Blue}}],
Plot[Evaluate[X[t]/.sol],{t,0,tmax},PlotRange->{0,All},PlotStyle->{{Red}}]
]

If you want to look under the hood, the trajectory is represented as a Mathematica TimeSeries object; you can access the values and times as shown below.  This will be more sensible if Vm=1...

trajectory

?TimeSeries

trajectory["Values"]

trajectory["Times"]

ListLinePlot[trajectory["Values"],AxesLabel->{X,Y}]

Predator-prey (Fox, Rabbit)

V=50;  (* volume is V * 1.66 fL when using nM in rate constants *)
tmax=50;
rsys={
rxn[R+F,2F,1.5/V],
rxn[R,2R,1],
rxn[F,0,1],
conc[{F,R},V/2]
};

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.

Plot[trajectory[t][[1]], {t,0,2}, PlotRange->All]

Digital Circuits with Stochastic Semantics

These definitions are taken directly from DigitalCircuits.nb, but with the following modifications for stochastic semantics:  (1) since the rate constants are interpreted as propensity rate constants in the default volume Vo = 1.66 fL, if we want to consider larger volumes V = m * Vo, we need to scale the given rate constants appropriately; (2) since conc[ ] statements are interpreted as counts, equal to nM concentrations only with the default volume, we ensure that initial values are integers and scale with m -- that is to say we are using ON = 1 nM and OFF = 0 nM regardless of volume, but integral multiples m of the default value will imply ON should have a count of m molecules; and (3) for non-input wires, rather than starting in an analog “half ON, half OFF” concentration, we start with an initially random number of molecules between 0 and m.

Even though our stochastic semantics is discrete, we retain the name “analogGate” for the four-reaction version without a restoring repeater, and the name “digitalGate” for the six-reaction version with a restoring repeater.

analogGate[input1_,input2_,output_,function_,m_]:= {
rxn[input1[0]+input2[0]+output[1-function[0,0]],
input1[0]+input2[0]+output[function[0,0]],1/m^2],
rxn[input1[0]+input2[1]+output[1-function[0,1]],
input1[0]+input2[1]+output[function[0,1]],1/m^2],
rxn[input1[1]+input2[0]+output[1-function[1,0]],
input1[1]+input2[0]+output[function[1,0]],1/m^2],
rxn[input1[1]+input2[1]+output[1-function[1,1]],
input1[1]+input2[1]+output[function[1,1]],1/m^2]
}

Wires[rxns_]:= SpeciesInRxnsys[rxns]/.{v_[_Integer]:>v} //Union

InitialConcentrations[rxns_,inputvalues_,m_]:=Module[{wires,inputwires,inputconcs,otherconcs},
wires=Wires[rxns];
inputwires=First/@inputvalues;
inputconcs=Join@@({conc[#[[1]][1],#[[2]]], conc[#[[1]][0],m-#[[2]]]}& /@ inputvalues);
otherconcs = Join@@(With[{c=RandomInteger[{0,m}]},{conc[#[0],c],conc[#[1],m-c]}]&/@Complement[wires,inputwires]);
Union[inputconcs,otherconcs]
]

repeaterGate[input_,output_,m_]:= {
rxn[input[0]+input[0]+output[1],input[0]+input[0]+output[0],1/m^2],
rxn[input[1]+input[1]+output[0],input[1]+input[1]+output[1],1/m^2]
}

digitalGate[input1_,input2_,output_,function_,m_]:= Module[{middle},
Join[analogGate[input1,input2,middle,function,m],
repeaterGate[middle,output,m]]]

fAND = Min;
fOR = Max;
fNOR =(1-Max[#1,#2])&;
fNAND = (1-Min[#1,#2])&;
fANDNOT=Min[#1,1-#2]&;

digitalSqrtCircuit[m_] := Join[
digitalGate[X1,X2,Z1,fNOR,m],
digitalGate[X3,X4,Z2,fANDNOT,m],
digitalGate[X3,X4,Y2,fOR,m],
digitalGate[Z1,X4,Z5,fAND,m],
digitalGate[Z5,X3,Z4,fNAND,m],
digitalGate[Z1,Z2,Z3,fOR,m],
digitalGate[Z4,Z3,Y1,fNAND,m]
];

analogSqrtCircuit[m_]:= Join[
analogGate[X1,X2,Z1,fNOR,m],
analogGate[X3,X4,Z2,fANDNOT,m],
analogGate[X3,X4,Y2,fOR,m],
analogGate[Z1,X4,Z5,fAND,m],
analogGate[Z5,X3,Z4,fNAND,m],
analogGate[Z1,Z2,Z3,fOR,m],
analogGate[Z4,Z3,Y1,fNAND,m]
];

The tuples listed in trajectory[“Values”] are the counts of species, in the same order as listed by the “SpeciesInRxnsys” function.  To figure out which index corresponds to a given species, you can use “Position”, which gives a list of each occurrence of a pattern.

SpeciesInRxnsys[digitalSqrtCircuit[1]]

Position[SpeciesInRxnsys[digitalSqrtCircuit[10]],Y1[1]]

To automate running the circuit in volume m, with max m counts per species, we create a script.  The plot is unintelligible when m=1 (so other means will be needed for making sense of that important case) but can be used to understand the circuit run with larger counts.

singleinputplot[circuit_,m_,{x4_,x3_,x2_,x1_}]:= Module[{rxns,rsys,trajectory},
rxns=circuit[m];
rsys= Join[rxns, 
InitialConcentrations[rxns,{{X1,x1},{X2,x2},{X3,x3},{X4,x4}},m]];
tmax=10;
trajectory=SimulateRxnsysSSA[rsys,tmax];
ListLinePlot[trajectory,
PlotRange->{0,All},
PlotLegends->SpeciesInRxnsys[rsys],
MaxPlotPoints->500 ]
]

singleinputplot[analogSqrtCircuit,1,{1,0,0,1}]

Here we see that, with perfect input, the output goes to perfect completion, unsurprisingly.

singleinputplot[analogSqrtCircuit,100,{100,0,0,0}]

With imperfect input, the output is worse -- although the digital circuit with restoration does manage the noise.

singleinputplot[analogSqrtCircuit,100,{80,20,20,20}]

singleinputplot[digitalSqrtCircuit,100,{80,20,20,20}]

There are so many species, it’s hard to tell what’s what, even with the legend.  Let’s just pull out the output variable Y1 and Y2, and run every input case.  “noise” is not a random variable here; it just indicated, deterministically, that the initial states are off by that many counts.

inputcases[circuit_,m_,noise_]:=Module[{rxns,rsys,trajectory,tmax,Y1i,Y2i},
Flatten[Table[
rxns=circuit[m];
rsys= Join[rxns, 
InitialConcentrations[rxns,{{X1,x1},{X2,x2},{X3,x3},{X4,x4}},m]];
Y1i=Position[SpeciesInRxnsys[rsys],Y1[1]][[1,1]];
Y2i=Position[SpeciesInRxnsys[rsys],Y2[1]][[1,1]];
tmax=10;
trajectory=SimulateRxnsysSSA[rsys,tmax];
tmax=Last[trajectory["Times"]];
{x4,x3,x2,x1,trajectory[tmax][[Y2i]],trajectory[tmax][[Y1i]]},
{x4,{noise,m-noise}},
{x3,{noise,m-noise}},
{x2,{noise,m-noise}},
{x1,{noise,m-noise}}],3]];

For unit counts (0 or 1), the circuit computes the (rounded-down) square root perfect.

inputcases[analogSqrtCircuit,1,0]//MatrixForm

When counts are increased, we see that the output is generally less perfect, but certain inputs get computed more noisily than others.  Every time you run it, the results will be a bit different, but certain inputs are fairly consistently poorly computed.  A 10-count deviation in the input becomes a 30-count deviation in some cases.

inputcases[analogSqrtCircuit,100,10]//MatrixForm

However, with restoration in every gate, the outputs remain no worse than the inputs (typically).

inputcases[digitalSqrtCircuit,100,10]//MatrixForm

Exploring reachability by simulation

Here we take a look at the five-reaction CRN discussed in class, for which we asked “Starting with #A=1 and #C=n, what is the highest #D that can be achieved?”  We argued that the answer is 2^n.   Below, we use lower-case species names because both C and D have defined meanings in Mathematica.

reach[n_]:={
rxn[a,b,1],
rxn[c,d,1],
rxn[b+c,a,1],
rxn[a+d,a+2e,1],
rxn[b+e,b+d,1],
conc[a,1], conc[c,n]
};

rsys=reach[5]

Because stochastic simulation explores the state space probabilistically based on propensities, it is likely to prefer certain areas of the state space.  For this system, because #C is strictly decreasing and a reaction doing so is always active, and because of the structure of the other reactions (can you see why?), the simulation is probabilistically guaranteed to terminate.  So we don’t need to set an end time.

trajectory=SimulateRxnsysSSA[rsys,Infinity];
ListLinePlot[trajectory,
PlotRange->{0,All},
PlotLegends->SpeciesInRxnsys[rsys],
MaxPlotPoints->500 ]

To process the trajectories, we note that #C and #D re the third and fourth entries in the counts list, respectively.  First taking a quick look at the first 5 steps of the trajectory, we can then plot what (#C,#D) values are visited.

SpeciesInRxnsys[rsys]

Take[trajectory["Values"],5]

Take[trajectory["Values"],5][[All,{3,4}]]

ListLinePlot[trajectory["Values"][[All,{3,4}]],AxesLabel->{C,D}]

We can combine the results from a large number of simulations. It is a bit slow, though, so more than 200 trials will take quite some time.

visited=Union@@Table[
trajectory=SimulateRxnsysSSA[rsys,Infinity];
Union[trajectory["Values"][[All,{3,4}]]],
{200}];
ListPlot[visited,PlotRange->{{-0.3,5},{-0.3,32}},AxesLabel->{C,D}]

You can see that the simulation doesn’t fully explore the space, as we know that for n=5, #D can reach 32.  Can you think of a more efficient way to exhaustively explore the reachable state space?