Needs CRNsimulator.m

Definitions for building genetic regulatory networks (GRNs) with (hetero)dimerization

It can be nice to keep parameters in symbolic form until the last moment, and only then convert them to numerical values.  This can help with debugging, among other things, because looking at numeric values leaves no trace of their origin.  In Mathematica, this last step can be done using a list of rewriting rules, which will be applied after the CRN is symbolically generated and is ready to be run.

DefaultParameters = {
kd->0.001 (* /s, e.g. dilution due to 15-minute cell doubling time *),
ktx0 -> 0.001 (* /nM/s, i.e. 10^6 /M/s as the "unit" transcription rate *),
ktl->0.001 (* /nM/s, i.e. 10^6 /M/s *),
Gtot->1.0 (* nM, i.e. one molecule per 1.66 fL cell *),
kbind->1.0 (* /nM/s, i.e. 10^9 /M/s, about as fast as one can go *)
}
(*
  steady state T[i] = Gtot * ktx / kd, so ktx = 0.001 /nM/s gives mean 1 nM RNA transcript
and similarly, steady state R[i] = Gtot ktx ktl / kd / kd, so ktx = 0.001 /nM/s gives mean 1 nM protein regulator.
in sum: the unit transcription rate results in 1 nM regulator, i.e. 1 copy per cell
*)

For each component of the GRN model, we define a function that returns the relevant CRN system specifications.  G[i] will be the gene coding for RNA transcript T[i] and regulatory protein R[i].  We use DD[i,j] for the (hetero)dimer between proteins R[i] and R[j] (avoiding use of D, which has built-in meaning to Mathematica).  To provide input to the network, we invoke an upstream always-ON gene that will produce the desired steady state regulator concentration.  Expression levels for genes are given in terms of the ultimate steady state concentration of the output regulator, via the unit transcription rate.

Gene[i_,constitutiveExpression_]:={
rxn[G[i],G[i]+T[i],ktx0 constitutiveExpression],
rxn[T[i],0,kd],
rxn[T[i],T[i]+R[i],ktl],
rxn[R[i],0,kd],
conc[G[i],Gtot]
}
Dimer[i_,j_,Km_]:={
revrxn[R[i]+R[j],DD[i,j],kbind,kbind Km],
rxn[DD[i,j],0,kd]
}
Regulation[j_,i_,Km_,regulatedExpression_]:={
revrxn[R[j]+G[i],RG[j,i],kbind,kbind Km],
rxn[RG[j,i],RG[j,i]+T[i],ktx0 regulatedExpression]
}
DimerRegulation[i_,j_,k_,Km_,regulatedExpression_]:={
revrxn[DD[i,j]+G[k],DDG[i,j,k],kbind,kbind Km],
rxn[DDG[i,j,k],DDG[i,j,k]+T[k],ktx0  regulatedExpression]
}
InputRegulator[i_,concentration_]:=Gene[i,concentration]

With these building block definitions in place, we can easily build and test functional genetic modules.  We’ll first try a simple “NOT” gate, and show that it doesn’t satisfy the digital abstraction.  We choose a Km for inhibitory regulation that is in the middle of the 0-100 nM range of interest; it is easy to try other values and see that they also don’t work.

WeakNOT[j_,i_]:=Join[
Gene[i,100],
Regulation[j,i,30,0]
]

You can play with inputs ranging from 0 to 100 nM.  Again, we set those values at the end, just for style.  Plotting the results allows us to examine each species.

rsys=Join[InputRegulator[1,input],WeakNOT[1,2]]/. DefaultParameters /. input->100

tmax=5*60*60; (* 5 hours *)
sol=SimulateRxnsys[rsys,tmax];plotter=#[t]/.sol&/@SpeciesInRxnsys[rsys];
Plot[Evaluate[plotter],{t,0,tmax},PlotLegends->SpeciesInRxnsys[rsys]]

We can run a bunch of different input concentrations, and extract the output concentrations.

InOutW=Table[
rsys=Join[InputRegulator[1,in],WeakNOT[1,2]] /. DefaultParameters;
tmax=5*60*60; (* 5 hours *)
sol=SimulateRxnsys[rsys,tmax];
R[2][tmax]/.sol,
{in,0,100}
]

Plotting the results along with boxes indicating a digital abstraction range, it is easy to see two violations:  For some inputs less than epsilon (i.e. valid OFF inputs), the output is below the box (i.e. not valid ON outputs); and for all outputs greater than 100-epsilon (i.e. valid ON inputs), the output is above the box (i.e. not valid OFF outputs).  It is an exercise to the reader to show that no choice of Km and no choice of epsilon and fix this problem.

With[{eps=10,m=100},Show[
ListLinePlot[InOutW,PlotRange->{{0,101},{0,101}},AxesLabel->{In,Out}],
Graphics[Line[{{m,0},{m-#,0},{m-#,#},{m,#},{m,0}}]&@eps],
Graphics[Line[{{0,m},{#,m},{#,m-#},{0,m-#},{0,m}}]&@eps]
]]

We now investigate a “NOT” gate based on dimerization of the input regulatory protein. In the class lecture, we showed that the geometric mean of the Km for dimerization and the Km for regulatory binding should be 30% of the maximum regulatory concentration, if we want a 10% digital abstraction.   That means for a 100 nM maximum signal, KmDjj * Kmjji = 900.  But there’s more: to get a good quadratic approximation for the dimer, its concentration can’t be too high; setting dimerization binding to be very weak, and regulation to be very tight, seems to provide a solution. By simulation, we seem to get a 12% digital abstraction.

DimerNOT[j_,i_]:=Join[
Gene[i,100],
Dimer[j,j,KmDjj],
DimerRegulation[j,j,i,Kmjji,0]
] /. {KmDjj -> 10*900, Kmjji->1/10}

rsys=Join[InputRegulator[1,100],DimerNOT[1,2]]/. DefaultParameters

tmax=5*60*60; (* 5 hours *)
sol=SimulateRxnsys[rsys,tmax];plotter=#[t]/.sol&/@SpeciesInRxnsys[rsys];
Plot[Evaluate[plotter],{t,0,tmax},PlotLegends->SpeciesInRxnsys[rsys]]

InOut=Table[
rsys=Join[InputRegulator[1,in],DimerNOT[1,2]] /. DefaultParameters;
tmax=5*60*60; (* 5 hours *)
sol=SimulateRxnsys[rsys,tmax];
R[2][tmax]/.sol,
{in,0,100}
]

With[{eps=12,m=100},Show[
ListLinePlot[InOut,PlotRange->{{0,101},{0,101}},AxesLabel->{In,Out}],
Graphics[Line[{{m,0},{m-#,0},{m-#,#},{m,#},{m,0}}]&@eps],
Graphics[Line[{{0,m},{#,m},{#,m-#},{0,m-#},{0,m}}]&@eps]
]]

OK, that is a working NOT gate, barely. We can make it stronger by cascading NOT gates.  Three genes in series allows us to obtain a 17% digital abstraction.

InOut3=Table[
rsys=Join[
InputRegulator[1,in],DimerNOT[1,2],DimerNOT[2,3],DimerNOT[3,4]
] /. DefaultParameters;
tmax=10*60*60; (* 10 hours *)
sol=SimulateRxnsys[rsys,tmax];
R[4][tmax]/.sol,
{in,0,100}
]

With[{eps=17,m=100},Show[
ListLinePlot[InOut3,PlotRange->{{0,101},{0,101}},AxesLabel->{In,Out}],
Graphics[Line[{{m,0},{m-#,0},{m-#,#},{m,#},{m,0}}]&@eps],
Graphics[Line[{{0,m},{#,m},{#,m-#},{0,m-#},{0,m}}]&@eps]
]]

Five genes in series gives us a NOT function satisfying a 25% digital abstraction.

InOut5=Table[
rsys=Join[
InputRegulator[1,in],DimerNOT[1,2],DimerNOT[2,3],DimerNOT[3,4],
DimerNOT[4,5],DimerNOT[5,6]
] /. DefaultParameters;
tmax=10*60*60; (* 10 hours *)
sol=SimulateRxnsys[rsys,tmax];
R[6][tmax]/.sol,
{in,0,100}
]

With[{eps=25,m=100},Show[
ListLinePlot[InOut5,PlotRange->{{0,101},{0,101}},AxesLabel->{In,Out}],
Graphics[Line[{{m,0},{m-#,0},{m-#,#},{m,#},{m,0}}]&@eps],
Graphics[Line[{{0,m},{#,m},{#,m-#},{0,m-#},{0,m}}]&@eps]
]]

It turns out that even the pathetic single-gene version of the dimerization nonlinearity is enough to make a valid digital NOR gate.  We just have two dimerizing inputs, either of which is sufficient to inhibit transcription.  With a 12% digital abstraction, we see that the output is ON (white) if both inputs are OFF, while the output is OFF (blue) in the other three valid input cases.

DimerNOR[i_,j_,k_]:=Join[
Gene[k,100],
Dimer[i,i,10*900],
DimerRegulation[i,i,k,1/10,0],
Dimer[j,j,10*900],
DimerRegulation[j,j,k,1/10,0]
]

InOutNOR=Flatten[Table[
rsys=Join[
InputRegulator[1,xinit],
InputRegulator[2,yinit],
DimerNOR[1,2,3]
] /. DefaultParameters;
tmax=5*60*60;
sol=SimulateRxnsys[rsys,tmax];
{xinit,yinit,R[3][tmax]/.sol},
{xinit,0,100,5},{yinit,0,100,5}
],1];

With[{eps=12,m=100},Show[
ListContourPlot[InOutNOR,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]
]]

Let us build a circuit that computes the 2-bit square root of a 4-bit number, i.e. if x = x4 x3 x2 x1 (in binary), then y = y2 y1 (in binary) will be the (rounded-down) square root of x.   In table form:
0 0 0 0 --> 0 0
0 0 0 1 --> 0 1
0 0 1 0 --> 0 1
0 0 1 1 --> 0 1
0 1 0 0 --> 1 0 
0 1 0 1 --> 1 0
0 1 1 0 --> 1 0
0 1 1 1 --> 1 0
1 0 0 0 --> 1 0
1 0 0 1 --> 1 1
1 0 1 0 --> 1 1
1 0 1 1 --> 1 1
1 1 0 0 --> 1 1
1 1 0 1 --> 1 1 
1 1 1 0 --> 1 1
1 1 1 1 --> 1 1
It should be obvious how to draw the digital circuit -- made exclusively of NOR gates and NOT gates -- from the specification below.

DimerSQRT=Join[
InputRegulator[1,x1in],
InputRegulator[2,x2in],
InputRegulator[3,x3in],
InputRegulator[4,x4in],
DimerNOR[1,2,5],
DimerNOT[3,6],DimerNOR[4,6,7],
DimerNOR[3,4,8],DimerNOT[8,18], (* output y2 *)
DimerNOT[3,9],DimerNOT[4,10],DimerNOR[9,10,11],
DimerNOR[5,7,12],
DimerNOT[5,13],DimerNOT[11,14],DimerNOR[13,14,15],
DimerNOR[12,15,16],DimerNOT[16,17] (* output y1 *)
]

InOutSQRT=With[{eps=12,m=100},Flatten[Table[
rsys=DimerSQRT /. DefaultParameters;
tmax=15*60*60;
sol=SimulateRxnsys[rsys,tmax];
y1out=R[17][tmax]/.sol;
y2out=R[18][tmax]/.sol;
Print[x4in," ",x3in, " ", x2in, " ", x1in, " -> ", y2out, " ", y1out]; 
{x4in,x3in,x2in,x1in,y2out,y1out},
{x4in,{eps,m-eps}},{x3in,{eps,m-eps}},
{x2in,{eps,m-eps}},{x1in,{eps,m-eps}}
],3]];

Well, that looks almost right, but many of the outputs aren’t in the 12% digital abstraction range that we thought we guaranteed!  What went wrong?

There are two things to consider.  They are both related to fan-out: when a single wire (i.e. regulator signal) drives more than one downstream gate.  In our earlier CRN construction for digital circuits, we made sure that signals acted catalytically on downstream gates, and thus they can serve as input to as many downstream gates as we like.  Thus, first, when a regulatory dimer binds to gate A and to gate B, its concentration is reduced more than if it were binding only to gate A.  In the extreme of having thousands of distinct downstream gates driven by the same signal, we can see that would correspond to a thousand-fold larger Gtot for a single gate, whereupon it is obvious that the gate won’t be fully suppressed (and our assumption that the gene concentrations are much lower than the signal concentration will fail).  A second consideration relates not to the GRN so much as to our modules for building them:  if two genes are regulated by the same signal i, then the reactions for Dimer[i,j,Km] will appear twice in the resulting CRN specification.  Since the CRNSimulator package adds rate constants when the same reaction appears multiple times, this results in faster equilibration of the dimerization.  Perhaps not what we intended, but unlikely to cause trouble.  So we’ll focus on the first issue.

If we have limited fan-out, e.g. each signal driving no more than three downstream gates as in this case, then the simplest solution is to accept that a single-gene NOT gate or NOR gate may no longer output a valid signal concentration, because its input is somewhat depleted -- but it will be close.  Since we already saw that a series of NOT gates provides greater restoration and enables a stronger digital abstraction, one thing we can do -- overkill, actually -- is to replace each original gate by [that same gate followed by two NOT gates].  

Below, we do that by automatically numbering the new gates using an offset that is larger than the number of gates in the original circuit -- this is just a hack to give each new gate a unique integer number.

offsethack=100;
DimerNOT3[i_,j_]:=Join[DimerNOT[i,j+offsethack],DimerNOT[j+offsethack,j+2offsethack],DimerNOT[j+2offsethack,j]]
DimerNOR3[i_,j_,k_]:=
Join[DimerNOR[i,j,k+offsethack],DimerNOT[k+offsethack,k+2offsethack],DimerNOT[k+2offsethack,k]]

DimerSQRT3=Join[
InputRegulator[1,x1in],
InputRegulator[2,x2in],
InputRegulator[3,x3in],
InputRegulator[4,x4in],
DimerNOR3[1,2,5],
DimerNOT3[3,6],DimerNOR3[4,6,7],
DimerNOR3[3,4,8],DimerNOT3[8,18], (* output y2 *)
DimerNOT3[3,9],DimerNOT3[4,10],DimerNOR3[9,10,11],
DimerNOR3[5,7,12],
DimerNOT3[5,13],DimerNOT3[11,14],DimerNOR3[13,14,15],
DimerNOR3[12,15,16],DimerNOT3[16,17] (* output y1 *)
];

InOutSQRT=With[{eps=12,m=100},Flatten[Table[
rsys=DimerSQRT3 /. DefaultParameters;
tmax=15*60*60;
sol=SimulateRxnsys[rsys,tmax];
y1out=R[17][tmax]/.sol;
y2out=R[18][tmax]/.sol;
Print[x4in," ",x3in, " ", x2in, " ", x1in, " -> ", y2out, " ", y1out]; 
{x4in,x3in,x2in,x1in,y2out,y1out},
{x4in,{eps,m-eps}},{x3in,{eps,m-eps}},
{x2in,{eps,m-eps}},{x1in,{eps,m-eps}}
],3]];

It works quite nicely.  However, it took 46 genes to robustly implement a digital circuit with 7 NOR gates (and a bunch of NOT gates).

Length[Union[Cases[DimerSQRT3,rxn[G[x_],_,_]:>x]]]

Finally, we show that a three-gene NOT gate based on titration provides the strongest threshold behavior yet -- a 45% digital abstraction!  The input regulates (by inhibiting) the production of an “inhibitor”; the inhibitor in turn is titrated by a competitor species; if any competitor remains, then it tightly turns off the output.  Double negatives all the way.  Look at the extremes:  no input, so inhibitor tries to go to 100 but is drained to 50 by the competitor; thus the competitor is shut off completely, and the output is free to go to 100.  Or, if the input is 100, then that exceeds the Km for regulating the inhibitor and we get only about 33 of the inhibitor; since this is swamped by the competitor, about 17 of competitor remains, which is plenty to shut off the output nearly completely. 

In between the extremes, the transition is sharp.  This can be argued intuitively as follows:  So long as the inhibitor is being produced faster than the competitor, the competitor will be pretty much wiped out by titration.  That will happen when the input is less than 50.  But, for larger input levels, the competitor will survive; it only has to be in excess by 9 in order to inhibit the output 10-fold.

TitrationNOT[input_,output_,competitor_,inhibitor_]:=Join[
Gene[inhibitor,100],Regulation[input,inhibitor,50,0],
Gene[competitor,50], Dimer[competitor,inhibitor,0.01],
Gene[output,100],
Regulation[competitor,output,1,0]
]

rsys=TitrationNOT[1,4,2,3]/.DefaultParameters

tmax=5*60*60; (* 5 hours *)
sol=SimulateRxnsys[rsys,tmax];plotter=#[t]/.sol&/@SpeciesInRxnsys[rsys];
Plot[Evaluate[plotter],{t,0,tmax},PlotLegends->SpeciesInRxnsys[rsys]]

InOutT=Table[
rsys=Join[InputRegulator[1,in],TitrationNOT[1,4,2,3]] /. DefaultParameters;
tmax=10*60*60; (* 10 hours *)
sol=SimulateRxnsys[rsys,tmax];
R[4][tmax]/.sol,
{in,0,100}
]

With[{eps=45,m=100},Show[
ListLinePlot[InOutT,PlotRange->{{0,101},{0,101}},AxesLabel->{In,Out}],
Graphics[Line[{{m,0},{m-#,0},{m-#,#},{m,#},{m,0}}]&@eps],
Graphics[Line[{{0,m},{#,m},{#,m-#},{0,m-#},{0,m}}]&@eps]
]]

It is easy to convert this mechanism into a NOR gate.  We get a 25% digital abstraction immediately.

TitrationNOR[input1_,input2_,output_,competitor_,inhibitor_]:=Join[
Gene[inhibitor,100],
Regulation[input1,inhibitor,50,0],
Regulation[input2,inhibitor,50,0],
Gene[competitor,50], Dimer[competitor,inhibitor,0.01],
Gene[output,100],
Regulation[competitor,output,1,0]
]

InOutNORT=Flatten[Table[
rsys=Join[
InputRegulator[1,xinit],
InputRegulator[2,yinit],
TitrationNOR[1,2,3,4,5]
] /. DefaultParameters;
tmax=5*60*60;
sol=SimulateRxnsys[rsys,tmax];
{xinit,yinit,R[3][tmax]/.sol},
{xinit,0,100,5},{yinit,0,100,5}
],1];

With[{eps=25,m=100},Show[
ListContourPlot[InOutNORT,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]
]]

To build a circuit, it’s convenient to automatically number the “internal” competitor and inhibitor species for each module.  We use Mathematica’s ability to have different rules for the same function with different arguments:

TitrationNOT[i_,k_]:=Module[{inh,comp},TitrationNOT[i,k,comp,inh]]
TitrationNOR[i_,j_,k_]:=Module[{inh,comp},TitrationNOR[i,j,k,comp,inh]]

TitrationSQRT=Join[
InputRegulator[1,x1in],
InputRegulator[2,x2in],
InputRegulator[3,x3in],
InputRegulator[4,x4in],
TitrationNOR[1,2,5],
TitrationNOT[3,6],TitrationNOR[4,6,7],
TitrationNOR[3,4,8],TitrationNOT[8,18], (* output y2 *)
TitrationNOT[3,9],TitrationNOT[4,10],TitrationNOR[9,10,11],
TitrationNOR[5,7,12],
TitrationNOT[5,13],TitrationNOT[11,14],TitrationNOR[13,14,15],
TitrationNOR[12,15,16],TitrationNOT[16,17] (* output y1 *)
]

InOutSQRTT=With[{eps=12,m=100},Flatten[Table[
rsys=TitrationSQRT /. DefaultParameters;
tmax=15*60*60;
sol=SimulateRxnsys[rsys,tmax];
y1out=R[17][tmax]/.sol;
y2out=R[18][tmax]/.sol;
Print[x4in," ",x3in, " ", x2in, " ", x1in, " -> ", y2out, " ", y1out]; 
{x4in,x3in,x2in,x1in,y2out,y1out},
{x4in,{eps,m-eps}},{x3in,{eps,m-eps}},
{x2in,{eps,m-eps}},{x1in,{eps,m-eps}}
],3]];

Gotta be happy with that!  However, it is still the same number of genes as the previous circuit; in both cases (dimerization with extra restoration, and the titration module above) each NOT takes 3 genes and each NOR takes 3 genes.

Length[Union[Cases[TitrationSQRT,rxn[G[x_],_,_]:>x]]]

We can also look at the somewhat smaller circuit found by LogicFriday / Espresso.  It comes to just 4 input genes + 27 circuit genes = 31 genes.

TitrationSQRT2=Join[
InputRegulator[1,x1in],
InputRegulator[2,x2in],
InputRegulator[3,x3in],
InputRegulator[4,x4in],
TitrationNOR[1,2,5],TitrationNOR[3,4,6],
TitrationNOT[3,7], TitrationNOT[4,8], TitrationNOT[6,9], (* output y2 *)
TitrationNOR[3,5,10],TitrationNOR[7,8,11],
TitrationNOR[10,11,12],TitrationNOT[12,13] (* output y1 *)
]

InOutSQRT2T=With[{eps=12,m=100},Flatten[Table[
rsys=TitrationSQRT2 /. DefaultParameters;
tmax=15*60*60;
sol=SimulateRxnsys[rsys,tmax];
y1out=R[13][tmax]/.sol;
y2out=R[9][tmax]/.sol;
Print[x4in," ",x3in, " ", x2in, " ", x1in, " -> ", y2out, " ", y1out]; 
{x4in,x3in,x2in,x1in,y2out,y1out},
{x4in,{eps,m-eps}},{x3in,{eps,m-eps}},
{x2in,{eps,m-eps}},{x1in,{eps,m-eps}}
],3]];

Length[Union[Cases[TitrationSQRT2,rxn[G[x_],_,_]:>x]]]

The next challenge would be to devise a four-gene titration module that produces both a NOT signal as well as a IDENTITY signal, which would obviate the need for all the NOT gates in the circuit, at a modest cost for the NORs. 

Better yet is the linear threshold unit (LTU) construction, using 3 genes. (If we want both a positive and a negative output, we can use 4 genes.  But in point of fact, this gives no extra power in an LTU circuit.) The LTU is defined by a vector of (positive or negative) weights and (positive or negative) threshold.  We’ll first examine the model by evaluating a three-LTU circuit that computes the same 4-bit square-root function that was implemented with NOR and NOT gates above (nine of them in the minimized circuit).

LTU[inputs_,weights_,threshold_]:= If[inputs.weights-threshold >0, 1, 0]

InOutSQRTbyLTU=Module[{x1,x2,x3,x4,z1,y1,y2},Flatten[Table[
z1 = LTU[{x1,x2,x3},{1,1,-2},0.5];
y1=LTU[{z1,x3,x4},{2,1,1},1.5];
y2=LTU[{x3,x4},{1,1},0.5];
Print[x4," ",x3, " ", x2, " ", x1, " -> ", y2, " ", y1]; 
{x4,x3,x2,x1,y2,y1},
{x4,0,1},{x3,0,1},{x2,0,1},{x1,0,1}
],3]];

OK, that does the job.  Now we’ll make our three-gene LTU using just inhibition and titration.

TitrationLTU[inputgates_,weights_,threshold_, genei_,geneip_,outputgate_]:=Module[{ci,cip,Km,Kmp},
(* 
inputgates is a list of indices for the inputs;
weights is a list of positive or negative real numbers, for the corresponding inputs;
threshold is a positive of negative real number;
genei and geneip are the indices for the two intermediate genes;
outputgate is the index of the output gene.
*)
(* fill in the code here! *)
]

(* 
For CRNSimulator, the indices don't have to be integers; they can also be symbols.  So the following definition allows one to not bother figuring out indices for the intermediate genes,
just producing a random symbol ("local variable") upon each invocation of the Module[ ] statement.
*)
TitrationLTU[inputgates_,weights_,threshold_,outputgate_]:=Module[{genei, geneip},
TitrationLTU[inputgates,weights,threshold, genei,geneip,outputgate]
]

It makes sense to first test a 1-input LTU, and look at the sharpness of the decision boundary for the equivalent of a NOT gate and a IDENTITY gate.  In fact, the parameters embedded in the above construction were tuned to make the threshold reasonably sharp and symmetric.  Note two things about these parameters: First, the weights and threshold are rescaled to our concentration range such that the threshold (cip - ci) is always either 50 nM or -50 nM.   Multiplying both weights and threshold simultaneously does not alter the decision boundary (although it could affect the sharpness of the boundary).  Second, the LTU is specified with input and output being 0 or 1, while the GRN has input and output in the range 0 nM to 100 nM, so for equivalent function, the weights must be divided by 100.  Since w_ji

TitrationLTU[{1},{1},0.5,2,3,4]

{G[2] (Overscript[→,  (50 ktx0) ]) G[2]+T[2] ,T[2] (Overscript[→,  kd ]) 0 ,T[2] (Overscript[→,  ktl ]) T[2]+R[2] ,R[2] (Overscript[→,  kd ]) 0 ,conc[G[2],Gtot],R[1]+G[2] (Overscript[→,  kbind ]) RG[1,2] ,RG[1,2] (Overscript[→,  (10000 kbind) ]) R[1]+G[2] ,RG[1,2] (Overscript[→,  0 ]) RG[1,2]+T[2] ,G[3] (Overscript[→,  (100 ktx0) ]) G[3]+T[3] ,T[3] (Overscript[→,  kd ]) 0 ,T[3] (Overscript[→,  ktl ]) T[3]+R[3] ,R[3] (Overscript[→,  kd ]) 0 ,conc[G[3],Gtot],R[1]+G[3] (Overscript[→,  kbind ]) RG[1,3] ,RG[1,3] (Overscript[→,  (49.504950495049506` kbind) ]) R[1]+G[3] ,RG[1,3] (Overscript[→,  0 ]) RG[1,3]+T[3] ,R[2]+R[3] (Overscript[→,  kbind ]) DD[2,3] ,DD[2,3] (Overscript[→,  0 ]) R[2]+R[3] ,DD[2,3] (Overscript[→,  kd ]) 0 ,G[4] (Overscript[→,  (100 ktx0) ]) G[4]+T[4] ,T[4] (Overscript[→,  kd ]) 0 ,T[4] (Overscript[→,  ktl ]) T[4]+R[4] ,R[4] (Overscript[→,  kd ]) 0 ,conc[G[4],Gtot],R[3]+G[4] (Overscript[→,  kbind ]) RG[3,4] ,RG[3,4] (Overscript[→,  (0.1` kbind) ]) R[3]+G[4] ,RG[3,4] (Overscript[→,  0 ]) RG[3,4]+T[4] }

InOutT=Table[
rsys=Join[InputRegulator[1,in],TitrationLTU[{1},{1},0.5,2,3,4]] /. DefaultParameters;
tmax=10*60*60; (* 10 hours *)
sol=SimulateRxnsys[rsys,tmax];
R[4][tmax]/.sol,
{in,0,100}
]

With[{eps=45,m=100},Show[
ListLinePlot[InOutT,PlotRange->{{0,101},{0,101}},AxesLabel->{In,Out}],
Graphics[Line[{{m,0},{m-#,0},{m-#,#},{m,#},{m,0}}]&@eps],
Graphics[Line[{{0,m},{#,m},{#,m-#},{0,m-#},{0,m}}]&@eps]
]]

TitrationLTU[{1},{-1},-0.5,2,3,4]

InOutT=Table[
rsys=Join[InputRegulator[1,in],TitrationLTU[{1},{-1},-0.5,2,3,4]] /. DefaultParameters;
tmax=10*60*60; (* 10 hours *)
sol=SimulateRxnsys[rsys,tmax];
R[4][tmax]/.sol,
{in,0,100}
]

With[{eps=45,m=100},Show[
ListLinePlot[InOutT,PlotRange->{{0,101},{0,101}},AxesLabel->{In,Out}],
Graphics[Line[{{m,0},{m-#,0},{m-#,#},{m,#},{m,0}}]&@eps],
Graphics[Line[{{0,m},{#,m},{#,m-#},{0,m-#},{0,m}}]&@eps]
]]

When we then test a 2-input LTU, note that the safe digital abstraction range must shrink, because of the way that the linear decision boundary cuts through the input plane.  Note also that a slightly shifted boundary would allow us to expand that range without affecting the boolean function being computed. And that with different weights, and hence decision boundaries at different angles, we can expect different digital abstraction ranges to work.  Adjusting the LTU so that the decision boundary leaves as wide a margin as possible will help.  Moving into higher dimensions, the margin is likely to shrink (e.g. for a n-input OR gate LTU with the boundary (∑_j)(Subscript[w, j,i])(Subscript[x, i])>0.5, the margin drops to (1/(2(Sqrt[n]))) prior to rescaling to concentrations).

TitrationLTU[{1,2},{1,-1},0.5,3,4,5]

InOutLTU=Flatten[Table[
rsys=Join[
InputRegulator[1,xinit],
InputRegulator[2,yinit],
TitrationLTU[{1,2},{1,-1},0.5,3,4,5]
] /. DefaultParameters;
tmax=15*60*60;
sol=SimulateRxnsys[rsys,tmax];
{xinit,yinit,R[5][tmax]/.sol},
{xinit,0,100,5},{yinit,0,100,5}
],1]

ListPlot3D[InOutLTU,PlotRange->All,AxesLabel->{x,y}]

With[{eps=20,m=100},Show[ (* note that the jaggeness of the boundary is an artifact of our discrete sampling *)
ListContourPlot[InOutLTU,Contours->{eps,m-eps},PlotLegends->Automatic,PlotRange->All],
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]
]]

Finally we can implement the SQRT function with just three LTU units and 9 computing genes.  The results are satisfactory.

TitrationSQRTLTU=Join[
InputRegulator[1,x1in],
InputRegulator[2,x2in],
InputRegulator[3,x3in],
InputRegulator[4,x4in],
TitrationLTU[{1,2,3},{1,1,-2},0.5,5],
TitrationLTU[{5,3,4},{2,1,1},1.5,6], (* output y1 *)
TitrationLTU[{3,4},{1,1},0.5,7] (* output y2 *)
];

InOutSQRTLTU=With[{eps=12,m=100},Flatten[Table[
rsys=TitrationSQRTLTU /. DefaultParameters;
tmax=15*60*60;
sol=SimulateRxnsys[rsys,tmax];
y1out=R[6][tmax]/.sol;
y2out=R[7][tmax]/.sol;
Print[x4in," ",x3in, " ", x2in, " ", x1in, " -> ", y2out, " ", y1out]; 
{x4in,x3in,x2in,x1in,y2out,y1out},
{x4in,{eps,m-eps}},{x3in,{eps,m-eps}},
{x2in,{eps,m-eps}},{x1in,{eps,m-eps}}
],3]];