Needs CRNSimulator.m

The idea here is to automatically transform a polynomial ODE system, expressed in the symbolic notation used for example by Mathematica’s NDSolve routine, into an equivalent CRN system.  This involves using a few Mathematica functions that extract parts of symbolic expressions (e.g. MonomialList) and using symbolic rewriting rules to modify and create variables, e.g. going from the ODE time-varying signal x[t] to the dual-rail CRN species x[pos] and x[neg] via “x[t] /. v_[t] -> v[pos]” and similar expressions.

The main intended lesson is that EVERYONE can write a “compiler” for building CRNs of a certain type, from a certain representation, and it’s “easy”.

The ODE-to-CRN “compiler” definitions

PolyFor[v_,ode_]:=Cases[ode,(v'[t]==e_):> (e/.x_[t]->x)] //First

DualRailPolyFor[v_,ode_,vars_]:=With[{r=(#->(#[pos]-#[neg])&)/@vars},
PolyFor[v,ode]/.r
]

PositiveMonomialsFor[v_,ode_,vars_]:=With[{r=(#[_]->1&)/@vars},
Select[ MonomialList[DualRailPolyFor[v,ode,vars]] , (#/.r)>0&]
]

NegativeMonomialsFor[v_,ode_,vars_]:=With[{r=(#[_]->1&)/@vars},
-Select[ MonomialList[DualRailPolyFor[v,ode,vars]] , (#/.r)<0&]
]

Clear[Catalysts]
Catalysts[_?NumericQ]:=0
Catalysts[Times[_?NumericQ,v__]]:=Plus[v]/.{Power->Times}
Catalysts[Times[v1_,v__]]:=Plus[v1,v]/.{Power->Times}
Catalysts[v_]:=v/.{Power->Times}

ReactionForMonomial[v_,mono_,sign_]:=With[{cat=Catalysts[mono],k=mono/.{_[_]->1}},
Evaluate/@rxn[cat,cat+v[sign],k]
]

ReactionsFor[v_,ode_,vars_]:=Join[
ReactionForMonomial[v,#,pos]& /@ PositiveMonomialsFor[v,ode,vars],
ReactionForMonomial[v,#,neg]& /@ NegativeMonomialsFor[v,ode,vars]
]

ReactionsForODE[ode_,vars_,kfast_]:=Join[
Union@@(ReactionsFor[#,ode,vars]&/@vars),
rxn[#[pos]+#[neg],0,kfast]&/@vars
]

Example runs of the ODE-to-CRN “compiler”.

(* Rossler attractor *)
Rossler={
x'[t]==-z[t]-y[t],
y'[t]==x[t]+0.2y[t], 
z'[t]==0.2+z[t](x[t]-5.7)
}

(* first simulate the Rossler attractor dynamics directly using numerical integration *)
tmax=200;
sol=NDSolve[Join[Rossler,{x[0]==0,y[0]==0,z[0]==0}],{x,y,z},{t,0,tmax}]

Plot[Evaluate[{x[t],y[t],z[t]}/.sol],{t,0,tmax},PlotRange->{{0,tmax},{-15,+25}},PlotStyle->{{Red},{Green},{Blue}}]

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol],{t,0,tmax},PlotRange->{{-15,+15},{-15,+15},{0,25}}]

rsys=ReactionsForODE[Rossler,{x,y,z},10]

tmax=200;
sol=SimulateRxnsys[rsys,tmax];
Plot[Evaluate[{x[pos][t]-x[neg][t],y[pos][t]-y[neg][t],z[pos][t]-z[neg][t]}/.sol],{t,0,tmax},PlotRange->{{0,tmax},{-15,+25}},PlotStyle->{{Red},{Green},{Blue}}]

ParametricPlot3D[Evaluate[{x[pos][t]-x[neg][t],y[pos][t]-y[neg][t],z[pos][t]-z[neg][t]}/.sol],{t,0,tmax},PlotRange->{{-15,+15},{-15,+15},{0,25}}]

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

rsys=ReactionsForODE[Brusselator,{x,y},100]

tmax=45;
sol=SimulateRxnsys[Join[rsys,{conc[x[pos],1],conc[y[pos],1]}],tmax];
Plot[Evaluate[{x[pos][t]-x[neg][t],y[pos][t]-y[neg][t]}/.sol],{t,0,tmax},PlotRange->{{0,tmax},{0,5}},PlotStyle->{{Red},{Green}}]

(* Lorenz attractor *)
Lorenz={
x'[t]==10(y[t]-x[t]),
y'[t]==28x[t]-y[t]-x[t]z[t],
z'[t]==x[t]y[t]-8/3 z[t]
}

rsys=ReactionsForODE[Lorenz,{x,y,z},100]

tmax=30;
sol=SimulateRxnsys[Join[rsys,{conc[x[pos],1],conc[y[pos],1],conc[z[pos],1]}],tmax];
Plot[Evaluate[{x[pos][t]-x[neg][t],y[pos][t]-y[neg][t],z[pos][t]-z[neg][t]}/.sol],{t,0,tmax},PlotRange->{{0,tmax},{-30,60}},PlotStyle->{{Red},{Green},{Blue}}]

ParametricPlot3D[Evaluate[{x[pos][t]-x[neg][t],y[pos][t]-y[neg][t],z[pos][t]-z[neg][t]}/.sol],{t,0,tmax},PlotRange->{{-30,+60},{-30,+60},{-30,+60}}]

Examples of testing the functions defined in the ODE-to-CRN “compiler”

MonomialList[0.2+2x+4.5(x+3)(z-2.3)+a]

PolyFor[z,Rossler]

DualRailPolyFor[z,Rossler,{x,y,z}] //Expand

PositiveMonomialsFor[z,Rossler,{x,y,z}]

Catalysts /@ NegativeMonomialsFor[z,Rossler,{x,y,z}]

ReactionForMonomial[z,3.3x[pos]y[neg],neg]