Difference between revisions of "Template:LMU-TUM Munich ModelCode"

import numpy as np
import matplotlib.pylab as plt
%matplotlib inline

from monomers import *

def init_monomers(n_cells, n_proteins, n_cell_binding = 30, n_protein_binding = 4):
    Cell.n_cells=0
    cells=[]
    for i in range(n_cells):
        cells.append(Cell(n_cell_binding))

    .n_proteins=0
    =[]
    for i in range(n_proteins):
        .append((n_protein_binding))
    return cells, 

def polymerize(cells, ):
    free_cell_bindings = []
    for c in cells:
        free_cell_bindings += c.free_bindings

    free_protein_bindings = []
    for p in :
        free_protein_bindings += p.free_bindings

    while len(free_protein_bindings)>0 and len(free_cell_bindings)>0:
        p_binding_idx = np.random.randint(len(free_protein_bindings))
        p_binding = free_protein_bindings.pop(p_binding_idx)
        c_binding_idx = np.random.randint(len(free_cell_bindings))
        c_binding = free_cell_bindings.pop(c_binding_idx)
        #print c_binding.parent.id, p_binding.parent.id
        p_binding.bind(c_binding)
    return len(free_protein_bindings) , len(free_cell_bindings)

def polymerize_spacial(cells, , n_protein_binding, concentration=0.01, reactivity=0.2, mobility=0.05, max_time=10.):
    time = 0.
    position_variance = 1.
    distance_binding_probability = lambda x: np.exp(-x ** 2 / position_variance ** 2) / (
    2 * np.pi * position_variance ** 2)

    free_protein_bindings = []
    for p in :
        free_protein_bindings += p.free_bindings

    while (time < max_time):
        for i in range(int(np.ceil(reactivity * len(free_protein_bindings)))):
            # choose a random PBS
            p_binding_idx = np.random.randint(len(free_protein_bindings))
            p_binding = free_protein_bindings[p_binding_idx]
            # do tests on parent 
            p = p_binding.parent
            # position: the liwer to concentration, the higher the distance
            p_pos = float(p.id) / (len() * concentration)
            # (float) idx of nearest cell
            c_id_mean = p_pos * len(cells) * concentration
            # max distance of reachable cells from nearest cell
            max_range = position_variance * 3.2
            c_id_range = max_range * len(cells) * concentration
            c_id_min = max(int(c_id_mean - c_id_range), 0)
            c_id_max = min(int(c_id_mean + c_id_range), len(cells))
            candidate_cells = cells[c_id_min:c_id_max]
            # each has a probability according to the distance
            if len(candidate_cells) == 0:
                continue
            min_distance = max(p_pos - position_variance * 3.2, 0) - p_pos
            max_distance = min(-p_pos + position_variance * 3.2,
                               len() - 1 / len() / concentration) + p_pos
            probabilities = np.arange(min_distance, max_distance,
                                      (max_distance - min_distance) / len(candidate_cells))[:len(candidate_cells)]
            probabilities = distance_binding_probability(probabilities)
            probabilities /= sum(probabilities)
            # pick one
            cell = np.random.choice(candidate_cells, p=probabilities)
            try:
                cell.free_bindings.next()
            except StopIteration:
                # no free bindings left in this cell
                continue
            # let react with probability respective to free bindings
            protein_partners = list(set([b.partner.parent for b in p.bindings if b.partner is not None]))
            if len(protein_partners) == 1:
                # if a  is already attatched to a cell it will most likely not bind to other cells
                if cell not in protein_partners and np.random.rand()<0.9:
                    continue
            if len(protein_partners)>1 and cell not in protein_partners:
                # if a  was binding to two cells already, it will not bind to any other
                continue
#            if np.random.rand() < float(4 - len(list(p.free_bindings))) / 3.:
#                # can only bind to cells it's already bound to
#                if cell not in [b.partner.parent for b in p.bindings if b.partner is not None]:
#                    continue
            cell.free_bindings.next().bind(p_binding)
            free_protein_bindings.pop(p_binding_idx)

        position_variance += mobility
        time += 0.01

    free_cell_bindings = []
    for c in cells:
        free_cell_bindings += c.free_bindings
    return len(free_protein_bindings), len(free_cell_bindings)

draw = True

_,_,_=random_graph(cell_graph, 30, 120, n_cell_binding=400, n_protein_binding=4, concentration=0.02, reactivity=0.1, mobility=0.05, max_time=5.0)

draw = False

from graph_tool.all import *

from graph_tool import clustering,stats

def cell_graph(cells,,draw=False):
    '''vertices are cells, cells connected via a  are connected in the graph with an edge'''
    g = Graph(directed=False)

    vertices = list(g.add_vertex(len(cells)))

    for c_id in range(len(cells)):
        for c_b in cells[c_id].bindings:
            try:
                for p_b in c_b.partner.parent.bindings:
                    try:
                        c2_id = p_b.partner.parent.id
                        if not c_id==c2_id:
                            g.add_edge(vertices[c_id],vertices[c2_id])
                    except AttributeError:
                        pass
            except AttributeError:
                pass
    if draw:
        graph_draw(g, vertex_text=g.vertex_index, edge_pen_width=2 , vertex_font_size=30, output="/tmp/graph.png",
            output_size=(2000, 2000))
    return g

#graph_draw(g, vertex_text=g.vertex_index, vertex_font_size=18, output="/tmp/graph.png",
#            output_size=(2000, 2000))

def random_graph(graph_method, n_cells, n_proteins, n_cell_binding = 30, n_protein_binding = 4,concentration=0.02,reactivity=0.01,mobility=0.1,max_time=2.):
#    cells,=init_monomers(n_cells, n_proteins, n_cell_binding, n_protein_binding)
#    polymerize(cells,)
    cells,=init_monomers(n_cells, n_proteins, n_cell_binding = 30, n_protein_binding = 4)
    polymerize_spacial(cells, , n_protein_binding=n_protein_binding, concentration=concentration, reactivity=reactivity, mobility=mobility, max_time=max_time)
    return graph_method(cells,,draw),cells,

results=[]
draw=False
for n_proteins in np.arange(-1.5,1.5,0.1):
    n_proteins = int(len(cells)*10**n_proteins)
    g,cells, = random_graph(cell_graph, 30, n_proteins, n_cell_binding=400, n_protein_binding=4, concentration=0.02, reactivity=0.1, mobility=0.05, max_time=5.0)
    
    global_clustering = clustering.global_clustering(g)[0]
    
    tmp=graph_tool.topology.max_cardinality_matching(g)
    n_max_matching = tmp.count(tmp)
    
    max_shortest_path = graph_tool.topology.pseudo_diameter(g)[0]
    
    tmp=graph_tool.topology.label_largest_component(g)
    biggest_cluster = tmp.get_array().shape[0]
    
    tmp = graph_tool.topology.label_components(g)[1]
    tmp = np.sort(tmp)[::-1]
    n_clusters = len(tmp)
    biggest_cluster = tmp[0]
    mean_cluster = np.mean(tmp)
    
    
    n_free_protein_bindings = 0
    n_free_proteins = 0
    for p in :
        p_free = len(list(p.free_bindings))
        n_free_protein_bindings += p_free
        
        if p_free == len(p.bindings):
            n_free_proteins += 1
    
    n_free_cell_bindings = 0
    n_free_cells = 0
    for c in cells:
        c_free = len(list(c.free_bindings))
        n_free_cell_bindings += c_free
        if c_free == len(c.bindings):
            n_free_cells += 1
    
    print n_proteins, n_free_cells, n_free_cell_bindings, n_free_proteins, n_free_protein_bindings, n_clusters, mean_cluster, biggest_cluster, global_clustering, max_shortest_path#, n_max_matching, biggest_cluster
    results.append([n_proteins, n_free_cells, n_free_cell_bindings, n_free_proteins, n_free_protein_bindings, n_clusters, mean_cluster, biggest_cluster, global_clustering, max_shortest_path])
    

tmp_res=results

results=np.asarray(results).T
n_proteins=results[0]

results=results[1:]

titles=["n_free_cells", "n_free_cell_bindings", "n_free_proteins", "n_free_protein_bindings", "n_clusters", "mean_cluster", "biggest_cluster", "global_clustering", "max_shortest_path"]
f, axes = plt.subplots(len(titles)/2+1,2, figsize=(16,8*(len(titles)/2+1)))
for i,result in enumerate(results):
    ax = axes[i/2,i%2]
    ax.set_title(titles[i])
    ax.plot(n_proteins,result.T,"-x")
    ax.set_xscale("log")

clustering.global_clustering(g)

stats.distance_histogram(g)

tmp=graph_tool.topology.max_cardinality_matching(g)
tmp.count(tmp)

graph_tool.topology.max_independent_vertex_set(g)

graph_tool.topology.pseudo_diameter(g)

tmp=graph_tool.topology.label_largest_component(g)
tmp.get_array().shape[0]

g = Graph(directed=False)

v_size = g.new_vertex_property("int",val=20)
v_group = g.new_vertex_property("int32_t")

#v_color = vprop_vint = g.new_vertex_property("<int>", vals=[])

#cell_vertices = list(g.add_vertex(n_cells))
for cell in cells:
    cell.vertex = g.add_vertex()
    v_size[cell.vertex] = 100

for p in :
    last_vertex = None
    for b in p.bindings:
        b.vertex=g.add_vertex()
        v_size[b.vertex] = 10
        v_group[cell.vertex] = p.id

        #intraprotein connections
        if last_vertex is not None:
            g.add_edge(b.vertex,last_vertex)
        last_vertex=b.vertex
        #cell connections
        if b.partner is not None:
            cell_vertex = b.partner.parent.vertex
            g.add_edge(b.vertex,cell_vertex)
            
        
pos=sfdp_layout(g,groups=v_group)
#graph_draw(g, vertex_text=g.vertex_index, vertex_font_size=18, output="/tmp/graph.png",
#            output_size=(2000, 2000))
graph_draw(g, pos=pos,edge_pen_width=2, vertex_text=g.vertex_index, vertex_size=v_size, vertex_font_size=0, output_size=(2000, 2000))
#graphviz_draw(g, elen=100, vsize=v_size, size=(2000, 2000))

def pc_graph(cells, , draw=False):
    g = Graph(directed=False)

    v_size = g.new_vertex_property("int")
    #v_group = g.new_vertex_property("i'nt32_t")
    #v_color = vprop_vint = g.new_vertex_property("<int>", vals=[])

    for cell in cells:
        cell.vertex = g.add_vertex()
        v_size[cell.vertex] = 120

    for p in :
        p.vertex=g.add_vertex()
        v_size[p.vertex] = 40
        for b in p.bindings:
            #cell connections
            if b.partner is not None:
                cell_vertex = b.partner.parent.vertex
                if g.edge(p.vertex,cell_vertex) is None: #only one connection per pair
                    g.add_edge(p.vertex,cell_vertex)

    pos=fruchterman_reingold_layout(g)
    #graph_draw(g, vertex_text=g.vertex_index, vertex_font_size=18, output="/tmp/graph.png",
    #            output_size=(2000, 2000))
    if draw:
        graph_draw(g, edge_pen_width=2, vertex_text=g.vertex_index, vertex_size=v_size, vertex_font_size=15, output_size=(2000, 2000))
    return g

g.vertex_index[vertices[1]]

cells[23].bindings[5].partner.parent.bindings[0].partner.parent