Modeling, Simulation and Visualization of Gene Networks

Newton's Method can be viewed as a Discrete Dynamical System in that at prescribed times (corresponding to iterations) we move from point to point in the complex plane under the action of a specific rule, z=F(z) where F(z)=z-f(z)/f'(z) is a rational function built from the complex polynomial f. In fact, what we have done is to consider 

  z   F(z)   F(F(z))   F(F(F(z)))  ...
either until a new iterate produces no change (because f=0) or until we reach a preordained maximum number of iterations. We shall, in this chapter, widen our scope beyond the complex rational F of the last chapter. Now z will be the state of a gene network and F will implement the inexorable logic of gene regulation.

Our first goal will be to model and study networks of genes along the lines of the work of Wuensche . For a glimpse at the complexity of gene regulation, please open (in a new window) the Encyclopedia of Escherichia coli K-12 Genes and Metabolism. Once there, 

   click on Genome Browser, 
   enter fuma in the Gene name box and click go, 
   mouse over fumA (blue hatched at center) and click,
   scroll down to Genetic Regulation Schematic and mouse about.

Please visit iGEM to learn about genetically engineered machines at Rice.

We shall assume that a gene is either on (1) or off (0) and that a gene network is a merely a Boolean Network, i.e., a List of nodes, with a wiring list and a logic rule for each node. It will be convenient to assume that each gene is regulated by exactly three of its neighbors. As such there are then 8 possible inputs at each gene and so 256 possible logic rules at each gene. For example, the 6 node net associated with Fig 12 in Wuensche is

  






Node   Wiring	Rule

 1     4 2 3	231
 2     5 3 2	64
 3     3 6 1	5
 4     5 4 6	108
 5     6 1 2	61
 6     3 5 6	62

The rule numbers correspond to the rule table 

   000	001	010	011	100	101	110	111

    1	 1	 1 	 0	 0	 1	 1	 1	231
    0	 0	 0	 0	 0	 0	 1	 0	64
    1	 0	 1	 0	 0	 0	 0	 0	5
    0	 0	 1	 1	 0	 1	 1	 0      108
    1	 0	 1	 1	 1	 1	 0	 0	61
    0	 1	 1	 1	 1	 1	 0	 0	62

Several of these logic functions are fairly complex while several are pretty simple. For example, gene 2, inhibits itself and requires activation from both genes 5 and 3. In other words 

    gene 2 = (gene 5 AND gene 3) AND NOT gene 2
Similarly, the logic at gene 6 follows 
    gene 6 = (gene 3 XOR gene 5) OR (gene 6 AND NOT (gene 3 OR gene 5))

When we run/execute/animate this net we find (box = on, blank = off)

Do you see the three attractors found by Wuensche? Do they have the correct periods? How was our boxy picture produced?

  1. Convert rule to ruletable by running down powers of two.
  2. Start from a random state, say, e.g., s = [1 0 1 0 1 0].
  3. Run through each gene in s to make the next state, say, ns.
  4. In order to find ns(1) we first go to wire and see that wire(1,:)=[4 2 3] and so we look up the state of genes 4, 2 and 3 in s and find the pattern 001.
  5. This now tells us what column in row 1 of the rule table to use. Namely, 001 corresponds to column 2. As ruletab(1,2)=1 we set ns(1)=1.
  6. Take care of the rest of the genes in precisely the same way.
  7. Now plot the new state by putting a box for each on gene. Here the x-axis is used for gene state and y-axis for time. If j is the "time" counter and ong is the list of genes that are on then something like
    plot(ong,-j*ones(size(ong)),'s')
    should do the trick. How do find which genes are on? Use find.

Regarding the more general class of cellular automata I recommend you take a look at the text of Wolfram and indulge in the game of life.

With your code you can now run forward from any state and so, eventually, land in an attractor. Once there, we'd like to work backwards and find every state that feels the pull.

We shall pursue a representation like that in Wuensche, our as on the cover of Origins of Order. For the coming week we will put last week's expertise into use in the creation of figures like

This is one attractor of the net whose wire, rule and seed are as stated. Let us break down its production steps

  1. Follow the seed into its attractor.
  2. Plot the attractor.
  3. Pick an attractor state and find all of its predecessors.
  4. Plot its fan of predecessors.
  5. Find and plot the predecessors of the predecessors until you reach the Garden of Eden.
  6. Repeat for all other attractor states.
You can see these steps in action when you examine the diary associated with the attractor above. You see that the seed and rule are random but that wire was generated by gring.

The most difficult step in this procedure is to find all predecessors of a chosen state. The brute force means, of advancing every state and then collecting the ones that hit the chosen, is neither intellectually satisfying nor practical on anything but the smallest of nets. It is, however, easy to code, and so we should all code it first before embarking on the extra credit bit.

Brute Force

For n genes
  1. make an (2^n)-by-n matrix, s1, of all possible states,
  2. use last week's code to make a (2^n)-by-n matrix, s2, of successors of s1,
  3. given a state s, use ismember to find each of its occurrences in s2. The predecessors of s are then each of the associated rows of s1.
Note that s1 and s2 can be made once and for all upon entry to boa.
For help on plotting please follow the hints on the assignment page.

Extra Credit

We will now work backwards through wire and rule to uncover the chosen one's ancestors. We fix our ideas let us return to last week's familiar net 
   000     001     010     011     100     101     110     111

    1       1       1       0       0       1       1       1      231
    0       0       0       0       0       0       1       0      64
    1       0       1       0       0       0       0       0      5
    0       0       1       1       0       1       1       0      108
    1       0       1       1       1       1       0       0      61
    0       1       1       1       1       1       0       0      62
and work backwards from the attractor 
   s = [1 0 0 0 0 1]
We see that s(1)=1 is on and so the rule table informs us that genes 4, 2 and 3 must have been either 
   [0 0 0], [0 0 1], [0 1 0], [1 0 1], [1 1 0] or [1 1 1]     (***)
Starting with the first we nominate the candidate 
  [NaN 0 0 0 NaN NaN]
and sequentially run it past the other genes. In particular, as s(2)=0 we know that genes 5, 3, and 2 must have been either 
   [0 0 0], [0 0 1], [0 1 0], [0 1 1], [1 0 0], [1 0 1], or [1 1 1]
We see that all but the first and fifth contradict our candidate, and so we are left with two candidates 
  [NaN 0 0 0 0 NaN]  and [NaN 0 0 0 1 NaN]
Proceeding, as s(3)=0, we know that genes 3, 6, and 1 must have been either 
   [0 0 1], [0 1 1], [1 0 0], [1 0 1], [1 1 0], or [1 1 1]
The last 4 contradict our two candidates and so our candidate list grows to 
  [1 0 0 0 0 0], [1 0 0 0 0 1], [1 0 0 0 1 0], [1 0 0 0 1 1].
Proceeding, as s(4)=0 we know that genes 5, 4, and 6 must have been either 
   [0 0 0], [0 0 1], [1 0 0] or [1 1 1]
The last contradicts each candidate while the second and third inform that genes 5 and 6 could not have both been on. As such our candidate list shrinks to 
  [1 0 0 0 0 0], [1 0 0 0 0 1], [1 0 0 0 1 0]
Proceeding, as s(5)=0 we know that genes 6, 1, and 2 must have been either 
   [0 0 1], [1 1 0] or [1 1 1]
The first and third contradict our candidates while the second informs us that gene 6 must have been on. Our candidate list is therefore 
   [1 0 0 0 0 1]
Finally, as s(6)=1 we know that genes 3, 5, and 6 must have been either 
   [0 0 1], [0 1 0], [0 1 1], [1 0 0] or [1 0 1]
As the first concurs with our candidate it is the sole predecessor - stemming from the first choice in (***). This is a bit underwhelming - for we have landed where we started. I imagine however, that if we start from the third (or 5th or 6th) choice in (***) that you will uncover 
[1 1 0 0 0 1].
After checking each of these steps by hand I recommend that you follow a second state (by hand) to its predecessors. Following that exercise you will yearn to teach your logic to Matlab.