Example of modular construction of digital logic circuits

Note: when going through the examples below, I encourage you to try out variations and test functions and code snippets in different ways yourself -- just start a new notebook cell above or below the given code.  Also note that hovering your mouse over a defined function should pull up an option to get more information about that function.  

Remember that CRNSimulator.m must be run before attempting to use this notebook.

For a wire signal X, the CRN will have species X[0] and X[1], with conserved total concentration m.
To construct a generic CRN subcircuit, we will define a function that is given the input and output wire signal symbols and the desired Boolean function, and that returns a list of the four reactions needed for an analog “approximate” implementation of that function.  (Later on, we will define a restoring gate CRN, and compose the two to create a CRN that satisfies the digital abstraction.)

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

Good practice for using Mathematica is to test small components and verify that they are doing what you want, then automate the process on larger examples.  Here we build an AND gate, using the function Times to indicate the output should be 1 if and only if both inputs are 1.

analogANDgate = analogGate[X,Y,Z,Times]

It will be helpful later to get a list of all the wire symbols, because these indicate which species will need to have initial concentrations specified.

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

Wires[analogANDgate]

To produce a list of the initial species concentrations, the user provides the circuit (from which the full list of wires is extracted), a list of {wire, value} pairs, and the maximal signal value m.  For each wire X, the function produces a concentration assignment for species X[0] and species X[1], either according to the input specification or “half-on” for the unspecified wires.

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@@({conc[#[0],m/2],conc[#[1],m/2]}&/@Complement[wires,inputwires]);
Union[inputconcs,otherconcs]
]

InitialConcentrations[analogANDgate,{{X,2.5},{Y,0.5}},3]

A simulatable reaction system consists of a list of reactions as well as initial concentrations.  So here is a test case.

rsys= Join[analogANDgate, InitialConcentrations[analogANDgate,{{X,2.5},{Y,0.5}},3]
]

Now we simulate the example AND gate, and plot the results, showing all six species for the three wires.  Note that the “Evaluate” statement is necessary because the “Plot” function uses “Hold” on its first argument -- that means that the first argument is passed as-is, symbolically, to the function.  So, without “Evaluate”, “Plot” just sees the symbol “plotter”, and the automatic routine for generating the legend doesn’t know what to do.

tmax=2;
sol=SimulateRxnsys[rsys,tmax];plotter=#[t]/.sol&/@SpeciesInRxnsys[rsys];
Plot[Evaluate[plotter],{t,0,tmax},PlotLegends->SpeciesInRxnsys[rsys]]

It is straightforward to run this simulation with many possible inputs.  “Table” produces a matrix of {X,Y,Z} triples, scanning across values for X and Y.  The “Flatten” statement turns the matrix into a list, as needed for ListContourPlot.

inputscan=With[{rxns=analogANDgate},
Flatten[Table[
rsys= Join[rxns, InitialConcentrations[rxns,{{X,xinit},{Y,yinit}},3]];
tmax=2;
sol=SimulateRxnsys[rsys,tmax];
{xinit,yinit,Z[1][tmax]/.sol},
{xinit,0,3,0.1},{yinit,0,3,0.1}
],1]];

ListContourPlot[inputscan,PlotLegends->Automatic]

It is easy to see that the above analog AND gate does not satisfy the digital abstraction for a margin of 0.5:  for inputs in the range [2.5,3.0]x[2.5,3.0], only about half of the outputs are greater than 2.5.

It is easier to see if we draw the regions relevant to the digital abstraction, and make contour lines exactly where the output is within valid ranges.  Note that it is better to compute the contours analytically; our purpose here is to demonstrate how to automate CRN simulations over many inputs to characterize the behavior of a circuit, which will be more valuable when the circuit is complex enough that no analytic solution is in hand.

With[{eps=0.5,m=3},Show[
ListContourPlot[inputscan,Contours->{eps,m-eps},PlotLegends->Automatic],
Graphics[Line[{{0,0},{#,0},{#,#},{0,#},{0,0}}]&@eps],
Graphics[Line[{{m,0},{m-#,0},{m-#,#},{m,#},{m,0}}]&@eps],
Graphics[Line[{{0,m},{#,m},{#,m-#},{0,m-#},{0,m}}]&@eps],
Graphics[Line[{{m,m},{m-#,m},{m-#,m-#},{m,m-#},{m,m}}]&@eps]
]]

Next, we will add a restoring repeater subcircuit, and show that a digital abstraction is obtained.  First just look at the repeater itself.

repeaterGate[input_,output_]:= {
rxn[input[0]+input[0]+output[1],input[0]+input[0]+output[0],1],
rxn[input[1]+input[1]+output[0],input[1]+input[1]+output[1],1]
}

inputscan=With[{rxns=repeaterGate[X,Z]},
Table[
rsys= Join[rxns, InitialConcentrations[rxns,{{X,xinit}},3]];
tmax=2;
sol=SimulateRxnsys[rsys,tmax];
{xinit,Z[1][tmax]/.sol},
{xinit,0,3,0.1}
]];

Show[ListPlot[inputscan,Joined->True],Graphics[{Dashed,Line[{{0,0},{3,3}}]}]]

To wire together an analog gate and a repeater, we need a new wire that is the output of the analog gate and the input to the repeater.  Assuming that the user doesn’t care about the name of that wire, we use Mathematica’s “Module” function, which produces a unique name every time the function is invoked.  (If you think that’s too ugly, a good exercise would be to redefine this function to let the user choose the name of the intermediate wire.)

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

digitalANDgate = digitalGate[X,Y,Z,Times]

inputscan=With[{rxns=digitalANDgate},
Flatten[Table[
rsys= Join[rxns, InitialConcentrations[rxns,{{X,xinit},{Y,yinit}},3]];
tmax=2;
sol=SimulateRxnsys[rsys,tmax];
{xinit,yinit,Z[1][tmax]/.sol},
{xinit,0,3,0.1},{yinit,0,3,0.1}
],1]];

Now we see that this six-reaction CRN satisfies the digital abstraction for margins less than about 0.5.

With[{eps=0.5,m=3},Show[
ListContourPlot[inputscan,Contours->{eps,m-eps},PlotLegends->Automatic],
Graphics[Line[{{0,0},{#,0},{#,#},{0,#},{0,0}}]&@eps],
Graphics[Line[{{m,0},{m-#,0},{m-#,#},{m,#},{m,0}}]&@eps],
Graphics[Line[{{0,m},{#,m},{#,m-#},{0,m-#},{0,m}}]&@eps],
Graphics[Line[{{m,m},{m-#,m},{m-#,m-#},{m,m-#},{m,m}}]&@eps]
]]

For building larger circuits, it will be convenient to name a few I/O functions. (Try adding a few others yourself.)

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

The following is a digital circuit that computes the (rounded-down) 2-bit square root of a 4-bit binary number, where input X4 is the most significant input bit and Y2 is the most significant output bit.

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

inputcases=With[{rxns=SqrtCircuit,m=3,eps=0.4},
Flatten[Table[
rsys= Join[rxns, 
InitialConcentrations[rxns,{{X1,x1},{X2,x2},{X3,x3},{X4,x4}},m]];
tmax=10;
sol=SimulateRxnsys[rsys,tmax];
{Round[x4/m],Round[x3/m],Round[x2/m],Round[x1/m],
Round[Y2[1][tmax]/m/.sol,0.01],Round[Y1[1][tmax]/m/.sol,0.01]},
{x4,{eps,m-eps}},
{x3,{eps,m-eps}},
{x2,{eps,m-eps}},
{x1,{eps,m-eps}}],3]];

inputcases//MatrixForm

Was all that fuss and bother worth it?  Let’s see how well the “same” circuit would work if we had not used the repeater.  When inputs are off by up to 0.4 out of 3.0, i.e 13%, the output is off by as much as 34%.  While one might argue that the correct True/False distinction is still discernible in the output of the analog circuit, but it is clear that this approach will not scale.

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

inputcases=With[{rxns=analogSqrtCircuit,m=3,eps=0.4},
Flatten[Table[
rsys= Join[rxns, 
InitialConcentrations[rxns,{{X1,x1},{X2,x2},{X3,x3},{X4,x4}},m]];
tmax=10;
sol=SimulateRxnsys[rsys,tmax];
{Round[x4/m],Round[x3/m],Round[x2/m],Round[x1/m],
Round[Y2[1][tmax]/m/.sol,0.01],Round[Y1[1][tmax]/m/.sol,0.01]},
{x4,{eps,m-eps}},
{x3,{eps,m-eps}},
{x2,{eps,m-eps}},
{x1,{eps,m-eps}}],3]];

inputcases//MatrixForm