Faure2006dag

<bibtex>@article{faure2006dag,
  title={{Dynamical analysis of a generic Boolean model for the control of the mammalian cell cycle}},
  author={Faure, A. and Naldi, A. and Chaouiya, C. and Thieffry, D.},
  journal={Bioinformatics},
  volume={22},
  number={14},
  year={2006},
  publisher={Oxford Univ Press}
}
</bibtex>
[http://scholar.google.com/scholar?q=dynamical+analysis+faure+naldi+boolean&hl=en&lr=&output=search link]



These folks too a biological regulatory network for the mammalian cell cycle and converted it to a boolean network.  Much of the paper is comparing the behavior of the Boolean model to that of the actual network.  The reduction of the literature to binary rules is great.  I coded it into Python with little trouble.

Link to [http://discover.nci.nih.gov/kohnk/fig6a.html Kohn's map]

Citation of Chaves et al (2005).  "One option is to consider the possibility that the realiszation of some transitions requires several updating steps.  [They] have explored this option to improve that analysis..."  Look to them for alternate updating schemes.

Their results with asynchronous updating were interesting.  They tiered molecules and than uupdated asynchronously within that and superimposed all the resulting transition graphs into one, that was more bioologically realistic than the synchronous scheme, which ended up as a subset of the asynchronous juxtaposition network.

Here is my code reproducing their results and drawing the transition graph for synch updating"
<pre>
from igraph import *
import random

g = Graph(directed=True)        # directed edges

g.add_vertices(5)
g.add_edges( [ (0, 0) ] );      # graphObject.add_edges( [ (start_vertex, end_vertex) ] )
g.add_edges( [ (0, 1) ] );
g.add_edges( [ (1, 2) ] );
g.add_edges( [ (2, 3) ] );
g.add_edges( [ (3, 4) ] );
g.add_edges( [ (4, 3) ] );
g.add_edges( [ (3, 5) ] );
g.add_edges( [ (4, 4) ] );
g.add_edges( [ (5, 2) ] );

# g.layout_graphopt() is just one layout algorithm,
# others are available @ http://cneurocvs.rmki.kfki.hu/igraph/doc/python/igraph.GraphBase-class.html
#plot(g, layout = g.layout_graphopt() )

class RBN:
    "Random boolean network implementation"

    
    
    def __init__(self):
        self.nodes=[]
        self.__previousSystemStates=[]

    def update(self):
        "update state of network"
        
        systemState=[]
        for i, n in enumerate(self.nodes):
            neighborStates = tuple([inNode.nodeState for inNode in n.inputs])
            systemState.insert(i, n.functions.get(neighborStates, 0) )
            #print [n.nodeName, n.functions.get(neighborStates, 0), neighborStates, n.nodeState, systemState[i] ]
        for i, n in enumerate(self.nodes):
            n.nodeState = systemState[i]
        self.__updatePreviousStates(  systemState )

        

    def initialize(self, initialSystemState=None):
        """ resets RBN.

        @param default 0 to randomize state.  otherwise it is a list of initial states
        """
        systemState=[]
        if (initialSystemState):
            for i, node in enumerate(self.nodes):
                node.nodeState = initialSystemState[i]
            systemState=initialSystemState
        else:
            for i, node in enumerate(self.nodes):
                node.nodeState = random.randint(0,1)
                systemState.insert(i, node.nodeState)
        self.__previousSystemStates = []
        self.__updatePreviousStates( systemState )

    def __updatePreviousStates(self, systemState):
        "tracks all the states of the network, for finding attractors"
        self.__previousSystemStates.append( tuple( systemState) )
        if len(self.__previousSystemStates) > len(systemState):
            self.__previousSystemStates.pop( 0) 

    def previousStates(self):
        for systemState in self.__previousSystemStates:
            print systemState

    def currentSystemState(self):
        "returns a tuple fo the current state"
        return self.__previousSystemStates[-1]
        
    def hasAttractor(self):
        """returns false if the RBN has not found an attractor.
        Returns the attractor (a list of tuples) if it has"""
        
        return not len(self.__previousSystemStates) == len(list(set(self.__previousSystemStates)))
    
    def attractor(self):
        """ if an attractor has been found, this function returns it, as a tuple of tuples.  """
        
        if self.hasAttractor():
            pStates = self.__previousSystemStates
            finalSystemState = pStates[-1]
            attractorPeriod = pStates[::-1].index( finalSystemState, 1 )
            attractor = pStates[len(pStates) - attractorPeriod - 1:-1]

            #sort the attractor, to easily spot limit cycles
            order = sorted(attractor)[0]
            while order != attractor[0]:
                t = attractor.pop(0)
                attractor.append(t)
            return tuple(attractor)

    def graph(self):
        """returns a directed igraph graph"""

        g = Graph(directed=True)        # directed edges

        g.add_vertices(len(self.nodes)-1)
        for node in self.nodes:
            for inNode in node.inputs:
                g.add_edges( [ (inNode.nodeNumber, node.nodeNumber) ] )      # graphObject.add_edges( [ (start_vertex, end_vertex) ] )
        return g

    def graphLabels(self):
        "returns a list of the nodeNames for graphing"
        return [n.nodeName for n in self.nodes]


    def stateToInt(self, state):
        "converts a tuple state of the RBN to a unique integer"
        return int(''.join(`x` for x in state), 2)
        
        
    def toPrint(self):
        print [n.nodeState for n in self.nodes]
    def toPrintNames(self):
        print [n.nodeName  for n in self.nodes]
    

class Node:
    "The node of an RBN.  constructor takes its name"
    
    population=0
    def __init__(self, name):
        self.inputs=[]
        self.functions={}
        self.nodeName=name
        self.nodeNumber=Node.population
        self.nodeState=0
        Node.population+=1

        
        
cycD = Node("CycD")
cycE = Node("CycE")
cycA = Node("CycA")
cycB = Node("CycB")
rb = Node("Rb")
p27 = Node("p27")
e2F = Node("E2F")
cdc20 = Node("Cdc20")
cdh1 = Node("Cdh1")
ubcH10 = Node("UbcH10")

cycD.inputs = [cycD]
cycD.functions[(1,)]=1

rb.inputs = [cycD, cycE, cycA, cycB, p27]
rb.functions[(0,   0,    0,    0,    0)] = 1
rb.functions[(0,   0,    0,    0,    1)] = 1
rb.functions[(0,   1,    1,    0,    1)] = 1
rb.functions[(0,   0,    1,    0,    1)] = 1
rb.functions[(0,   1,    0,    0,    1)] = 1

e2F.inputs = [rb, cycA, cycB, p27]
e2F.functions[(0, 0,    0,    0)] = 1 
e2F.functions[(0, 0,    0,    1)] = 1 
e2F.functions[(0, 1,    0,    1)] = 1 
e2F.functions[(0, 0,    0,    1)] = 1 

cycE.inputs = [e2F, rb]
cycE.functions[(1,  0)] = 1

cycA.inputs = [rb, cycA, e2F, cdc20, cdh1, ubcH10]
cycA.functions[(0, 1,    1,   0,     0,    0)] = 1
cycA.functions[(0, 1,    1,   0,     1,    0)] = 1
cycA.functions[(0, 1,    1,   0,     0,    1)] = 1
cycA.functions[(0, 0,    1,   0,     0,    0)] = 1
cycA.functions[(0, 0,    1,   0,     1,    0)] = 1
cycA.functions[(0, 0,    1,   0,     0,    1)] = 1

cycA.functions[(0, 1,    1,   0,     0,    0)] = 1
cycA.functions[(0, 1,    1,   0,     1,    0)] = 1
cycA.functions[(0, 1,    1,   0,     0,    1)] = 1
cycA.functions[(0, 1,    0,   0,     0,    0)] = 1
cycA.functions[(0, 1,    0,   0,     1,    0)] = 1
cycA.functions[(0, 1,    0,   0,     0,    1)] = 1

p27.inputs = [cycD, cycE, cycA, cycB, p27]
p27.functions[(0,   0,    0,    0,    1)] = 1
p27.functions[(0,   0,    0,    0,    0)] = 1
p27.functions[(0,   0,    1,    0,    1)] = 1
p27.functions[(0,   1,    0,    0,    1)] = 1

cdc20.inputs = [cycB]
cdc20.functions[(1,)] = 1

cdh1.inputs = [cycA, cycB, cdc20, p27]
cdh1.functions[(0,   0,    1,     1)] = 1
cdh1.functions[(0,   0,    0,     1)] = 1
cdh1.functions[(0,   0,    1,     0)] = 1
cdh1.functions[(0,   0,    0,     0)] = 1

cdh1.functions[(1,   1,    1,     1)] = 1
cdh1.functions[(1,   1,    1,     0)] = 1
cdh1.functions[(0,   1,    1,     1)] = 1
cdh1.functions[(0,   1,    1,     0)] = 1
cdh1.functions[(1,   0,    1,     1)] = 1
cdh1.functions[(1,   0,    1,     0)] = 1
cdh1.functions[(0,   0,    1,     1)] = 1
cdh1.functions[(0,   0,    1,     0)] = 1

cdh1.functions[(1,   0,    1,     1)] = 1
cdh1.functions[(1,   0,    0,     1)] = 1

ubcH10.inputs = [cdh1, ubcH10, cdc20, cycA, cycB]
ubcH10.functions[(0,  1,      1,     1,    1)] = 1
ubcH10.functions[(0,  1,      1,     1,    0)] = 1
ubcH10.functions[(0,  1,      1,     0,    1)] = 1
ubcH10.functions[(0,  1,      1,     0,    0)] = 1
ubcH10.functions[(0,  1,      0,     1,    1)] = 1
ubcH10.functions[(0,  1,      0,     1,    0)] = 1
ubcH10.functions[(0,  1,      0,     0,    1)] = 1
ubcH10.functions[(0,  1,      0,     0,    0)] = 1
ubcH10.functions[(0,  0,      1,     1,    1)] = 1
ubcH10.functions[(0,  0,      1,     1,    0)] = 1
ubcH10.functions[(0,  0,      1,     0,    1)] = 1
ubcH10.functions[(0,  0,      1,     0,    0)] = 1
ubcH10.functions[(0,  0,      0,     1,    1)] = 1
ubcH10.functions[(0,  0,      0,     1,    0)] = 1
ubcH10.functions[(0,  0,      0,     0,    1)] = 1
ubcH10.functions[(0,  0,      0,     0,    0)] = 1

ubcH10.functions[(1,  1,      1,     1,    1)] = 1
ubcH10.functions[(1,  1,      1,     1,    0)] = 1
ubcH10.functions[(1,  1,      1,     0,    1)] = 1
ubcH10.functions[(1,  1,      0,     1,    1)] = 1
ubcH10.functions[(1,  1,      0,     0,    1)] = 1
ubcH10.functions[(1,  1,      0,     1,    0)] = 1
ubcH10.functions[(1,  1,      1,     0,    0)] = 1

cycB.inputs = [cdc20, cdh1]
cycB.functions[(0,    0)] = 1


print
random.seed( 1 )
attractors = []
graph1 = Graph(directed=True)        # directed edges
counter=0
dSystemStates={}

for i in range(30):
    network1 = RBN()
    network1.nodes = [cycD, cycE ,cycA ,cycB ,rb ,p27 ,e2F, cdc20 ,cdh1, ubcH10 ]
    network1.initialize()
    prevState = network1.stateToInt( network1.currentSystemState() )
    dSystemStates[prevState]=counter
    graph1.add_vertices(1)
    counter+=1
    
    #network1.toPrint()
    while not network1.hasAttractor():
        network1.update()        
        currentState = network1.stateToInt( network1.currentSystemState() )
        
        if  not currentState in dSystemStates:
            dSystemStates[currentState]=counter
            graph1.add_vertices(1)
            counter+=1            
        graph1.add_edges( [ (dSystemStates[prevState], dSystemStates[currentState] ) ] ) # graphObject.add_edges( [ (start_vertex, end_vertex) ] )
        prevState = currentState
    attractors.append( network1.attractor() )
graph1.delete_vertices( counter )
for attractor in list(set(attractors)):
    print attractor
    
network2 = RBN()
network2.nodes = [cycD, cycE ,cycA ,cycB ,rb ,p27 ,e2F, cdc20 ,cdh1, ubcH10 ]
graph2 = network2.graph()
# g.layout_graphopt() is just one layout algorithm,
# others are available @ http://cneurocvs.rmki.kfki.hu/igraph/doc/python/igraph.GraphBase-class.html

#plot(graph2, layout = graph2.layout_graphopt(), vertex_label=network2.graphLabels(),
#     vertex_label_dist=3, vertex_label_angle=-math.pi/6, margin=20 )

plot(graph1, layout = graph1.layout_fruchterman_reingold(), 
     vertex_label_dist=3, vertex_label_angle=-math.pi/6, margin=20 )
</pre>