Note: requires seesawSimulatorDM.wl to be loaded.

Functions

PlotSimInput function plots a simulation including multiple trajectories of the same output, each has a distinct input.
SIMcircuit is a simulation. 
input is a list of input values, relative to a standard concentration 1×.
time defines the range of time plotted, starting from 0. The unit is hours.
InputLabel is a label for the legend of inputs. The default is none. 
CircuitLabel is a label for the plot. It can be the function of the circuit. The default is none. 
size is the size of image. The default is 400.

PlotSimInput[SIMcircuit_,input_,time_,InputLabel_:"",CircuitLabel_:"",size_:400]:=
Plot[SIMcircuit,{t,0,time},
PlotLabel->Style[CircuitLabel,24],
Frame->True,FrameLabel->{"Time (hours)","Output"},
PlotStyle->Table[{AbsoluteThickness[3],ColorData["Rainbow"][i]},{i,If[Length[input]<8,0.12,0],0.96,0.96/Length[input]}],
PlotLegends->SwatchLegend[Automatic,Table[Row[Table[input[[i,j]],{j,1,Length[input[[i]]]}]],{i,1,Length[input]}],LegendLabel->InputLabel,LegendMarkerSize->14],
LabelStyle->Directive[Gray,FontSize->20,FontFamily->"Helvetica"],
GridLines->Automatic,
PlotRange->{All,{-0.05,1.05}},AspectRatio->1/1.3,ImageSize->size]

PlotSimOutput function plots a simulation including multiple trajectories of distinct outputs, with the same input.
SIMcircuit is a simulation. 
output is a list of output labels.
time defines the range of time plotted, starting from 0. The unit is hours.
OutputLabel is a label for the legend of outputs. The default is none. 
CircuitLabel is a label for the plot. It can be the function of the circuit. The default is none. 
size is the size of image. The default is 400.

PlotSimOutput[SIMcircuit_,output_,time_,OutputLabel_:"",CircuitLabel_:"",size_:400]:=
Plot[SIMcircuit,{t,0,time},
PlotLabel->Style[CircuitLabel,20],
Frame->True,FrameLabel->{"Time (hours)","Output"},
PlotStyle->Table[{AbsoluteThickness[3],ColorData["FallColors"][Floor[i/2]/(Length[output]/2-1)],If[EvenQ[i],Dashed]},{i,0,Length[output]-1}],
PlotLegends->LineLegend[Automatic,output,LegendLabel->OutputLabel,LegendMarkerSize->14],
LabelStyle->Directive[Gray,FontSize->20,FontFamily->"Helvetica"],
GridLines->Automatic,
PlotRange->{All,{-0.05,1.05}},AspectRatio->1/1.3,ImageSize->size]

A single logic gate

input = Flatten[Table[{x1,x2},{x1,0,1},{x2,0,1}],1];
time = 10;  (* unit: hour *)

OR=Table[
gatesys={
seesawOR[12,13,{5,9},{36},1],
reporter[36,13],
concd[w[5,12],x[[1]]*c],
concd[w[9,12],x[[2]]*c],
conc[W,2*c],
conc[G,4.5*c],
conc[TH,0.6*c]
};
sol=SimulateRxnsys[gatesys,time*60*60];
Fluor[36][t*60*60]/maxSignal[srb]/c/.sol,{x,input}];

PlotSimInput[OR,input,time,"   (Subscript[x, 1])(Subscript[x, 2])","(Subscript[x, 1]) OR (Subscript[x, 2])"]

AND=Table[
gatesys={
seesawAND[12,13,{5,9},{36},1],
reporter[36,13],
concd[w[5,12],x[[1]]*c],
concd[w[9,12],x[[2]]*c],
conc[W,2*c],
conc[G,4.5*c],
conc[TH,1.2*c]
};
sol=SimulateRxnsys[gatesys,time*60*60];
Fluor[36][t*60*60]/maxSignal[srb]/c/.sol,{x,input}];

PlotSimInput[AND,input,time,"   (Subscript[x, 1])(Subscript[x, 2])","(Subscript[x, 1]) AND (Subscript[x, 2])"]

A rule 110-124 circuit

input=Flatten[Table[{L,C,R},{L,0,1},{C,0,1},{R,0,1}],2];
output={"(R110^0)","(R110^1)","(R124^0)","(R124^1)"};
time = 24;  (* unit: hour *)

R110124circuit=Table[
gatesys={
seesawAND[4,5,{27,31,35},{12,20},0],
seesawOR[6,7,{25,29,33},{14,22},0],
seesawAND[8,9,{29,33},{12},0],
seesawOR[10,11,{31,35},{14},0],
seesawOR[12,13,{5,9},{36},1],
seesawAND[14,15,{7,11},{38},1],
seesawAND[16,17,{25,29},{20},0],
seesawOR[18,19,{27,31},{22},0],
seesawOR[20,21,{5,17},{40},1],
seesawAND[22,23,{7,19},{42},1],

inputfanout[25,24,{6,16}],
inputfanout[27,26,{4,18}],
inputfanout[29,28,{6,8,16}],
inputfanout[31,30,{4,10,18}],
inputfanout[33,32,{6,8}],
inputfanout[35,34,{4,10}],

reporter[36,13],
reporter[38,15],
reporter[40,21],
reporter[42,23],

concd[w[24,25],(1-x[[1]])*c],
concd[w[26,27],x[[1]]*c],
concd[w[28,29],(1-x[[2]])*c],
concd[w[30,31],x[[2]]*c],
concd[w[32,33],(1-x[[3]])*c],
concd[w[34,35],x[[3]]*c],

conc[W,55*c],
conc[G,54*c],
conc[TH,11.2*c]
};

sol=SimulateRxnsys[gatesys,time*60*60];
{Fluor[36][t*60*60]/maxSignal[srb]/c/.sol,
Fluor[38][t*60*60]/maxSignal[srb]/c/.sol,
Fluor[40][t*60*60]/maxSignal[srb]/c/.sol,
Fluor[42][t*60*60]/maxSignal[srb]/c/.sol},
{x,input}];

PlotPerRow=2; (* number of plots per row in a grid of multiple plots *)

Grid[Partition[Table[PlotSimInput[Table[R110124circuit[[i,j]],{i,1,Length[input]}],input,time,"LCR",output[[j]],250],{j,1,Length[output]}],PlotPerRow]]

PlotPerRow=2; (* number of plots per row in a grid of multiple plots *)

Grid[Partition[Table[PlotSimOutput[Table[R110124circuit[[i,j]],{j,1,Length[output]}],output,time,"",Row[{"LCR=",Row[input[[i]]]}],250],{i,1,Length[input]}],PlotPerRow]]

A four-bit square root circuit (worse case)

with default additional threshold adjustment in logic gates at layer two and above

input=Flatten[Table[{x4,x3,x2,x1},{x4,0,1},{x3,0,1},{x2,0,1},{x1,0,1}],3];
output={"((Subscript[y, 2])^0)","((Subscript[y, 2])^1)","((Subscript[y, 1])^0)","((Subscript[y, 1])^1)"};
time = 24;  (* unit: hour *)

squareroot=Table[
gatesys={
inputfanout[13,12,{28,32,38}],
inputfanout[15,14,{30,34,36}],
inputfanout[17,16,{30,32,38}],
inputfanout[19,18,{28,34,36}],

seesawOR[20,21,{7,11},{38,40},0],
seesawAND[22,23,{5,9},{36,42},0],
seesawOR[28,29,{13,19},{40},0],
seesawAND[30,31,{15,17},{42},0],
seesawAND[32,33,{13,17},{48},0],
seesawOR[34,35,{15,19},{50},0],
seesawAND[36,37,{23,15,19},{46},1],
seesawOR[38,39,{21,13,17},{44},1],
seesawAND[40,41,{21,29},{46},1],
seesawOR[42,43,{23,31},{44},1],
seesawAND[44,45,{39,43},{52},2],
seesawOR[46,47,{37,41},{54},2],

reporter[48,33],
reporter[50,35],
reporter[52,45],
reporter[54,47],

concd[w[5,22],(1-x[[4]])*c],
concd[w[7,20],x[[4]]*c],
concd[w[9,22],(1-x[[3]])*c],
concd[w[11,20],x[[3]]*c],
concd[w[12,13],(1-x[[2]])*c],
concd[w[14,15],x[[2]]*c],
concd[w[16,17],(1-x[[1]])*c],
concd[w[18,19],x[[1]]*c],

conc[W,56*c],
conc[G,58*c],
conc[TH,12.6*c]
};
sol=SimulateRxnsys[gatesys,time*60*60];
{Fluor[48][t*60*60]/maxSignal[srb]/c/.sol,
Fluor[50][t*60*60]/maxSignal[srb]/c/.sol,
Fluor[52][t*60*60]/maxSignal[srb]/c/.sol,
Fluor[54][t*60*60]/maxSignal[srb]/c/.sol},
{x,input}];

PlotPerRow=2; (* number of plots per row in a grid of multiple plots *)

Grid[Partition[Table[PlotSimOutput[Table[squareroot[[i,j]],{j,1,Length[output]}],output,time,"",Row[{"(Subscript[x, 4])(Subscript[x, 3])(Subscript[x, 2])(Subscript[x, 1])=",Row[input[[i]]]}],250],{i,1,Length[input]}],PlotPerRow]]

A four-bit square root circuit (better case)

without additional threshold adjustment in logic gates at layer two and above

ThInc = 0;

input=Flatten[Table[{x4,x3,x2,x1},{x4,0,1},{x3,0,1},{x2,0,1},{x1,0,1}],3];
output={"((Subscript[y, 2])^0)","((Subscript[y, 2])^1)","((Subscript[y, 1])^0)","((Subscript[y, 1])^1)"};
time = 24;  (* unit: hour *)

squareroot=Table[
gatesys={
inputfanout[13,12,{28,32,38}],
inputfanout[15,14,{30,34,36}],
inputfanout[17,16,{30,32,38}],
inputfanout[19,18,{28,34,36}],

seesawOR[20,21,{7,11},{38,40},0],
seesawAND[22,23,{5,9},{36,42},0],
seesawOR[28,29,{13,19},{40},0],
seesawAND[30,31,{15,17},{42},0],
seesawAND[32,33,{13,17},{48},0],
seesawOR[34,35,{15,19},{50},0],
seesawAND[36,37,{23,15,19},{46},1],
seesawOR[38,39,{21,13,17},{44},1],
seesawAND[40,41,{21,29},{46},1],
seesawOR[42,43,{23,31},{44},1],
seesawAND[44,45,{39,43},{52},2],
seesawOR[46,47,{37,41},{54},2],

reporter[48,33],
reporter[50,35],
reporter[52,45],
reporter[54,47],

concd[w[5,22],(1-x[[4]])*c],
concd[w[7,20],x[[4]]*c],
concd[w[9,22],(1-x[[3]])*c],
concd[w[11,20],x[[3]]*c],
concd[w[12,13],(1-x[[2]])*c],
concd[w[14,15],x[[2]]*c],
concd[w[16,17],(1-x[[1]])*c],
concd[w[18,19],x[[1]]*c],

conc[W,56*c],
conc[G,58*c],
conc[TH,12.6*c]
};
sol=SimulateRxnsys[gatesys,time*60*60];
{Fluor[48][t*60*60]/maxSignal[srb]/c/.sol,
Fluor[50][t*60*60]/maxSignal[srb]/c/.sol,
Fluor[52][t*60*60]/maxSignal[srb]/c/.sol,
Fluor[54][t*60*60]/maxSignal[srb]/c/.sol},
{x,input}];

PlotPerRow=2; (* number of plots per row in a grid of multiple plots *)

Grid[Partition[Table[PlotSimOutput[Table[squareroot[[i,j]],{j,1,Length[output]}],output,10,"",Row[{"(Subscript[x, 4])(Subscript[x, 3])(Subscript[x, 2])(Subscript[x, 1])=",Row[input[[i]]]}],250],{i,1,Length[input]}],2]]