Here is the steady-state approximation for genetic repression without a dimerizing repressor.

inv[x_,a_] := a/(1+x)

Plot[{inv[inv[x,2],2],x},{x,0,2}]

Two inverters don’t make a restoring digital gate; it doesn't satisfy the digital abstraction.  The output of a cascade will be nearly independent of the input.

Plot[{inv[inv[x,10],10],x},{x,0,10}]

Playing with alpha doesn’t help.  We need to get a sigmoid.  Let’s try a dimerizing repressor.

inv2[x_,a_] := a/(1+x^2)

Plot[inv2[x,2],{x,0,2}]

Plot[{inv2[inv2[x,2],2],x},{x,0,2}]

Using alpha = 2 isn’t sufficient.  There is only one fixed point where output=input, and in a long cascade that is what the output will approach, no matter what the input is.

Plot[{inv2[inv2[inv2[inv2[x,2],2],2],2],x},{x,0,2}]

Therefore we can try a larger value of alpha, and cross our fingers.

Plot[inv2[x,4],{x,0,4}]

Plot[{inv2[inv2[x,4],4],x},{x,0,4}]

This is a decent restoring gate.  And now, it only gets better with longer cascades.

Plot[{inv2[inv2[inv2[inv2[x,4],4],4],4],x},{x,0,4}]

To build a NAND gate and keep the max output in the range 0 to 4, we must use alpha=2.

weakNAND[x_,y_,a_] := a/(1+x^2) + a/(1+y^2)

ContourPlot[weakNAND[x,y,2],{x,0,4},{y,0,4}, ContourLabels->True, Contours->10]

Digital abstraction:  ON = 3.6 to 4.0, OFF = 0 to 0.4.... fails.  Perhaps changing the range to 0 to 10, allowing alpha=5, will work?

ContourPlot[weakNAND[x,y,5],{x,0,10},{y,0,10}, ContourLabels->True, Contours->{0,.5,1.0,1.5,2.0,8.0,8.5,9.0,9.5,10.0}]

Digital abstraction: ON = 9 to 10, OFF = 0 to 1... fails.  Again.  Let’s keep the 0 to 4 range, but add a restoring gate after the NAND.  Now, we can use alpha=4 on the integrating gate itself, because the restoring gate will bring signal levels back within the standard.

weakNOT[x_,a_] := a/(1+x^2)

RESTORE[x_,a_]:= weakNOT[weakNOT[x,a],a]

ContourPlot[RESTORE[RESTORE[weakNAND[x,y,4],4],4],{x,0,4},{y,0,4}, ContourLabels->True, Contours->{0.5,1,2,3,3.5}]

Digital abstraction: ON = 3.5 to 4.0, OFF = 0 to 0.5 -- very clean, should have a nice noise boundary too.

NAND[x_,y_] := RESTORE[weakNAND[x,y,4],4]

NOT[x_] := RESTORE[weakNOT[x,4],4]

Plot[NOT[x],{x,0,4}]

Clearly satisfies the above digital abstraction.

XOR[x_,y_] := NAND[NAND[NOT[x],y],NAND[x,NOT[y]]]

ContourPlot[XOR[x,y],{x,0,4},{y,0,4}, ContourLabels->True, Contours->{0.5,1,2,3,3.5}]

Nice!

One can go ahead and make other gates that work at steady-state... but since NAND is universal, we really don’t need to.  Skip to the section on dynamics.

ANDNOT[x_,y_] := NOT[NAND[x,NOT[y]]]

ContourPlot[ANDNOT[x,y],{x,0,4},{y,0,4}, ContourLabels->True, Contours->{0.5,1,2,3,3.5}]

OR[x_,y_]  := RESTORE[RESTORE[x,4]+RESTORE[y,4],4]

ContourPlot[OR[x,y],{x,0,4},{y,0,4}, ContourLabels->True, Contours->{0.5,1,2,3,3.5}]

XOR2[x_,y_] := OR[ANDNOT[x,y],ANDNOT[y,x]]

ContourPlot[XOR2[x,y],{x,0,4},{y,0,4}, ContourLabels->True, Contours->{0.5,1,2,3,3.5}]

Also works ... XOR2 uses 26 genes rather than 18 for XOR ... why bother?

Now let' s define a framework for dynamic circuits, and simulate them.
i --> j means there is a gene repressed by i that produces j.  We use a=4, i.e. we consider every gene to have identical biophysical parameters, which simplifies the design despite perhaps being unrealistic.  So, a network is just a list of ordered pairs.  We’ll use symbolic variable names (and
Unique[ ]) to indicate repressors, so that they can be double-used as concentrations.

Clear[f,g,h,i,j,k,l,x,y,z,t, X,Y,Z]

Warning: if symbols have definitions, they can’t be used as nodes in the network. When results are bizarre, this is usually the reason.  Clearing the relevant variables is often a solution.  Watch out for Mathematica-defined symbols, such as c[ ], C[ ], and N[ ].  The later two take single arguments, and thus will definitely cause trouble.  The most common problem is to accidentally use a single equals (=, assignment) rather than double equals (==, asserting equality) when setting up initial conditions for a simulation.

To automate the construction and simulation of GRNs, we will first try to do things by hand, and then work our way to automated functions.  We start with a NAND gate simulated by hand...

s=NDSolve[{z'[t]==-z[t]+ 4/(1+x[t]^2) + 4/(1+y[t]^2),x'[t]==0, y'[t]==0,z[0]==0, x[0] ==0.2, y[0]==3.7},{x,y,z},{t,0,30}]

Plot[Evaluate[{x[t],y[t],z[t]}/.s],{t,0,30},PlotRange->All]

Here we introduce “GRN[ list of pairs of symbols ]” to specify a genetic regulatory network of repressors, where each pair {x,y} indicates that there is a gene repressed by x that otherwise produces y.   We try out NAND, and then explore how to automatically construct the differential equations.

NANDpairs = GRN[{{x,i},{y,i},{i,j},{j,z}}]

First[NANDpairs]

Cases[First[NANDpairs],{x_,i}->x]

Cases[First[NANDpairs],{x_,i}->4/(1+x^2)]

We find the equation of motion for a single repressor’s concentration, with respect to the full GRN network.  Nodes that only have outward arrows are considered to be “inputs” and thus their derivate is set to zero, rather than having the exponential decay term.

GRNeqn[node_,GRN_] := If[Length[Cases[First[GRN],{_,node}]]>0,(node'[t] == -node[t] + Plus @@ Cases[First[GRN],{x_,node}->4/(1+x[t]^2)]), (node'[t] == 0)]

GRNeqn[j,NANDpairs]

GRNeqn[x, NANDpairs]

To generate the full set of equations, we first collect the full set of repressor names, and then build the equations for each one.

GRNnodes[GRN_] := Union[Flatten[First[GRN]]]

GRNeqns[GRN_] := Map[ GRNeqn[#,GRN]&, GRNnodes[GRN] ]

GRNeqns[NANDpairs]

The next step in a simulation involves setting the initial conditions.  It will be convenient to only specify non-zero initial concentrations, so we need a way to find which nodes (repressors) were not specified by the user, so we can define their initial condition to be zero.

Cases[{x[0]==3.4,b[0]==2.2},v_[0]== _ -> v]

(#[t] == 0)& /@ {b,x,z}

Given a GRN network and a few specified initial conditions, we create the full list of initial conditions for all repressors in the system.

GRNinitial[GRN_,init_] := With[{nodes= Complement[Flatten[First[GRN]],Cases[init,v_[0]== _ -> v]]}, Join[init,(#[0]==0)& /@ nodes ]]

GRNinitial[NANDpairs,{x[0]==3.2}]

Since Mathematica’s NDSolve function takes both differential equations and initial conditions as a single list, it is convenient to put them together as a “system”.

GRNsystem[GRN_,init_] := Join[GRNeqns[GRN],GRNinitial[GRN,init]]

s=NDSolve[GRNsystem[NANDpairs,{x[0]==3.2,y[0]==0.2}],{x,y,z},{t,0,30}]

Plot[Evaluate[{x[t],y[t],z[t]}/.s],{t,0,30},PlotRange->All,PlotLegends->{x,y,z}]

Now, our final function:  given a GRN, some specified initial conditions, and the variables that we want to plot, GRNsimulate[ ] will perform the simulation and plot the results.

GRNsimulate[GRN_,init_,watch_,time_] := With[{ s=NDSolve[GRNsystem[GRN,init],watch,{t,0,time}]},
Plot[Evaluate[(#[t]& /@ watch)/.s],{t,0,time},PlotRange->All, PlotLegends->watch] ]

GRNsimulate[NANDpairs,{x[0]==3.2,y[0]==0.3},{x,y,z},5]

Here is a NOT gate: three repressors in a row.

GRNsimulate[GRN[{{i,j},{j,k},{k,l}}],{i[0]==0.2},{i,j,l},15]

Here is a (putative) oscillator: three repressors in a ring.

GRNsimulate[GRN[{{X,Y},{Y,Z},{Z,X}}],{X[0]==3.2},{X,Y,Z},15]

It doesn’t work.  The GRN is essentially one digital NOT gate with its output hooked up to its input; we know that it cannot reach a steady state in which X (or Y or Z) is in the digital range (0 to 0.5 or 3.5 to 4.0) because the NOT gate satisfies the digital abstraction: if input X is a valid 0, then the output X must be a valid 1.  The ring finds a way out: it converges on a solution in which the digital abstraction is not obeyed in either the “input” nor the “output”.

However, make the ring a bit longer, and we get sustained oscillations.

GRNsimulate[GRN[{{V,W},{W,X},{X,Y},{Y,Z},{Z,V}}],{X[0]==3.2},{V,W,X,Y,Z},100]

Can we use this oscillator to drive other circuitry?  Here, we plug oscillating Z into a NOT gate that outputs k.

GRNsimulate[GRN[{{V,W},{W,X},{X,Y},{Y,Z},{Z,V},{Z,i},{i,j},{j,k}}],{X[0]==3.2},{Z,k},100]

k should be NOT Z, but it takes a while to compute, and by that time, the oscillation has moved on. The oscillator is too fast.  If we look closer, we also notice that Z varies between 0.6 and 2.4, which is not in the range of our digital abstraction.  Thus the output of the NOT gate is not guaranteed to be digital either... and it isn’t.

GRNsimulate[GRN[{{T,U},{U,V},{V,W},{W,X},{X,Y},{Y,Z},{Z,T},{Z,i},{i,j},{j,k}}],{X[0]==3.2},{Z,k},100]

Making a ring of 7 repressors is a bit better, but not good enough.

GRNsimulate[GRN[{{R,S},{S,T},{T,U},{U,V},{V,W},{W,X},{X,Y},{Y,Z},{Z,R},{Z,i},{i,j},{j,k}}],{X[0]==3.5},{X,k},100]

This 9-repressor oscillator barely makes it to touch the digital abstraction range, but it does make it.  Note that it is exactly equivalent to a ring of three digital NOT gates, since each NOT gate is three repressors.

Our next task will be to come up with a convenient way to build large GRNs using simple modular building blocks.  It’s actually easy, so long as the GRNs don’t use the same symbol name to represent different repressors:   we just take the union of all the arrows (pairs).  Here, we make use of a Mathematica notion called “upvalues”.  The /: notation defines an upvalue rule for GRN[ ] whenever it is found inside a Plus[ ] expression.

CLOCKpairs = GRN[{{R,S},{S,T},{T,U},{U,V},{V,W},{W,X},{X,Y},{Y,Z},{Z,R}}]

GRN /: GRN[pairs1_]+GRN[pairs2_] := GRN[Union[pairs1,pairs2]]

CLOCKpairs + NANDpairs

Wow, we made a GRN consisting of an oscillator and a NAND gate that aren’t connected to each other.  Useless.

Below, we will define modules that build subcircuits using input and output names that the user provides.  The subcircuits also have internal nodes, and we ask Mathematica to provide unique symbol names for those repressors.  They don’t look pretty, but they work.  Yes, $10344 is a symbol, just like X and Y and Z.

GRNnand[x_,y_,z_] := With[{i=Unique[],j=Unique[]},GRN[{{x,i},{y,i},{i,j},{j,z}}]]

GRNnand[X,Y,Z]

GRNsimulate[GRNnand[X,Y,Z],{X[0]==3.5,Y[0]==0.3},{X,Y,Z},10]

GRNnot[x_,y_] := With[{i=Unique[],j=Unique[]},GRN[{{x,i},{i,j},{j,y}}]]

GRNand[x_,y_,z_] := With[{i=Unique[]},GRNnand[x,y,i]+GRNnot[i,z]]

GRNxor[x_,y_,z_] := With[{x1=Unique[],y1=Unique[],x2=Unique[],y2=Unique[]}, GRNnot[x,x1]+GRNnot[y,y1]+GRNnand[x1,y,x2]+GRNnand[y1,x,y2]+GRNnand[x2,y2,z]]

GRNsimulate[GRNxor[X,Y,Z],{X[0]==0.5,Y[0]==4.3},{X,Y,Z},30]

The wiggles that you see prior to t=5 are due to the initial conditions:  every repressor concentration other than X and Y starts at zero, including the automatically generated internal nodes.  It takes a while for those “arbitrary” signals to flush through the system.

Here is a circuit that takes a 6-bit input x=(X0,X1,X2,X3,X4,X5) and calculates y=(Y0,Y1,Y2,Y3,Y4,Y5)=x+1.  It increments a 6-bit binary number.  The circuit uses a ripple carry, using carry bits C1,C2,C3,C4.

GRNsimulate[GRNnot[x0,y0]+GRNxor[x0,x1,y1]+ GRNand[x0,x1,c1]+GRNxor[c1,x2,y2]+GRNand[c1,x2,c2]+GRNxor[c2,x3,y3]+GRNand[c2,x3,c3]+GRNxor[c3,x4,y4]+GRNand[c3,x4,c4]+GRNxor[c4,x5,y5],{x0[0]==3.5,x1[0]== 3.5,x2[0]==3.5,x3[0]==0.5,x4[0]==0.5,x5[0]==3.5 },{y0,y1,y2,y3,y4,y5},20]

For x = 100111, y computes to 101000, which is the correct increment.  

Can we build a binary counter, by hooking the output y up to the inputs x?  Thus, once y is computed, it feeds back to x, and gets incremented again, and this feeds back and gets incremented again, and so forth.  A working circuit should see x0 flipping between ON and OFF, x1 flipping at every other cycle, x2 flipping every fourth cycle, etc, and x5 flipping every 64th cycle. Same for the y bits, of course. Let’s see...

GRNrestore[x_,y_] := With[{i=Unique[]},GRNnot[x,i]+GRNnot[i,y]]

GRNsimulate[GRNnot[x0,y0]+GRNxor[x0,x1,y1]+ GRNand[x0,x1,c1]+GRNxor[c1,x2,y2]+GRNand[c1,x2,c2]+GRNxor[c2,x3,y3]+GRNand[c2,x3,c3]+GRNxor[c3,x4,y4]+GRNand[c3,x4,c4]+GRNxor[c4,x5,y5] + GRNrestore[y0,x0]+GRNrestore[y1,x1]+GRNrestore[y2,x2]+GRNrestore[y3,x3]+GRNrestore[y4,x4]+GRNrestore[y5,x5],{x0[0]==3.5,x1[0]== 3.5,x2[0]==3.5,x3[0]==0.5,x4[0]==0.5,x5[0]==3.5 },{y0,y1,y2,y3,y4,y5},200]

Something interesting is happening, but the circuit hit a limit cycle with just y0 and y1 flipping... and then, only in a weird way, with the same frequency (rather than decreasing frequency for higher order bits).  This is bad.  What’s happening?  When we close the feedback loop, the shorter paths in the circuit start processing the “next round of incrementing” before the longer paths through the circuit are done processing the previous round.  And everything gets garbled.  We need a smarter approach.  

The standard approach in digital circuitry is to use a clock, and have the clock gate the feedback step:  only when the clock goes high, only then does the output y get transferred back to x.  So it happens for all bits synchronously.  (A more elegant “asyncronous” approach, which uses self-timing circuits and handshaking, has been developed by Caltech’s Alain Martin and others, and is probably quite relevant for biochemical circuitry -- but the clocked approach is simpler to explain, so that’s what we do here.)

GRNclock[clock_] := With[{i=Unique[],j=Unique[]},GRNnot[clock,i]+GRNnot[i,j]+GRNnot[j,clock]]

This is exactly the same 9-repressor ring oscillator as before, but now we build it explicitly from NOT gates.

Do not use C for the clock species, as it appears to be a reserved word, like N.

GRNsimulate[GRNclock[c],{c[0]==3.5},{c},100]

We will want oscillators that have sharper ON/OFF transitions and have longer duty cycles.  Since we tied our hands behind our backs by making all genes have the same biophysical parameters, we can’t simply slow down the production and degration rates to slow down the oscillator, or play other games.  Instead, we make longer cascades.  Happily, this gives us nice sharp digital clocks.  (Again, this is not a good design, but it is easy to explain.)

GRNcycle[clock_,N_] := Module[{nodes,i},nodes=Unique[]& /@ Range[N-1];GRNnot[clock,nodes[[1]]]+Sum[GRNnot[nodes[[i]],nodes[[i+1]]],{i,1,N-2}] + GRNnot[nodes[[N-1]],clock]]

GRNclock[ _ ] is equivalent to GRNcycle[ _ , 3 ].   For even N, it’s not a clock...

GRNsimulate[GRNcycle[c,11],{c[0]==3.5},{c},200]

We will later want both the original clock signal (c) as well as a complementary clock (cc).  This is the same kind of thing that we tried earlier, hooking up a NOT gate to an oscillatory species, but now with the slower and sharper clock, it works better.  There is still a transition region where neither signal is a valid digital signal, and this is the danger zone.

GRNsimulate[GRNcycle[c,11]+GRNnot[c,notc],{c[0]==3.5},{c,notc},200]

Our second dynamical subcircuit is a bistable memory that holds one bit:  either x is on and y is off, or vice versa.  It is easy to make by just adding two NOT gates.

GRNbistable[x_,y_]:= GRNnot[x,y]+GRNnot[y,x]

GRNsimulate[GRNbistable[X,Y],{X[0]==1.8,Y[0]==1.2},{X,Y},20]

Note the wiggling signals at the beginning -- this is because we started with initial signals that don’t satisfy the digital abstraction.  (Additionally, the internal nodes in each NOT gate start at zero, and these need some time to settle into proper values.)

Enhancing the bistable memory to contain ways to set it ON or reset it OFF, we obtain the simple latch subcircuit.  When Sb and Rb are both ON, the circuit remembers the current value of Q, and its logical complement Qb.  If input Sb is pulled down to OFF, then Q is set to ON (and Qb set to OFF).  Conversely, if input Rb is pulled down to OFF, then Q is reset to OFF (and Qb is reset to ON).  Don’t pull both Rb and Sb down at the same time.

GRNlatch[Sb_,Rb_,Q_,Qb_]:= GRNnand[Sb,Qb,Q]+GRNnand[Q,Rb,Qb]

GRNsimulate[GRNlatch[Sb,Rb,Q,Qb],{Qb[0]==3.3,Q[0]==0.5,Sb[0]==0.2,Rb[0]==3.9},{Sb,Rb,Q,Qb},20]

This isn't quite what we want for iterative computations.  We would like a latch that memorizes a particular input signal when the clock is ON, but when the clock is OFF it ignores the input signal and just outputs the last value that it memorized.  So, we feed the input into Sb, clocked with a NAND gate; and we feed the negated input into Rb, clocked with another NAND gate.  Voila!

GRNclockedlatch[clock_,in_,out_, notout_]:= With[{reset=Unique[],mid1=Unique[],mid2=Unique[]},GRNnot[in,reset]+GRNnand[clock,in,mid1]+GRNnand[clock,reset,mid2]+GRNnand[mid1,notout,out]+GRNnand[mid2,out,notout]]

GRNsimulate[GRNclockedlatch[c,y,x,xb],{y[0]==0.2,x[0]==3.5,xb[0]==0.5,c[0]==4.1},{c,y,x},20]

Here we look at what happens when the clock is ON:  whatever y was, x will be too (regardless of x’s initial values).  You can play with initial conditions and see that if the clock is OFF, then x will retain its original value (after wiggling a bit because of how we set up these circuits with all-zero initial concentrations for intermediate species... which is an undesirable feature of this simple simulator).  Be careful to set xb[0] to the opposite of x[0]; with invalid initial conditions, the circuit won’t hold the bit correctly.

Let’s try this out with a two-bit binary counter.  The output bits y0 and y1 are fed into the clocked latch, and when the clock goes high, the latch loads the y0 and y1 values onto x0 and x1 respectively.

GRNsimulate[GRNcycle[c,35]+GRNnot[x0,y0]+GRNclockedlatch[c,y0,x0]+GRNxor[x0,x1,y1]+GRNclockedlatch[c,y1,x1],{c[0]==3.5 },{c,x0,x1},400]

This isn’t what we want.  When the clock goes high, rather than x0 taking on the new value from y0, it oscillates wildly.  So does x1, but a little less wildly.  The problem is that when the clock is high, the latch just passes everything through, so x0 connects to y0 and y0 connects to x0 and we get an oscillation.  We need a pair of latches.  Now, no matter what the clock cycle is, the feedback loop is broken.

GRNsimulate[GRNcycle[c,11]+GRNnot[x0,y0]+GRNclockedlatch[c,y0,z0]+GRNxor[x0,x1,y1]+GRNclockedlatch[c,y1,z1]+GRNnot[c,cc]+GRNclockedlatch[cc,z0,x0]+GRNclockedlatch[cc,z1,x1],{c[0]==3.5 },{c,x0,x1},400]

And it works!

Let’s be ambitious, and try the full 6-bit circuit.  We need to make the clock slower, to allow time for the longer paths in the full circuit.

GRNlatchpair[c_,cc_,x_,y_]:=With[{i=Unique[]},GRNclockedlatch[c,x,i]+GRNclockedlatch[cc,i,y]]

BinaryCounter6 = GRNcycle[c,15]+GRNnot[c,cc]+GRNnot[x0,y0]+GRNxor[x0,x1,y1]+ GRNand[x0,x1,c1]+GRNxor[c1,x2,y2]+GRNand[c1,x2,c2]+GRNxor[c2,x3,y3]+GRNand[c2,x3,c3]+GRNxor[c3,x4,y4]+GRNand[c3,x4,c4]+GRNxor[c4,x5,y5] + GRNlatchpair[c,cc,y0,x0]+GRNlatchpair[c,cc,y1,x1]+GRNlatchpair[c,cc,y2,x2]+GRNlatchpair[c,cc,y3,x3]+GRNlatchpair[c,cc,y4,x4]+GRNlatchpair[c,cc,y5,x5];

GRNsimulate[BinaryCounter6,{c[0]==3.5 },{x0,x1,x2,x3,x4,x5},2000]

It is hard to see what’s going on.  Let’s stack the plots, by adding meaningless offsets.

GRNsimulate2[GRN_,init_,watch_,time_] := With[{ s=NDSolve[GRNsystem[GRN,init],watch,{t,0,time}]},
Plot[Evaluate[MapThread[#1[t]+5*#2&, {watch,Range[Length[watch]]}]/.s],{t,0,time},PlotRange->All, PlotLegends->watch] ]

GRNsimulate2[BinaryCounter6,{c[0]==3.5 },{x0,x1,x2,x3,x4,x5},4000]

Hey, that looks like a binary counter!

Conclusion: The “simplest possible” genetic regulatory mechanism, repression with dimerization, is sufficient to construct arbitrary sequential (iterated) digital logic circuits.  If we were to allow more complex regulation mechanisms (varying the biophysical parameters from gene to gene, activation in addition to repression, combinatorial regulation such as two repressors both targetting the same promoter, alternative thresholding mechanisms such as competitively binding species, etc etc) then it should be possible to build much more compact circuits.  A single gene with complex regulation might be able to perform the same task as 10 genes with simple repression.  The point of our exercise was to understand the logical principles for constructing large circuits from simple components, and to drive home the remarkable potential of genetic regulatory function.  This is the “grappling hook” approach to scientific and engineering:  throw a lightweight hook to snare the conclusion, and then reel it in with more careful optimization later.

But just how bad was our construction, in terms of unnecessary complexity?

Length[First[BinaryCounter6]]

We have 397 genes.   I would wager that with optimization and only slightly more complicated regulatory mechanisms, it could be done with on the order of 30 genes.  Anyway, just for fun, here’s what Mathematica thinks the GRN looks like:

GRNgraph[GRN_]:= (#[[1]] -> #[[2]])& /@ First[GRN]

GRNgraph[NANDpairs]

GraphPlot[GRNgraph[NANDpairs],DirectedEdges->True,VertexLabeling->True]

GraphPlot[GRNgraph[BinaryCounter6],DirectedEdges->True, Method->"SpringEmbedding"]

Oh well.  So much for that.