This Mathematica notebook requires the CRNSimulator package.  It can be downloaded from: http://dna.caltech.edu/SupplementaryMaterial/DNA_for_CRNs/

We try to load it from the same directory as this notebook:

nbdir=DirectoryName[ToFileName["FileName"/.NotebookInformation[SelectedNotebook[]]]];

Get["CRNSimulator.m",Path->nbdir]

Code

Code to convert formal CRN to DNA implementation

Code assumes cmax and qmax globals.

Produces reactions in srxn[[reaction index],{x1,x2},{x3},k] format ("structured reaction system"):

rxnsysToSrxnsys[rsys_]:=
Module[{i=1},
rsys/.
rxn[rs_,ps_,k_]:>
srxn[i++,
{Unevaluated[rs]/.{Plus->Seq,c_Integer*s_:>Seq@@Table[s,{c}] }},{Unevaluated[ps]/.{Plus->Seq,c_Integer*s_:>Seq@@Table[s,{c}] }},
k]]

Totals formal rate constants of bimolecular reactions with s as the first reactant:

computeSigmaj[s_,rsys_]:=Total[Cases[rxnsysToSrxnsys[rsys],srxn[_,{s,_},_,k_]:>k]]

Maximum over s:

computeSigma[rsys_]:=Max[computeSigmaj[#,rsys]&/@SpeciesInRxnsys[rsys]]

Compute γ as discussed in the text:

computeGamma[rsys_]:=(qmax^(-1))(qmax-computeSigma[rsys])

Find the lower-bound on qmax as discussed in the text. Note that the generateBufferingModules -> False option can be used when buffering modules are not wanted.

Options[computeQmaxLowerBound]={generateBufferingModules->True};

computeQmaxLowerBound[rsys_,OptionsPattern[]]:=
Module[{σ=computeSigma[rsys],srsys=rxnsysToSrxnsys[rsys]},
Max[
Cases[srsys,srxn[_,{_},_,k_]:>k/cmax+σ],
Cases[srsys,srxn[_,{_,_},_,k_]:>k+σ],
If[OptionValue[generateBufferingModules],
(2σ-computeSigmaj[#,rsys])&  /@  SpeciesInRxnsys[rsys],
0]
]]

Definitions of reaction and buffering modules:

Bimolecular reaction module:

dnarxn[i_,{r1_,r2_},{ps__},k_,γ_]:=
Seq[
revrxn[r1+gateL[i],gateH[i]+strandB[i],(γ^(-1))k ,qmax],
rxn[r2+gateH[i],strandO[i],qmax],
rxn[strandO[i]+gateT[i],Plus[ps],qmax],

conc[{gateL[i],strandB[i],gateT[i]},cmax]
]

Unimolecular reaction module:

dnarxn[i_,{r1_},{ps___},k_,γ_]:=
Seq[
rxn[r1+gateG[i],strandO[i],(γ^(-1)) k /cmax],
rxn[strandO[i]+gateT[i],Plus[ps],qmax],

conc[{gateG[i],gateT[i]},cmax]
]

Buffering reaction module:

dnabalrxn[s_,γ_,σ_,σj_]:=
Module[{j=ToString[s]},
Seq[
revrxn[s+gateLS[j],gateHS[j]+strandBS[j],(γ^(-1))(σ-σj),qmax],

conc[{gateLS[j],strandBS[j]},cmax]
]]

Main function to convert a formal CRN to the DNA implementation.  Note that the generateBufferingModules -> False option can be used when buffering modules are not wanted.

Options[rxnsysToDnasys]={generateBufferingModules->True};

rxnsysToDnasys[rsys_,OptionsPattern[]]:=
Module[
{γ=computeGamma[rsys],
σ=computeSigma[rsys]},
Join[
rxnsysToSrxnsys[rsys]/.
{srxn[i_,rs_,ps_,k_]:>dnarxn[i,rs,ps,k,γ],
conc[s_,c_]:>conc[s,(γ^(-1))c]},
If[computeSigmaj[#,rsys]<σ && OptionValue[generateBufferingModules],
dnabalrxn[#,γ,σ,computeSigmaj[#,rsys]],Seq[]]&  
/@  SpeciesInRxnsys[rsys]]
]

Computes the concentration and time scaling factors for the DNA system so that largest concentration cmax is maxconc (usually 10^-5), and max rate constant qmax is maxrateconst (usually 10^6):

computeConcAndTimeScaling[maxconc_,maxrateconst_]:=
{maxconc/cmax, (((maxconc/cmax))^(-1))(((maxrateconst/qmax))^(-1))}

Code for controlling simulations

Given a solution, generates a new reaction system with the initial concentrations given by the final concentrations of the solution.

rsysForRestart[rsys_,sol_,tmax_]:=
Join[
rsys/.conc[___]->Seq[],
conc[#,(#[tmax]/.sol)]&/@SpeciesInRxnsys[rsys]]

For formal functions of one argument, whose definitions are valid starting at 0 ending at end1 and end2.  Make new function which uses the second one when first one ends:

joinFunctions[f1_,f2_,end1_,end2_]:=
Function[Evaluate[Which@@{0≤#<end1, f1[#],end1≤#<end1+end2, f2[#-end1]}]]

Main function for simulating a reaction system using a schedule of addition/removal of species, rather than just starting in one initial configuration.  (In rxnsys, x1 must not have initial conc):

simulateRxnsysWithSchedule[rxnsys_,schedule_,x1_]:=
simulateRxnsysWithSchedule[rxnsys,schedule,{x1}]

simulateRxnsysWithSchedule::undershoot="Scheduled decrease (schedule step `1`: `2`) resulted in negative concentration. Going to 0.";

simulateRxnsysWithSchedule[rxnsys_,schedule_,xs_List]:=
Module[
{n=Length[schedule],rsys,lastsol,sol},

rsys=Join[rxnsys, MapIndexed[conc[#1,schedule[[1,1+#2[[1]]]]]&,xs]];
lastsol=SimulateRxnsys[rsys,schedule[[1,1]]];
sol=lastsol;

Do[
rsys=
rsysForRestart[rsys,lastsol,schedule[[i-1,1]]]/.
MapIndexed[
conc[#1,s_]:>
(If[s+schedule[[i,1+#2[[1]]]]<0,
Message[simulateRxnsysWithSchedule::undershoot,i,schedule[[i,1+#2[[1]]]]]];
conc[#1,Max[0,s+schedule[[i,1+#2[[1]]]]]])&,
xs];
lastsol=SimulateRxnsys[rsys,schedule[[i,1]]];
sol=MapThread[#1/.(s_->f1_):>(s->joinFunctions[f1,#2/.(_->f2_):>f2,Total[schedule[[;;i-1,1]]],Total[schedule[[;;i,1]]]])&,
{sol,lastsol}],

{i,2,n}];

sol
]

Since some species are buffered, we may not be able to instantaneously decrease concentration because this may result in going negative.  We do 1/3 of the distance, then 1/3 of the remaining distance, etc (here 3 is the stepfactor).  Rescaled so that after all iterations we cover the whole distance.

steppizeSchedule[schedule_,iterations_,stepfactor_:3,stepduration_:0.01]:=
Module[{sum=1-((stepfactor-1)/(stepfactor))^(iterations)},
Map[
Seq@@Table[Prepend[N[Rest[#]]* (stepfactor-1)^(i-1)(stepfactor)^(-i) / sum,If[i==iterations,First[#],stepduration]],{i,1,iterations}]&,
schedule]]

Code for modifying reaction system

Adds waste product w[i] to reaction i.  This is useful for monitoring the total flow through reactions.

addWasteProducts[rsys_,w_]:=
Module[{i=1},
rsys/.  
rxn[rs_,ps_,k_]:>
With[{pse=Replace[Plus[ps,w[i++]],1->0,1]},rxn[rs,pse,k]]
]

Examples

Predator-Prey

Ideal system

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=
#[t]/.sol&/@{x1,x2};

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

DNA implementation

$HistoryLength=0;ClearSystemCache[];

cmax=1000;

qmax=computeQmaxLowerBound[rsys]

γ=computeGamma[rsys]

dnasys=rxnsysToDnasys[rsys];

sol=SimulateRxnsys[dnasys,tmax];

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

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

What fraction of auxiliary species was used up by the end?

1-(#[tmax]/.sol&/@SpeciesInRxnsys[dnasys,gateT[_]])/cmax

Scale into realistic regime:

{cscale,tscale}=computeConcAndTimeScaling[10^-5,10^6];

plotterscaled=cscale(#/.t->(tscale^(-1))t)&/@plotter;

Plot[plotterscaled,{t,0,tscale tmax},PlotPoints->300,PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue}}]

Oregonator

Ideal system

rsys={
rxn[x2,x1,ka*k1],
rxn[x2+x1,1,k2],
rxn[x1,2x1+x3,kb*k3],
rxn[2 x1,1,k4],
rxn[x3,x2,f*k5],
rxn[x3,1,(1-f)*k5],

conc[x1,1.001*5.2408379728095877`*^-8],
conc[x2,0.999*3.2271479050679757`*^-7],
conc[x3,0.999*6.110394049015663`*^-8]
}/.{ka->0.065,kb->0.065,k1->1.34,k2->1.6 10^9,k3->8 10^3,k4->4 10^7,k5->446,f->0.994}

Find equilibrium:

Solve[RxnsysToOdesys[rsys]/._'[_]->0,{x1[t],x2[t],x3[t]}]

tmax=4.6;

sol=SimulateRxnsys[rsys,tmax,AccuracyGoal->13,PrecisionGoal->13,WorkingPrecision->15];

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

Plot[plotter,{t,0,tmax},PlotPoints->300,PlotRange->{0,5 10^-7},PlotStyle->{{Red},{Green},{Blue}}]

DNA implementation

$HistoryLength=0;ClearSystemCache[];

cmax=0.0007;

qmax=computeQmaxLowerBound[rsys]

γ=computeGamma[rsys]

dnasys=rxnsysToDnasys[rsys];

sol=SimulateRxnsys[dnasys,tmax];

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

Plot[plotter,{t,0,tmax},PlotPoints->300,PlotRange->{0,5 10^-7},PlotStyle->{{Red},{Green},{Blue}}]

What fraction of auxiliary species was used up by the end?

1-(#[tmax]/.sol&/@SpeciesInRxnsys[dnasys,gateT[_]])/cmax

Scale system into realistic regime:

{cscale,tscale}=computeConcAndTimeScaling[10^-5,10^6];

plotterscaled=cscale(#/.t->(tscale^(-1))t)&/@plotter;

Plot[plotterscaled,{t,0,tscale tmax},PlotPoints->300,PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue}}]

Rossler

Ideal system

rsys={
rxn[x1,2x1,30],
rxn[2x1,x1,0.5],
rxn[x2+x1,2x2,1],
rxn[x2,1,10],
rxn[x1+x3,1,1],
rxn[x3,2x3,16.5],
rxn[2x3,x3,0.5],

conc[{x1,x2,x3},10]
};

tmax=40/4.1;

sol=SimulateRxnsys[rsys,tmax];

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

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

DNA implementation

$HistoryLength=0;ClearSystemCache[];

cmax=60000;

qmax=computeQmaxLowerBound[rsys]

γ=computeGamma[rsys]

dnasys=rxnsysToDnasys[rsys];

sol=SimulateRxnsys[dnasys,tmax];

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

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

What fraction of auxiliary species was used up by the end?

1-(#[tmax]/.sol&/@SpeciesInRxnsys[dnasys,gateT[_]])/cmax

Scale system into realistic regime:

{cscale,tscale}=computeConcAndTimeScaling[10^-5,10^6];

plotterscaled=cscale(#/.t->(tscale^(-1))t)&/@plotter;

Plot[plotterscaled,{t,0,tscale tmax},PlotPoints->300,PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue}}]

Pulse counter

Module definitions

repgate[x_,y_]:=
Module[{xx0,xx1},
Seq[
revrxn[2x[1],xx1,0.1/2,2/2],
rxn[y[0]+xx1,y[1]+xx1,20],

revrxn[2x[0],xx0,0.1/2,2/2],
rxn[y[1]+xx0,y[0]+xx0,20]
]]

andgate[x_[xv_],y_[yv_],z_]:=
Module[{bin,o},
Seq[
rxn[x[1-xv]+o[1],x[1-xv]+o[0],1],
rxn[y[1-yv]+o[1],y[1-yv]+o[0],1],
revrxn[x[xv]+y[yv],bin,0.1/2,2/2],
rxn[bin+o[0],bin+o[1],20],

repgate[o,z],

conc[o[0],1]
]]

norgate[x_,y_,z_]:=andgate[x[0],y[0],z]

latchgate[s_,q_]:=
Module[{r,qn},
Seq[
norgate[r,qn,q],
norgate[q,s,qn],

conc[{qn[1],r[0]},1]
]]

Edge-triggered master-slave flip-flop:

edgeflipflop[p_,q_]:=
Module[{a1,a2,a3,a4,a5,a6,a7},
Seq[
andgate[a7[0],p[1],a1],
andgate[a1[0],a5[0],a2],
andgate[a2[0],p[0],a3],
andgate[a3[0],a7[0],q],
andgate[p[1],q[0],a4],
andgate[a2[0],a4[0],a5],
andgate[p[0],a5[0],a6],
andgate[q[0],a6[0],a7],
conc[{a1[0],a2[0],a3[1],a4[0],a5[1],a6[0],a7[1]},1]
]]

4-Pulse Counter Ideal System

rsys={
repgate[in,p],
edgeflipflop[p,o1],
edgeflipflop[o1,o2],

conc[{p[0],o1[0],o2[0]},1]
};

schedule={
{15,1,0},
{5,-1,1},
{75-15,1,-1},
{5,-1,1},
{135-75,1,-1},
{5,-1,1},
{300-135,1,-1},
{5,-1,1},
{340-300,1,-1},
{5,-1,1},
{420-340,1,-1},
{5,-1,1},
{500-420,1,-1}
};

schedule=steppizeSchedule[schedule,10,4,0.2];

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

sol=simulateRxnsysWithSchedule[rsys,schedule,{in[0],in[1]}];

plotter=
#[t]/.sol&/@Join[{in[1],o1[1],o2[1]}];

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

DNA implementation, no buffering modules

$HistoryLength=0;ClearSystemCache[];

cmax=150;

qmax=computeQmaxLowerBound[rsys,generateBufferingModules->False]

γ=computeGamma[rsys]

dnasys=rxnsysToDnasys[rsys,generateBufferingModules->False];

dnaschedule=Prepend[#[[2;;]] (γ^(-1)),#[[1]]]&/@schedule;

sol=simulateRxnsysWithSchedule[dnasys,dnaschedule,{in[0],in[1]}];

plotter=
#[t]/.sol&/@{in[1],o1[1],o2[1]};

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

What fraction of auxiliary species was used up by the end?

1-(#[tmax]/.sol&/@SpeciesInRxnsys[dnasys,gateT[_]])/cmax

Scale system into realistic regime:

{cscale,tscale}=computeConcAndTimeScaling[10^-5,10^6];

plotterscaled=cscale(#/.t->(tscale^(-1))t)&/@plotter;

Plot[plotterscaled,{t,0,tscale tmax},PlotPoints->300,PlotRange->{0,All}]

DNA implementation, with buffering modules

$HistoryLength=0;ClearSystemCache[];

computeQmaxLowerBound[rsys]

cmax=150; qmax=50;

γ=computeGamma[rsys]

dnasys=rxnsysToDnasys[rsys];

dnaschedule=Prepend[#[[2;;]] (γ^(-1)),#[[1]]]&/@schedule;

sol=simulateRxnsysWithSchedule[dnasys,dnaschedule,{in[0],in[1]}];

plotter=
#[t]/.sol&/@{in[1],o1[1],o2[1]};

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

What fraction of auxiliary species was used up by the end?

1-(#[tmax]/.sol&/@SpeciesInRxnsys[dnasys,gateT[_]])/cmax

Scale system into realistic regime:

{cscale,tscale}=computeConcAndTimeScaling[10^-5,10^6];

plotterscaled=cscale(#/.t->(tscale^(-1))t)&/@plotter;

Plot[plotterscaled,{t,0,tscale tmax},PlotPoints->300,PlotRange->{0,All}]

Register machine

Ideal system

rsys={
(* clock reactions *)
rxn[s[1]+s[2],2s[2],2],
rxn[s[2]+s[3],2s[3],2],
rxn[s[3]+s[4],2s[4],2],
rxn[s[4]+s[5],2s[5],2],
rxn[s[5]+s[6],2s[6],2],
rxn[s[6]+s[7],2s[7],2],
rxn[s[7]+s[8],2s[8],2],
rxn[s[8]+s[9],2s[9],2],
rxn[s[9]+s[10],2s[10],2],
rxn[s[10]+s[11],2s[11],2],
rxn[s[11]+s[12],2s[12],2],
rxn[s[12]+s[1],2s[1],2],

(* clock initial concentrations *)
conc[{s[1],s[2]},6*5],
conc[{s[3],s[4],s[5],s[6],s[7],s[8],s[9],s[10],s[11],s[12]},6*1/500 10^-4],



(* logic reactions. w: stores the count in the current red-red cycle, v: used to store how many green peaks are left in a red-red cycle, r: backups the value 1.
*)
(*   logically occur in first separated clock cycle: *)
rxn[v2+c2,1,1],
rxn[v1+s[1],s[1]+v2+v3,1],
rxn[c1+s[1],s[1]+c2+c3,1],
rxn[w2+s[1],s[1]+w1,1],

(*   logically occur in second separated clock cycle: *)
rxn[r1+s[5],s[5]+r2+r3,1],
rxn[r2+v2,v2+d,1],
rxn[c2+r2,c2+i+w1,1],
rxn[i+w1,i+v1+w2,1],

(*   logically occur in third separated clock cycle: *)
rxn[v2+s[9],s[9],1],
rxn[c2+s[9],s[9],1],
rxn[i+s[9],s[9],1],
rxn[v3+s[9],s[9]+v1,1],
rxn[c3+s[9],s[9]+c1,1],
rxn[r3+s[9],s[9]+r1,1],
rxn[d+v1,1,1],


(* initial concentration of logic species *)
conc[{v3,r1},6*10],
conc[{w1},0],
conc[{c2},6*1/2*10],
conc[{c3},6*(1+1/2)*10]
};

tmax=52;

sol=SimulateRxnsys[rsys,tmax];

plotter=
#[t]/.sol&/@{i,d};

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

DNA implementation, no buffering modules

$HistoryLength=0;ClearSystemCache[];

cmax=20000;

qmax=computeQmaxLowerBound[rsys,generateBufferingModules->False]

dnasys=rxnsysToDnasys[rsys,generateBufferingModules->False];

sol=SimulateRxnsys[dnasys,tmax];

plotter=
#[t]/.sol&/@{i,d};

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

What fraction of auxiliary species was used up by the end?

1-(#[tmax]/.sol&/@SpeciesInRxnsys[dnasys,gateT[_]])/cmax

Scale system into realistic regime:

{cscale,tscale}=computeConcAndTimeScaling[10^-5,10^6];

plotterscaled=cscale(#/.t->(tscale^(-1))t)&/@plotter;

Plot[plotterscaled,{t,0,tscale tmax},PlotPoints->300,PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue}}]

DNA implementation, with buffering modules

$HistoryLength=0;ClearSystemCache[];

cmax=20000;

qmax=computeQmaxLowerBound[rsys]

dnasys=rxnsysToDnasys[rsys];

sol=SimulateRxnsys[dnasys,tmax];

plotter=
#[t]/.sol&/@{i,d};

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

What fraction of auxiliary species was used up by the end?

1-(#[tmax]/.sol&/@SpeciesInRxnsys[dnasys,gateT[_]])/cmax

Scale system into realistic regime:

{cscale,tscale}=computeConcAndTimeScaling[10^-5,10^6];

plotterscaled=cscale(#/.t->(tscale^(-1))t)&/@plotter;

Plot[plotterscaled,{t,0,tscale tmax},PlotPoints->300,PlotRange->{0,All},PlotStyle->{{Red},{Green},{Blue}}]