(* ::Package:: *) (* ::Section:: *) (*CRNSimulator*) (* ::Text:: *) (*Chemical Reaction Network (CRN) Simulator package is developed by David Soloveichik. Copyright 2009-2012. *) (*http://www.dna.caltech.edu/~davids/*) (* ::Subsection::Closed:: *) (*Public interface specification*) Notation`AutoLoadNotationPalette = False; Needs["Notation`"] BeginPackage["CRNSimulator`", {"Notation`"}]; rxn::usage="Represents an irreversible reaction. eg. rxn[a+b,c,1]"; revrxn::usage="Represents a reversible reaction. eg. revrxn[a+b,c,1,1]"; conc::usage="Initial concentration: conc[x,10] or conc[{x,y},10]."; term::usage="Represents an additive term in the ODE for species x. \ Species concentrations must be expressed in x[t] form. eg. term[x, -2 x[t]*y[t]]"; SimulateRxnsys::usage= "SimulateRxnsys[rxnsys,endtime] simulates the reaction system rxnsys for time 0 \ to endtime. In rxnsys, initial concentrations are specified by conc statements. \ If no initial condition is set for a species, its initial concentration is set to 0. \ Rxnsys can include term[] statements (eg. [x, -2 x[t]]) which are additively combined \ together with term[]s derived from rxn[] statements. \ Rxnsys can also include ODEs directly for some species (Eg x'[t]==...) that are \ passed on to NDSolve without modification. Any options specified (eg WorkingPrecision->30) \ are passed to NDSolve."; SpeciesInRxnsys::usage= "SpeciesInRxnsys[rxnsys] returns the species in reaction system rxnsys. \ SpeciesInRxnsys[rxnsys,pttrn] returns the species in reaction system rxnsys \ matching Mathematica pattern pttrn (eg x[1,_])."; SpeciesInRxnsysStringPattern::usage= "SpeciesInRxnsysPattern[rxnsys,pttrn] returns the species in reaction system rxnsys \ matching Mathematica string pattern pttrn. \ (Eg \"g$*\" matches all species names starting with \"g$\" ; \ can also do RegularExpression[\"o..d.\$.*\"].)"; RxnsysToOdesys::usage= "RxnsysToOdesys[rxnsys,t] returns the ODEs corresponding to reaction system rxnsys, \ with initial conditions. If no initial condition is set for a species, its initial \ concentration is set to 0. \ The time variable is given as the second argument; if omitted it is set to Global`t."; (*To use instead of Sequence in functions with Hold attribute but not HoldSequence, like Module, If, etc*) Seq:=Sequence (* ::Subsection::Closed:: *) (*Private*) Begin["`Private`"]; Notation[ParsedBoxWrapper[RowBox[{"r_", " ", OverscriptBox["\[RightArrow]", RowBox[{" ", "k_", " "}]], " ", "p_", " "}]] \[DoubleLongLeftRightArrow] ParsedBoxWrapper[RowBox[{"rxn", "[", RowBox[{"r_", ",", "p_", ",", "k_"}], "]"}]]] Notation[ParsedBoxWrapper[RowBox[{"r_", " ", UnderoverscriptBox["\[RightArrowLeftArrow]", "k2_", "k1_"], " ", "p_", " "}]] \[DoubleLongLeftRightArrow] ParsedBoxWrapper[RowBox[{"revrxn", "[", RowBox[{"r_", ",", "p_", ",", "k1_", ",", "k2_"}], "]"}]]] AddInputAlias[ParsedBoxWrapper[RowBox[{"\[Placeholder]", " ", OverscriptBox["\[RightArrow]", RowBox[{" ", "\[Placeholder]", " "}]], " ", "\[Placeholder]", " "}]],"rxn"] AddInputAlias[ParsedBoxWrapper[RowBox[{"\[Placeholder]", " ", UnderoverscriptBox["\[RightArrowLeftArrow]", "\[Placeholder]", "\[Placeholder]"], " ", "\[Placeholder]", " "}]],"revrxn"] (*We want rxn[a+b,c,1] to be different from rxn[b+a,c,1], so we have to set attribute HoldAll. But we also want to evaluate if any variables can be evaluated.*) SetAttributes[{rxn,revrxn}, HoldAll] rxn[rs_Plus,ps_,k_]:= ReleaseHold[ReplacePart[rxn[1,ps,k],1->Hold[Plus]@@List@@Unevaluated[rs]]]/; Hold@@Unevaluated[rs] =!= Hold@@List@@Unevaluated[rs] rxn[rs_,ps_Plus,k_]:= ReleaseHold[ReplacePart[rxn[rs,1,k],2->Hold[Plus]@@List@@Unevaluated[ps]]]/; Hold@@Unevaluated[ps] =!= Hold@@List@@Unevaluated[ps] rxn[rs_,ps_,k_]:= (With[{rse=rs},rxn[rse,ps,k]])/;Head[rs]=!=Plus&&Unevaluated[rs]=!=rs rxn[rs_,ps_,k_]:= (With[{pse=ps},rxn[rs,pse,k]])/;Head[ps]=!=Plus&&Unevaluated[ps]=!=ps rxn[rs_,ps_,k_]:= (With[{ke=k},rxn[rs,ps,ke]])/;Unevaluated[k]=!=k revrxn[r_,p_,k1_,k2_]:=Sequence[rxn[r,p,k1],rxn[p,r,k2]] conc[xs_List,c_]:=Seq@@(conc[#,c]&/@xs) (* Species as products or reactants in rxn[] statements, as well as defined in x'[t]== or x[t]== statements \ or term statements, or conc statements *) SpeciesInRxnsys[rxnsys_]:= Union[ Cases[Cases[rxnsys,rxn[r_,p_,_]:>Seq[r,p]]/.Times|Plus->Seq,s_Symbol|s_Symbol[__]], Cases[rxnsys, x_'[_]==_ | x_[_]==_ | term[x_,__] | conc[x_,_] :> x]] SpeciesInRxnsys[rsys_,pattern_]:=Cases[SpeciesInRxnsys[rsys],pattern] SpeciesInRxnsysStringPattern[rsys_,pattern_]:=Select[SpeciesInRxnsys[rsys],StringMatchQ[ToString[#],pattern]&] (* Check if a species' initial value is set in a odesys *) InitialValueSetQ[odesys_,x_]:= MemberQ[odesys,x[_]==_] (* Check if a species is missing an ODE or a direct definition (x[t]=_) in odesys. *) MissingODEQ[odesys_,x_,t_Symbol]:= !MemberQ[odesys, D[x[t],t]==_ | x[t]==_] RxnsysToOdesys[rxnsys_,t_Symbol:Global`t]:= Module[ {spcs=SpeciesInRxnsys[rxnsys], termssys, odesys, eqsFromTerms, eqsFromConcs}, (* Convert rxn[] to term[] statements *) termssys=rxnsys /. rxn:rxn[__]:>Seq@@ProcessRxnToTerms[rxn,t]; (* ODEs from parsing terms *) eqsFromTerms = ProcessTermsToOdes[Cases[termssys,term[__]],t]; (* initial values from parsing conc statements *) eqsFromConcs = Cases[termssys,conc[x_,c_]:>x[0]==c]; (* Remove term and conc statements from rxnsys and add eqs generated from them. If there is a conflict, use pass-through equations *) odesys = DeleteCases[termssys, term[__]|conc[__]]; odesys = Join[odesys, DeleteCases[eqsFromTerms, Alternatives@@(#'[t]==_& /@ Cases[odesys,(x_'[t]|x_[t])==_:>x])], eqsFromConcs]; (* For species still without initial values, add zeros *) odesys = Join[odesys, #[0]==0& /@ Select[spcs, !InitialValueSetQ[odesys,#]&]]; (* For species still without ODE or direct definition, add zero time derivative *) (* This can happen for example if conc is defined, but nothing else *) Join[odesys, D[#[t],t]==0& /@ Select[spcs, MissingODEQ[odesys,#,t]&]] ] (* Create list of ODEs from parsing term statements. terms should be list of term[] statements. *) ProcessTermsToOdes[terms_,t_Symbol]:= Module[{spcs=Union[Cases[terms,term[s_,_]:>s]]}, #'[t]==Total[Cases[terms,term[#,rate_]:>rate]] & /@ spcs]; (* Create list of term[] statements from parsing a rxn statement *) ProcessRxnToTerms[reaction:rxn[r_,p_,k_],t_Symbol]:= Module[{spcs=SpeciesInRxnsys[{reaction}], rrate, spccoeffs,terms}, (* compute rate of this reaction *) rrate = k (r/.{Times[b_,s_]:>s^b,Plus->Times}); (*for each species, get a net coefficient*) spccoeffs=Coefficient[p-r,#]& /@ spcs; (*create term for each species*) terms=MapThread[term[#1,#2*rrate]&,{spcs, spccoeffs}]; (*change all species variables in the second arg in term[] to be functions of t*) terms/.term[spc_,rate_]:>term[spc,rate/.s_/;MemberQ[spcs,s]:>s[t]]]; SimulateRxnsys[rxnsys_,endtime_,opts:OptionsPattern[NDSolve]]:= Module[{spcs=SpeciesInRxnsys[rxnsys],odesys=RxnsysToOdesys[rxnsys,Global`t]}, Quiet[NDSolve[odesys, spcs, {Global`t,0,endtime},opts,MaxSteps->Infinity,AccuracyGoal->MachinePrecision],{NDSolve::"precw"}][[1]]] End[]; EndPackage[]; (* ::Section:: *) (*seesawSimulatorDM*) (* ::Text:: *) (*Seesaw Simulator DM (i.e. with Defective Molecules) package is developed by Chris Thachuk and Lulu Qian. *) (* ::Subsection:: *) (*Define seesaw parameters*) (* ::Text:: *) (*Rate constants:*) kf = 2*10^6; (* fast strand displacement rate with extended toehold, unit: /M/s *) ks = 5*10^4; (* slow strand displacement rate with universal toehold, unit: /M/s *) krf = 26; (* fast strand dissociation rate with universal toehold, unit: /s *) krs = 1.3; (* slow strand dissociation rate with extended toehold, unit: /s *) kl = 10; (* strand displacement leak rate, unit: /M/s *) (* ::Text:: *) (*Standard concentration 1x:*) c = 100*10^-9; (* unit: M *) (* ::Text:: *) (*Synthesis and stoichiometry errors:*) srb = 0.01; (* synthesis error per base *) cr = 0; (* stoichiometry error per gate *) (* ::Text:: *) (*Global concentration adjustment:*) ThtoSig = 1.4; (* \[Beta]/\[Alpha] concentration adjustment for thresholds in all logic gates *) ThInc = 0.2; (* \[Delta]\[Times]\[Alpha]/\[Beta] additional concentration adjustment for thresholds in logic gates at layer two and above *) (* ::Subsection:: *) (*Define seesaw functions*) (* ::Text:: *) (*Some helper functions to calculate the probability of generating molecules with errors:*) (* Calculates the probability of synthesizing a strand of n number of errors, of length l, and synthesis error per base r. *) probStrand[r_, l_, n_]:=Binomial[l, n]*(1-r)^(l-n)*r^n; (* Calculates the probability of synthesizing a signal strand with s number of errors in any long domain, and t number of errors in the middle toehold. The parameter r is the synthesis error rate per base. *) probSignal[r_, s_, t_] := probStrand[r, 15, s] * probStrand[r, 5, t] * probStrand[r, 15, 0]; (* Calculates the probability of synthesizing a gate. We care about every domain except the bottom long domain -- we assume it is always synthesized correctly. There are therefore 3 toehold domains of interest, and 2 long domains of interest. The parameter r is the synthesis error rate per base, t is number of errors in any toehold domain, and s is the number of errors in any long domain. We make the simplifying assumption that all errors of a certain type are within the same domain. (This is fine when at most 1 error can occur, as in our model. *) probGate[r_, s_, t_]:=probStrand[r, 15, s]*probStrand[r, 5, t]*probStrand[r, 25, 0]; (* Calculates the probability for a threshold gate. We only care if there is 0 or more errors in the toehold region. We assume errors elsewhere make no difference, so they are ignored in the calculation. *) probThreshold[r_, t_]:=probStrand[r, 10, t]; (* Calculates the maximal signal that can be released from a gate, with r being the synthesis error rate per base and cr being the stoichiometry error per gate. The default value of cr is 0. *) maxSignal[r_,cr_:0]:=probGate[r,0,0]+2*probGate[r,1,0]+2*probGate[r,0,1]+cr; (* ::Text:: *) (*Converts a seesaw signal, gate or threshold molecule into a list of molecules including defected ones:*) concd[name_,c_]:= Switch[name[[0]], w,{ conc[name,probSignal[srb,0,0]*c], conc[ReplacePart[name,0->wd3S],probSignal[srb,1,0]*c], conc[ReplacePart[name,0->wdT],probSignal[srb,0,1]*c], conc[ReplacePart[name,0->wd5S],probSignal[srb,1,0]*c]}, g,{ conc[name,probGate[srb,0,0]*c], conc[ReplacePart[name,{2,0}->wd3S],probGate[srb,1,0]*c], conc[ReplacePart[name,{2,0}->wdT],probGate[srb,0,1]*c], conc[ReplacePart[name,{2,0}->wd5S],probGate[srb,1,0]*c], conc[ReplacePart[name,0->gd5T],probGate[srb,0,1]*c], conc[ReplacePart[name,0->gd3T],probGate[srb,0,1]*c], (* We assume that if there is an excess of signal strands in the gate molecules, due to competition, the free-floating signal strands would have synthesis errors either in the toehold or in the 3' end long domain *) conc[ReplacePart[name[[2]],0->wd3S],cr*probSignal[srb,1,0]/(probSignal[srb,1,0]+probSignal[srb,0,1])*c], conc[ReplacePart[name[[2]],0->wdT],cr*probSignal[srb,0,1]/(probSignal[srb,1,0]+probSignal[srb,0,1])*c]}, th,{ conc[name,probThreshold[srb,0]*c], conc[ReplacePart[name,0->thdT],probThreshold[srb,1]*c]} ]/.List->Sequence (* ::Text:: *) (*Translates a seesaw gate into a list of reactions:*) seesaw[x_,l_List,r_List]:={ (* Toehold exchange reactions *) Outer[revrxn[#3[#1,x]+#4[x,#5[x,#2]],#4[#3[#1,x],x]+#5[x,#2],ks,ks]&,l,r,{w,wd3S},{g},{w,wd5S}], Outer[revrxn[#3[#1,x]+#4[x,#5[x,#2]],#4[#3[#1,x],x]+#5[x,#2],ks/100,ks]&,l,r,{wdT,wd5S},{g},{w,wd5S}], Outer[revrxn[#3[#1,x]+#4[x,#5[x,#2]],#4[#3[#1,x],x]+#5[x,#2],ks/100,ks]&,l,r,{w,wd3S},{gd5T},{w,wd5S}], Outer[revrxn[#3[#1,x]+#4[x,#5[x,#2]],#4[#3[#1,x],x]+#5[x,#2],ks,ks/100]&,l,r,{w,wd3S},{g},{wd3S,wdT}], Outer[revrxn[#3[#1,x]+#4[x,#5[x,#2]],#4[#3[#1,x],x]+#5[x,#2],ks,ks/100]&,l,r,{w,wd3S},{gd3T},{w,wd5S}], (* Thresholding reactions *) Outer[rxn[#2[#1,x]+#3[w[#1,x],x],waste,kf]&,l,{w,wd3S,wd5S},{th}], Outer[rxn[#2[#1,x]+#3[w[#1,x],x],waste,kf/100]&,l,{wdT},{th}], Outer[rxn[#2[#1,x]+#3[w[#1,x],x],waste,kf/100]&,l,{w,wd3S,wd5S},{thdT}], Outer[rxn[#2[x,#1]+#3[x,w[x,#1]],waste,kf]&,r,{w,wd3S,wd5S},{th}], Outer[rxn[#2[x,#1]+#3[x,w[x,#1]],waste,kf/100]&,r,{wdT},{th}], Outer[rxn[#2[x,#1]+#3[x,w[x,#1]],waste,kf/100]&,r,{w,wd3S,wd5S},{thdT}], (* Universal toehold binding reactions *) Outer[revrxn[#2[#3[#1,x],x]+W,#2[#3[#1,x],x,W],kf,krf]&,l,{g,gd5T},{w,wd3S,wdT,wd5S}], Outer[revrxn[#2[#3[#1,x],x]+W,#2[#3[#1,x],x,W],kf,krf*10]&,l,{gd3T},{w,wd3S,wdT,wd5S}], Outer[revrxn[#2[x,#3[x,#1]]+W,#2[W,x,#3[x,#1]],kf,krf]&,r,{g,gd3T},{w,wd3S,wdT,wd5S}], Outer[revrxn[#2[x,#3[x,#1]]+W,#2[W,x,#3[x,#1]],kf,krf*10]&,r,{gd5T},{w,wd3S,wdT,wd5S}], Outer[revrxn[#2[w[#1,x],x]+W,#2[W,w[#1,x],x],kf,krs]&,l,{th}], Outer[revrxn[#2[w[#1,x],x]+W,#2[W,w[#1,x],x],kf,krs*10]&,l,{thdT}], Outer[revrxn[#2[x,w[x,#1]]+W,#2[x,w[x,#1],W],kf,krs]&,r,{th}], Outer[revrxn[#2[x,w[x,#1]]+W,#2[x,w[x,#1],W],kf,krs*10]&,r,{thdT}], Outer[revrxn[#2[#1,x]+G,#2[#1,x,G],kf,krf]&,l,{w,wd3S,wd5S}], Outer[revrxn[#2[#1,x]+G,#2[#1,x,G],kf,krf*10]&,l,{wdT}], Outer[revrxn[#2[#1,x]+TH,#2[#1,x,TH],kf,krs]&,l,{w,wd3S,wd5S}], Outer[revrxn[#2[#1,x]+TH,#2[#1,x,TH],kf,krs*10]&,l,{wdT}], (* Leak reactions *) Outer[rxn[#3[#4[#1,x],x]+#5[#2,x],#3[#5[#2,x],x]+#4[#1,x],kl]&,l,l,{g,gd5T,gd3T},{w,wd3S},{w,wd3S}], Outer[rxn[#3[#4[#1,x],x]+#5[#2,x],#3[#5[#2,x],x]+#4[#1,x],kl*2]&,l,l,{g,gd5T,gd3T},{wdT,wd5S},{w,wd3S}], Outer[rxn[#3[x,#4[x,#1]]+#5[x,#2],#3[x,#5[x,#2]]+#4[x,#1],kl]&,r,r,{g,gd5T,gd3T},{w,wd5S},{w,wd5S}], Outer[rxn[#3[x,#4[x,#1]]+#5[x,#2],#3[x,#5[x,#2]]+#4[x,#1],kl*2]&,r,r,{g,gd5T,gd3T},{wd3S,wdT},{w,wd5S}] }/.List->Sequence (* ::Text:: *) (*Translates a reporter into a list of reactions and the initial concentration of the reporter:*) reporter[x_,l_]:={ (* reporting reactions *) Outer[rxn[#[l,x]+Rep[x],Fluor[x],2*ks]&,{w,wd3S,wd5S}], Outer[rxn[#[l,x]+Rep[x],Fluor[x],2*ks/100]&,{wdT}], (* Universal toehold binding reactions *) revrxn[Rep[x]+W, Rep[x,W],kf,krf], Outer[revrxn[#[l,x]+G,#[l,x,G],kf,krf]&,{w,wd3S,wd5S}], Outer[revrxn[#[l,x]+G,#[l,x,G],kf,krf*10]&,{wdT}], Outer[revrxn[#[l,x]+TH,#[l,x,TH],kf,krs]&,{w,wd3S,wd5S}], Outer[revrxn[#[l,x]+TH,#[l,x,TH],kf,krs*10]&,{wdT}], conc[Rep[x],1.5*c] }/.List->Sequence (* ::Text:: *) (*Translates logic OR operation into a list of seesaw gates and the concentrations of initial species:*) (* m is the number of logic circuit layers upstream of the logic gate. *) seesawOR[x1_,x2_,l_List,r_List,m_]:= Module[{f}, {seesaw[x1,l,{x2}], seesaw[x2,{x1},{r,f}], concd[g[x1,w[x1,x2]],Length[l]*c], Outer[concd[g[x2,w[x2,#]],1*c]&,r], concd[w[x2,f],2*Length[r]*c], concd[th[w[x1,x2],x2],ThtoSig*(0.35+(Length[l]-2)*0.05+If[m>0,ThInc,0])*c] }]/.List->Sequence (* ::Text:: *) (*Translates logic AND operation into a list of seesaw gates and the concentrations of initial species:*) (* m is the number of logic circuit layers upstream of the logic gate. *) seesawAND[x1_,x2_,l_List,r_List,m_]:= Module[{f}, {seesaw[x1,l,{x2}], seesaw[x2,{x1},{r,f}], concd[g[x1,w[x1,x2]],Length[l]*c], Outer[concd[g[x2,w[x2,#]],1*c]&,r], concd[w[x2,f],2*Length[r]*c], concd[th[w[x1,x2],x2],ThtoSig*((Length[l]-1+0.2)/ThtoSig+If[m>0,ThInc,0])*c] }]/.List->Sequence (* ::Text:: *) (*Translates inputs with fan-out into a list of seesaw gates and the concentrations of initial species:*) inputfanout[x_,l_,r_List]:= Module[{f}, {seesaw[x,{l},{r,f}], Outer[concd[g[x,w[x,#]],1*c]&,r], concd[w[x,f],2*Length[r]*c], concd[th[w[l,x],x],ThtoSig*0.2*c] }]/.List->Sequence